Mechanistic species distribution modelling
This vignette describes the available options of incorporating mechanistic modelling approaches in the ibis.iSDM package. These approaches can broadly be separated as those that are “added” to existing modelling routines, or how ibis.iSDM outputs can be used as input in mechanistic modelling. For both this package provides some basic wrappers.
Before we delve into the options to do mechanistic SDMs in the ibis.iSDM package, it is useful to remind us what the term ‘mechanism’ actually means. In the literature there are a range of different definitions, sometimes referring to mechanistic SDMs as those that incorporate ecological processes (e.g. demography, dispersal, eco-evolutionary principles). Yet often, correlative SDMs are also declared “mechanistic” if they somehow incorporate a specific constrain or response function towards an environmental variable. For example, the micro-climatic limits to the persistence of a species Briscoe et al. 2023, or the presence of biotic interactions (estimated through a separate SDM of a different species) are also sometimes referred to as limiting mechanisms (Ohlmann et al. 2023).
The latter approaches - which are largely a fine-tuning of specific
response function - can to some extent be emulated through creating
specific derivates or by adding covariate priors
add_priors()) or model predictors
add_predictors_model()) to a SDM. Further methods will be
added to the package as they become available or can be readily
incorporated into the modelling framework. Some types of integration can
also be directly modelled through integration. More details are provided
in other vignettes such as the data preparation for creating
derivates or the vignette on data
integration. Users of the package are also directed to the various
add_constraint() functions, many of which enable
corrections of projected scenarios.
There are a range of wrappers implemented in ibis.iSDM that
allow a convenient passing of outputs or parameters to other mechanistic
modelling packages. These wrappers support the convenient addition of
ecological processes such as dispersal to scenarios data integration. Or they
enable ibis.iSDM outputs to directly become inputs to further
simulations. In case key parameters are not available, package users are
encouraged to check out the various options in the
Most of these mechanistic approaches require quite extensive model understanding and in many cases additional training. Furthermore a range of parameters are usually required so that outputs are meaningful. It is beyond the scope of this vignette to provide an introduction to the various models. Rather, it is only demonstrated how linkages between ibis.iSDM and other models can be made, and the reader is referred to the original publication underlying each approach (see help page for references).
The KISSMig model provide a simple model to estimate and limit dispersal in species distribution models (Nobis & Normand, 2014). It does not include any ecological mechanism related to recruitment as such, but instead works as a simple stochastic migration estimator that allows the inclusion of time-lagged dispersal to local neighborhoods. In the ibis.iSDM package the KISSMig simulator can be added as a dispersal constraint (among others) to scenario objects.
library(ibis.iSDM) library(terra) #> terra 1.7.55 library(ggplot2) # Don't print out as many messages options("ibis.setupmessages" = FALSE) # Background and biodiversity data background <- terra::rast(system.file('extdata/europegrid_50km.tif', package='ibis.iSDM')) virtual_points <- sf::st_read(system.file('extdata/input_data.gpkg', package='ibis.iSDM'), 'points', quiet = TRUE) # Add some pseudo-absence information for later poa <- virtual_points |> add_pseudoabsence(field_occurrence = 'Observed', template = background) # Note we are loading different predictors than in previous examples # These are in netcdf4 format, a format specific for storing spatial-temporal data including metadata. ll <- list.files(system.file("extdata/predictors_presfuture/", package = "ibis.iSDM", mustWork = TRUE), "*.nc",full.names = TRUE) # From those list of predictors are first loading the current ones as raster data # We are loading only data from the very first, contemporary time step for model fitting pred_current <- terra::rast() for(i in ll) suppressWarnings( pred_current <- c(pred_current, terra::rast(i, lyrs = 1) ) ) names(pred_current) <- tools::file_path_sans_ext( basename(ll) ) # Get future predictors # These we will load in some time steps using the stars package and ignoring the first time step suppressWarnings( pred_future <- stars::read_stars(ll) |> stars:::slice.stars('Time', seq(2,86,by=10)) ) sf::st_crs(pred_future) <- sf::st_crs(4326) # Set projection # Rename future predictors to those of current names(pred_future) <- names(pred_current) # ------ # # Fit a model fit <- distribution(background) |> add_biodiversity_poipa(poa, field_occurrence = 'Observed', name = 'Virtual points') |> # Note that we scale the predictors here add_predictors(pred_current, transform = 'scale',derivates = 'none') |> engine_glmnet(alpha = 0) |> # Train the model train(verbose = FALSE) |> # Add simple percentile thresholds threshold(method = 'percentile', value = .33) # Show the threshold fit$plot_threshold()
Now lets add KISSMig as dispersal constraints. The constrain
directly used the fitted suitability estimates from each projected
timestep and also makes use of the created thresholded layer. Per
time-step dispersal events are stochastically simulated and constraint
range expansions in the next modelling steps. See
?kissmig::kissmig and the help-page
for further help and explanations of parameters.
# Create a scenario object sc <- scenario(fit) |> # Apply the same variable transformations as above. add_predictors(pred_future, transform = 'scale') |> # Calculate thresholds at each time step. The threshold estimate is taken from # the fitted model object. threshold() # Add KISSMig constraint sc1 <- sc |> add_constraint_dispersal(method = "kissmig", type = "DIS", # Final distribution result value = 10, # Number of iteration steps # These parameters are for KISSMig and get passed on # Probablitiy of local extinction between iterations pext = 0.5, # Probability corner cells are colonized. pcor = 0.2 ) sc2 <- sc |> add_constraint_dispersal(method = "kissmig", type = "DIS", # Final distribution result value = 10, # Number of iteration steps # These parameters are for KISSMig and get passed on # Probablitiy of local extinction between iterations pext = 0.9, # Probability corner cells are colonized. pcor = 0.1 ) # Project two scenarios with varying local extinction probability df1 <- project(sc1, verbose = FALSE) |> summary() #> Linking to GEOS 3.10.2, GDAL 3.4.1, PROJ 8.2.1; sf_use_s2() is FALSE df2 <- project(sc2, verbose = FALSE) |> summary() df <- dplyr::bind_rows(df1 |> dplyr::mutate(scenario = "low"), df2 |> dplyr::mutate(scenario = "high") ) # ------ # ggplot(df, aes(x = band, y= area_km2/1e6, group = scenario, color = scenario)) + theme_bw(base_size = 16) + geom_line() + labs(x = "Time", y = "Occupied area (Million km2)")
Another dispersal simulator is MIGCLIM (Engler et al. 2014), a stochastic simulator which innovatively allows to differentiate between short and long-distance dispersal events as well as varying propagule pressure. Unfortunately is currently not available on CRAN anymore (stand September 2023), possibly because of lack of maintenance or missing dependency. The package can still be downloaded from github however (https://github.com/robinengler/MigClim/).
Assuming the user is able to install the MigClim package with all it’s dependencies (some of which have also disappared from CRAN), it can be run in ibis.iSDM as follows:
prj <- scenario(fit) |> # Apply the same variable transformations as above. add_predictors(pred_future, transform = 'scale') |> # Calculate thresholds at each time step. The threshold estimate is taken from # the fitted model object. threshold() |> # Check the help files for the function for an explanation of the parameters. add_constraint_MigClim(rcThresholdMode = 'continuous', dispSteps = 1, dispKernel = c(1.0, 0.4, 0.16, 0.06, 0.03), barrierType = "strong", lddFreq = 0, lddRange = c(1000, 10000), iniMatAge = 1, propaguleProdProb = c(0.2, 0.6,0.8, 0.95), replicateNb = 10) |> # Project the model project() # MIGCLIM outputs are provided a single updated layer and can be plotted through # a customized plotting function. prj$plot_migclim()
This example will be updated as an update for current R versions (>3.0) becomes available.
The steps package implements a spatial-temporal explicit metapopulation simulator (Visintin et al. 2021) that is able to account for varying vital rates, dispersal with barriers and density dependence. steps is a simulator, thus makes use of a range of parameters and critical correlative assumptions to estimate for example abundance of a given time step.
In the ibis.iSDM package the linkage to steps can be established directly in the scenario projections by simply adding it as a separate module. When added, steps is used to make spatial-temporal abundance estimates that are aligned with each projection time step, eventual specified barriers and the provided parameters with regards to vital rates and density-dependence.
Note: The wrapper functionality implemented in the ibis.iSDM package is based on the assumption that higher habitat suitability (as estimated by a correlative SDM) is linearly correlated to higher population abundance. It should be noted that such assumptions have been questioned before and should be interpreted with caution (Lee-Yaw et al. 2021).
Users should always clearly understand the rationale behind parameter choices!
require("steps") #> Loading required package: steps # Define some arbitrary vital rates for the transition for this purpose # Define vital rates vt <- matrix(c(0.0,0.52,0.75, 0.52,0.28,0.0, 0.0,0.52,0.75), nrow = 3, ncol = 3, byrow = TRUE) colnames(vt) <- rownames(vt) <- c('juvenile','subadult','adult') # We again specify a scenario as before using the fitted model prj <- scenario(fit) |> # Apply the same variable transformations as above. add_predictors(pred_future, transform = 'scale') |> # Calculate thresholds at each time step. The threshold estimate is taken from # the fitted model object. threshold() |> # We then specify that we we simulate_population_steps(vital_rates = vt) # Notice how we have added steps as additional simulation outcome prj #> Spatial-temporal scenario: #> Used model: GLMNET-Model #> --------- #> Predictors: bio01, bio12, crops, ... (9 predictors) #> Time period: 2016-01-01 -- 2096-01-01 (79.9 years) #> --------- #> Threshold: 0.514 (percentile) #> Simulations: simulation (steps) #> --------- #> Scenarios fitted: None # Now project scenario1 <- project(prj) plot(scenario1, "population")
# Also see a different one where we add a dispersal constraint and density dependence dispersal <- steps::fast_dispersal(dispersal_kernel = steps::exponential_dispersal_kernel(distance_decay = 1)) scenario2 <- project(prj |> simulate_population_steps(vt, dispersal = dispersal, density_dependence = steps::ceiling_density(3) ) ) # We can see that the dispersal constraint and higher density dependence cleary # results in a population abundance that tends to be concentrated in central Europe. plot(scenario2, "population")