| Title: | Optimize Resistance Surfaces Using Genetic Algorithms |
|---|---|
| Description: | Optimizes resistance surfaces using genetic algorithms. Continuous and categorical surfaces can be optimized individually or jointly to create composite resistance surfaces. Resistance distances can be calculated using 'Circuitscape' through 'Julia' or with the 'gdistance' package. The framework extends Peterman (2018) <doi:10.1111/2041-210X.12984>. |
| Authors: | Bill Peterman [aut, cre] |
| Maintainer: | Bill Peterman <[email protected]> |
| License: | GPL (>= 3) |
| Version: | 1.2.0 |
| Built: | 2026-11-12 16:47:08 UTC |
| Source: | https://github.com/wpeterman/ResistanceGA2 |
A wrapper function to run both SS_optim and
MS_optim to optimize all combinations of
resistance surfaces with the Genetic Algorithms package GA. Following
optimization, Resist.boot is run to conduct a
bootstrap analysis.
all_comb( gdist.inputs = NULL, jl.inputs = NULL, GA.inputs, results.dir, max.combination = 4, iters = 1000, replicate = 1, sample.prop = 0.75, nlm = FALSE, dist_mod = TRUE, null_mod = TRUE, n.cores = NULL, ... )all_comb( gdist.inputs = NULL, jl.inputs = NULL, GA.inputs, results.dir, max.combination = 4, iters = 1000, replicate = 1, sample.prop = 0.75, nlm = FALSE, dist_mod = TRUE, null_mod = TRUE, n.cores = NULL, ... )
gdist.inputs |
Object created from running
|
jl.inputs |
Object created from running
|
GA.inputs |
Object created from running
|
results.dir |
Directory in which to write the analysis. A nonempty directory is used only after interactive confirmation; the user may instead choose a new time-stamped subdirectory. |
max.combination |
Maximum number of surfaces per candidate model, or a two-element vector giving the minimum and maximum. A minimum greater than one skips single-surface optimization. Values cannot exceed the number of raster layers. |
iters |
Positive integer giving the number of population-subsampling iterations passed to [Resist.boot()]. Default is 1000. |
replicate |
Positive integer giving the number of independent GA runs for each candidate set. Default is 1. |
sample.prop |
Proportion of populations used to fit each bootstrap model. Default is 0.75; see [Resist.boot()]. |
nlm |
Reserved for a possible final nonlinear optimization step. This
feature is not implemented and must remain |
dist_mod |
Logical. Include a geographic-distance baseline in the
candidate set? Default is |
null_mod |
Logical. Include an intercept-only baseline? Default is
|
n.cores |
Number of parallel workers used by [Resist.boot()].
|
... |
Reserved for backward compatibility. Additional arguments are not currently used. |
Candidate-set comparisons are valid only for models fitted to the same response and retained population pairs. GA replication evaluates numerical stability of the optimizer. Population subsampling evaluates how candidate rankings change across sampled population sets. Neither is a substitute for inspecting fitted coefficients, diagnostic plots, or the scientific plausibility of each surface.
An object with class resga_all_comb. For one replicate, the
object contains summary.table, boot.results, model metadata
(all.k), resistance-distance matrices (all.cd), the genetic
distance matrix, and the underlying ss.results and
ms.results. With multiple replicates, these elements are nested in
named rep_1, rep_2, and subsequent lists. CSV summaries and
serialized optimization outputs are also written beneath
results.dir.
Bill Peterman <[email protected]>
[SS_optim()], [MS_optim()], [Resist.boot()], [summary.resga_all_comb()]
## Not run: pts <- terra::vect(sample_pops[[1]], type = "points") gd <- lower(Dc_list[[1]]) gdist.inputs <- gdist.prep( n.Pops = nrow(sample_pops[[1]]), response = gd, samples = pts ) GA.inputs <- GA.prep( raster = raster_orig, Results.dir = "all.comb", quiet = TRUE ) results_dir <- file.path(tempdir(), "ResistanceGA2-all-comb") combo.out <- all_comb( gdist.inputs = gdist.inputs, GA.inputs = GA.inputs, results.dir = results_dir, max.combination = 2, iters = 10, replicate = 1, n.cores = 2 ) ## End(Not run)## Not run: pts <- terra::vect(sample_pops[[1]], type = "points") gd <- lower(Dc_list[[1]]) gdist.inputs <- gdist.prep( n.Pops = nrow(sample_pops[[1]]), response = gd, samples = pts ) GA.inputs <- GA.prep( raster = raster_orig, Results.dir = "all.comb", quiet = TRUE ) results_dir <- file.path(tempdir(), "ResistanceGA2-all-comb") combo.out <- all_comb( gdist.inputs = gdist.inputs, GA.inputs = GA.inputs, results.dir = results_dir, max.combination = 2, iters = 10, replicate = 1, n.cores = 2 ) ## End(Not run)
Combines multiple resistance surfaces by transforming each layer with specified parameters and summing the results.
Combine_Surfaces( PARM, gdist.inputs = NULL, jl.inputs = NULL, GA.inputs, out = NULL, File.name = paste(GA.inputs$parm.type$name, collapse = "."), rescale = TRUE, p.contribution = FALSE )Combine_Surfaces( PARM, gdist.inputs = NULL, jl.inputs = NULL, GA.inputs, out = NULL, File.name = paste(GA.inputs$parm.type$name, collapse = "."), rescale = TRUE, p.contribution = FALSE )
PARM |
Parameters to transform continuous or categorical resistance
surfaces. Provide a vector with parameters in the order of resistance
surfaces as output by |
gdist.inputs |
Object from |
jl.inputs |
Object from |
GA.inputs |
Object from |
out |
Directory to write the combined |
File.name |
Name for the output |
rescale |
Logical. If |
p.contribution |
Logical. If |
For continuous surfaces, PARM should specify, per layer, three
values: (1) transformation code (1–9), (2) shape, and (3) maximum value.
Transformation codes:
| 1 = "Inverse-Reverse Monomolecular" | |
| 2 = "Inverse-Reverse Ricker" | |
| 3 = "Monomolecular" | |
| 4 = "Ricker" | |
| 5 = "Reverse Monomolecular" | |
| 6 = "Reverse Ricker" | |
| 7 = "Inverse Monomolecular" | |
| 8 = "Inverse Ricker" | |
| 9 = "Distance" |
A SpatRaster of the combined resistance surface (sum of all
transformed/reclassified layers), or a list when
p.contribution = TRUE.
Bill Peterman <[email protected]>
## Not run: # After running MS_optim(), use the best-fit PARM vector: combined <- Combine_Surfaces( PARM = ms.out$PARM, gdist.inputs = gdist.inputs, GA.inputs = GA.inputs ) ## End(Not run)## Not run: # After running MS_optim(), use the best-fit PARM vector: combined <- Combine_Surfaces( PARM = ms.out$PARM, gdist.inputs = gdist.inputs, GA.inputs = GA.inputs ) ## End(Not run)
Sample file to be used with examples in the vignette
data(Dc_list)data(Dc_list)
A list of length 3
Each matrix depicts the pairwise genetic distance (measured as chord distance) between sample locations
Dc_cat. Pairwise chord distance matrix simulated across the categorical resistance surface
Dc_cont. Pairwise chord distance matrix simulated across the continuous resistance surface
Dc_multi. Pairwise chord distance matrix simulated across the multivariate resistance surface
Fits a mixed-effects model with pairwise effects (MLPE) to an optimized resistance-distance predictor and saves a four-panel diagnostic figure.
Diagnostic.Plots( resistance.mat, genetic.dist, XLAB = "Estimated resistance", YLAB = "Genetic distance", plot.dir, type = c("categorical", "continuous"), name = NULL, ID = NULL, ZZ = NULL )Diagnostic.Plots( resistance.mat, genetic.dist, XLAB = "Estimated resistance", YLAB = "Genetic distance", plot.dir, type = c("categorical", "continuous"), name = NULL, ID = NULL, ZZ = NULL )
resistance.mat |
A square resistance-distance matrix, a
lower-triangle distance vector, a |
genetic.dist |
Numeric vector of pairwise genetic distances in the same
lower-triangle order as |
XLAB |
Character string used for the resistance-distance axis label. |
YLAB |
Character string used for the response axis label. |
plot.dir |
Directory in which to write the TIFF figure. The directory is created when necessary. |
type |
Surface type, either |
name |
Output-file stem. A stem is derived from a Circuitscape file name
when possible; otherwise |
ID |
Optional pair-identifier data frame from [gdist.prep()] or [jl.prep()]. When omitted, identifiers are constructed from the inferred population count. |
ZZ |
Optional sparse MLPE design matrix from [gdist.prep()] or [jl.prep()]. This argument is retained for compatibility with legacy workflows. Current pair-identifier metadata is used when available. |
The response panel includes an ordinary least-squares line as a visual trend guide only. Inference should come from the fitted MLPE model, which accounts for the nonindependence of pairwise observations. Strong curvature, changing residual spread, or marked departures in the quantile plot indicate that the fitted mean or variance structure may be inadequate. These plots assess in-sample fit and do not demonstrate transferability to new populations or landscapes.
Invisibly returns a list with the fitted model, the model
data, normalized output file path, and surface type.
The saved TIFF contains observed response against resistance distance,
residuals against marginal fitted values, a residual histogram, and a
residual normal quantile plot.
Bill Peterman <[email protected]>
[plot.resga_ss_optim()], [plot.resga_ms_optim()], [MLPE.lmm()]
## Not run: set.seed(4) xy <- matrix(stats::runif(12), ncol = 2) genetic <- lower(as.matrix(stats::dist(xy))) resistance <- lower(as.matrix(stats::dist(cbind(xy[, 1], xy[, 2]^2)))) diagnostics <- Diagnostic.Plots( resistance.mat = resistance, genetic.dist = genetic, plot.dir = tempdir(), type = "continuous", name = "example" ) diagnostics$model ## End(Not run)## Not run: set.seed(4) xy <- matrix(stats::runif(12), ncol = 2) genetic <- lower(as.matrix(stats::dist(xy))) resistance <- lower(as.matrix(stats::dist(cbind(xy[, 1], xy[, 2]^2)))) diagnostics <- Diagnostic.Plots( resistance.mat = resistance, genetic.dist = genetic, plot.dir = tempdir(), type = "continuous", name = "example" ) diagnostics$model ## End(Not run)
Validates a multilayer resistance raster, classifies its layers as continuous or categorical, defines parameter bounds and starting values, and stores the settings used by [SS_optim()], [MS_optim()], and [all_comb()].
GA.prep( raster, Results.dir = NULL, min.cat = NULL, max.cat = 1000, max.cont = 1000, shape.min = NULL, shape.max = NULL, cont.shape = NULL, select.trans = NULL, cat.levels = 15, method = "LL", k.value = 2, pop.mult = 15, percent.elite = 5, type = "real-valued", pcrossover = 0.85, pmutation = 0.125, maxiter = 1000, run = NULL, keepBest = TRUE, population = gaControl(type)$population, selection = gaControl(type)$selection, crossover = "gareal_blxCrossover", mutation = gaControl(type)$mutation, pop.size = NULL, parallel = FALSE, gaisl = FALSE, island.pop = 20, numIslands = NULL, migrationRate = NULL, migrationInterval = NULL, optim = FALSE, optim.method = "L-BFGS-B", poptim = 0, pressel = 1, control = list(fnscale = -1, maxit = 100), hessian = FALSE, opt.digits = NULL, seed = NULL, monitor = TRUE, quiet = FALSE )GA.prep( raster, Results.dir = NULL, min.cat = NULL, max.cat = 1000, max.cont = 1000, shape.min = NULL, shape.max = NULL, cont.shape = NULL, select.trans = NULL, cat.levels = 15, method = "LL", k.value = 2, pop.mult = 15, percent.elite = 5, type = "real-valued", pcrossover = 0.85, pmutation = 0.125, maxiter = 1000, run = NULL, keepBest = TRUE, population = gaControl(type)$population, selection = gaControl(type)$selection, crossover = "gareal_blxCrossover", mutation = gaControl(type)$mutation, pop.size = NULL, parallel = FALSE, gaisl = FALSE, island.pop = 20, numIslands = NULL, migrationRate = NULL, migrationInterval = NULL, optim = FALSE, optim.method = "L-BFGS-B", poptim = 0, pressel = 1, control = list(fnscale = -1, maxit = 100), hessian = FALSE, opt.digits = NULL, seed = NULL, monitor = TRUE, quiet = FALSE )
raster |
A |
Results.dir |
Directory to export optimization results. If the directory
does not exist, it will be created automatically. If it exists and is not
empty, a message will note that existing results may be overwritten. It is
critical that there are NO SPACES in the directory path. If using the
|
min.cat |
The minimum value to be assessed during optimization of categorical resistance surfaces (Default = 1 / max.cat) |
max.cat |
The maximum value to be assessed during optimization of categorical resistance surfaces (Default = 1000) |
max.cont |
The maximum value to be assessed during optimization of continuous resistance surfaces (Default = 1000) |
shape.min |
The minimum value for the shape parameter used for transforming resistance surfaces. If unspecified, used 0.5 |
shape.max |
The maximum value for the shape parameter used for transforming resistance surfaces. If unspecified, used 14.5 |
cont.shape |
A vector of hypothesized relationships that each continuous resistance surface will have in relation to the genetic distance response (Default = NULL; see details) |
select.trans |
Option to specify which transformations are applied to continuous surfaces. Must be provided as a list. "A" = All, "M" = Monomolecular only, "R" = Ricker only. Default = "M"; see Details. |
cat.levels |
Number of unique levels to permit in categorical surface (Default = 15). See Details |
method |
Objective function to be optimized. Select "AIC", "R2", or "LL" to optimize resistance surfaces based on AIC, variance explained (R2), or log-likelihood. (Default = "LL") |
k.value |
Specification of how k, the number of parameters in the mixed effects model, is determined. Specify 1, 2, 3, or 4 (Default = 2; see details). 1 –> k = 2; 2 –> k = number of parameters optimized plus the intercept; 3 –> k = the number of parameters optimized plus the intercept and the number of layers optimized; 4 –> k = the number of layers optimized plus the intercept |
pop.mult |
Value will be multiplied with number of parameters in surface to determine 'popSize' in GA. By default this is set to 15. |
percent.elite |
An integer percent used to determine the number of best fitness individuals to survive at each generation ('elitism' in GA). By default the top 5% individuals will survive at each iteration. |
type |
Default is "real-valued" |
pcrossover |
Probability of crossover. Default = 0.85 |
pmutation |
Probability of mutation. Default = 0.125 |
maxiter |
Maximum number of iterations to run before the GA search is halted. If using standard |
run |
Number of consecutive generations or epochs without any improvement in objective function before the GA is stopped. If using standard |
keepBest |
A logical argument specifying if best solutions at each iteration should be saved (Default = TRUE) |
population |
Default is "gareal_Population" from GA |
selection |
Default is "gareal_lsSelection" from GA |
crossover |
Default = "gareal_blxCrossover". This crossover method greatly improved optimization during preliminary testing |
mutation |
Default is "gareal_raMutation" from GA |
pop.size |
Number of individuals to create each generation. If |
parallel |
Controls parallel fitness evaluation in the GA package. Use
|
gaisl |
Should the genetic algorithm use the islands parallel optimization? (Default = FALSE) |
island.pop |
The number of individuals to populate each island. (Default = 20) |
numIslands |
If |
migrationRate |
If |
migrationInterval |
If |
optim |
A logical defaulting to |
optim.method |
The method to be used among those available in |
poptim |
A value in the range [0,1] specifying the probability of performing a local search at each iteration of GA (default 0.0). Only change if your optimization is relatively fast. |
pressel |
A value in the range [0,1] specifying the pressure selection (default 1.00). The local search is started from a random solution selected with probability proportional to fitness. High values of pressel tend to select the solutions with the largest fitness, whereas low values of pressel assign quasi-uniform probabilities to any solution. |
control |
A list of control parameters. See 'Details' section in |
hessian |
Logical. Should a numerically differentiated Hessian matrix be returned? This will allow for the calculation of standard errors on parameter estimates (not yet implemented). Default = FALSE |
opt.digits |
The number of significant digits that the objective function will be assessed at. By default, no rounding occurs. |
seed |
Integer random number seed to replicate |
monitor |
Default = TRUE, which prints the average and best fitness values at each iteration. |
quiet |
Logical. If TRUE, the objective function and step run time will not be printed to the screen after each step. Only |
Provide the resistance surfaces directly as a SpatRaster. If
multiple surfaces are to be optimized, include them as layers in that
object. File-path and directory-based raster inputs are no longer accepted.
The Default for k.value is 2, which sets k equal to the number of parameters optimized, plus 1 for the intercept term. Prior to version 3.0-0, k.value could not be specified by the user and followed setting 2, such that k was equal to the number of parameters optimized plus the intercept term.
cont.shape can take values of "Increase", "Decrease", or "Peaked". If you believe a resistance surface is related to your response in a particular way, specifying this here may decrease the time to optimization. cont.shape is used to generate an initial set of parameter values to test during optimization. If specified, a greater proportion of the starting values will include your believed relationship. If unspecified (the Default), a completely random set of starting values will be generated.
If it is desired that only certain transformations be assessed for
continuous surfaces, then this can be specified using
select.trans. By default, only monomolecular transformations
will be assessed for continuous surfaces unless otherwise specified.
This default is also the recommended starting point for most analyses:
monomolecular families are usually easier to interpret and generally
reduce the risk of overfitting relative to Ricker families. Ricker
transformations should usually be explored only when there is a strong
biological or ecological justification for a unimodal or quadratic-type
relationship. Specific transformations can be specified by providing a
vector of values (e.g., c(1,3,5)), with values corresponding to the
equation numbers as detailed in Resistance.tran.
If multiple rasters are to be optimized from the same SpatRaster
object, then a list of transformations must be provided in the order that
the raster surfaces will be assessed. For example:select.trans = list("M", "A", "R", c(5,6))
will result in surface one only being optimized with Monomolecular transformations, surface two with all transformations, surface three with only Ricker transformations, and surface four with Reverse Ricker and Reverse Monomolecular only. If a categorical surface is among the rasters to be optimized, it is necessary to specify NA to accommodate this.
cat.levels defaults to 15. This means that when a raster surface has <= 15 unique levels, it will be treated as a categorical surface in the analysis. This value can be increased, but optimization of surfaces with many levels may take more time. Additionally, depending upon the prevalence and configuration of categorical features and spatial sample locations, some levels are likely to be poorly estimated. This may be evident if estimated resistance values vary substantially between runs of ResistanceGA.
Setting gaisl = TRUE has the potential greatly reduce the optimization run time, potentially with greater accuracy. This is a distributed multiple-population GA, where the population is partitioned into several subpopulations and assigned to separated islands. Independent GAs are executed in each island, and only occasionally sparse exchanges of individuals are performed among the islands.
It is recommended to first run GA optimization with the default settings.
A named list used as GA.inputs. It contains the validated
resistance raster (Resistance.stack and raster), layer names
and inferred surface types, parameter indices and bounds, starting
suggestions, output directories, the objective and parameter-count rules,
random seed, and all standard or island-GA control settings. When
Results.dir is "all.comb" or "all_comb", directory
fields are left NULL and an inputs list is retained so that
[all_comb()] can create replicate-specific settings.
Bill Peterman <[email protected]>
[SS_optim()], [MS_optim()], [all_comb()]
ga.inputs <- GA.prep( raster = raster_orig, Results.dir = file.path(tempdir(), "ResistanceGA2-example"), quiet = TRUE, monitor = FALSE ) ga.inputs$layer.namesga.inputs <- GA.prep( raster = raster_orig, Results.dir = file.path(tempdir(), "ResistanceGA2-example"), quiet = TRUE, monitor = FALSE ) ga.inputs$layer.names
gdistance
Creates a necessary input for optimizing resistance surfaces based on
pairwise cost distances, implemented using the gdistance library.
gdist.prep( n.Pops, response = NULL, samples, covariates = NULL, formula = NULL, transitionFunction = function(x) 1/mean(x), directions = 8, longlat = FALSE, method = "commuteDistance", min.max_dist = NULL, keep = NULL, commute.approx = c("none", "aggregate"), approx.factor = 4L, approx.scale = TRUE )gdist.prep( n.Pops, response = NULL, samples, covariates = NULL, formula = NULL, transitionFunction = function(x) 1/mean(x), directions = 8, longlat = FALSE, method = "commuteDistance", min.max_dist = NULL, keep = NULL, commute.approx = c("none", "aggregate"), approx.factor = 4L, approx.scale = TRUE )
n.Pops |
The number of populations that are being assessed. |
response |
Vector of pairwise genetic distances (lower half of pairwise matrix). |
samples |
A |
covariates |
Data frame of additional covariates to include in the MLPE model during optimization. |
formula |
If covariates are included in the model, specify the R formula
for the fixed effects portion of the MLPE model (e.g.,
|
transitionFunction |
Function to calculate the |
directions |
Directions in which cells are connected (4, 8, 16, or other). Default = 8. |
longlat |
Logical. If |
method |
Pairwise distance method: |
min.max_dist |
Optional two-element vector |
keep |
Optional vector of length equal to the number of pairwise
observations, with 1 to keep an observation and 0 to drop it. Use
|
commute.approx |
Approximation backend for iterative
|
approx.factor |
Integer aggregation factor used when
|
approx.scale |
Logical. If |
Sample locations must be supplied as a point SpatVector. When
Run_gdistance is called, the resistance surface is converted
internally to the raster class required by gdistance.
Use recommend_commute_settings if you want a quick
benchmark-based suggestion for commute.approx and
approx.factor before calling gdist.prep().
Use graph_to_keep if you want to restrict optimization to a
graph-pruned subset of population pairs.
A covariate formula is expanded to include the resistance-distance predictor and MLPE pair effect. Covariates therefore adjust the resistance association; they do not become part of the resistance surface itself.
A named list used as gdist.inputs. It contains the response,
point locations, optional covariates and expanded MLPE formula, transition
and correction settings, pair identifiers (ID), MLPE design matrix
(ZZ), retained-pair indicator (keep), population count,
distance method, commute-approximation settings, and model data
(df). The object can be passed directly to [Run_gdistance()],
[SS_optim()], or [MS_optim()].
Bill Peterman <[email protected]>
[Run_gdistance()], [GA.prep()], [recommend_commute_settings()]
pts <- terra::vect(sample_pops[[1]], type = "points") commute_plan <- recommend_commute_settings( raster = terra::subset(resistance_surfaces, "continuous"), n.Pops = nrow(sample_pops[[1]]) ) gdist.inputs <- gdist.prep( n.Pops = nrow(sample_pops[[1]]), response = lower(Dc_list[[1]]), samples = pts, method = "commuteDistance", commute.approx = commute_plan$commute.approx, approx.factor = commute_plan$approx.factor ) head(gdist.inputs$df)pts <- terra::vect(sample_pops[[1]], type = "points") commute_plan <- recommend_commute_settings( raster = terra::subset(resistance_surfaces, "continuous"), n.Pops = nrow(sample_pops[[1]]) ) gdist.inputs <- gdist.prep( n.Pops = nrow(sample_pops[[1]]), response = lower(Dc_list[[1]]), samples = pts, method = "commuteDistance", commute.approx = commute_plan$commute.approx, approx.factor = commute_plan$approx.factor ) head(gdist.inputs$df)
Translate a pruned population graph, adjacency matrix, or edge list into the
integer keep vector used by gdist.prep and
jl.prep to retain only a subset of pairwise observations.
graph_to_keep(graph, sampled_pops = NULL, pop_n = NULL, labels = NULL)graph_to_keep(graph, sampled_pops = NULL, pop_n = NULL, labels = NULL)
graph |
Graph-like object describing which population pairs should be retained. Supported inputs are:
|
sampled_pops |
Integer number of sampled populations. Required when it
cannot be inferred from |
pop_n |
Optional integer vector giving the number of sampled
individuals per population. When supplied, the returned keep vector is
expanded to the individual-pair scale used by
|
labels |
Optional character vector giving population labels in the same
order as the genetic distance matrix and sample coordinates supplied to
|
ResistanceGA2 assumes pairwise distances are symmetric, so directed
graph inputs are treated as undirected. For square matrices without
dimnames, row order is assumed to already match the order of populations in
your pairwise genetic distance data.
This helper is intended to bridge graph-based pruning workflows from tools
such as popgraph or graph4lg into ResistanceGA2 without
requiring those packages as hard dependencies.
Integer vector of 1s and 0s aligned to lower /
To.From.ID ordering. Pass this to gdist.prep(keep = )
or jl.prep(pairs_to_include = ).
adj <- matrix( c(0, 1, 0, 1, 0, 1, 0, 1, 0), nrow = 3, byrow = TRUE, dimnames = list(c("A", "B", "C"), c("A", "B", "C")) ) keep <- graph_to_keep(adj, labels = c("A", "B", "C")) keep graph_to_keep(adj, labels = c("A", "B", "C"), pop_n = c(1L, 2L, 1L))adj <- matrix( c(0, 1, 0, 1, 0, 1, 0, 1, 0), nrow = 3, byrow = TRUE, dimnames = list(c("A", "B", "C"), c("A", "B", "C")) ) keep <- graph_to_keep(adj, labels = c("A", "B", "C")) keep graph_to_keep(adj, labels = c("A", "B", "C"), pop_n = c(1L, 2L, 1L))
A SpatRaster produced by a Circuitscape Julia run.
data(jl_out)data(jl_out)
SpatRaster with 1 layer
Creates all inputs required by the optimization functions when using Circuitscape run through Julia (via JuliaConnectoR).
jl.prep( n.Pops, response = NULL, CS_Point.File, covariates = NULL, formula = NULL, JULIA_HOME = NULL, Neighbor.Connect = 8, pairs_to_include = NULL, pop2ind = NULL, nb = NULL, parallel = FALSE, cores = NULL, cholmod = TRUE, precision = FALSE, run_test = TRUE, write.files = NULL, write.criteria = NULL, silent = TRUE, scratch = NULL, rm.files = TRUE )jl.prep( n.Pops, response = NULL, CS_Point.File, covariates = NULL, formula = NULL, JULIA_HOME = NULL, Neighbor.Connect = 8, pairs_to_include = NULL, pop2ind = NULL, nb = NULL, parallel = FALSE, cores = NULL, cholmod = TRUE, precision = FALSE, run_test = TRUE, write.files = NULL, write.criteria = NULL, silent = TRUE, scratch = NULL, rm.files = TRUE )
n.Pops |
The number of populations being assessed. |
response |
Vector of pairwise genetic distances (lower half of pairwise
matrix). If using |
CS_Point.File |
A |
covariates |
Data frame of additional covariates for the MLPE model. |
formula |
R formula for the fixed effects of the MLPE model (e.g.,
|
JULIA_HOME |
Path to the Julia |
Neighbor.Connect |
Connectivity scheme for Circuitscape: 4 or 8 (Default = 8). |
pairs_to_include |
Integer vector of 1 (keep) or 0 (drop) for each
pairwise observation. Use |
pop2ind |
Integer vector giving the number of sampled individuals per
population. Default = |
nb |
Maximum distance (in coordinate units) for treating locations as
the same neighborhood. Default = |
parallel |
Logical. Run Circuitscape in parallel? Default = |
cores |
Number of cores for parallel processing (only used when
|
cholmod |
Logical. Use the CHOLMOD direct solver? Default = |
precision |
Logical. Use experimental single-precision mode?
Default = |
run_test |
Logical. Test the Julia/Circuitscape connection before
continuing? Default = |
write.files |
Directory path. If provided, the |
write.criteria |
Minimum run time (seconds) threshold for writing files
when |
silent |
Logical. Suppress Circuitscape progress output?
Default = |
scratch |
Scratch directory for temporary files when the default
|
rm.files |
Logical. Remove temporary files after each Julia run?
Default = |
Julia and Circuitscape.jl must be installed in the selected Julia environment (see https://github.com/Circuitscape/Circuitscape.jl).
Set JULIA_HOME to the bin/ subdirectory of your Julia
installation (e.g., "/usr/local/julia/bin").
The CHOLMOD solver (cholmod = TRUE) uses a direct Cholesky
decomposition that can be much faster than the default algebraic multigrid
solver for moderate problem sizes, but may consume large amounts of memory
for very large problems. It cannot be used with single precision.
Point locations are supplied as a SpatVector and written to a
temporary Circuitscape point file internally.
The requested Circuitscape parallel settings are stored in the returned object. They control Circuitscape execution and are separate from GA fitness parallelism in [GA.prep()] and bootstrap parallelism in [Resist.boot()].
A named list used as jl.inputs. It contains retained and full
pair identifiers (ID and ID.all), MLPE design matrices,
response and covariate data, fitted-model formula, sample points,
connectivity and solver settings, Julia location, file controls, pair
inclusion metadata, and population-to-individual mapping. The object can be
passed directly to [Run_CS.jl()], [SS_optim()], or [MS_optim()].
Bill Peterman <[email protected]>
[Run_CS.jl()], [Run_CS.ini()], [GA.prep()]
## Not run: pts <- terra::vect(sample_pops[[1]], type = "points") jl.inputs <- jl.prep( n.Pops = nrow(sample_pops[[1]]), response = lower(Dc_list[[1]]), CS_Point.File = pts, JULIA_HOME = Sys.getenv("JULIA_BINDIR"), run_test = FALSE ) ## End(Not run)## Not run: pts <- terra::vect(sample_pops[[1]], type = "points") jl.inputs <- jl.prep( n.Pops = nrow(sample_pops[[1]]), response = lower(Dc_list[[1]]), CS_Point.File = pts, JULIA_HOME = Sys.getenv("JULIA_BINDIR"), run_test = FALSE ) ## End(Not run)
Returns matrix entries below the diagonal in the column-major order used by ResistanceGA2 pairwise responses, predictors, and keep vectors.
lower(matrix)lower(matrix)
matrix |
Square matrix. Dimension names are permitted but are not retained in the output. |
Pairwise vectors supplied to the same analysis must all be created with the same population order. A warning is issued for a vector or a nonsquare input.
A vector of length choose(nrow(matrix), 2) containing entries
for pairs (2,1), (3,1), ..., (n,n - 1). This is the same ordering as
matrix[lower.tri(matrix)] and [To.From.ID()].
Bill Peterman <[email protected]>
[To.From.ID()], [mlpe_data()]
m <- matrix(1:16, nrow = 4) lower(m) # returns 6-element vectorm <- matrix(1:16, nrow = 4) lower(m) # returns 6-element vector
Fit maximum-likelihood population effects (MLPE) models using standard
lme4-style formulas while defining dyadic random-intercept terms from
explicit endpoint columns.
mlpe(formula, data, pairs, REML = FALSE, keep = NULL, ...)mlpe(formula, data, pairs, REML = FALSE, keep = NULL, ...)
formula |
Mixed-effects model formula. |
data |
Data frame containing the response, covariates, and endpoint columns. |
pairs |
Definition of one or more MLPE dyad terms. Supply either a
character vector of length two, which is interpreted as
|
REML |
Logical. If TRUE, fit by restricted maximum likelihood for linear mixed models. |
keep |
Optional logical/integer vector used to retain a subset of rows
from |
... |
Additional arguments passed to |
Dyadic random effects are specified in formula using standard random
intercept syntax, for example (1 | pair). The corresponding endpoint
columns are declared through pairs, e.g.
pairs = c("from", "to") or
pairs = list(pair = c("from", "to")).
Additional ordinary lme4 random effects may be included in the same
model. The dyadic covariance contribution is rebuilt after any row subsetting
so omitted pairwise observations are handled correctly.
Fixed-effect coefficients are conditional associations given the declared dyadic and ordinary random effects. Compare likelihood-based criteria only across models fitted to the same rows and by the same likelihood method.
A fitted merMod object. Linear models return an
lmerMod; generalized models return a glmerMod.
[mlpe_data()], [MLPE.lmm()]
y <- matrix(rnorm(25), 5) x <- matrix(rnorm(25), 5) dat <- mlpe_data(response = y, resistance = x) fit <- mlpe( response ~ resistance + (1 | pair), data = dat, pairs = c("from", "to") )y <- matrix(rnorm(25), 5) x <- matrix(rnorm(25), 5) dat <- mlpe_data(response = y, resistance = x) fit <- mlpe( response ~ resistance + (1 | pair), data = dat, pairs = c("from", "to") )
Convert pairwise responses and covariates into a long-form data frame that is
ready for mlpe. Inputs can be supplied either as
symmetric matrices/dist objects, or as long vectors paired with an
explicit two-column pair index.
mlpe_data( response, ..., pairs = NULL, labels = NULL, pair_names = c("from", "to"), keep = NULL )mlpe_data( response, ..., pairs = NULL, labels = NULL, pair_names = c("from", "to"), keep = NULL )
response |
Pairwise response as a lower-triangle vector, symmetric
matrix, |
... |
Named pairwise covariates supplied in the same format as
|
pairs |
Optional data frame or matrix containing at least two columns
that identify the endpoints for each observation. When supplied,
|
labels |
Optional vector of endpoint labels used when converting a
symmetric matrix, |
pair_names |
Names to use for the generated endpoint columns when
|
keep |
Optional logical/integer vector used to retain a subset of pair observations. |
Matrix and dist inputs are converted in the ordering used by
[lower()]. Every response, covariate, pair row, and keep value must
describe the same ordered pair. This function restructures data only; it
does not test whether a distance response is scientifically appropriate for
the intended model.
A data frame with one row per retained pairwise observation. The
first two columns identify pair endpoints, named by pair_names when
generated internally. Supplied covariates retain their names, and the
response is stored in response.
y <- matrix(rnorm(25), 5) x <- matrix(rnorm(25), 5) dat <- mlpe_data(response = y, resistance = x) head(dat) id <- data.frame(from = c(1, 1, 2), to = c(2, 3, 3)) mlpe_data(response = c(0.2, 0.4, 0.6), x = c(1, 2, 3), pairs = id)y <- matrix(rnorm(25), 5) x <- matrix(rnorm(25), 5) dat <- mlpe_data(response = y, resistance = x) head(dat) id <- data.frame(from = c(1, 1, 2), to = c(2, 3, 3)) mlpe_data(response = c(0.2, 0.4, 0.6), x = c(1, 2, 3), pairs = id)
This function is retained for backward compatibility and internal
optimization workflows that still pass a precomputed ZZ matrix. For
new standalone analyses, prefer mlpe_data and
mlpe, which construct dyadic random effects from
explicit endpoint columns and avoid manual To.From.ID() setup.
mlpe_rga(formula, data, REML = FALSE, ZZ = NULL, keep = NULL, ...)mlpe_rga(formula, data, REML = FALSE, ZZ = NULL, keep = NULL, ...)
formula |
An |
data |
Data frame containing aligned pairwise response, predictor, and endpoint variables. |
REML |
Logical. Fit a linear model by restricted maximum likelihood?
Default is |
ZZ |
Optional sparse MLPE design matrix. It is inferred for a complete
pairwise data set when neither |
keep |
Optional binary or logical vector selecting rows of |
... |
Additional arguments passed to |
Runs MLPE as detailed by Clarke et al. (2002) using a standard lme4
formula and either inferred or supplied dyadic random-effect structure.
An AIC value will only be returned if REML = FALSE. In the legacy
one-dyad workflow, the random effect is usually the population vector
pop1 generated by To.From.ID.
A generalized MLPE model can be fit if an appropriate family. It is not possible to use REML when fitting a generalized model.
A fitted lmerMod or glmerMod object, depending on
whether a generalized-model family is supplied through ....
Bill Peterman <[email protected]>
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. doi:10.1198/108571102320
[mlpe()], [mlpe_data()], [MLPE.lmm()]
# Prefer mlpe_data() + mlpe() for new analyses. # Create square 'distance' matrices y <- matrix(rnorm(25), 5) y.pois <- matrix(rpois(25, 5), 5) x <- matrix(rnorm(25), 5) # Create to-from object (5 populations sampled) id <- To.From.ID(5) # Create data frame df <- data.frame(y = lower(y), y.pois = lower(y.pois), x = lower(x), pop = id$pop1) # Fit MLPE model out <- mlpe_rga(formula = y ~ x + (1 | pop), REML = TRUE, data = df) # Fit generalized MLPE model out.pois <- mlpe_rga(formula = y.pois ~ x + (1 | pop), family = poisson, data = df) # Fit model with only select pairs keep <- c(1,1,1,1,1,1,0,1,1,0) out.select <- mlpe_rga(formula = y ~ x + (1 | pop), data = df, keep = keep)# Prefer mlpe_data() + mlpe() for new analyses. # Create square 'distance' matrices y <- matrix(rnorm(25), 5) y.pois <- matrix(rpois(25, 5), 5) x <- matrix(rnorm(25), 5) # Create to-from object (5 populations sampled) id <- To.From.ID(5) # Create data frame df <- data.frame(y = lower(y), y.pois = lower(y.pois), x = lower(x), pop = id$pop1) # Fit MLPE model out <- mlpe_rga(formula = y ~ x + (1 | pop), REML = TRUE, data = df) # Fit generalized MLPE model out.pois <- mlpe_rga(formula = y.pois ~ x + (1 | pop), family = poisson, data = df) # Fit model with only select pairs keep <- c(1,1,1,1,1,1,0,1,1,0) out.select <- mlpe_rga(formula = y ~ x + (1 | pop), data = df, keep = keep)
This function accepts a vector, matrix, dist object, or file path for a
single pairwise resistance-distance predictor and fits a one-predictor MLPE
model. It is retained for backward compatibility. For new analyses, prefer
mlpe_data and mlpe,
which support explicit endpoint-pair columns, omitted observations, and
multiple dyadic random-intercept terms.
MLPE.lmm( resistance, pairwise.genetic, REML = FALSE, ID = NULL, ZZ = NULL, scale = TRUE )MLPE.lmm( resistance, pairwise.genetic, REML = FALSE, ID = NULL, ZZ = NULL, scale = TRUE )
resistance |
Path to a Circuitscape resistance output file, or a
vector, matrix, or |
pairwise.genetic |
Numeric genetic-distance vector in the same
lower-triangle order as |
REML |
Logical. Fit using restricted maximum likelihood? Use
|
ID |
Optional pair-identifier data frame from [gdist.prep()] or [jl.prep()]. It is constructed automatically for a complete set of pairs. |
ZZ |
Optional sparse MLPE design matrix corresponding to |
scale |
Logical. Center and standardize resistance distances before
fitting? Default is |
The resistance coefficient estimates the conditional association
between genetic and resistance distance after accounting for the MLPE
covariance structure. It is not, by itself, a causal effect or evidence of
predictive transfer to new landscapes. When scale = TRUE, the
coefficient is the expected response change for a one-standard-deviation
increase in resistance distance.
A fitted lmerMod object. Use summary() for fixed-effect
estimates and uncertainty, fitted() and residuals() for
diagnostics, and AIC() only for maximum-likelihood fits compared on
the same observations.
Bill Peterman <[email protected]>
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. doi:10.1198/108571102320
[mlpe_data()], [mlpe()], [Diagnostic.Plots()]
set.seed(10) xy <- matrix(stats::runif(12), ncol = 2) genetic <- lower(as.matrix(stats::dist(xy))) resistance <- lower(as.matrix(stats::dist(cbind(xy[, 1], xy[, 2]^2)))) fit <- MLPE.lmm( resistance = resistance, pairwise.genetic = genetic, REML = TRUE ) summary(fit)set.seed(10) xy <- matrix(stats::runif(12), ncol = 2) genetic <- lower(as.matrix(stats::dist(xy))) resistance <- lower(as.matrix(stats::dist(cbind(xy[, 1], xy[, 2]^2)))) fit <- MLPE.lmm( resistance = resistance, pairwise.genetic = genetic, REML = TRUE ) summary(fit)
Optimize multiple resistance surfaces simultaneously using a genetic algorithm.
MS_optim( gdist.inputs = NULL, jl.inputs = NULL, GA.inputs, diagnostic_plots = TRUE )MS_optim( gdist.inputs = NULL, jl.inputs = NULL, GA.inputs, diagnostic_plots = TRUE )
gdist.inputs |
Object from |
jl.inputs |
Object from |
GA.inputs |
Object from |
diagnostic_plots |
Logical. Generate and save diagnostic plots? Default = |
All supplied layers are transformed and combined during one GA search. Percent contributions are spatial averages of each transformed layer's fraction of the composite resistance surface. They describe the fitted composite and are not independent regression effect sizes. Inspect the final coefficient, response plot, and contribution plot together. Do not multiply the output raster by the contribution percentages again.
A list with class resga_ms_optim. It contains the fitted
GA.summary, the maximum-likelihood and REML MLPE fits
(MLPE.model and MLPE.model_REML), a one-row model-fit table
(AICc.tab), optimized resistance-distance matrix (cd), mean
proportional surface contributions (percent.contribution), and
parameter metadata (k). Optimized surfaces and tabular results are
also written beneath GA.inputs$Results.dir.
Bill Peterman <[email protected]>
[GA.prep()], [SS_optim()], [Combine_Surfaces()], [summary.resga_ms_optim()]
## Not run: pts <- terra::vect(sample_pops[[1]], type = "points") gdist.inputs <- gdist.prep( n.Pops = nrow(sample_pops[[1]]), response = lower(Dc_list[[1]]), samples = pts ) GA.inputs <- GA.prep( raster = raster_orig, Results.dir = file.path(tempdir(), "ResistanceGA2-ms"), pop.size = 20, maxiter = 10, run = 5, quiet = TRUE, monitor = FALSE ) ms.out <- MS_optim( gdist.inputs = gdist.inputs, GA.inputs = GA.inputs, diagnostic_plots = FALSE ) ## End(Not run)## Not run: pts <- terra::vect(sample_pops[[1]], type = "points") gdist.inputs <- gdist.prep( n.Pops = nrow(sample_pops[[1]]), response = lower(Dc_list[[1]]), samples = pts ) GA.inputs <- GA.prep( raster = raster_orig, Results.dir = file.path(tempdir(), "ResistanceGA2-ms"), pop.size = 20, maxiter = 10, run = 5, quiet = TRUE, monitor = FALSE ) ms.out <- MS_optim( gdist.inputs = gdist.inputs, GA.inputs = GA.inputs, diagnostic_plots = FALSE ) ## End(Not run)
Plots a transformed continuous resistance surface against the original resistance values
Plot.trans( PARM, Resistance, transformation, print.dir = NULL, marginal.plot = TRUE, marg.type = "histogram", Name = "layer" )Plot.trans( PARM, Resistance, transformation, print.dir = NULL, marginal.plot = TRUE, marg.type = "histogram", Name = "layer" )
PARM |
Parameters to transform a continuous surface or resistance values of a categorical surface. A vector of two parameters is required. The first term is the value of shape parameter (c), and the second term is the value of maximum scale parameter (b) |
Resistance |
A |
transformation |
Transformation equation to apply. Can be provided as the name of the transformation or its numeric equivalent (see details) |
print.dir |
Specify the directory where a .tiff of the transformation will be written (Default = NULL) |
marginal.plot |
Logical. Should distribution plots be added to the margins of the response plot (Default = TRUE). Requires that the resistance surface be specified for 'Resistance'. |
marg.type |
Type of marginal plot to add. One of: [density, histogram, boxplot]. See |
Name |
Name of resistance surface being transformed (optional). This will be added to the output file name. |
This function will create a ggplot object and plot, so it requires ggplot2 to be installed.
Equation names can be:
| 1 = "Inverse-Reverse Monomolecular" | |
| 2 = "Inverse-Reverse Ricker" | |
| 3 = "Monomolecular" | |
| 4 = "Ricker" | |
| 5 = "Reverse Monomolecular" | |
| 6 = 'Reverse Ricker" | |
| 7 = "Inverse Monomolecular" | |
| 8 = "Inverse Ricker" | |
| 9 = "Distance" |
The "Distance" equation sets all cell values equal to 1.
Because of the flexibility of the Ricker function to take a monomolecular
shape (try Plot.trans(PARM=c(10,100), Resistance=c(1,10),
transformation="Ricker") to see this), whenever a shape parameter >6 is
selected in combination with a Ricker family transformation, the
transformation reverts to a Distance transformation. In most applications,
monomolecular families should be preferred and Ricker families should be
explored only when there is a strong biological or ecological rationale for
a unimodal or quadratic-type relationship. Otherwise, the extra flexibility
of Ricker curves can increase the risk of overfitting. This constraint has
not been implemented in the Plot.trans function.
A ggplot object, or a ggExtraPlot object when a
marginal distribution is requested. If print.dir is supplied, a
TIFF copy is written there and the plot object is returned invisibly.
Bill Peterman <[email protected]>
[Resistance.tran()], [GA.prep()]
Plot.trans( PARM = c(2.5, 100), Resistance = raster_orig[["cont_orig"]], transformation = "Monomolecular", marginal.plot = FALSE )Plot.trans( PARM = c(2.5, 100), Resistance = raster_orig[["cont_orig"]], transformation = "Monomolecular", marginal.plot = FALSE )
A SpatRaster with two layers generated from the RandomFields
package. These are the raw (untransformed) surfaces used as optimization inputs.
data(raster_orig)data(raster_orig)
SpatRaster with 2 layers
cat_orig. 5-class categorical resistance surface.
cont_orig. Continuous resistance surface.
A SpatRaster with three layers generated from the RandomFields
package. These surfaces were used with individual-based genetic simulations
conducted using PopGenReport.
data(raster_true)data(raster_true)
SpatRaster with 3 layers
cat_true. 4-class categorical resistance surface. Reclass: 1=1, 2=75, 3=200, 4=25, 5=1.
cont_true. Continuous resistance surface transformed via inverse monomolecular (shape = 2.75, max = 100).
multi_true. Composite surface: sum of transformed categorical and continuous surfaces, rescaled to minimum of 1.
Recommend gdistance commute-distance settings for optimization
workflows based on raster size, number of populations, and whether a
selected-pair workflow is being used.
recommend_commute_settings( raster, n.Pops, keep = NULL, pairs_to_include = NULL )recommend_commute_settings( raster, n.Pops, keep = NULL, pairs_to_include = NULL )
raster |
A |
n.Pops |
Integer number of populations or sample locations. |
keep |
Optional 1/0 or logical keep vector from a
|
pairs_to_include |
Optional 1/0 or logical keep vector from a
|
This helper is intentionally lightweight: it does not change package
defaults, and it does not run any benchmarks itself. Instead, it maps a new
workflow onto the benchmark regimes already represented under
inst/benchmarks/:
benchmark_gdistance_vs_julia.R compares exact
gdistance::commuteDistance() runs with Julia/Circuitscape on
50, 100, and 250 cell-edge rasters using 20 populations.
benchmark_gdistance_approx.R compares aggregate
approximations on a 1000 cell-edge raster using 100 populations and
aggregation factors 2 and 4.
The recommendation path is:
Workflows at or below the 250 x 250 / 20-population benchmark are
treated as "small" and keep exact gdistance settings
(commute.approx = "none").
Workflows between the exact and large aggregate benchmark regimes are
treated as "intermediate" and use
commute.approx = "aggregate" with approx.factor = 2.
Workflows at or above the 1000 x 1000 / 100-population benchmark are
treated as "large" and use approx.factor = 4.
Selected-pair workflows shift the recommendation one band toward more exact settings because the package already reduces commute work when only a subset of pairs is retained.
When an aggregate approximation is recommended, rerun the final best surface once without approximation for final reporting.
A named list containing the recommended
commute.approx/approx.factor settings, the inferred benchmark
workload band, and a short rationale.
r <- terra::rast( nrows = 250, ncols = 250, xmin = 0, xmax = 250, ymin = 0, ymax = 250 ) recommend_commute_settings(r, n.Pops = 20) keep <- c(rep(1L, 25), rep(0L, 20)) recommend_commute_settings(r, n.Pops = 10, keep = keep)r <- terra::rast( nrows = 250, ncols = 250, xmin = 0, xmax = 250, ymin = 0, ymax = 250 ) recommend_commute_settings(r, n.Pops = 20) keep <- c(rep(1L, 25), rep(0L, 20)) recommend_commute_settings(r, n.Pops = 10, keep = keep)
Sample file to be used with examples in the vignette
data(resist_list)data(resist_list)
A list of length 3
Each matrix depicts the pairwise effective resistance distance (measured using 'gdistance') between sample locations across 'resist_true' surfaces
resist_cat. Pairwise effective distance matrix calculated across the categorical resistance surface
resist_cont. Pairwise effective distance matrix calculated across the continuous resistance surface
resist_cont. Pairwise effective distance matrix calculated across the multivariate resistance surface
Repeatedly partitions populations into fitting and holdout sets, refits an MLPE model for each candidate resistance-distance matrix, and summarizes how consistently each candidate ranks first.
Resist.boot( mod.names, dist.mat, n.parameters, sample.prop = 0.75, iters, obs, rank.method = "AIC", genetic.mat, keep = NULL, n.cores = NULL )Resist.boot( mod.names, dist.mat, n.parameters, sample.prop = 0.75, iters, obs, rank.method = "AIC", genetic.mat, keep = NULL, n.cores = NULL )
mod.names |
Character vector naming the candidate surfaces. Names must
correspond by position to |
dist.mat |
List of symmetric, |
n.parameters |
Numeric vector giving the number of fitted parameters in each candidate model, including the intercept and variance parameters used for information-criterion calculations. |
sample.prop |
Proportion of populations sampled without replacement for
model fitting in each iteration. Must leave at least two populations in
at least four populations in the fitting set and two in the holdout set.
Default is |
iters |
Positive integer giving the number of subsampling iterations. |
obs |
Positive integer giving the number of populations represented by each matrix. At least six populations are required. |
rank.method |
Criterion used to rank candidates within each iteration.
One of |
genetic.mat |
Symmetric, |
keep |
Optional binary vector in lower-triangle order. A value of
|
n.cores |
Number of parallel workers. |
This procedure is a population-subsampling stability assessment, not a nonparametric bootstrap of independent observations. Pairwise distances are dependent, so whole populations are assigned to fitting or holdout sets before pairwise vectors are constructed.
AIC, AICc, marginal R-squared, and log likelihood describe fits to the
sampled populations. RMSE is calculated from marginal predictions for the
held-out populations. A candidate with a low average rank and a high
Percent.top is comparatively stable across the sampled population
sets. These summaries do not establish that a resistance surface is the
true ecological process, and they should be interpreted alongside effect
estimates, uncertainty, diagnostics, and a geographic-distance baseline.
Set a random seed before calling this function when reproducible population partitions are required. The same seed produces the same partitions for serial and parallel execution.
A data frame with class resga_bootstrap. Each row represents a
candidate surface and includes avg.AIC, avg.AICc, the mean
AICc weight (avg.weight), mean within-iteration rank
(avg.rank), mean marginal R-squared (avg.R2m), mean log
likelihood (avg.LL), mean holdout RMSE (avg.RMSE), the count
and percentage of first-place finishes (n and
Percent.top), and the supplied parameter count (k).
Bill Peterman <[email protected]>
[MLPE.lmm()], [SS_optim()], [MS_optim()]
## Not run: set.seed(42) xy <- matrix(stats::runif(12), ncol = 2) genetic <- as.matrix(stats::dist(xy)) distance <- as.matrix(stats::dist(xy * cbind(1.2, 0.8))) alternative <- as.matrix(stats::dist(cbind(xy[, 1], xy[, 2]^2))) boot.out <- Resist.boot( mod.names = c("distance", "alternative"), dist.mat = list(distance, alternative), n.parameters = c(4, 4), sample.prop = 0.67, iters = 20, obs = nrow(genetic), rank.method = "RMSE", genetic.mat = genetic, n.cores = 2 ) boot.out ## End(Not run)## Not run: set.seed(42) xy <- matrix(stats::runif(12), ncol = 2) genetic <- as.matrix(stats::dist(xy)) distance <- as.matrix(stats::dist(xy * cbind(1.2, 0.8))) alternative <- as.matrix(stats::dist(cbind(xy[, 1], xy[, 2]^2))) boot.out <- Resist.boot( mod.names = c("distance", "alternative"), dist.mat = list(distance, alternative), n.parameters = c(4, 4), sample.prop = 0.67, iters = 20, obs = nrow(genetic), rank.method = "RMSE", genetic.mat = genetic, n.cores = 2 ) boot.out ## End(Not run)
A SpatRaster with three layers for use in vignette examples.
data(resistance_surfaces)data(resistance_surfaces)
SpatRaster with 3 layers
categorical. 3-class categorical resistance surface.
continuous. Continuous resistance surface.
feature. 2-class categorical (feature) resistance surface.
Applies one of eight resistance transformations (or a distance/flat transformation) to a continuous resistance surface.
Resistance.tran(transformation, shape, max, r, out = NULL)Resistance.tran(transformation, shape, max, r, out = NULL)
transformation |
Transformation to apply. Either a numeric code (1–9) or the name of the transformation (see Details). |
shape |
Shape parameter value. |
max |
Maximum value parameter. |
r |
Resistance surface as a single-layer |
out |
Directory path for exporting the transformed surface as an
|
Valid values for transformation:
| 1 = "Inverse-Reverse Monomolecular" | |
| 2 = "Inverse-Reverse Ricker" | |
| 3 = "Monomolecular" | |
| 4 = "Ricker" | |
| 5 = "Reverse Monomolecular" | |
| 6 = "Reverse Ricker" | |
| 7 = "Inverse Monomolecular" | |
| 8 = "Inverse Ricker" | |
| 9 = "Distance" |
The Distance transformation sets all cell values equal to 1. Because the Ricker function can approximate a monomolecular shape at high shape parameter values, whenever a shape parameter > 6 is selected with a Ricker family transformation, the transformation automatically reverts to Distance. In most optimization workflows, monomolecular families should be treated as the recommended default, and Ricker families should only be explored when there is strong biological or ecological support for a unimodal or quadratic-type relationship. Otherwise, the additional flexibility may increase the risk of overfitting.
A SpatRaster with the transformation applied.
Bill Peterman <[email protected]>
r_trans <- Resistance.tran( transformation = "Monomolecular", shape = 2.5, max = 100, r = raster_orig[["cont_orig"]] ) r_transr_trans <- Resistance.tran( transformation = "Monomolecular", shape = 2.5, max = 100, r = raster_orig[["cont_orig"]] ) r_trans
S3 methods that provide concise printing, structured summaries, fitted
response plots, and model-support plots for outputs returned by
Resist.boot, SS_optim, MS_optim,
and all_comb.
## S3 method for class 'resga_bootstrap' print(x, n = 6L, ...) ## S3 method for class 'resga_bootstrap' summary(object, n = 6L, ...) ## S3 method for class 'summary.resga_bootstrap' print(x, ...) ## S3 method for class 'resga_bootstrap' plot( x, metric = c("avg.weight", "Percent.top", "avg.rank"), top_n = 10L, main = NULL, ... ) ## S3 method for class 'resga_ss_optim' print(x, n = 6L, ...) ## S3 method for class 'resga_ss_optim' summary(object, n = 6L, ...) ## S3 method for class 'summary.resga_ss_optim' print(x, ...) ## S3 method for class 'resga_ss_optim' plot( x, type = c("response", "ranking"), metric = c("AICc", "delta.AICc", "R2m", "LL"), top_n = 10L, model = NULL, main = NULL, ... ) ## S3 method for class 'resga_ms_optim' print(x, n = 6L, ...) ## S3 method for class 'resga_ms_optim' summary(object, n = 6L, ...) ## S3 method for class 'summary.resga_ms_optim' print(x, ...) ## S3 method for class 'resga_ms_optim' plot( x, type = c("response", "contribution", "fit"), metric = c("mean", "AICc", "R2m", "LL"), top_n = 10L, main = NULL, ... ) ## S3 method for class 'resga_all_comb' print(x, n = 6L, ...) ## S3 method for class 'resga_all_comb' summary(object, n = 6L, ...) ## S3 method for class 'summary.resga_all_comb' print(x, ...) ## S3 method for class 'resga_all_comb' plot( x, metric = c("avg.weight", "Percent.top", "AICc"), top_n = 10L, main = NULL, ... )## S3 method for class 'resga_bootstrap' print(x, n = 6L, ...) ## S3 method for class 'resga_bootstrap' summary(object, n = 6L, ...) ## S3 method for class 'summary.resga_bootstrap' print(x, ...) ## S3 method for class 'resga_bootstrap' plot( x, metric = c("avg.weight", "Percent.top", "avg.rank"), top_n = 10L, main = NULL, ... ) ## S3 method for class 'resga_ss_optim' print(x, n = 6L, ...) ## S3 method for class 'resga_ss_optim' summary(object, n = 6L, ...) ## S3 method for class 'summary.resga_ss_optim' print(x, ...) ## S3 method for class 'resga_ss_optim' plot( x, type = c("response", "ranking"), metric = c("AICc", "delta.AICc", "R2m", "LL"), top_n = 10L, model = NULL, main = NULL, ... ) ## S3 method for class 'resga_ms_optim' print(x, n = 6L, ...) ## S3 method for class 'resga_ms_optim' summary(object, n = 6L, ...) ## S3 method for class 'summary.resga_ms_optim' print(x, ...) ## S3 method for class 'resga_ms_optim' plot( x, type = c("response", "contribution", "fit"), metric = c("mean", "AICc", "R2m", "LL"), top_n = 10L, main = NULL, ... ) ## S3 method for class 'resga_all_comb' print(x, n = 6L, ...) ## S3 method for class 'resga_all_comb' summary(object, n = 6L, ...) ## S3 method for class 'summary.resga_all_comb' print(x, ...) ## S3 method for class 'resga_all_comb' plot( x, metric = c("avg.weight", "Percent.top", "AICc"), top_n = 10L, main = NULL, ... )
x |
A ResistanceGA2 result object. |
n |
A positive integer giving the maximum number of table rows to print. |
... |
Additional graphical arguments passed to |
object |
A ResistanceGA2 result object. |
metric |
Character string naming the model-support metric to plot. Lower values are better for AIC, AICc, delta AICc, rank, and RMSE. Higher values are better for model weight, percent top-ranked, R-squared, and log-likelihood. |
top_n |
A positive integer giving the maximum number of models or surfaces to display. |
main |
Optional character string used as the plot title. |
type |
Plot to draw. For single-surface results, use
|
model |
Optional model name for a single-surface response plot. The top-ranked model is used by default. |
The response plot compares observed values with predictions from the fixed effects of the fitted MLPE model. The dashed one-to-one line marks exact agreement. Systematic curvature or changing vertical spread indicates that the fitted mean or variance structure may be inadequate. This plot is a diagnostic, not evidence that an optimized resistance surface will predict new populations or landscapes.
Model-support plots retain the scale of the selected metric. AIC and AICc values are meaningful only within the candidate set fitted to the same response and observations. Percent contributions describe the fitted composite surface; they are not independent effect sizes and should not be applied to the component rasters a second time.
The print() methods return x invisibly. The
summary() methods return a named list with the class shown below:
summary.resga_bootstrap: n_models,
best_surface, metric, ranked table, and print
limit n.
summary.resga_ss_optim: model counts, top surface and AICc,
run time, ranked table, the top fitted model, its
fixed_effects table, and print limit n.
summary.resga_ms_optim: top surface, fit_table,
contribution_table, fitted model,
fixed_effects, and print limit n.
summary.resga_all_comb: whether the analysis contains
replicated runs, the number of runs, top models by run, replicate
consensus, and the single-run model and bootstrap tables when
applicable.
The plot() methods return the data shown in the plot invisibly.
Response plots return columns observed and fitted, where
fitted excludes MLPE random effects and is on the response scale.
SS_optim, MS_optim,
all_comb, Resist.boot,
Diagnostic.Plots
ss <- readRDS(system.file( "extdata", "vignette-results", "ss_optim_results.rds", package = "ResistanceGA2" )) summary(ss) plot(ss, type = "response") plot(ss, type = "ranking", metric = "AICc") ms <- readRDS(system.file( "extdata", "vignette-results", "ms_optim_results.rds", package = "ResistanceGA2" )) summary(ms) plot(ms, type = "contribution")ss <- readRDS(system.file( "extdata", "vignette-results", "ss_optim_results.rds", package = "ResistanceGA2" )) summary(ss) plot(ss, type = "response") plot(ss, type = "ranking", metric = "AICc") ms <- readRDS(system.file( "extdata", "vignette-results", "ms_optim_results.rds", package = "ResistanceGA2" )) summary(ms) plot(ms, type = "contribution")
A named list bundling a small Columbia spotted frog example derived from the GeNetIt and LandGenCourse teaching materials.
data(rga2_demo)data(rga2_demo)
A named list of length 5
sites. A 12-row data frame of selected sites from
GeNetIt::ralu.site, including demo_id, the original site
attributes, and projected UTM coordinates (x, y).
sample_coords. A 12 x 2 numeric matrix of sample coordinates in the same order as the genetic matrix.
site_environment. A data frame of raster values extracted at the 12
sites from rga2_demo_covariates.
genetic_matrix. A 12 x 12 symmetric chord-distance matrix imported from the LandGenCourse radish example.
genetic_vector. The lower-triangle vector of
genetic_matrix in package order.
See data-raw/build_rga2_demo.R for the pinned source revisions used
to construct this demo.
A SpatRaster containing the six Idaho raster layers used to construct
rga2_demo. The layers are derived from GeNetIt and include
cti, err27, ffp, gsp, hli, and
nlcd.
data(rga2_demo_covariates)data(rga2_demo_covariates)
SpatRaster with 6 layers
A 12 x 12 symmetric chord-distance matrix used in rga2_demo. Row and
column names match the demo_id values in rga2_demo$sites.
data(rga2_demo_genetic)data(rga2_demo_genetic)
A 12 x 12 numeric matrix
A low-level wrapper around Circuitscape.jl that writes a complete
Circuitscape .ini file from R, runs compute(), and returns the
raw result plus any generated output files.
Run_CS.ini( r = NULL, CS_Point.File = NULL, JULIA_HOME = NULL, EXPORT.dir = NULL, mask_file = NULL, source_file = NULL, ground_file = NULL, included_pairs_file = NULL, reclass_file = NULL, cs_options = list(), return = c("all", "compute", "matrix", "current_map"), rm.files = FALSE, quiet = TRUE, scratch = NULL, is_resistance = TRUE, sanitize_raster = FALSE, ... )Run_CS.ini( r = NULL, CS_Point.File = NULL, JULIA_HOME = NULL, EXPORT.dir = NULL, mask_file = NULL, source_file = NULL, ground_file = NULL, included_pairs_file = NULL, reclass_file = NULL, cs_options = list(), return = c("all", "compute", "matrix", "current_map"), rm.files = FALSE, quiet = TRUE, scratch = NULL, is_resistance = TRUE, sanitize_raster = FALSE, ... )
r |
A single-layer |
CS_Point.File |
A |
JULIA_HOME |
Path to the Julia |
EXPORT.dir |
Directory where the temporary |
mask_file |
Optional Circuitscape |
source_file |
Optional Circuitscape |
ground_file |
Optional Circuitscape |
included_pairs_file |
Optional Circuitscape
|
reclass_file |
Optional Circuitscape |
cs_options |
Named list of Circuitscape |
return |
What to return: |
rm.files |
Logical. Remove generated |
quiet |
Logical. Suppress Circuitscape progress output emitted through
Julia? Default = |
scratch |
Optional scratch directory for temporary raster / point files. |
is_resistance |
Logical. Should the input habitat map be treated as
resistances ( |
sanitize_raster |
Logical. If |
... |
Additional Circuitscape |
This function complements Run_CS.jl. Whereas
Run_CS.jl is tailored to the package optimization workflow,
Run_CS.ini is intended for direct Circuitscape runs where the user
wants to specify Circuitscape settings in a way that closely mirrors the
.ini file.
Option names in cs_options and ... should use the same keys
that appear in a Circuitscape .ini file. The function currently
exposes the full raster-mode option template written by the package for
Circuitscape 4.0.5 and validates that overridden keys are known.
The file-oriented arguments (mask_file, source_file,
ground_file, included_pairs_file, and reclass_file)
are convenience layers over the corresponding Circuitscape options. They are
especially useful when working directly with terra objects or small R tables.
For any more specialized file layout, users can still supply a file path
directly through these arguments or through cs_options.
Depending on return, either the raw Julia result, a parsed
pairwise matrix, a cumulative current-map raster, or a list containing all
of those plus the written .ini file path, generated output files,
and the final Circuitscape configuration used for the run.
Bill Peterman <[email protected]>
## Not run: pts <- terra::vect(sample_pops[[1]], type = "points") cs.out <- Run_CS.ini( r = raster_orig[["cont_orig"]], CS_Point.File = pts, JULIA_HOME = Sys.getenv("JULIA_BINDIR"), scenario = "pairwise", connect_four_neighbors_only = TRUE, write_cur_maps = TRUE, write_cum_cur_map_only = TRUE, rm.files = FALSE ) ## End(Not run)## Not run: pts <- terra::vect(sample_pops[[1]], type = "points") cs.out <- Run_CS.ini( r = raster_orig[["cont_orig"]], CS_Point.File = pts, JULIA_HOME = Sys.getenv("JULIA_BINDIR"), scenario = "pairwise", connect_four_neighbors_only = TRUE, write_cur_maps = TRUE, write_cum_cur_map_only = TRUE, rm.files = FALSE ) ## End(Not run)
Executes Circuitscape through Julia (via JuliaConnectoR) and returns pairwise resistance distances.
Run_CS.jl( jl.inputs = NULL, r, CurrentMap = FALSE, full.mat = FALSE, EXPORT.dir = NULL, output = "matrix", CS_Point.File = NULL, pairs_to_include = NULL, pop2ind = NULL, parallel = FALSE, cores = NULL, cholmod = TRUE, precision = FALSE, JULIA_HOME = NULL, rm.files = TRUE, scratch = NULL, is_resistance = TRUE )Run_CS.jl( jl.inputs = NULL, r, CurrentMap = FALSE, full.mat = FALSE, EXPORT.dir = NULL, output = "matrix", CS_Point.File = NULL, pairs_to_include = NULL, pop2ind = NULL, parallel = FALSE, cores = NULL, cholmod = TRUE, precision = FALSE, JULIA_HOME = NULL, rm.files = TRUE, scratch = NULL, is_resistance = TRUE )
jl.inputs |
Object from |
r |
A single-layer |
CurrentMap |
Logical. Generate the cumulative current map?
Default = |
full.mat |
Logical. Return the full square distance matrix instead of
only the lower triangle? Default = |
EXPORT.dir |
Directory where Circuitscape results are written. Required
when |
output |
|
CS_Point.File |
Sample locations. Only needed when |
pairs_to_include |
1/0 vector selecting pairs. Only needed when
|
pop2ind |
Per-population sample-size vector. Only needed when
|
parallel |
Logical. Run Circuitscape in parallel? Default = |
cores |
Number of cores for parallel runs. |
cholmod |
Logical. Use CHOLMOD solver? Default = |
precision |
Logical. Use single precision? Default = |
JULIA_HOME |
Path to the Julia |
rm.files |
Logical. Remove temporary files after the run?
Default = |
scratch |
Scratch directory for temporary files. |
is_resistance |
Logical. Is |
Users should supply terra objects directly. The raster and point files needed by Circuitscape are written to temporary files internally for each run.
Lower-half pairwise resistance vector, a full square matrix
(full.mat = TRUE), or a SpatRaster current map
(output = "raster").
Bill Peterman <[email protected]>
## Not run: pts <- terra::vect(sample_pops[[1]], type = "points") jl.inputs <- jl.prep( n.Pops = nrow(sample_pops[[1]]), response = lower(Dc_list[[1]]), CS_Point.File = pts, JULIA_HOME = Sys.getenv("JULIA_BINDIR"), run_test = FALSE ) cs.result <- Run_CS.jl(jl.inputs = jl.inputs, r = raster_orig[["cont_orig"]]) ## End(Not run)## Not run: pts <- terra::vect(sample_pops[[1]], type = "points") jl.inputs <- jl.prep( n.Pops = nrow(sample_pops[[1]]), response = lower(Dc_list[[1]]), CS_Point.File = pts, JULIA_HOME = Sys.getenv("JULIA_BINDIR"), run_test = FALSE ) cs.result <- Run_CS.jl(jl.inputs = jl.inputs, r = raster_orig[["cont_orig"]]) ## End(Not run)
Calculates pairwise cost or commute distances across a resistance surface
using the gdistance package.
Run_gdistance( gdist.inputs, r, scl = TRUE, return.error.value = FALSE, commute.approx = NULL, approx.factor = NULL, approx.scale = NULL )Run_gdistance( gdist.inputs, r, scl = TRUE, return.error.value = FALSE, commute.approx = NULL, approx.factor = NULL, approx.scale = NULL )
gdist.inputs |
Object created by |
r |
A single-layer |
scl |
Logical. Scale the correction values (default |
return.error.value |
Logical. If |
commute.approx |
Optional override for the commute-distance
approximation stored in |
approx.factor |
Optional integer aggregation factor override used when
|
approx.scale |
Optional logical override. If |
gdistance still operates on raster objects. This function keeps
the public API terra-based by coercing the supplied SpatRaster
internally before building the transition object. The default transition
function uses vectorized sparse builders for projected 4- and 8-neighbor
transition matrices and projected geographic corrections. Selected-pair runs
reuse a single distance calculation over the required endpoints.
Commute-distance runs solve the reduced Laplacian for all requested endpoints
in one batch.
A numeric distance vector (or matrix) of pairwise cost/commute
distances. If return.error.value = TRUE, returns -99999 on
error or warning.
Bill Peterman <[email protected]>
pts <- terra::vect(samples[, 2:3], type = "points") gdist.inputs <- gdist.prep( n.Pops = nrow(samples), samples = pts, method = "costDistance" ) r_surface <- terra::subset(terra::unwrap(resistance_surfaces), "continuous") cd <- Run_gdistance(gdist.inputs, r_surface) length(cd)pts <- terra::vect(samples[, 2:3], type = "points") gdist.inputs <- gdist.prep( n.Pops = nrow(samples), samples = pts, method = "costDistance" ) r_surface <- terra::subset(terra::unwrap(resistance_surfaces), "continuous") cd <- Run_gdistance(gdist.inputs, r_surface) length(cd)
50 sample locations per list element corresponding to populations used in genetic simulations. Each element is a two-column numeric matrix (x, y).
data(sample_pops)data(sample_pops)
A list of length 3; each element is a 50 x 2 coordinate matrix.
sample_cat. Sample locations for categorical surface simulation.
sample_cont. Sample locations for continuous surface simulation.
sample_multi. Sample locations for multi-surface simulation.
Sample file to be used with examples in the vignette
data(samples)data(samples)
A 25 x 3 data frame
A data frame that can be saved as a .txt file for running examples in the vignette
Calculates the pairwise di and dhm covariates described by
Prunier et al. (2017) from site-level effective population sizes.
SHNe(n.samples, pop.size, output = c("df", "matrix"))SHNe(n.samples, pop.size, output = c("df", "matrix"))
n.samples |
Positive integer giving the number of sampled populations. |
pop.size |
Positive numeric vector of effective population sizes, ordered to match the populations used in all other pairwise inputs. |
output |
Output representation: |
For populations and , the function calculates
and
. These are demographic covariates. Their
coefficients in an MLPE model describe conditional associations and should
not be interpreted as causal effects without an appropriate study design.
With output = "df", a data frame with choose(n.samples,
2) rows and columns dhm and di, in [lower()] order. With
output = "matrix", a named list of two n.samples square
matrices whose lower triangles contain those values.
Bill Peterman <[email protected]>
Prunier, J. G., V. Dubut, L. Chikhi, and S. Blanchet. 2017. Contribution of spatial heterogeneity in effective population sizes to the variance in pairwise measures of genetic differentiation. Methods in Ecology and Evolution 8:1866-1877. doi:10.1111/2041-210X.12820
[site_pairwise_diff()], [mlpe_data()], [mlpe()]
Ne <- rep(c(5, 25, 100), each = 3) SHNe(n.samples = length(Ne), pop.size = Ne) # Return as list of matrices SHNe(n.samples = length(Ne), pop.size = Ne, output = "matrix")Ne <- rep(c(5, 25, 100), each = 3) SHNe(n.samples = length(Ne), pop.size = Ne) # Return as list of matrices SHNe(n.samples = length(Ne), pop.size = Ne, output = "matrix")
Convert site-level environmental variables into pairwise covariates that are aligned with the lower-triangle ordering used throughout ResistanceGA2.
site_pairwise_diff( data, vars = NULL, scale = FALSE, transform = c("absolute", "squared"), keep = NULL, pop_n = NULL, suffix = NULL )site_pairwise_diff( data, vars = NULL, scale = FALSE, transform = c("absolute", "squared"), keep = NULL, pop_n = NULL, suffix = NULL )
data |
Numeric vector, matrix, or data frame of site-level environmental variables. |
vars |
Optional character or integer index selecting columns from
|
scale |
Logical. If |
transform |
Transformation applied to pairwise differences:
|
keep |
Optional logical/integer vector used to retain a subset of the final pairwise observations. |
pop_n |
Optional integer vector giving the number of sampled
individuals per population. When supplied, population-level environmental
differences are expanded to the individual-level pair structure used by
|
suffix |
Optional suffix appended to output column names. Defaults to
|
These helpers are intended for isolation-by-environment workflows, where
environmental mismatch between sample locations is used as one or more fixed
effects in an MLPE model. The output can be supplied directly to the
covariates argument of gdist.prep or
jl.prep.
A data frame of pairwise environmental-difference covariates.
env <- data.frame( temp = c(10, 12, 15, 11), precip = c(100, 120, 80, 110) ) site_pairwise_diff(env) site_pairwise_diff(env, vars = "temp", scale = TRUE)env <- data.frame( temp = c(10, 12, 15, 11), precip = c(100, 120, 80, 110) ) site_pairwise_diff(env) site_pairwise_diff(env, vars = "temp", scale = TRUE)
Convert one or more site-level environmental variables into a single pairwise environmental-distance covariate aligned with the lower-triangle ordering used throughout ResistanceGA2.
site_pairwise_dist( data, vars = NULL, scale = TRUE, method = "euclidean", keep = NULL, pop_n = NULL, name = "env_dist", ... )site_pairwise_dist( data, vars = NULL, scale = TRUE, method = "euclidean", keep = NULL, pop_n = NULL, name = "env_dist", ... )
data |
Numeric vector, matrix, or data frame of site-level environmental variables. |
vars |
Optional character or integer index selecting columns from
|
scale |
Logical. If |
method |
Distance metric passed to |
keep |
Optional logical/integer vector used to retain a subset of the final pairwise observations. |
pop_n |
Optional integer vector giving the number of sampled
individuals per population. When supplied, population-level environmental
distances are expanded to the individual-level pair structure used by
|
name |
Name of the output covariate column. Default = |
... |
Additional arguments passed to |
A one-column data frame containing pairwise environmental distance.
env <- data.frame( temp = c(10, 12, 15, 11), precip = c(100, 120, 80, 110) ) site_pairwise_dist(env) site_pairwise_dist(env, method = "manhattan", name = "climate_dist")env <- data.frame( temp = c(10, 12, 15, 11), precip = c(100, 120, 80, 110) ) site_pairwise_dist(env) site_pairwise_dist(env, method = "manhattan", name = "climate_dist")
Optimize each surface in a SpatRaster stack using a genetic algorithm
executed with the ga function in GA.
SS_optim( gdist.inputs = NULL, jl.inputs = NULL, GA.inputs, dist_mod = TRUE, null_mod = TRUE, diagnostic_plots = TRUE )SS_optim( gdist.inputs = NULL, jl.inputs = NULL, GA.inputs, dist_mod = TRUE, null_mod = TRUE, diagnostic_plots = TRUE )
gdist.inputs |
Input list returned by [gdist.prep()]. Supply exactly one distance backend. |
jl.inputs |
Input list returned by [jl.prep()]. Supply exactly one distance backend. |
GA.inputs |
Settings and resistance surfaces returned by [GA.prep()]. |
dist_mod |
Logical. Add a geographic-distance baseline to the common
model comparison? Default is |
null_mod |
Logical. Add an intercept-only baseline? Default is
|
diagnostic_plots |
Logical. Save [Diagnostic.Plots()] for optimized
candidates? Default is |
Each raster layer is optimized independently. AICc comparisons are
meaningful only within the returned table, where candidates use the same
response and retained observations. The smallest AICc identifies the
best-supported candidate among those evaluated, not a uniquely true
resistance process. Use summary() to inspect the leading model and
its fixed effects, and plot(..., type = "response") to assess its
response relationship.
A list with class resga_ss_optim. It contains continuous and
categorical optimization tables, the common model-comparison table
(AICc), coefficient table (MLPE), elapsed Run.Time,
fitted REML models by candidate (MLPE.list), optimized pairwise
distance matrices (cd), parameter metadata (k), and fitted
GA objects (ga). CSV tables, fitted objects, optimized surfaces, and
optional diagnostics are written beneath GA.inputs$Results.dir.
Bill Peterman <[email protected]>
[GA.prep()], [MS_optim()], [Resist.boot()], [summary.resga_ss_optim()]
## Not run: pts <- terra::vect(sample_pops[[1]], type = "points") gdist.inputs <- gdist.prep( n.Pops = nrow(sample_pops[[1]]), response = lower(Dc_list[[1]]), samples = pts ) GA.inputs <- GA.prep( raster = raster_orig, Results.dir = file.path(tempdir(), "ResistanceGA2-ss"), pop.size = 20, maxiter = 10, run = 5, quiet = TRUE, monitor = FALSE ) ss.out <- SS_optim( gdist.inputs = gdist.inputs, GA.inputs = GA.inputs, dist_mod = FALSE, null_mod = FALSE, diagnostic_plots = FALSE ) ## End(Not run)## Not run: pts <- terra::vect(sample_pops[[1]], type = "points") gdist.inputs <- gdist.prep( n.Pops = nrow(sample_pops[[1]]), response = lower(Dc_list[[1]]), samples = pts ) GA.inputs <- GA.prep( raster = raster_orig, Results.dir = file.path(tempdir(), "ResistanceGA2-ss"), pop.size = 20, maxiter = 10, run = 5, quiet = TRUE, monitor = FALSE ) ss.out <- SS_optim( gdist.inputs = gdist.inputs, GA.inputs = GA.inputs, dist_mod = FALSE, null_mod = FALSE, diagnostic_plots = FALSE ) ## End(Not run)
A SpatRaster cumulative current map produced by
Run_CS.jl for the "continuous" layer of
resistance_surfaces (after terra::unwrap() when
needed) and samples.
data(vignette_current_map)data(vignette_current_map)
SpatRaster with 1 layer
A named list containing lightweight, precomputed outputs from the current terra-based workflow.
data(vignette_pairwise_outputs)data(vignette_pairwise_outputs)
A named list of length 2
gdistance_commute. Full square commute-distance matrix returned by
Run_gdistance for the "continuous" layer of
resistance_surfaces (after terra::unwrap() when
needed) and samples.
circuitscape_matrix. Lower-half pairwise resistance vector returned
by Run_CS.jl for the same surface and sample locations.