Resistance Modeling in Landscape Genetics

The central question in landscape genetics is deceptively simple: why do some pairs of populations exchange genes more readily than others? Geography matters (nearby populations tend to be more similar) but geography is rarely the whole story. The landscape between populations shapes movement, and movement shapes gene flow. A river, a road, a ridge line, or simply an expanse of unsuitable habitat can reduce effective connectivity far below what geographic distance alone would predict.

Quantifying that connectivity requires a resistance surface: a raster in which each cell value represents the cost of moving through that location. Low resistance values allow easy passage; high resistance values impede movement. By accumulating resistance across paths between locations, we obtain pairwise effective distances that can be compared against observed genetic differentiation.

The problem is that resistance values are not directly observable. Ecologists and geneticists have traditionally assigned them from expert opinion or tested a small discrete set of candidate values. ResistanceGA2 replaces that approach with stochastic optimization: a genetic algorithm searches the parameter space and identifies the resistance parameterization that best explains the observed pattern of pairwise genetic distances, as judged by the fit of a Maximum Likelihood Population Effects (MLPE) mixed model.

The Genetic Algorithm Approach

The genetic algorithm (GA) is a metaheuristic that mimics natural selection over a population of candidate solutions. In each generation:

1. Every candidate parameter set is used to transform the input raster into a resistance surface.
2. Pairwise effective distances are computed from that surface (via gdistance or Julia/Circuitscape).
3. An MLPE model is fit with those distances as the predictor of pairwise genetic distance.
4. The objective function (by default, log-likelihood) is returned as the “fitness” of that candidate.
5. High-fitness candidates are selected, recombined, and mutated to form the next generation.

This process repeats until the population converges or a stopping criterion is met (maximum iterations, or a run of generations with no improvement). The result is the parameter set (transformation shape, maximum resistance, and/or class-specific values) that maximizes the explanatory power of the resistance surface for the genetic data at hand.

The MLPE Framework

Standard regression models cannot be applied directly to pairwise genetic and landscape distance data because observations sharing a common sampling location are not independent (any pair involving population A is correlated with any other pair involving population A). The MLPE model (Clarke et al. 2002) handles this by including a random effect structured to absorb that dependency:

genetic_distance ~ resistance_distance + (1 | pair)

In the standalone mlpe() interface, the dyadic random effect is declared as a named term such as (1 | pair) and constructed automatically from the endpoint columns in the pairwise data. Inside the optimization workflow, the package still builds the legacy internal pop/grp terms behind the scenes so existing result objects and helper functions continue to work. Users do not need to build a ZZ matrix or call To.From.ID() directly when using the main optimization functions or the modern MLPE interface.

ResistanceGA2 provides a dyad-aware interface for fitting MLPE models directly:

• mlpe_data(response, ..., pairs, labels) constructs a pairwise data frame by combining a response vector with named covariate vectors (IBD, IBE, etc.) and from/to endpoint columns.
• mlpe(formula, data, pairs, ...) then fits the MLPE model using lme4-style syntax while building the dyadic random-effect covariance from those endpoint columns automatically.

The older mlpe_rga(), MLPE.lmm(), and To.From.ID() functions remain in the package for backward compatibility and internal optimization workflows, but they are no longer the recommended API for new standalone analyses.

For a worked example of the mlpe()/mlpe_data() interface in a joint IBD/IBE/IBR context, see the IBE_IBR_Framework vignette.

Choosing an Objective Function

Three objective functions are available via GA.prep(method = ...):

• "LL" (default): Maximize the log-likelihood of the MLPE model. This is the recommended starting point. LL is a direct measure of model fit and provides a clean, sensitive objective surface for the GA.
• "AIC": Minimize AIC. Penalizes model complexity, but the penalty depends on k.value (see below).
• "R2": Maximize marginal R² of the MLPE model. Intuitive, but R² tends to be a less sensitive and less reliable guide for the GA than LL, which is why log-likelihood is the default objective. This option is retained for backward compatibility but is not recommended as a primary search objective.

Regardless of which objective drives the GA, the output always includes AIC, AICc, R²_m, R²_c, and LL for all fitted models so that you can compare results on multiple criteria after the search.

Counting Parameters: k.value

AIC and AICc require knowing the number of free parameters k in the model. For a standard linear model, this is straightforward. For a pairwise MLPE model fitted to distances that were themselves produced by optimizing several underlying parameters, the answer is less clear.

ResistanceGA2 exposes four counting rules via k.value in GA.prep():

There is no universally correct choice. k.value = 2 is a practical starting point, but you should report which option you used and be consistent across all models in a comparison. Using different k.value settings for different models in the same study is not valid.

Transformation Families

For continuous surfaces, the GA searches a two-dimensional space: the shape parameter and the maximum resistance parameter of a chosen transformation family. The choice of transformation families to include in the search is controlled by select.trans in GA.prep().

The default behavior includes all Monomolecular families ("M") for continuous layers and assigns NA (no transformation, class-specific values) for categorical layers. Adding Ricker families ("A" for all, or "R" for Ricker only) increases search space and flexibility but also increases the risk of overfitting; Ricker curves can approximate almost any monotone shape when the shape parameter is extreme, at the cost of less interpretable parameters.

Recommendation: Start with select.trans = list("M") for continuous layers. Only include Ricker families if you have a biological hypothesis about a peaked or threshold-type resistance response. See the Surface_Transformations vignette for a detailed discussion.

The Dataset

This vignette uses the raster_orig surface and the sample_pops/Dc_list data for code demonstrations, and the precomputed results stored in inst/extdata/vignette-results/ for all displayed outputs. This keeps render time short while showing realistic results from a properly parameterized run.

# Larger dataset suitable for full optimization
pts  <- vect(sample_pops[[1]], type = "points")
gd   <- lower(Dc_list[[1]])

pts
#> class       : SpatVector
#> geometry    : points
#> dimensions  : 50, 0  (geometries, attributes)
#> extent      : 3.5, 20.5, 3.5, 20.5  (xmin, xmax, ymin, ymax)
#> coord. ref. :
raster_orig
#> class       : SpatRaster
#> size        : 25, 25, 2  (nrow, ncol, nlyr)
#> resolution  : 1, 1  (x, y)
#> extent      : 0, 25, 0, 25  (xmin, xmax, ymin, ymax)
#> coord. ref. : 
#> source(s)   : memory
#> names       : cat_orig, cont_orig
#> min values  :        1, -3.533382
#> max values  :        5,   4.03984

# Summary statistics of the genetic distance vector
summary(gd)
#>    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
#>  0.2417  0.4059  0.5142  0.5284  0.6186  0.9528

plot(raster_orig[["cont_orig"]], main = "Continuous resistance surface")
points(pts, pch = 19, col = "white", cex = 0.8)
points(pts, pch = 1,  col = "black", cex = 0.8)

Continuous resistance surface with sample population locations (white dots). This is the surface to be optimized.

Preparing GA Settings

ga_inputs <- GA.prep(
  raster      = raster_orig,
  Results.dir = file.path(tempdir(), "ResistanceGA2-tutorial"),
  method      = "LL",
  pop.size    = 25,
  maxiter     = 200,
  run         = 25,
  monitor     = FALSE,
  quiet       = TRUE
)

ga_inputs$layer.names
#> [1] "cat_orig"  "cont_orig"
ga_inputs$surface.type
#> [1] "cat"  "cont"
ga_inputs$select.trans
#> [[1]]
#> [1] NA
#> 
#> [[2]]
#> [1] 1 3 5 7 9

The select.trans component shows which transformation families will be searched for each layer. For continuous layers, the default is "M" (Monomolecular families). You can override this:

# Include all transformation families for a more exhaustive search
ga_inputs_full <- GA.prep(
  raster       = raster_orig,
  Results.dir  = file.path(tempdir(), "ResistanceGA2-tutorial-full"),
  select.trans = list("A"),   # "A" = all 8 families
  pop.size     = 50,
  maxiter      = 500,
  run          = 50,
  monitor      = FALSE,
  quiet        = TRUE
)

# Restrict to Monomolecular only (recommended default)
ga_inputs_mono <- GA.prep(
  raster       = raster_orig,
  Results.dir  = file.path(tempdir(), "ResistanceGA2-tutorial-mono"),
  select.trans = list("M"),   # "M" = Monomolecular families only
  pop.size     = 50,
  maxiter      = 500,
  run          = 50,
  monitor      = FALSE,
  quiet        = TRUE
)

Preparing gdistance Inputs

gdist_inputs <- gdist.prep(
  n.Pops   = nrow(sample_pops[[1]]),
  response = gd,
  samples  = pts,
  method   = "costDistance"
)

names(gdist_inputs)
#>  [1] "response"           "samples"            "covariates"        
#>  [4] "formula"            "transitionFunction" "directions"        
#>  [7] "ID"                 "ZZ"                 "keep"              
#> [10] "n.Pops"             "longlat"            "method"            
#> [13] "commute.approx"     "approx.factor"      "approx.scale"      
#> [16] "df"

The response argument is the vector of pairwise genetic distances obtained from lower() applied to a square distance matrix. Including it here links the MLPE model to the genetic data during optimization.

For commuteDistance workflows, you do not need to guess when the aggregate approximation is likely to be worthwhile. recommend_commute_settings() summarizes the benchmark guidance built into the package and returns a ready-to-pass recommendation:

commute_plan <- recommend_commute_settings(
  raster = raster_orig,
  n.Pops = nrow(sample_pops[[1]])
)

commute_plan[c("commute.approx", "approx.factor", "benchmark_band")]
#> $commute.approx
#> [1] "aggregate"
#> 
#> $approx.factor
#> [1] 2
#> 
#> $benchmark_band
#> [1] "intermediate"

gdist_commute_inputs <- gdist.prep(
  n.Pops         = nrow(sample_pops[[1]]),
  response       = gd,
  samples        = pts,
  method         = "commuteDistance",
  commute.approx = commute_plan$commute.approx,
  approx.factor  = commute_plan$approx.factor
)

When the helper recommends an aggregate approximation, use it for the iterative GA search and then rerun the final best surface once with exact commuteDistance() output for reporting.

Note: If you want to use Circuitscape instead of gdistance, replace gdist.prep() with jl.prep() and pass jl.inputs instead of gdist.inputs to all downstream functions. The interface is otherwise identical. See the Julia_Guide vignette.

Running SS_optim

SS_optim() fits one GA optimization per surface layer, plus optional null (intercept-only) and distance (IBD) reference models. Each surface is optimized independently.

ss_out <- SS_optim(
  gdist.inputs     = gdist_inputs,
  GA.inputs        = ga_inputs,
  dist_mod         = TRUE,    # include IBD model in comparison
  null_mod         = TRUE,    # include null model in comparison
  diagnostic_plots = TRUE     # generate residual and QQ plots
)

This is computationally expensive for large rasters or many sample locations. The code above is shown for reference; the results below are loaded from a precomputed file.

The optimization result is still a named list, but it now carries an S3 class with package methods for print(), summary(), and plot().

class(ss_out)
summary(ss_out)
plot(ss_out, type = "response")
plot(ss_out, type = "ranking")

Multi-Surface Optimization

Loading and Displaying Results

ss_rds <- system.file(
  "extdata/vignette-results/ss_optim_results.rds",
  package = "ResistanceGA2"
)

if (file.exists(ss_rds)) {
  ss_out <- readRDS(ss_rds)
} else {
  message("Precomputed SS results not found. Run data-raw/generate_optimization_results.R first.")
  ss_out <- NULL
}

knitr::kable(
  ss_out$AICc,
  digits  = 3,
  caption = "Single-surface model comparison. Lower AICc = better fit. LL = log-likelihood. R2m = marginal R-squared. k = number of parameters."
)

The AICc table compares all surfaces tested plus the null and distance reference models. Columns:

• AIC / AICc: Akaike Information Criteria (corrected for small samples). Lower = better.
• R2m: Marginal R², variance explained by fixed effects (resistance distance) alone.
• R2c: Conditional R², variance explained by fixed and random effects combined.
• LL: Log-likelihood. Higher = better fit.
• k: Number of parameters counted (governed by k.value in GA.prep()).

Recovered Transformation Parameters

if (!is.null(ss_out$ContinuousResults)) {
  knitr::kable(
    ss_out$ContinuousResults,
    digits  = 3,
    caption = "Optimized transformation parameters for continuous surfaces."
  )
}

For each continuous surface, the GA returns:

• Transformation: Which of the 8 transformation families was selected.
• shape: The curvature parameter. Controls how steeply the curve accelerates.
• max: The maximum resistance value. Controls the amplitude of the transformation.
• LL / AICc: Model fit metrics for this surface’s best parameter set.

For categorical surfaces, class-specific resistance values are returned in $CategoricalResults.

Interpreting SS_optim Output

The single-surface results tell you, for each landscape variable tested in isolation, how much explanatory power it contributes to pairwise genetic distances. A surface with a substantially lower AICc than the null model is worth investigating further. A surface that does not improve on the null, or that is outperformed by the simple geographic distance (IBD) model, may not be an important driver of genetic structure for this system.

Key interpretation points:

1. Compare against the null model first. If the best resistance surface barely improves on the null, either the landscape features tested are not driving connectivity or the genetic data lack sufficient signal.

2. Compare against the IBD model. If a resistance surface outperforms geographic distance, that surface captures structure beyond simple isolation by distance.

3. Look at R²_m in absolute terms. An improvement in AICc can be small in absolute fit terms. R²_m gives you a sense of how much genetic variance is actually accounted for.

4. Do not over-interpret single-surface results. In real systems, multiple features usually contribute simultaneously. Multi-surface optimization (below) tests whether combining surfaces improves on the best individual surface.

Warning: A surface that ranks first in single-surface analysis is not guaranteed to be the most important driver in a multi-surface analysis. Variables can have correlated effects, and their relative importance may shift when optimized jointly.

Running MS_optim

MS_optim() simultaneously optimizes the transformation parameters for all surface layers and estimates their relative contributions to a single composite resistance surface. The composite is a linear combination of functional transformations defining the maximum resistance and pattern of change:

composite_resistance = f(surface1_transformed) + f(surface2_transformed) + ...

where the transformations are also estimated by the GA.

ms_out <- MS_optim(
  gdist.inputs     = gdist_inputs,
  GA.inputs        = ga_inputs,
  diagnostic_plots = TRUE
)

Loading Multi-Surface Results

class(ms_out)
summary(ms_out)
plot(ms_out, type = "response")
plot(ms_out, type = "contribution")
plot(ms_out, type = "fit")

If you want to reconstruct the optimized composite surface after a multisurface run:

ms_rds <- system.file(
  "extdata/vignette-results/ms_optim_results.rds",
  package = "ResistanceGA2"
)

if (file.exists(ms_rds)) {
  ms_out <- readRDS(ms_rds)
} else {
  message("Precomputed MS results not found. Run data-raw/generate_optimization_results.R first.")
  ms_out <- NULL
}

knitr::kable(
  ms_out$AICc,
  digits  = 3,
  caption = "Multi-surface model comparison."
)

Percent Contribution

if (!is.null(ms_out$percent.contribution)) {
  # percent.contribution is a data frame with columns: surface, mean
  pct <- ms_out$percent.contribution
  pct$mean <- round(pct$mean * 100, 1)  # convert to percentage
  names(pct) <- c("Surface", "Contribution (%)")
  knitr::kable(
    pct,
    caption = "Percent contribution of each surface to the composite resistance model.",
    row.names = FALSE
  )
}

Percent contribution reflects the relative weight of each surface in the optimized composite. A surface with 90% contribution dominates the composite; a surface with 10% contribution has little influence. However, note that these percentages are sensitive to scaling and the correlation structure among surfaces. They should be interpreted as a rough guide to relative importance, not as precise effect sizes.

Reconstructing the Composite Surface

After a multi-surface optimization, you can reconstruct the optimized composite surface for visualization or downstream analysis:

composite <- Combine_Surfaces(
  PARM         = ms_out$GA.summary@solution[1, ],
  gdist.inputs = gdist_inputs,
  GA.inputs    = ga_inputs,
  rescale      = TRUE
)

plot(composite, main = "Optimized composite resistance surface")

Combine_Surfaces() applies the optimal transformations to each layer and combines them with the optimal weights. Setting rescale = TRUE scales the composite to the range [1, 10].

When to Use all_comb

all_comb() is useful when you have multiple candidate surfaces and want to systematically evaluate all possible combinations up to a specified size. Rather than choosing between single-surface and multi-surface analyses, you test every subset: surface A alone, surface B alone, A+B together, A+C together, and so on.

This is more thorough than running SS_optim() and MS_optim() separately, because it evaluates partial combinations that neither function covers.

Setting Up all_comb

all_comb() requires Results.dir = "all.comb" in the GA.prep() call. A separate results.dir argument specifies where all-combinations outputs are written.

# GA settings must use Results.dir = "all.comb"
ga_allcomb <- GA.prep(
  raster      = raster_orig,
  Results.dir = "all.comb",
  pop.size    = 25,
  maxiter     = 200,
  run         = 25,
  monitor     = FALSE,
  quiet       = TRUE
)

results_dir <- file.path(tempdir(), "ResistanceGA2-all-comb")

combo_out <- all_comb(
  gdist.inputs     = gdist_inputs,
  GA.inputs        = ga_allcomb,
  results.dir      = results_dir,
  max.combination  = 2,   # test 1- and 2-surface combinations
  iters            = 500, # bootstrap iterations
  replicate        = 1    # GA replicates per combination
)

Loading all_comb Results

class(combo_out)
summary(combo_out)
plot(combo_out)

ac_rds <- system.file(
  "extdata/vignette-results/all_comb_results.rds",
  package = "ResistanceGA2"
)

if (file.exists(ac_rds)) {
  combo_out <- readRDS(ac_rds)
} else {
  message("Precomputed all_comb results not found.")
  combo_out <- NULL
}

if (!is.null(combo_out$summary.table)) {
  knitr::kable(
    combo_out$summary.table,
    digits  = 3,
    caption = "All-combinations summary table. Rows are surface combinations ranked by AICc."
  )
}

Why Bootstrap?

Single-run optimization results depend on the stochastic nature of the GA and the composition of the specific sample. Bootstrap resampling (Resist.boot()) asks: how consistently does each surface combination rank highly across many subsamples of the pairwise data? A surface that ranks first in 90% of bootstrap replicates provides much stronger support than one that ranks first only 55% of the time.

Running Resist.boot

Resist.boot() resamples the pairwise observations to assess how consistently each surface ranks highly across subsampled datasets. When you use all_comb(), bootstrap is run automatically and the results are stored in combo_out$boot.results; you do not need to call Resist.boot() separately.

If you want to run bootstrap directly after SS_optim() and MS_optim() (without all_comb()), you must build the distance matrix list and k vector from those outputs:

# Collect distance matrices and k values from SS and MS runs
mat_list <- c(ss_out$cd, list(ms_out$cd))
k_vec    <- c(ss_out$k[, 2], nrow(ms_out$k))  # k per model
names(mat_list) <- c(names(ss_out$cd), "multi_surface")

# Build square genetic distance matrix (required by Resist.boot)
n_pops <- gdist_inputs$n.Pops
gd_mat <- matrix(0, n_pops, n_pops)
gd_mat[lower.tri(gd_mat)] <- gdist_inputs$response

boot_out <- Resist.boot(
  mod.names    = names(mat_list),
  dist.mat     = mat_list,
  n.parameters = k_vec,
  sample.prop  = 0.75,
  iters        = 1000,
  obs          = n_pops,
  genetic.mat  = gd_mat
)

Note: dist.mat must be a named list of square pairwise distance matrices (one per model). n.parameters is a numeric vector with the parameter count for each model (the k column from $k in SS/MS output). obs is the number of sampling locations as an integer. genetic.mat is the square genetic distance matrix without row or column names.

Loading Bootstrap Results

boot_rds <- system.file(
  "extdata/vignette-results/boot_results.rds",
  package = "ResistanceGA2"
)

if (file.exists(boot_rds)) {
  boot_out <- readRDS(boot_rds)
} else {
  message("Precomputed bootstrap results not found.")
  boot_out <- NULL
}

knitr::kable(
  boot_out,
  digits  = 3,
  caption = "Bootstrap model support summary from Resist.boot()."
)

Interpreting Bootstrap Results

The Resist.boot() output is a tibble with one row per surface combination. Key columns:

Percent.top is usually the most informative column. A surface with Percent.top > 50% is the modal top model across bootstrap replicates. A surface with avg.weight substantially higher than all alternatives has consistently stronger support.

Tip: If no single surface achieves Percent.top > 50%, the data may not have sufficient power to distinguish among the candidate surfaces, or the top surfaces may all be plausible models (model uncertainty). In that case, multi-model inference (model averaging using AIC weights) is more appropriate than selecting a single best model.