--- title: "Getting Started with ResistanceGA2" author: "Bill Peterman" date: "`r format(Sys.Date(), '%d %B %Y')`" output: rmarkdown::html_vignette: toc: true toc_depth: 3 number_sections: true vignette: | %\VignetteIndexEntry{Getting Started with ResistanceGA2} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r setup, include=FALSE} knitr::opts_chunk$set( collapse = TRUE, comment = "#>", fig.width = 7, fig.height = 5 ) ``` ```{r libraries, message=FALSE, warning=FALSE} library(ResistanceGA2) library(terra) surfaces <- resistance_surfaces if (inherits(surfaces, "PackedSpatRaster")) { surfaces <- terra::unwrap(surfaces) } ``` --- # 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. ## Spatial classes: terra only All spatial inputs must now use the `terra` package. There is no fallback to `raster` or `sp` classes. | Old | New | |-----|-----| | `RasterLayer`, `RasterStack` | `SpatRaster`; use `terra::rast()` | | `SpatialPoints`, `SpatialPointsDataFrame` | `SpatVector`; use `terra::vect()` | | `raster::stack()`, `dropLayer()` | `terra::rast()`, `terra::subset()` | Bundled package datasets (`resistance_surfaces`, `raster_orig`, etc.) are stored as wrapped `SpatRaster` objects and unwrapped automatically when the package loads, so no manual conversion is needed. ## Removed features | Removed | Replacement | |---------|-------------| | `CS.prep()`, `Run_CS()` | `jl.prep()`, `Run_CS.jl()` for Julia/Circuitscape workflows | | `JuliaCall` dependency | `JuliaConnectoR` | | Standalone Circuitscape `.exe` | Julia + Circuitscape.jl | | `Grid.Search()` | No equivalent; use GA with restricted families | | `k.smooth` parameter | Removed; use `select.trans` in `GA.prep()` | | `SS_optim.scale()`, `MS_optim.scale()` | Merged into `SS_optim()`, `MS_optim()` | ## New features in ResistanceGA2 - **Conditional IBR optimization**: `gdist.prep()` and `jl.prep()` accept `covariates` and `formula` arguments to condition resistance search on IBD or IBE fixed effects. - **`mlpe()` and `mlpe_data()`**: New dyad-aware interface for fitting MLPE models directly, with `lme4`-style formula syntax. - **`site_pairwise_diff()` / `site_pairwise_dist()`**: Helper functions for constructing pairwise IBE covariates from site-level environmental measurements. - **S3 methods on results**: `SS_optim()`, `MS_optim()`, and `all_comb()` return objects with `print()`, `summary()`, and `plot()` methods. - **`Run_CS.ini()`**: Low-level Circuitscape `.ini` wrapper for advanced scenarios (source/ground mode, mask files, included-pairs files). - **Idaho demo data**: `rga2_demo` and `rga2_demo_covariates` for the `IBE_IBR_Framework` vignette. - **Parallel improvements**: `SpatRaster` objects are serialized via `terra::wrap()` before being dispatched to parallel workers, fixing Windows-specific crashes. --- # 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. ## Isolation by Distance, Resistance, and Environment Three hypotheses are commonly evaluated in landscape genetics: **Isolation by Distance (IBD):** Genetic differentiation increases with geographic distance because dispersal is limited. This is the null expectation under random dispersal across a homogeneous landscape. **Isolation by Resistance (IBR):** Genetic differentiation is better explained by *effective* distance (the cost of moving through a heterogeneous landscape) than by straight-line distance. Different habitat types impede movement to different degrees, so a resistance surface is parameterized and optimized. **Isolation by Environment (IBE):** Genetic differentiation reflects differences in local environmental conditions (temperature, elevation, soil type) that reduce the probability of successful migration between ecologically dissimilar sites, even if those sites are geographically close. These three hypotheses are not mutually exclusive. In many systems, genetic structure reflects a combination of all three. `ResistanceGA2` supports workflows that test these hypotheses individually and in combination, including the `mlpe()` formula interface and `mlpe_data()` dyad constructor that let you include environmental covariates alongside resistance distances in a single MLPE model. See the `IBE_IBR_Framework` vignette for a full worked example of concurrent IBD/IBE/IBR testing. ## The Conceptual Workflow ``` 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 ``` The genetic algorithm (GA) is the engine that drives the search. At each generation, a population of candidate parameter sets is evaluated: each candidate defines a specific transformation of the input raster, a distance matrix is calculated, an MLPE model is fit, and the objective function (default: log-likelihood) is returned. The GA then selects, recombines, and mutates the best candidates to produce the next generation, repeating until convergence. Inside the optimization workflow, `ResistanceGA2` still prepares the legacy `pop`/`grp` random-effect terms needed for the package internals. For standalone analysis, however, the recommended interface is the dyad-aware `mlpe_data()`/`mlpe()` workflow shown in the `IBE_IBR_Framework` vignette. --- # Installation Install the development version from GitHub: ```{r install, eval=FALSE} 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: ```{r install-no-vignettes, eval=FALSE} 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 ```{r resistance-surfaces} resistance_surfaces ``` `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). ```{r plot-surfaces, fig.cap="The three bundled resistance surface layers. Categorical (left) has discrete class values. Continuous (center) varies smoothly. Feature (right) represents a linear landscape element."} plot(surfaces, main = c("Categorical", "Continuous", "Feature")) ``` ## Sample Locations ```{r sample-data} # Small dataset: 25 populations pts_small <- vect(samples[, 2:3], type = "points") pts_small # Larger dataset for full optimization runs pts_large <- vect(sample_pops[[1]], type = "points") pts_large ``` `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. ```{r plot-locations, fig.cap="Sample locations overlaid on the continuous resistance surface. The 25-point small dataset (filled circles) is useful for quick examples."} 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) ``` ## The `raster_orig` Dataset ```{r raster-orig} raster_orig ``` `raster_orig` contains the `cont_orig` layer, a continuous surface intended for use with the larger `sample_pops` dataset in full optimization demonstrations. Use this layer when you want to run `SS_optim()` or `MS_optim()` with the larger sample sets. ## Precomputed Pairwise Outputs Expensive computations (Circuitscape runs, full GA optimizations) are precomputed and stored as package data so that vignettes render quickly: ```{r precomputed-outputs} # Precomputed pairwise distance matrices names(vignette_pairwise_outputs) cd_matrix <- vignette_pairwise_outputs$gdistance_commute dim(cd_matrix) cd_matrix[1:4, 1:4] cs_vector <- vignette_pairwise_outputs$circuitscape_matrix length(cs_vector) head(cs_vector) ``` ```{r current-map} # Precomputed Circuitscape current map (stored as a PackedSpatRaster) vignette_current_map ``` --- # 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. ## Step 1: Prepare Genetic Algorithm Settings with `GA.prep()` `GA.prep()` configures the genetic algorithm and records metadata about each resistance surface layer (type, transformation families to consider, parameter bounds). ```{r ga-prep, message=FALSE, warning=FALSE} 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 ga_inputs$surface.type ``` Key arguments to know: | Argument | Default | Purpose | |---|---|---| | `raster` | required | `SpatRaster` with one layer per surface hypothesis | | `Results.dir` | required | Directory where tables and plots will be saved | | `method` | `"LL"` | Objective function: `"LL"`, `"AIC"`, or `"R2"` | | `select.trans` | auto | List of transformation families per layer (`"M"`, `"R"`, `"A"`, or `NA` for categorical) | | `pop.size` | auto | GA population size (individuals per generation); computed from the number of parameters when not set | | `maxiter` | 1000 | Maximum GA generations | | `run` | 25 | Generations of no improvement before stopping | | `parallel` | FALSE | Use parallel processing | | `k.value` | 2 | How to count model complexity for AIC/AICc | > **Tip:** For continuous surfaces, restrict `select.trans` to Monomolecular families (`"M"`) unless you have a specific biological reason to expect a unimodal response. Ricker families are more flexible but increase overfitting risk. For example: `select.trans = list("M", "M")` for a two-layer continuous raster. ## Step 2: Prepare Distance Engine Inputs ### Using gdistance `gdist.prep()` sets up everything needed to compute pairwise least-cost or commute-time distances inside R. No external software is required. ```{r gdist-prep, warning=FALSE} gdist_inputs <- gdist.prep( n.Pops = nrow(samples), samples = pts_small, method = "costDistance" ) names(gdist_inputs) ``` The `method` argument selects between: - **`"costDistance"`**: Dijkstra's least-cost path algorithm. Faster; analogous to the least-cost corridor concept. - **`"commuteDistance"`**: Random-walk commute time. Slower but conceptually closer to Circuitscape's effective resistance. If you want a quick benchmark-based recommendation for commute settings before you call `gdist.prep()`, use `recommend_commute_settings()`: ```{r commute-recommendation} commute_plan <- recommend_commute_settings( raster = terra::subset(surfaces, "continuous"), n.Pops = nrow(samples) ) commute_plan[c("commute.approx", "approx.factor", "benchmark_band")] commute_plan$rationale ``` You can then pass the recommendation directly into `gdist.prep()`: ```{r gdist-commute-prep, eval=FALSE} 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 ) ``` ### Using Julia / Circuitscape For Circuitscape-based distances, use `jl.prep()`. This requires a working Julia installation with `Circuitscape.jl`. See the `Julia_Guide` vignette for full setup instructions. ```{r jl-prep, eval=FALSE} 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 ) ``` > **Note:** You only need to prepare *one* distance engine (gdistance or Julia, not both). For iterative optimization over many surfaces, `gdistance` is usually faster. For final Circuitscape-matched outputs or current maps, use Julia. ## Step 3: Run Optimization With both `ga_inputs` and a distance engine object ready, call the appropriate optimization function: - `SS_optim()`: optimize each surface individually (single-surface) - `MS_optim()`: optimize all surfaces simultaneously as a composite (multi-surface) - `all_comb()`: systematically test all combinations of surfaces up to a specified size ```{r ss-optim-eval, eval=FALSE} # 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 ) ``` The returned object has an S3 class with built-in `print()`, `summary()`, and `plot()` methods: ```{r ss-s3, eval=FALSE} class(ss_out) print(ss_out) summary(ss_out) plot(ss_out, type = "response") plot(ss_out, type = "ranking") ``` ## Multisurface optimization ```{r ms-workflow, eval=FALSE} ms_out <- MS_optim( gdist.inputs = gdist_inputs, GA.inputs = ga_inputs, diagnostic_plots = FALSE ) ``` Multisurface results support the same style of inspection. ```{r ms-s3, eval=FALSE} class(ms_out) summary(ms_out) plot(ms_out, type = "response") plot(ms_out, type = "contribution") plot(ms_out, type = "fit") ``` ## All-combinations analysis For a full worked example with precomputed results, see the `ResistanceGA` vignette. ## Step 4: Interpret and Compare Results Optimization functions return structured lists. The most important components are: - `$AICc`: model comparison table with AIC, AICc, R^2^~m~, R^2^~c~, LL, and k values - `$ContinuousResults` / `$CategoricalResults`: optimized parameter estimates - `$cd`: list of pairwise distance matrices for the winning surfaces - `$MLPE.list`: fitted MLPE model objects for downstream inference ```{r results-preview, eval=FALSE} # View model comparison table knitr::kable(ss_out$AICc, digits = 2, caption = "Single-surface model comparison") # View optimized transformation parameters ss_out$ContinuousResults ``` `all_comb()` results also have S3 methods. ```{r all-comb-s3, eval=FALSE} 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. ```{r plot-trans, fig.cap="A Monomolecular transformation applied to the continuous resistance surface. Low values in the original raster receive low resistance; the curve saturates at high values."} Plot.trans( PARM = c(2.5, 100), Resistance = terra::subset(surfaces, "continuous"), transformation = "Monomolecular", marginal.plot = FALSE ) ``` You can also apply a transformation directly to obtain a new `SpatRaster`: ```{r resistance-tran} transformed <- Resistance.tran( transformation = "Monomolecular", shape = 2.5, max = 100, r = terra::subset(surfaces, "continuous") ) transformed ``` 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 ```{r getting-started-quick-reference, eval=FALSE} 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 | Function | Purpose | |---|---| | `GA.prep()` | Prepare raster layers, transformation bounds, output directories, and GA control settings | | `gdist.prep()` | Prepare sample points and pairwise data for R-based distance calculations | | `jl.prep()` | Prepare sample points and pairwise data for Julia/Circuitscape distance calculations | | `SS_optim()` | Optimize each candidate resistance surface separately | | `MS_optim()` | Optimize a composite resistance surface from multiple layers | | `all_comb()` | Compare candidate surface combinations and optional replicate runs | | `Resist.boot()` | Summarize model support across pseudo-bootstrap subsamples | # 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 | Vignette | Description | |---|---| | `ResistanceGA` | Full workflow tutorial with single-surface, multi-surface, all-combinations, and bootstrap analyses using precomputed results | | `Surface_Transformations` | Detailed explanation of all eight transformation families, when to use each, and how to constrain the GA search | | `IBE_IBR_Framework` | Concurrent testing of IBD, IBE, and IBR using `mlpe()`/`mlpe_data()` with pairwise environmental covariates and the Idaho demo dataset | | `Julia_Guide` | Installing Julia, configuring `JuliaConnectoR`, running Circuitscape, and troubleshooting | Access any vignette with: ```{r vignette-access, eval=FALSE} 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. 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. 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. 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.