--- title: "ResistanceGA2: Full Optimization Workflow" 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{ResistanceGA2: Full Optimization Workflow} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r setup, include=FALSE} knitr::opts_chunk$set( collapse = TRUE, comment = "#>", fig.width = 7, fig.height = 5, message = FALSE, warning = FALSE ) ``` ```{r libraries} library(ResistanceGA2) library(terra) library(ggplot2) ``` --- > **Migrating from ResistanceGA?** See the *Key Changes from ResistanceGA* section in the `GettingStarted` vignette for a complete list of breaking changes, removed functions, and new features. --- # Background ## 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. --- # Thoughts on Optimization Before launching a GA search, it is worth thinking carefully about the choices that determine what the search is actually optimizing. ## 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^2^ of the MLPE model. Intuitive, but R^2^ 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^2^~m~, R^2^~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()`: | k.value | Rule | |---|---| | 1 | Parameters from the fitted MLPE model only (intercept + variance components) | | 2 | MLPE parameters + intercept (default) | | 3 | MLPE parameters + number of surface layers + intercept | | 4 | Number of surfaces + intercept (ignores MLPE variance components) | 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. --- # Example Data and Preparation ## 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. ```{r example-data} # Larger dataset suitable for full optimization pts <- vect(sample_pops[[1]], type = "points") gd <- lower(Dc_list[[1]]) pts raster_orig # Summary statistics of the genetic distance vector summary(gd) ``` ```{r plot-data, fig.cap="Continuous resistance surface with sample population locations (white dots). This is the surface to be optimized."} 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) ``` ## Preparing GA Settings ```{r ga-prep} 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 ga_inputs$surface.type ga_inputs$select.trans ``` 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: ```{r ga-prep-custom, eval=FALSE} # 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 ```{r gdist-prep} gdist_inputs <- gdist.prep( n.Pops = nrow(sample_pops[[1]]), response = gd, samples = pts, method = "costDistance" ) names(gdist_inputs) ``` 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: ```{r gdist-commute-recommend} commute_plan <- recommend_commute_settings( raster = raster_orig, n.Pops = nrow(sample_pops[[1]]) ) commute_plan[c("commute.approx", "approx.factor", "benchmark_band")] ``` ```{r gdist-commute-recommend-use, eval=FALSE} 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. --- # Single-Surface Optimization ## 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. ```{r ss-optim-code, eval=FALSE} 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()`. ```{r ss-optim-s3, eval=FALSE} class(ss_out) summary(ss_out) plot(ss_out, type = "response") plot(ss_out, type = "ranking") ``` ## Multi-Surface Optimization ## Loading and Displaying Results ```{r ss-results-load} 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 } ``` ```{r ss-aic-table, eval=exists("ss_out") && !is.null(ss_out)} 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^2^, variance explained by fixed effects (resistance distance) alone. - **R2c**: Conditional R^2^, 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 ```{r ss-cont-results, eval=exists("ss_out") && !is.null(ss_out)} 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^2^~m~ in absolute terms.** An improvement in AICc can be small in absolute fit terms. R^2^~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. --- # Multi-Surface Optimization ## 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. ```{r ms-optim-code, eval=FALSE} ms_out <- MS_optim( gdist.inputs = gdist_inputs, GA.inputs = ga_inputs, diagnostic_plots = TRUE ) ``` ## Loading Multi-Surface Results ```{r ms-optim-s3, eval=FALSE} 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: ```{r ms-results-load} 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 } ``` ```{r ms-aic-table, eval=exists("ms_out") && !is.null(ms_out)} knitr::kable( ms_out$AICc, digits = 3, caption = "Multi-surface model comparison." ) ``` ## Percent Contribution ```{r ms-contribution, eval=exists("ms_out") && !is.null(ms_out)} 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: ```{r combine-surfaces-code, eval=FALSE} 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]. --- # All-Combinations Analysis ## 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. ```{r all-comb-code, eval=FALSE} # 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 ```{r all-comb-s3, eval=FALSE} class(combo_out) summary(combo_out) plot(combo_out) ``` ```{r allcomb-results-load} 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 } ``` ```{r allcomb-table, eval=exists("combo_out") && !is.null(combo_out)} 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." ) } ``` --- # Bootstrap Analysis ## 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: ```{r boot-code, eval=FALSE} # 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 ```{r boot-results-load} 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 } ``` ```{r boot-table, eval=exists("boot_out") && !is.null(boot_out)} 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: | Column | Interpretation | |---|---| | `surface` | Surface combination label | | `avg.AICc` | Mean AICc across bootstrap replicates | | `avg.weight` | Mean AIC weight (0 to 1); higher values indicate more support | | `avg.rank` | Mean rank across replicates (1 = best) | | `avg.R2m` | Mean marginal R^2^ | | `avg.LL` | Mean log-likelihood | | `Percent.top` | % of replicates where this surface ranked first | | `n` | Number of replicates included | | `k` | Parameter count | **`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. --- # Comparing Single-Surface, Multi-Surface, and Null Models After running `SS_optim()`, `MS_optim()`, and optionally baseline IBD/null models, you can combine the AICc tables to compare the full candidate set. The precomputed results below illustrate this workflow. ```{r combined-comparison, eval=exists("ss_out") && !is.null(ss_out) && exists("ms_out") && !is.null(ms_out)} # Combine SS and MS AICc tables (both have the same column structure) if (!is.null(ss_out) && !is.null(ms_out)) { all_results <- rbind( ss_out$AICc, ms_out$AICc[!rownames(ms_out$AICc) %in% rownames(ss_out$AICc), ] ) # Sort by AICc all_results <- all_results[order(all_results$AICc), ] knitr::kable( all_results, digits = 3, caption = "Combined model comparison: single-surface and multi-surface results ranked by AICc." ) } ``` When interpreting the combined table: - **AICc differences (ΔAICc):** Models within 2 AICc units of the best model have substantial support. Models 4–7 units above the best have considerably less support. Models > 10 units above can generally be ruled out. - **Null model as a baseline:** If the best resistance model barely improves on the null, reconsider whether the landscape variables tested are meaningful drivers of connectivity. - **Multi-surface vs. single-surface:** If the multi-surface model substantially outperforms the best single-surface model, the composite landscape story is more informative than any single driver. --- # Diagnostic Plots When `diagnostic_plots = TRUE`, `SS_optim()` and `MS_optim()` automatically generate and save: 1. **Observed vs. predicted plot:** Genetic distance on the y-axis, resistance distance on the x-axis. A tight scatter around the regression line indicates good fit. 2. **Residual plots:** Standardized residuals vs. fitted values. Patterns in residuals (curvature, funnel shapes) suggest model misspecification. 3. **QQ plot:** Quantile–quantile plot of residuals against a normal distribution. Deviations suggest non-normality. 4. **GA convergence plot:** Best fitness value per generation. A well-converged run shows a plateau at the end. You can also call `Diagnostic.Plots()` directly on any MLPE model object: ```{r diag-plots-code, eval=FALSE} # Access fitted MLPE model from SS results mlpe_mod <- ss_out$MLPE.list[[1]] # Generate diagnostic plots Diagnostic.Plots( resistance.mat = ss_out$cd[[1]], genetic.mat = as.matrix(dist(sample_pops[[1]])), MLPE.list = ss_out$MLPE.list, plot.dir = file.path(tempdir(), "diag-plots") ) ``` --- # Quick reference: full optimization workflow ```{r resistancega-quick-reference, eval=FALSE} library(ResistanceGA2) library(terra) # 1. Prepare surfaces, sample points, and pairwise response. 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. Configure the genetic algorithm and distance engine. ga_inputs <- GA.prep( raster = surfaces, Results.dir = file.path(tempdir(), "rga2-full-workflow") ) gdist_inputs <- gdist.prep( n.Pops = nrow(samples), samples = pts, response = response, method = "commuteDistance" ) # 3. Fit single-surface and multisurface optimizations. ss_out <- SS_optim(GA.inputs = ga_inputs, gdist.inputs = gdist_inputs) ms_out <- MS_optim(GA.inputs = ga_inputs, gdist.inputs = gdist_inputs) # 4. Compare candidate models and inspect diagnostics. summary(ss_out) # single-surface model support and fitted transformations summary(ms_out) # multisurface support and percent contribution plot(ss_out, type = "response") # top single-surface response plot(ss_out, type = "ranking") # single-surface candidate support plot(ms_out, type = "response") # multisurface response plot(ms_out, type = "contribution") # component contributions ``` **How to read the output:** `SS_optim()` asks which individual surface has the strongest support. `MS_optim()` asks whether a weighted composite of surfaces explains the response better than any one layer alone. Use AICc, delta AICc, model weights, and consistency across replicate runs together. A model with a low AICc but unstable transformation parameters should be treated as a signal to run more replicates or simplify the candidate transformation set. # Summary of key functions | Function | Purpose | |---|---| | `GA.prep()` | Define raster inputs, transformation families, parameter bounds, and GA control settings | | `gdist.prep()` | Prepare R-based cost or commute-distance calculations | | `recommend_commute_settings()` | Compare exact and aggregate commute-distance settings for a dataset | | `SS_optim()` | Optimize candidate resistance surfaces one at a time | | `MS_optim()` | Optimize a weighted composite resistance surface | | `all_comb()` | Evaluate surface combinations and optional replicate runs | | `Resist.boot()` | Assess model support across pseudo-bootstrap subsamples | | `Diagnostic.Plots()` | Save diagnostic plots for fitted MLPE models | # See also - `vignette("GettingStarted", package = "ResistanceGA2")`, the shortest path through the package workflow. - `vignette("Surface_Transformations", package = "ResistanceGA2")`, interpreting the transformation families searched by the GA. - `vignette("IBE_IBR_Framework", package = "ResistanceGA2")`, conditioning resistance optimization on IBD and IBE covariates. - `vignette("Julia_Guide", package = "ResistanceGA2")`, switching the distance engine from `gdistance` to Julia/Circuitscape. --- # Tips for Running Full Analyses **Set `parallel = TRUE`** in `GA.prep()` for any run with `pop.size > 20`. Parallel evaluation of the GA population provides substantial speedups on multi-core machines. **Start with small `pop.size` and `maxiter`** (e.g., 25 and 100) for workflow testing. Only increase to production settings (50–100 and 500–1000) for final analyses. **Run multiple replicates.** The `replicate` argument in `all_comb()` runs the GA multiple times per surface combination. If the best-fit parameters are consistent across replicates, the optimization has converged reliably. **Use `run` to control early stopping.** The `run` argument stops the GA after `run` consecutive generations with no improvement greater than `run.diff`. Setting `run = 50` and `maxiter = 1000` allows long searches but stops early when the solution has plateaued. **Save your Results.dir.** All tables and plots are written to `Results.dir` automatically. Keep this directory so you can review optimization trajectories and compare across parameter settings without re-running. **Consider resolution carefully.** Very high-resolution rasters dramatically increase computation time for each distance calculation. Consider aggregating to a coarser resolution for initial optimization, then refining at higher resolution for final analyses. > **Warning:** Do not use `maxiter = 10` or `run = 5` for publication-quality results. These settings are only appropriate for workflow testing. An inadequately converged GA can return suboptimal parameter estimates that appear definitive but are actually artifacts of premature stopping. --- # References 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. 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. McRae, B. H. 2006. Isolation by resistance. *Evolution* 60:1551-1561. 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. Row, J. R., G. Blouin-Demers, and S. C. Lougheed. 2010. Habitat distribution influences dispersal and fine-scale genetic population structure of eastern foxsnakes (*Mintonius gloydi*) across a fragmented landscape. *Molecular Ecology* 19:5157–5171.