library(ResistanceGA2)
library(terra)

Overview

This vignette sketches a practical workflow for fitting isolation by environment (IBE) and isolation by resistance (IBR) concurrently with the new MLPE interface in ResistanceGA2.

The central idea is simple:

IBE enters the model as one or more pairwise endpoint-difference covariates.
IBR enters the model as the optimized resistance-distance term.
both are estimated in the same MLPE framework so we can ask whether one still explains genetic differentiation after controlling for the other.

This example uses the package’s Idaho demo data:

rga2_demo for sites, extracted site environments, and the chord-distance matrix
rga2_demo_covariates for the raster layers used to generate resistance distances

The analysis is structured so that the expensive optimization steps can be run outside the vignette, saved as package extdata, and then displayed here for later exploration.

The document is organized in two layers:

a separate builder script that generates the full analysis objects and output files
this vignette, which then loads those saved results by default for display and interpretation

The optimization examples below use an explicit seed so the worked example is reproducible. The external builder script also supports named runtime profiles and an optional core-count override so you can choose between quick smoke tests and deeper optimization runs without editing the script itself.

if (!have_saved_outputs) {
  cat(
    "Saved analysis outputs were not found under ",
    output_dir,
    ".\n\nRun `data-raw/build_ibe_ibr_vignette_outputs.R` first if you want ",
    "the full set of model tables, optimization summaries, and plots to be ",
    "stored beneath `inst/extdata/ibe-ibr-demo/` and appear throughout this ",
    "vignette.\n",
    sep = ""
  )
}

Design Choices

For a joint IBE/IBR analysis, it is important not to blur the mechanism of endpoint mismatch with the mechanism of landscape permeability. In this example:

IBE is represented by site-level mismatch in cti and hli
IBR is optimized over the separate raster hypotheses err27 and ffp
geographic distance is included as an additional covariate

This is only one defensible partition. The key point is to avoid using the same variable in both roles unless that dual interpretation is explicitly part of the biological hypothesis.

Demo Data

data(rga2_demo)
data(rga2_demo_covariates)

if (inherits(rga2_demo_covariates, "PackedSpatRaster")) {
  rga2_demo_covariates <- unwrap(rga2_demo_covariates)
}

r_covariate <- subset(rga2_demo_covariates, c("err27", "ffp", "cti", "hli"))

sites <- rga2_demo$sites
site_env <- rga2_demo$site_environment

sites[, c("demo_id", "SiteName", "x", "y")]
#>               demo_id       SiteName        x       y
#> Airplane     Airplane   AirplaneLake 688816.6 5003207
#> Bachelor     Bachelor BachelorMeadow 688494.4 4999093
#> BarkingFox BarkingFox BarkingFoxLake 687938.4 5000223
#> Bob               Bob        BobLake 690104.0 4999355
#> Cache           Cache      CacheLake 688742.5 4997481
#> Egg               Egg   EggWhiteLake 688539.3 4998146
#> Frog             Frog   FrogPondLake 688845.2 4998427
#> GentianL     GentianL    GentianLake 689781.6 5002972
#> ParagonL     ParagonL    ParagonLake 687421.5 4994948
#> Pothole       Pothole    PotholeLake 690481.4 5004352
#> ShipIsland ShipIsland ShipIslandLake 686908.5 5004064
#> Skyhigh       Skyhigh    SkyhighLake 688159.7 4999027
site_env
#>               demo_id       cti     err27       ffp      gsp  hli nlcd
#> Airplane     Airplane  6.551698 0.2211860  8.000000 286.0000 1438   42
#> Bachelor     Bachelor  7.257386 0.3399709  8.000000 286.0000 1718   71
#> BarkingFox BarkingFox  7.587409 0.2636720  9.000000 283.0000 1308   42
#> Bob               Bob 10.962278 0.5830295  6.000000 290.0000 1685   42
#> Cache           Cache  7.817918 0.4895256 11.000000 274.0000 1277   42
#> Egg               Egg  6.560551 0.4499388 10.000000 279.0000 1244   42
#> Frog             Frog  6.675704 0.3340760 10.155787 276.1548 1170   42
#> GentianL     GentianL  5.595447 0.2247885  6.000000 292.0000 1370   42
#> ParagonL     ParagonL  9.865539 0.4723573  7.000000 285.3084 1298   11
#> Pothole       Pothole 10.490610 0.4085309  6.000000 292.0000 1368   42
#> ShipIsland ShipIsland  8.181678 0.1197917 13.000000 273.0000 1300   11
#> Skyhigh       Skyhigh 12.254863 0.2937623  6.829502 290.0000 1316   11
r_covariate
#> class       : SpatRaster
#> size        : 426, 358, 4  (nrow, ncol, nlyr)
#> resolution  : 30, 30  (x, y)
#> extent      : 683282.5, 694022.5, 4992833, 5005613  (xmin, xmax, ymin, ymax)
#> coord. ref. : NAD83 / UTM zone 11N (EPSG:26911)
#> source(s)   : memory
#> varname     : covariates
#> names       :    err27, ffp,      cti,  hli
#> min values  : 0.039066,   0, 0.842985, 1014
#> max values  : 0.763764,  51, 23.71476, 9263

Step 1: Build Pairwise Geographic and Environmental Covariates

The pairwise response is the lower triangle of the 12 x 12 chord-distance matrix. We pair that response with:

Euclidean distance among sites (IBD)
pairwise endpoint differences in cti and hli (IBE)

There are two common ways to encode IBE, and they answer slightly different questions.

Absolute environmental differences:

use one pairwise covariate per environmental variable, such as |cti_i - cti_j| or |hli_i - hli_j|
keep interpretation simple, because each coefficient refers to a specific environmental mismatch axis
are the best starting point when you have a small set of clear, named hypotheses

Multivariate environmental distance:

compress several site-level variables into a single pairwise distance term
is useful when you want an overall environmental dissimilarity index rather than separate coefficients
is harder to interpret mechanistically, because one coefficient summarizes several variables at once

In practice, I would usually start with absolute differences for a small number of biologically motivated variables, then compare that model to a simpler multivariate-distance version as a sensitivity check.

That strategy is especially useful here because cti and hli have interpretable ecological meanings, while env_dist is intentionally more agnostic. Separate absolute-difference terms tell us which environmental axis is carrying signal; a multivariate distance tells us whether there is overall environmental mismatch even when we do not want to interpret individual axes too literally.

site_ids <- sites$demo_id
pair_index <- t(combn(site_ids, 2))
pair_index <- data.frame(
  from = pair_index[, 1],
  to = pair_index[, 2],
  stringsAsFactors = FALSE
)

coords <- rga2_demo$sample_coords
geo <- as.numeric(scale(lower(as.matrix(dist(coords)))))

ibe_diff <- site_pairwise_diff(
  site_env[, c("cti", "hli")],
  scale = TRUE
)

env_dist <- site_pairwise_dist(
  site_env[, c("cti", "hli")],
  scale = TRUE
)

head(ibe_diff)
#>      cti_diff  hli_diff
#> 1 0.338346773 1.6785034
#> 2 0.496578284 0.7793052
#> 3 2.114680594 1.4806798
#> 4 0.607097178 0.9651395
#> 5 0.004244835 1.1629631
#> 6 0.059455814 1.6065676
head(env_dist)
#>    env_dist
#> 1 1.7122653
#> 2 0.9240706
#> 3 2.5815280
#> 4 1.1402023
#> 5 1.1629708
#> 6 1.6076674

ibe_diff and env_dist are both valid IBE representations, but they should not usually be included in the same model because they encode overlapping information.

Step 1a: Check Predictor Overlap Before Fitting Models

Before interpreting coefficients, it helps to inspect how strongly the pairwise predictors overlap. If IBD, IBE, and later IBR are all highly correlated, then model ranking may still be informative but individual fixed-effect estimates will be harder to interpret mechanistically.

predictor_cor <- cor(
  data.frame(
    geo = geo,
    cti_diff = ibe_diff$cti_diff,
    hli_diff = ibe_diff$hli_diff,
    env_dist = env_dist$env_dist
  )
)

round(predictor_cor, 3)
#>             geo cti_diff hli_diff env_dist
#> geo       1.000   -0.095   -0.211   -0.269
#> cti_diff -0.095    1.000   -0.095    0.606
#> hli_diff -0.211   -0.095    1.000    0.701
#> env_dist -0.269    0.606    0.701    1.000

If you have already run the external builder script, the saved output bundle stores this same matrix for later review:

if (have_saved_outputs) {
  round(summary_bundle$predictor_correlation, 3)
}
#>             geo cti_diff hli_diff env_dist
#> geo       1.000   -0.095   -0.211   -0.269
#> cti_diff -0.095    1.000   -0.095    0.606
#> hli_diff -0.211   -0.095    1.000    0.701
#> env_dist -0.269    0.606    0.701    1.000

Step 2: Fit Baseline MLPE Models Before Optimization

Before optimizing any resistance surfaces, it is useful to ask whether there is support for IBD, IBE, or both.

pair_data <- mlpe_data(
  response = rga2_demo$genetic_vector,
  pairs = pair_index,
  geo = geo,
  cti_diff = ibe_diff$cti_diff,
  hli_diff = ibe_diff$hli_diff,
  env_dist = env_dist$env_dist
)

head(pair_data)
#>       from         to        geo    cti_diff  hli_diff  env_dist  response
#> 1 Airplane   Bachelor  0.1491438 0.338346773 1.6785034 1.7122653 0.4827800
#> 2 Airplane BarkingFox -0.3013249 0.496578284 0.7793052 0.9240706 0.4133083
#> 3 Airplane        Bob  0.1204355 2.114680594 1.4806798 2.5815280 0.4093260
#> 4 Airplane      Cache  0.8583245 0.607097178 0.9651395 1.1402023 0.4977017
#> 5 Airplane        Egg  0.5669958 0.004244835 1.1629631 1.1629708 0.4873166
#> 6 Airplane       Frog  0.4387835 0.059455814 1.6065676 1.6076674 0.4530400

Optional: Multiple Dyadic Random Effects

The examples above use one MLPE dyad term, (1 | pair), where pair is built from the from and to site columns. The current mlpe() implementation can also include more than one separate dyadic random-intercept term in the same model, as long as each term is linked to its own pair of endpoint columns through the pairs argument.

This can be useful when pairwise observations share endpoints at multiple sampling scales. For example, we may want one dyadic random effect for repeated use of the same sampled sites and a second dyadic random effect for repeated use of broader region labels attached to those sites.

One current constraint is that each dyadic term still expects the two endpoint labels to differ within every row for that term. A coarse region-level dyadic effect therefore works naturally for between-region pairs, but not for within-region rows where from_region == to_region. The example below keeps only between-region observations so the two-term fit is well-defined.

site_region <- setNames(
  ifelse(sites$x <= stats::median(sites$x), "west", "east"),
  as.character(sites$demo_id)
)

pair_data_region <- pair_data
pair_data_region$from_region <- site_region[as.character(pair_data_region$from)]
pair_data_region$to_region <- site_region[as.character(pair_data_region$to)]
pair_data_region <- pair_data_region[
  pair_data_region$from_region != pair_data_region$to_region,
  ,
  drop = FALSE
]
rownames(pair_data_region) <- NULL

m_site_region <- mlpe(
  response ~ geo + cti_diff + hli_diff +
    (1 | site_pair) + (1 | region_pair),
  data = pair_data_region,
  pairs = list(
    site_pair = c("from", "to"),
    region_pair = c("from_region", "to_region")
  )
)

summary(m_site_region)
#> Linear mixed model fit by maximum likelihood  ['lmerMod']
#> 
#>       AIC       BIC    logLik -2*log(L)  df.resid 
#>     -77.7     -66.6      45.9     -91.7        29 
#> 
#> Scaled residuals: 
#>      Min       1Q   Median       3Q      Max 
#> -2.30629 -0.53670 -0.03469  0.82154  1.99303 
#> 
#> Random effects:
#>  Groups      Name        Variance  Std.Dev.
#>  site_pair   (Intercept) 0.0002541 0.01594 
#>  region_pair (Intercept) 0.0000000 0.00000 
#>  Residual                0.0041339 0.06430 
#> Number of obs: 36, groups:  site_pair, 12; region_pair, 2
#> 
#> Fixed effects:
#>              Estimate Std. Error t value
#> (Intercept)  0.393343   0.026202  15.012
#> geo          0.073851   0.011699   6.312
#> cti_diff    -0.003229   0.014014  -0.230
#> hli_diff    -0.009549   0.013405  -0.712
#> 
#> Correlation of Fixed Effects:
#>          (Intr) geo    ct_dff
#> geo      -0.098              
#> cti_diff -0.656  0.003       
#> hli_diff -0.606  0.142  0.127

The important distinction is that these are parallel dyadic terms, not a nested dyadic syntax. Formula shorthand such as (1 | site_pair / region_pair) is not currently translated into a nested MLPE covariance structure, and dyadic terms with same-label endpoints are still rejected. If you need both site-level and region-level dyadic effects right now, write them as separate terms and ensure the retained rows satisfy that endpoint-label constraint.

m_null <- mlpe(
  response ~ 1 + (1 | pair),
  data = pair_data,
  pairs = c("from", "to")
)

m_ibd <- mlpe(
  response ~ geo + (1 | pair),
  data = pair_data,
  pairs = c("from", "to")
)

m_ibe_cti <- mlpe(
  response ~ cti_diff + (1 | pair),
  data = pair_data,
  pairs = c("from", "to")
)

m_ibe_hli <- mlpe(
  response ~ hli_diff + (1 | pair),
  data = pair_data,
  pairs = c("from", "to")
)

m_ibe_both <- mlpe(
  response ~ cti_diff + hli_diff + (1 | pair),
  data = pair_data,
  pairs = c("from", "to")
)

m_ibd_ibe <- mlpe(
  response ~ geo + cti_diff + hli_diff + (1 | pair),
  data = pair_data,
  pairs = c("from", "to")
)

m_ibd_envdist <- mlpe(
  response ~ geo + env_dist + (1 | pair),
  data = pair_data,
  pairs = c("from", "to")
)

baseline_table <- data.frame(
  model = c(
    "Null",
    "IBD",
    "IBE: cti",
    "IBE: hli",
    "IBE: cti + hli",
    "IBD + IBE",
    "IBD + env_dist"
  ),
  AIC = c(
    AIC(m_null),
    AIC(m_ibd),
    AIC(m_ibe_cti),
    AIC(m_ibe_hli),
    AIC(m_ibe_both),
    AIC(m_ibd_ibe),
    AIC(m_ibd_envdist)
  ),
  logLik = c(
    logLik(m_null),
    logLik(m_ibd),
    logLik(m_ibe_cti),
    logLik(m_ibe_hli),
    logLik(m_ibe_both),
    logLik(m_ibd_ibe),
    logLik(m_ibd_envdist)
  )
)

baseline_table[order(baseline_table$AIC), ]

This baseline comparison tells us whether the environmental mismatch terms have signal before we add the path-based IBR component.

The different baseline models answer different questions:

Null tells us how much structure is captured by the MLPE random effect alone
IBD asks whether simple geographic separation explains genetic distance
IBE: cti and IBE: hli isolate each environmental mismatch axis
IBE: cti + hli asks whether a small set of interpretable mismatch variables improves fit together
IBD + IBE tests whether environmental mismatch adds anything after accounting for geography
IBD + env_dist is a compact alternative when you prefer one overall environmental dissimilarity term

If saved outputs are available, the table below shows an actual model ranking from the example analysis.

if (have_saved_outputs) {
  summary_bundle$baseline_table
}
#>         model       AIC   logLik      sigma  deltaAIC       weight
#> 2         ibd -170.8078 89.40392 0.05695819  0.000000 6.273086e-01
#> 7 ibd_envdist -169.0868 89.54341 0.05686835  1.721012 2.653190e-01
#> 6     ibd_ibe -167.2776 89.63878 0.05674323  3.530271 1.073724e-01
#> 1        null -122.7880 64.39402 0.08169552 48.019802 2.344842e-11
#> 3     ibe_cti -121.4812 64.74061 0.08123091 49.326610 1.219955e-11
#> 4     ibe_hli -121.3974 64.69872 0.08219886 49.410390 1.169907e-11
#> 5    ibe_both -120.2543 65.12715 0.08173548 50.553528 6.605749e-12

This fixed-effect table is often even more informative than the rank table, because it shows how much the environmental coefficients shrink once geography enters the model.

if (have_saved_outputs) {
  summary_bundle$baseline_fixef
}
#>          model        term     estimate
#> 1         null (Intercept)  0.388676042
#> 2          ibd (Intercept)  0.388676042
#> 3          ibd         geo  0.073303219
#> 4      ibe_cti (Intercept)  0.402725418
#> 5      ibe_cti    cti_diff -0.012131977
#> 6      ibe_hli (Intercept)  0.401402489
#> 7      ibe_hli    hli_diff -0.011979799
#> 8     ibe_both (Intercept)  0.418607528
#> 9     ibe_both    cti_diff -0.013491980
#> 10    ibe_both    hli_diff -0.013467778
#> 11     ibd_ibe (Intercept)  0.399096675
#> 12     ibd_ibe         geo  0.072585107
#> 13     ibd_ibe    cti_diff -0.003117198
#> 14     ibd_ibe    hli_diff -0.006411196
#> 15 ibd_envdist (Intercept)  0.398006270
#> 16 ibd_envdist         geo  0.072466306
#> 17 ibd_envdist    env_dist -0.005231526

For this example, the fast-run saved outputs suggest a simple but instructive pattern:

IBD is clearly supported relative to Null
IBE by itself is weaker than IBD
adding cti_diff and hli_diff to IBD changes the coefficients only modestly in the fast run, suggesting that most of the baseline signal is geographic rather than environmental
the multivariate env_dist model behaves similarly, which makes it a useful sensitivity check rather than a replacement for the more interpretable univariate terms

Step 3: Prepare the Conditional IBR Workflow

Once the fixed covariates are chosen, they can be supplied directly to gdist.prep() so every resistance optimization is fit conditionally on those same IBD and IBE terms.

pts <- vect(
  rga2_demo$sample_coords,
  type = "points",
  crs = crs(r_covariate, proj = TRUE)
)

workflow_covariates <- data.frame(
  geo = geo,
  cti_diff = ibe_diff$cti_diff,
  hli_diff = ibe_diff$hli_diff
)

gdist_inputs <- gdist.prep(
  n.Pops = nrow(rga2_demo$sample_coords),
  response = rga2_demo$genetic_vector,
  samples = pts,
  covariates = workflow_covariates,
  formula = gd ~ geo + cti_diff + hli_diff,
  method = "costDistance"
)

gdist_inputs$formula
#> gd ~ geo + cti_diff + hli_diff + cd + (1 | pop)
head(gdist_inputs$df)
#>          gd        geo    cti_diff  hli_diff pop .mlpe_pop1 .mlpe_pop2
#> 1 0.4827800  0.1491438 0.338346773 1.6785034   1          1          2
#> 2 0.4133083 -0.3013249 0.496578284 0.7793052   1          1          3
#> 3 0.4093260  0.1204355 2.114680594 1.4806798   1          1          4
#> 4 0.4977017  0.8583245 0.607097178 0.9651395   1          1          5
#> 5 0.4873166  0.5669958 0.004244835 1.1629631   1          1          6
#> 6 0.4530400  0.4387835 0.059455814 1.6065676   1          1          7

The resulting optimization models are of the form:

gd ~ geo + cti_diff + hli_diff + cd + (MLPE random effects)

where cd is the resistance-distance term generated from the candidate raster(s).

This conditional formulation is the key to the joint framework. We are not asking whether IBR is better than IBE in isolation; we are asking whether a path-based resistance hypothesis still improves model fit after geography and endpoint environmental mismatch have already been accounted for.

Step 4: Run Single-Surface IBR Models Conditional on IBE

This example treats err27 and ffp as the candidate permeability hypotheses.

ibr_layers <- subset(r_covariate, c("err27", "ffp"))

GA_inputs_ss <- GA.prep(
  raster = ibr_layers,
  Results.dir = file.path(tempdir(), "rga2-ibe-ibr-single"),
  method = "LL",
  seed = 123,
  pop.size = 25,
  maxiter = 15,
  run = 5, 
  parallel = F, 
  keepBest = FALSE,
  monitor = FALSE,
  quiet = TRUE
)

ss_out <- SS_optim(
  gdist.inputs = gdist_inputs,
  GA.inputs = GA_inputs_ss,
  dist_mod = TRUE,
  null_mod = TRUE,
  diagnostic_plots = FALSE
)

ss_out$AICc
summary(ss_out)

That table compares:

the optimized IBR + IBD + IBE model for each candidate raster
a pure geographic Distance model
a Null model with no fixed predictors beyond the intercept

If saved outputs are available, we can inspect the actual optimization ranking:

if (have_saved_outputs) {
  summary_bundle$ss_aicc
}
#>     Surface                  Fixed.Effects obj.func_LL k       AIC      AICc
#> 12     Null      geo + cti_diff + hli_diff    89.63878 4 -171.2776 -165.5633
#> 11 Distance geo + cti_diff + hli_diff + cd    89.91863 5 -169.8373 -159.8373
#> 2       ffp geo + cti_diff + hli_diff + cd   100.62203 7 -187.2441 -159.2441
#> 1     err27 geo + cti_diff + hli_diff + cd    98.51117 7 -183.0223 -155.0223
#>          R2m       R2c        LL
#> 12 0.5961967 0.6527893  89.63878
#> 11 0.5988908 0.6564737  89.91863
#> 2  0.7003150 0.7169291 100.62203
#> 1  0.6926207 0.7458538  98.51117

The single-surface summary is also helpful because it packages the model ranking with a concise narrative about the best-supported surface.

if (have_saved_outputs) {
  summary(ss_out)
}
#> ResistanceGA2 single-surface optimization
#>   models evaluated: 4 
#>   categorical surfaces: 0 
#>   continuous surfaces: 2 
#>   best surface: Null 
#>   best AICc: -165.563 
#>   run time (min): 9.3415 
#> 
#>   Surface                  Fixed.Effects obj.func_LL k       AIC      AICc
#>      Null      geo + cti_diff + hli_diff    89.63878 4 -171.2776 -165.5633
#>  Distance geo + cti_diff + hli_diff + cd    89.91863 5 -169.8373 -159.8373
#>       ffp geo + cti_diff + hli_diff + cd   100.62203 7 -187.2441 -159.2441
#>     err27 geo + cti_diff + hli_diff + cd    98.51117 7 -183.0223 -155.0223
#>        R2m       R2c        LL
#>  0.5961967 0.6527893  89.63878
#>  0.5988908 0.6564737  89.91863
#>  0.7003150 0.7169291 100.62203
#>  0.6926207 0.7458538  98.51117
#> 
#> Fixed effects for the top model
#>         term     Estimate  Std. Error    t value
#>  (Intercept)  0.399096675 0.022359888 17.8487782
#>          geo  0.072585107 0.008509314  8.5300775
#>     cti_diff -0.003117198 0.010086539 -0.3090454
#>     hli_diff -0.006411196 0.009951220 -0.6442623

And the corresponding plot provides a quick visual ranking of the candidate surfaces:

if (have_saved_outputs) {
  plot(ss_out, type = "ranking")
}

The exact top-ranked surface can vary across systems, and with shallow search settings it can also vary across exploratory runs if no seed is set. The key interpretive question is whether at least one candidate resistance surface outranks the pure Distance model while IBD + IBE are already in the fixed effects. If it does, the sampled data support an additional conditional association with path-based landscape distance.

Step 5: Run a Multi-Surface IBR Model Conditional on IBE

If the single-surface results suggest complementary hypotheses, a multisurface fit lets us ask whether a composite IBR term improves support further.

GA_inputs_ms <- GA.prep(
  raster = ibr_layers,
  Results.dir = file.path(tempdir(), "rga2-ibe-ibr-multi"),
  method = "LL",
  seed = 123,
  pop.size = 25,
  maxiter = 15,
  run = 5,
  parallel = FALSE,
  keepBest = FALSE,
  monitor = FALSE,
  quiet = TRUE
)

ms_out <- MS_optim(
  gdist.inputs = gdist_inputs,
  GA.inputs = GA_inputs_ms,
  diagnostic_plots = FALSE
)

ms_out$AICc.tab
ms_out$percent.contribution
summary(ms_out)

Multisurface optimization is especially helpful when two rasters are plausible for different ecological reasons. It lets us ask whether the best explanation is truly composite or whether one layer dominates once they are optimized together.

if (have_saved_outputs) {
  summary(ms_out)
}
#> ResistanceGA2 multisurface optimization
#>   best surface: err27.ffp 
#> 
#> Best-fit table
#>    Surface                  Fixed.Effects obj.func_LL  k       AIC     AICc
#>  err27.ffp geo + cti_diff + hli_diff + cd     98.3008 10 -176.6016 43.39841
#>        R2m       R2c      LL
#>  0.6934124 0.7489791 98.3008
#> 
#> Percent contribution
#>  surface       mean
#>    err27 0.98589649
#>      ffp 0.01410351
#> 
#> Fixed effects
#>         term      Estimate  Std. Error     t value
#>  (Intercept)  0.3875158596 0.023269855 16.65312747
#>          geo -0.1165464117 0.044079401 -2.64401080
#>     cti_diff  0.0001726458 0.009223669  0.01871769
#>     hli_diff  0.0009039138 0.009675719  0.09342084
#>           cd  0.1910451025 0.043676807  4.37406289

if (have_saved_outputs) {
  plot(ms_out, type = "contribution")
}

The percent-contribution table is often more informative than the raw multisurface fit by itself. It tells us whether two plausible resistance hypotheses are both contributing meaningfully or whether one layer dominates once they are optimized together. That is scientifically useful because it shows that not every plausible raster necessarily adds unique explanatory power in the joint model.

Step 6: Generate Full Analysis Outputs Outside the Vignette

The intended workflow is to run a separate builder script once, then let this vignette pull those saved outputs back in for a fully populated worked example.

The package includes a builder script that runs the baseline MLPE models plus single- and multi-surface optimization and saves the results as .rds and .csv files beneath inst/extdata/ibe-ibr-demo/ by default so the completed worked example can ship with the package.

Run from the package root:

source("data-raw/build_ibe_ibr_vignette_outputs.R")
build_ibe_ibr_demo_outputs()

From an RStudio session, you can specify the parameters directly:

source("data-raw/build_ibe_ibr_vignette_outputs.R")

build_ibe_ibr_demo_outputs(
  profile = "default",
  cores = NULL,
  seed = 123
)

For a quick smoke test:

source("data-raw/build_ibe_ibr_vignette_outputs.R")

build_ibe_ibr_demo_outputs(
  profile = "fast",
  cores = NULL,
  seed = 123
)

For a deeper run:

source("data-raw/build_ibe_ibr_vignette_outputs.R")

build_ibe_ibr_demo_outputs(
  profile = "full",
  cores = NULL,
  seed = 123
)

The function arguments are:

profile = "fast" for quick smoke tests
profile = "default" for a moderate reproducible example
profile = "full" for a deeper optimization run
cores as an optional positive integer passed to GA.prep(parallel = ...)
seed for reproducible GA initialization

The cores argument is now passed directly through to GA.prep(parallel = ...), including on Windows. For large optimization runs, this is the main knob for scaling the workflow across available hardware.

The command-line interface still works if you prefer Rscript, but the function-based RStudio workflow is the recommended path for interactive use.

The script writes:

pair_data.rds
workflow_covariates.rds
baseline_models.rds
baseline_model_table.csv
baseline_fixed_effects.csv
single_surface/Results/ with optimizer summaries, rasters, and distance matrices
single_surface_optim.rds
single_surface_aicc.csv
multi_surface/Results/ with multisurface summaries and percent contribution
multi_surface_optim.rds
multi_surface_aicc.csv
analysis_summary.rds

Step 7: Load and Explore Saved Outputs

output_dir <- system.file("extdata", "ibe-ibr-demo", package = "ResistanceGA2")

baseline_models <- readRDS(file.path(output_dir, "baseline_models.rds"))
ss_out <- readRDS(file.path(output_dir, "single_surface_optim.rds"))
ms_out <- readRDS(file.path(output_dir, "multi_surface_optim.rds"))
summary_bundle <- readRDS(file.path(output_dir, "analysis_summary.rds"))

summary_bundle$baseline_table
summary_bundle$ss_aicc
summary_bundle$ms_aicc

plot(ss_out, type = "ranking")
plot(ms_out, type = "contribution")

Interpretation

This framework lets us separate three related questions:

Does genetic differentiation increase with geographic distance alone?
Does endpoint environmental mismatch explain additional variation?
Does a resistance-distance hypothesis still improve fit after controlling for both geography and environment?

Possible outcomes are:

IBE improves fit but IBR does not: endpoint mismatch matters more than permeability.
IBR improves fit after conditioning on IBE: there is support for a genuine landscape-connectivity effect.
both improve fit: connectivity and environmental mismatch likely contribute independently.
coefficients become unstable or change sign across models: the predictors are probably too correlated for strong mechanistic interpretation.

The example here also shows why it is useful to explore several model classes instead of going straight to optimization:

baseline MLPE models tell us whether IBE is promising enough to include at all
conditional single-surface optimization tells us whether specific IBR hypotheses improve fit beyond IBD + IBE
multisurface optimization tells us whether apparently useful resistance layers still contribute once they compete in the same model

That sequence turns the vignette into more than a software tutorial. It becomes an analysis template for asking:

What is the best-supported baseline explanation of pairwise genetic structure?
Does environmental mismatch add anything to geography?
Does landscape permeability add anything after both of those are already in the model?
If multiple resistance hypotheses are plausible, do they all matter, or does one dominate?

The same workflow can be extended by:

replacing cti and hli with other site-level environmental variables
adding a single multivariate env_dist term instead of several univariate mismatch terms
moving from gdistance to jl.prep() if Circuitscape-based pair filtering or current maps are required

Quick reference: conditional IBE and IBR workflow

library(ResistanceGA2)
library(terra)

# 1. Build pairwise geographic and environmental covariates.
pts <- terra::vect(samples[, 2:3], type = "points")
pair_data <- mlpe_data(
  response = lower(as.matrix(dist(samples[, 2:3]))),
  distance = site_pairwise_dist(terra::crds(pts), name = "distance")$distance
)

site_env <- data.frame(
  hli = sample_pops[[1]][, 1],
  cti = sample_pops[[1]][, 2]
)
env_diff <- site_pairwise_diff(site_env)
pair_data <- cbind(pair_data, env_diff)

# 2. Fit baseline MLPE models before optimizing resistance.
ibd_mod <- mlpe(
  response ~ distance + (1 | pair),
  data = pair_data,
  pairs = c("from", "to")
)
ibe_mod <- mlpe(
  response ~ distance + hli_diff + cti_diff + (1 | pair),
  data = pair_data,
  pairs = c("from", "to")
)

# 3. Prepare a conditional resistance optimization.
gdist_inputs <- gdist.prep(
  n.Pops = nrow(samples),
  samples = pts,
  response = pair_data$response,
  covariates = pair_data[, c("distance", "hli_diff", "cti_diff")],
  formula = response ~ distance + hli_diff + cti_diff,
  method = "commuteDistance"
)

# 4. Optimize resistance while holding IBD and IBE terms in the model.
ss_out <- SS_optim(
  GA.inputs = GA.prep(
    raster = {
      x <- resistance_surfaces
      if (inherits(x, "PackedSpatRaster")) x <- terra::unwrap(x)
      x
    },
    Results.dir = file.path(tempdir(), "ibe-ibr-results")
  ),
  gdist.inputs = gdist_inputs
)

How to read the output: The baseline models test whether geographic distance and endpoint environmental mismatch explain pairwise genetic structure before resistance is optimized. The conditional SS_optim() call asks whether the resistance distance term cd improves fit after those baseline terms are already included. A stable, positive cd coefficient and improved AICc support a conditional association between landscape permeability and the sampled genetic distances beyond geographic and environmental mismatch. They do not, by themselves, establish a causal landscape effect or out-of-sample prediction.

Summary of key functions

Function	Purpose
`mlpe_data()`	Build the pairwise data frame and shared-site metadata used by MLPE models
`site_pairwise_diff()`	Convert site-level environmental covariates into absolute pairwise differences
`site_pairwise_dist()`	Convert site-level covariates into multivariate pairwise environmental distances
`mlpe()`	Fit baseline MLPE models before resistance optimization
`gdist.prep()`	Prepare conditional resistance-distance inputs for optimization
`SS_optim()`	Optimize resistance surfaces while conditioning on IBD or IBE covariates
`MS_optim()`	Optimize multisurface resistance models under the same conditional framework

Joint Isolation by Environment and Isolation by Resistance