Key Changes from ResistanceGA

ResistanceGA2 is a major revision. If you are migrating scripts from the original ResistanceGA package, the summary below covers every breaking change.

What Is ResistanceGA2?

Landscape genetics asks how the physical environment shapes the movement of individuals; and by extension, how genes flow across the landscape. Two populations might be far apart in straight-line geographic distance, but if a highway, a mountain range, or a patch of inhospitable habitat lies between them, their genetic exchange may be far lower than geography alone would predict.

ResistanceGA2 is an R package for optimizing the resistance values assigned to those landscape features. Rather than deciding in advance that “forest = low resistance, pavement = high resistance,” the package uses a genetic algorithm to search the space of possible parameter values and find the resistance surface that best explains the observed pattern of genetic differentiation among sampled populations.

  Landscape rasters               Sample locations
  (categorical / continuous)      (SpatVector points)
         |                               |
         v                               v
  +----------------------------------------------+
  |  Resistance surface parameterization          |
  |  (8 transformation families or class values)  |
  +----------------------------------------------+
         |
         v
  +----------------------------------------------+
  |  Distance engine                              |
  |  gdistance (costDistance / commuteDistance)   |
  |  OR Julia / Circuitscape                      |
  +----------------------------------------------+
         |  pairwise effective distances
         v
  +----------------------------------------------+
  |  MLPE model                                   |
  |  genetic distance ~ resistance distance       |
  |           + dyadic random intercept           |
  +----------------------------------------------+
         |  log-likelihood (objective)
         v
  +----------------------------------------------+
  |  Genetic Algorithm                            |
  |  searches transformation parameter space      |
  |  maximizes LL / minimizes AIC                 |
  +----------------------------------------------+
         |
         v
  Optimized resistance surface + model comparison table

Installation

Install the development version from GitHub:

install.packages("devtools")
devtools::install_github("wpeterman/ResistanceGA2", build_vignettes = TRUE)

If Pandoc is not available on your machine, install without vignettes first and render the vignette sources later:

devtools::install_github("wpeterman/ResistanceGA2", build_vignettes = FALSE)

The core dependencies are terra, gdistance, lme4, GA, and ggplot2. For Julia/Circuitscape support, you also need JuliaConnectoR and a working Julia installation; see the Julia_Guide vignette for setup instructions.

Package Data Overview

ResistanceGA2 ships several datasets that you can use immediately for exploration and testing.

Resistance Surfaces

resistance_surfaces
#> class       : SpatRaster
#> size        : 50, 50, 3  (nrow, ncol, nlyr)
#> resolution  : 0.5, 0.5  (x, y)
#> extent      : 0, 25, 0, 25  (xmin, xmax, ymin, ymax)
#> coord. ref. : 
#> source(s)   : memory
#> names       : categorical, continuous, feature
#> min values  :           0,   6.514823,       0
#> max values  :           2,  13.376758,       1

resistance_surfaces is a three-layer SpatRaster:

• categorical: a categorically coded surface with class values representing distinct habitat types (e.g., land cover classes). The GA assigns a separate resistance value to each class.
• continuous: a continuously varying surface (e.g., canopy cover, elevation). The GA fits a parameterized transformation curve mapping raw values to resistance values.
• feature: a binary or near-binary surface representing a linear landscape feature (e.g., a road or river).

plot(surfaces, main = c("Categorical", "Continuous", "Feature"))

The three bundled resistance surface layers. Categorical (left) has discrete class values. Continuous (center) varies smoothly. Feature (right) represents a linear landscape element.

Sample Locations

# Small dataset: 25 populations
pts_small <- vect(samples[, 2:3], type = "points")
pts_small
#> class       : SpatVector
#> geometry    : points
#> dimensions  : 25, 0  (geometries, attributes)
#> extent      : 6.25, 19.25, 6.25, 19.25  (xmin, xmax, ymin, ymax)
#> coord. ref. :

# Larger dataset for full optimization runs
pts_large <- vect(sample_pops[[1]], type = "points")
pts_large
#> class       : SpatVector
#> geometry    : points
#> dimensions  : 50, 0  (geometries, attributes)
#> extent      : 3.5, 20.5, 3.5, 20.5  (xmin, xmax, ymin, ymax)
#> coord. ref. :

samples contains 25 sampling locations as a matrix with x/y coordinates. sample_pops is a list of coordinate matrices representing multiple replicate populations; Dc_list contains corresponding genetic distance matrices.

plot(terra::subset(surfaces, "continuous"), main = "Sample locations")
points(pts_small, pch = 19, col = "white", cex = 0.9)
points(pts_small, pch = 1,  col = "black", cex = 0.9)

Sample locations overlaid on the continuous resistance surface. The 25-point small dataset (filled circles) is useful for quick examples.

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

# Precomputed pairwise distance matrices
names(vignette_pairwise_outputs)
#> [1] "gdistance_commute"   "circuitscape_matrix"

cd_matrix <- vignette_pairwise_outputs$gdistance_commute
dim(cd_matrix)
#> [1] 25 25
cd_matrix[1:4, 1:4]
#>           1         2         3         4
#> 1   0.00000  79.71312 138.35098 200.72560
#> 2  79.71312   0.00000  58.63786 121.01249
#> 3 138.35098  58.63786   0.00000  62.37462
#> 4 200.72560 121.01249  62.37462   0.00000

cs_vector <- vignette_pairwise_outputs$circuitscape_matrix
length(cs_vector)
#> [1] 300
head(cs_vector)
#> [1] 5.538005 6.486677 7.276222 7.802246 5.089549 5.828889

# Precomputed Circuitscape current map (stored as a PackedSpatRaster)
vignette_current_map
#> class       : SpatRaster
#> size        : 50, 50, 1  (nrow, ncol, nlyr)
#> resolution  : 0.5, 0.5  (x, y)
#> extent      : 0, 25, 0, 25  (xmin, xmax, ymin, ymax)
#> coord. ref. : lon/lat WGS 84 (CRS84) (OGC:CRS84)
#> source(s)   : memory
#> varname     : continuous_cum_curmap
#> name        : continuous
#> min value   :   0.028609
#> max value   :   31.56336

The Four-Step Workflow

Every analysis in ResistanceGA2 follows the same four-step structure, regardless of whether you are testing one surface or many, using gdistance or Julia, or running a single replicate or a full bootstrap.

ga_inputs <- GA.prep(
  raster      = terra::subset(surfaces, c("categorical", "continuous")),
  Results.dir = file.path(tempdir(), "ResistanceGA2-getting-started"),
  monitor     = FALSE,
  quiet       = TRUE
)

# Layer names and types detected automatically
ga_inputs$layer.names
#> [1] "categorical" "continuous"
ga_inputs$surface.type
#> [1] "cat"  "cont"

gdist_inputs <- gdist.prep(
  n.Pops  = nrow(samples),
  samples = pts_small,
  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"

commute_plan <- recommend_commute_settings(
  raster = terra::subset(surfaces, "continuous"),
  n.Pops = nrow(samples)
)

commute_plan[c("commute.approx", "approx.factor", "benchmark_band")]
#> $commute.approx
#> [1] "aggregate"
#> 
#> $approx.factor
#> [1] 2
#> 
#> $benchmark_band
#> [1] "intermediate"
commute_plan$rationale
#> [1] "Raster size 50 x 50 (2,500 cells) with 25 populations maps to the 'intermediate' benchmark band. Recommendation: use commute.approx = 'aggregate' with approx.factor = 2 during optimization, then rerun the final best surface exactly once without approximation."

commute_inputs <- gdist.prep(
  n.Pops         = nrow(samples),
  response       = lower(as.matrix(dist(samples[, 2:3]))),
  samples        = pts_small,
  method         = "commuteDistance",
  commute.approx = commute_plan$commute.approx,
  approx.factor  = commute_plan$approx.factor
)

jl_inputs <- jl.prep(
  n.Pops        = nrow(samples),
  response      = lower(as.matrix(dist(samples[, 2:3]))),
  CS_Point.File = pts_small,
  JULIA_HOME    = Sys.getenv("JULIA_BINDIR"),
  run_test      = TRUE,
  silent        = TRUE
)

# Single-surface optimization (one run per layer)
ss_out <- SS_optim(
  gdist.inputs    = gdist_inputs,
  GA.inputs       = ga_inputs,
  dist_mod        = TRUE,   # also fit an IBD distance model
  null_mod        = TRUE,   # also fit a null (intercept-only) model
  diagnostic_plots = FALSE
)

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

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

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

# View model comparison table
knitr::kable(ss_out$AICc, digits = 2, caption = "Single-surface model comparison")

# View optimized transformation parameters
ss_out$ContinuousResults

class(all_out)
summary(all_out)
plot(all_out)

Transformations and Inspection

Visualizing Transformations

Before running a full optimization, it is useful to explore how the transformation families reshape a resistance surface. Plot.trans() shows the relationship between raw raster values and transformed resistance values for a given parameter set.

Plot.trans(
  PARM           = c(2.5, 100),
  Resistance     = terra::subset(surfaces, "continuous"),
  transformation = "Monomolecular",
  marginal.plot  = FALSE
)

A Monomolecular transformation applied to the continuous resistance surface. Low values in the original raster receive low resistance; the curve saturates at high values.

You can also apply a transformation directly to obtain a new SpatRaster:

transformed <- Resistance.tran(
  transformation = "Monomolecular",
  shape          = 2.5,
  max            = 100,
  r              = terra::subset(surfaces, "continuous")
)
transformed
#> class       : SpatRaster
#> size        : 50, 50, 1  (nrow, ncol, nlyr)
#> resolution  : 0.5, 0.5  (x, y)
#> extent      : 0, 25, 0, 25  (xmin, xmax, ymin, ymax)
#> coord. ref. : 
#> source(s)   : memory
#> name        : continuous
#> min value   :          1
#> max value   :  99.168436

The two parameters that define every continuous transformation are:

• shape: Controls the curvature of the response. Small values produce gradual curves; large values produce steep thresholds.
• max: The maximum resistance value assigned to the most resistant cells. The minimum is always scaled to 1.

For a complete description of all eight transformation families and guidance on choosing among them, see the Surface_Transformations vignette.

These objects still keep their original list structure, so existing code using $AICc, $AICc.tab, $cd, $k, and related components continues to work. The new S3 classes add a cleaner interactive layer on top of that structure.

For larger or slower analyses, it is often worth running multiple replicates and checking for consistency in the top-supported surfaces and parameter estimates.

Quick reference: first complete workflow

library(ResistanceGA2)
library(terra)

# 1. Load bundled surfaces and sample locations.
surfaces <- resistance_surfaces
if (inherits(surfaces, "PackedSpatRaster")) {
  surfaces <- terra::unwrap(surfaces)
}
pts <- terra::vect(samples[, 2:3], type = "points")
response <- lower(as.matrix(dist(samples[, 2:3])))

# 2. Prepare GA settings and the distance engine.
ga_inputs <- GA.prep(
  raster = surfaces,
  Results.dir = file.path(tempdir(), "rga2-results")
)

gdist_inputs <- gdist.prep(
  n.Pops = nrow(samples),
  samples = pts,
  response = response,
  method = "commuteDistance"
)

# 3. Run a small single-surface optimization for workflow testing.
ss_out <- SS_optim(
  GA.inputs = ga_inputs,
  gdist.inputs = gdist_inputs
)

# 4. Interpret and compare fitted surfaces.
summary(ss_out)                         # ranking and top-model coefficients
plot(ss_out, type = "response")         # fitted response relationship
plot(ss_out, type = "ranking")          # candidate model support

How to read the workflow: GA.prep() defines what the genetic algorithm is allowed to search. gdist.prep() or jl.prep() defines how resistance distances are calculated. SS_optim() fits one surface at a time and returns a ranked model object. Start with a small run to verify paths, data, and output interpretation, then increase pop.size, maxiter, and the number of replicates for production analyses.

Summary of key functions

Key Terms Glossary

Resistance surface: A raster in which each cell value represents the cost of moving through that cell. High values mean high resistance (impeded movement).

Effective distance: A pairwise measure of the total movement cost between two locations, computed by accumulating resistance across the landscape path.

MLPE (Maximum Likelihood Population Effects): A mixed-effects model for pairwise distance data that accounts for the non-independence of observations sharing a common sample location (Clarke et al. 2002).

Genetic algorithm (GA): A stochastic optimization method that maintains a population of candidate solutions, applies selection, crossover, and mutation operators, and converges toward high-fitness parameter regions.

Transformation family: One of eight parameterized curve shapes (Monomolecular, Reverse Monomolecular, Inverse Monomolecular, Inverse-Reverse Monomolecular, and four Ricker equivalents) that map raw raster values to resistance values.

Log-likelihood (LL): The default objective function for optimization. The GA maximizes LL of the MLPE model, which corresponds to finding the resistance surface that best explains pairwise genetic distances.

AICc: Corrected Akaike Information Criterion. Used for model comparison after optimization. Lower values indicate better-fitting models penalized for complexity.

k.value: A user-selected rule for counting the number of parameters in the AIC/AICc calculation. Different defensible choices exist; the default (k.value = 2) should be reported in publications.

IBD / IBR / IBE: Isolation by Distance / Resistance / Environment. Three competing (and complementary) hypotheses about the drivers of spatial genetic structure.

See also

Access any vignette with:

vignette("ResistanceGA",         package = "ResistanceGA2")
vignette("Surface_Transformations", package = "ResistanceGA2")
vignette("IBE_IBR_Framework",    package = "ResistanceGA2")
vignette("Julia_Guide",          package = "ResistanceGA2")

Tips for New Users

Start small. Use resistance_surfaces and samples for your first runs. Keep pop.size at 25 and maxiter at 50-100 to get fast results and verify your workflow before scaling up.

Use Monomolecular first. The Monomolecular family covers most ecologically realistic resistance–covariate relationships. It is harder to overfit and easier to interpret than Ricker families. Only include Ricker families if you have a biological hypothesis that predicts a peaked or threshold-type response.

Run replicates. Genetic algorithms are stochastic. If the same surface wins consistently across replicate runs with similar parameter estimates, that is a more trustworthy result than a single run.

Check diagnostic plots. Set diagnostic_plots = TRUE in SS_optim() or MS_optim(). Residual plots and the observed-vs-predicted scatter are among the fastest ways to detect poor fits.

Report k.value. Because there is no universally accepted way to count parameters in pairwise MLPE models, you should always report which k.value option you used and justify the choice in your methods.

Save your Results.dir outputs. The functions write tables and plots to Results.dir. These are useful for comparing runs and documenting the optimization trajectory.

Warning: Do not interpret optimization as hypothesis testing. Finding that a particular resistance surface has the lowest AICc does not mean the alternative surfaces are false. It means the data, as analyzed, are more consistent with that parameterization. Replication, cross-validation, and biological interpretation are essential complements to the optimization output.

References

Clarke, R. T., P. Rothery, and A. F. Raybould. 2002. Confidence limits for regression relationships between distance matrices: Estimating gene flow with distance. Journal of Agricultural, Biological, and Environmental Statistics 7:361-372. https://doi.org/10.1198/108571102320

Peterman, W. E. 2018. ResistanceGA: An R package for the optimization of resistance surfaces using genetic algorithms. Methods in Ecology and Evolution 9:1638-1647. https://doi.org/10.1111/2041-210X.12984

Peterman, W. E., G. M. Connette, R. D. Semlitsch, and L. S. Eggert. 2014. Ecological resistance surfaces predict fine-scale genetic differentiation in a terrestrial woodland salamander. Molecular Ecology 23:2402–2413.

McRae, B. H. 2006. Isolation by resistance. Evolution 60:1551-1561. https://doi.org/10.1111/j.0014-3820.2006.tb00500.x

Beninde, J., J. Wittische, and A. C. Frantz. 2024. Quantifying uncertainty in inferences of landscape genetic resistance due to choice of individual-based genetic distance metric. Molecular Ecology Resources 24:e13831.