Advances in global sensitivity analyses of demographic-based species distribution models to address uncertainties in dynamic landscapes

Whitebark pine coupled SDM-population dynamics model

Whitebark pine is classified as an Endangered species (Mahalovich & Stritch, 2013) in decline primarily due to blister rust infections (Cronartium ribicola J.C. Fisch.), although the mountain pine beetle (Dendroctonus ponderosae) epidemic and fire suppression are also major threats. The original stochastic spatially-explicit metapopulation dynamics model was undertaken to evaluate the influence of blister rust infection on whitebark pine viability and persistence in Mount Rainier National Park, Washington, USA, where it is in danger of local extinction (Ettl & Cottone, 2004). As such, it is also referred to as a population viability analysis (PVA).

The whitebark pine model (Ettl & Cottone, 2004) was implemented using RAMAS GIS software (Akçakaya, 2002) and relied on the RAMAS Spatial module to identify populations (i.e., habitat patches) using a raster-based habitat suitability map of a 420.3 km² region in the park. This population map provides the basis on which metapopulation dynamics are simulated in the RAMAS Metapop module. A total of 46 populations were identified by the RAMAS patch detection algorithm, which was consistent with the actual distribution of trees (Ettl & Cottone, 2004). The authors used aerial photography to estimate initial abundances within each of the populations. Vegetation plots were used to estimate the proportion of individuals in each of the 13 stage classes and the final stage matrix included 24 transition probabilities that modeled both healthy and trees infected with blister rust. The stage classes included four seedling stages, saplings, infected saplings, non-reproductive adults, infected non-reproductive adults, healthy adult trees (Class 1), and three adult stages with various degrees of blister rust infection (Classes 2–4; Table 1). Adults in Class 1 through 4, non-reproductive, and infected non-reproductive adults were affected by a ceiling model of density dependence.

Seed dispersal among pairs of populations was distance-dependent whereby dispersal rates are calculated as a function of distance between populations and the parameters of the dispersal-distance function of the original model. The general form of the model is: ${\displaystyle m_{ij}=a\times \exp \left(-D_{ij}^c∕b\right)}$ where m_ij is the dispersal rate between the ith and jth populations, a (scaling parameter), b, and c are function parameters, and D_ij is the distance between a pair of populations. The dispersal distance function parameters for the whitebark pine model were a = 0.1, b = 0.5, c = 1, and D_max = 2.5 km (maximum dispersal distance).

Blister rust disease was modeled as a catastrophe where the probability of invasion of a site corresponded to the average invasion time of 8 years. When blister rust occurs, healthy individuals become infected, leading to declines in survival and fecundity rates. Specifically, saplings, non-reproductive adults, and healthy adults (Class 1) transition to infected sapling, infected non-reproductive adults, and Class 2 adult stage, respectively. Once sapling or non-reproductive adult stages become infected, they continue on to death, while reproductive adult stages cycle through to Class 4, where fecundities are lower. The proportion of individuals that transition into an infected class was based on empirical data and calibration tests in the original model. The model includes both demographic and environmental stochasticity, with the latter based on a lognormal distribution. A full description of the model is detailed in Ettl & Cottone (2004) and stage classes are defined in Table 1.

RAMAS Spatial uses spatial data representing habitat requirements for a species, a user-defined habitat suitability function, and habitat-specific demographic relationships to derive metapopulation structure (Akçakaya & Root, 2005). The original study included a habitat suitability map at a resolution of 0.02 km with suitability values ranging from 0 to 5 and a threshold value of 1 was applied to distinguish between non-suitable and suitable habitat. We estimated habitat-specific demographic relationships because these were not provided by the authors. Using data available in the original whitebark pine PVA, we fit linear regression models to define the relationships of total habitat suitability (ths, sum of habitat suitability values of all raster cells belonging to a patch), to patch-specific initial abundances and carrying capacity, separately. We used ths instead of habitat area as this incorporates a measure of habitat area weighted by habitat quality. Our estimated relationship between ths and initial abundance was; ${\displaystyle \text{Initial abundance}=227.32\ast ths}$ with an r² = 0.780, p < 0.01. The relationship between ths and carrying capacity, K was; ${\displaystyle K=0.8449\ast ths}$ with an r² = 0.765, p < 0.01. Apart from standardizing the time horizon to 100 years and extinction threshold to 0, and formatting the input files for use with GRIP 2.0 (see Appendix S1) the model was not further modified. We assumed that the remaining parameter values and model structure synthesized the best available information.

Overview of GRIP 2.0

We created GRIP 2.0 in R (Generating Random Input Parameters; developed and tested in v. 2.7.1 through 3.2.1 (R Core Team, 2015), which interacts with RAMAS Spatial and Metapop, two modules of the software RAMAS GIS (Akçakaya & Root, 2005), and uses R-spatial packages to facilitate spatial and geostatistical analyses. There are two versions available, each compatible with a one of RAMAS GIS v 5.0 or v 6.0. GRIP 2.0 builds on GRIP 1.0 (Curtis & Naujokaitis-Lewis, 2008) and uses Monte Carlo methods to vary parameters (1) directly on the grid-based habitat suitability map used as an input to RAMAS Spatial, (2) specified in the RAMAS Spatial module, and (3) specified in the RAMAS Metapop module (Fig. 1). Unique sets of input parameters for the Spatial and Metapop modules of RAMAS GIS are drawn from user-defined random distributions reflecting parameter uncertainty and the range of plausible values.

Once the replicate landscapes and corresponding metapopulation models are defined, GRIP 2.0 runs the replicate simulations in batch mode, and collates the input parameter values and predictions into a comma delimited file for subsequent analyses. Thus, GRIP 2.0 automates the generation of unique sets of replicate stochastic simulation files and manages the simulations with batch files, which would otherwise be prohibitive if undertaken manually (McCarthy, Burgman & Ferson, 1995). The code of this global sensitivity analysis program is annotated and easily customized to reflect a particular species’ biology or address related research questions. In the following sections, we briefly describe how habitat attributes and demographic parameters are varied in GRIP 2.0, but details are also annotated in Appendix S1 and the landscape generator and routine for varying demographic parameters is described in more detail in Appendix S2.

Generation of alternative landscapes

GRIP 2.0 generates alternative landscapes by varying the number of patches (i.e., populations; herein we use the terms patch and population interchangeably), patch size, and the habitat suitability value associated with each raster spatial data cell directly on the original habitat suitability and patch maps (Table 1). The current version of GRIP 2.0 includes four options for varying habitat suitability maps: ‘random.normal’, ‘spatially autocorrelated’, ‘ensemble’, and ‘no variation’ (Table 1). In the ‘random.normal’ option, habitat suitability (HS) values are simulated by drawing a value for each grid cell from a normal distribution based on the mean and standard deviation of HS values within the reference landscape. There is no spatial autocorrelation amongst habitat suitability values. The ‘spatially autocorrelated’ option provides flexibility to vary the degree of autocorrelation among HS values simulated using the RandomFields R package (Schlather et al., 2015). Base settings include a Gaussian model of spatial autocorrelation with a mean of 0, variance of 5, nugget value of 1, and a scale of 10, creating a generally highly correlated surface. All HS values are set to be equal to or greater than the newly sampled HS threshold value (see below). The ‘ensemble’ approach can evaluate SDM-based uncertainties associated with alternative models (e.g., based on different model algorithms) whereby a number of models are fit (an ensemble) and then combined into a consensus prediction. In GRIP 2.0, users must include a raster-based ensemble SDM prediction layer reflecting a consensus estimate of habitat suitability and a raster reflecting uncertainties in that prediction. The sensitivity analysis of habitat suitability thus reflects cell-based uncertainty in the prediction and captures error propagation across multiple SDMs. Finally, the ‘no variation’ option ensures that HS is not included in the GSA. Given information contained in the original model, we varied HS for the whitebark pine based on the ‘random.normal’ option.

Additionally, two user-specified settings in RAMAS Spatial that inform its patch detection algorithm were varied: the habitat suitability threshold, which is used to distinguish between unsuitable and suitable cells, and the neighbourhood distance value, used to identify spatially discrete patches of suitable habitat. For further details see Table 1.

Variation of demographic parameters

GRIP 2.0 generates unique replicate metapopulation models (one for each replicate landscape model) by varying parameters specified in RAMAS Metapop. A total of 23 demographic parameters specified in the whitebark pine metapopulation model were varied in the GSA (Table 1). These included stage-specific survival and fecundity rates, dispersal rates, dispersal survival, correlation of vital rates among populations, and catastrophe parameters. All demographic parameters were varied using sampling distributions and parameter ranges as specified in GRIP 1.0, a version of this freeware developed for spatial PVAs that are not based on spatially-explicit maps of habitat or habitat suitability (Curtis & Naujokaitis-Lewis, 2008). We did not vary the model of density dependence to retain the original model structure to the extent possible. Further details governing the sampling distributions, mean and measure of variation for each parameter varied in the GSA are described in Table 1 and Appendix S2. As with other parameters not varied in this version of GRIP 2.0, the code is extensively annotated and customizable. To the extent possible, we used information from the original model to select biologically relevant uncertainty estimates, and in the case where these were not specified we either applied a 10% coefficient of variation or sampled from a uniform distribution (Table 1). While this 10% value is somewhat arbitrary, we used the best available information to select realistic uncertainty estimates. Users are able to modify these ranges, and even sampling distributions to better reflect the research and management context of their study system.

Simulations and data analysis

Using the whitebark pine map as the original reference landscape, we created a total of 10,000 replicates, where refers to the total number of final model configurations that includes variation in landscape/habitat suitability parameters (Fig. 2), as well as demographic parameters. Each of the replicate metapopulation dynamics models consisted of 1,000 stochastic runs using a 100 year time period. For each replicate landscape file, we calculated landscape composition and configuration metrics using RAMAS GIS. The landscape metrics included number of patches, total amount of suitable habitat (an integrated measure of both habitat amount and suitability), mean patch area, edge to area ratio, and connectivity, calculated as the proportion of all population pairs that are linked through dispersal. A full description of the landscape metrics calculated by RAMAS GIS and associated mathematical formulae is available in Akçakaya & Root (2005).

We applied a boosted regression tree (BRT) to rank the relative influence of habitat-based measures of landscape pattern and demographic parameters on the binary response variable, conservation status. Conservation status was calculated based on the probability of extinction over a 100-year time period with 0.1 defined as the threshold for distinguishing metapopulations expected to be not at risk from those expected to be at risk of extinction. This benchmark corresponds to international criterion for listing species as Vulnerable (Criterion E; IUCN, 2001). We used the machine-learning method of BRTs based on our interest in understanding the relative importance of different parameter uncertainties on model outcomes and additional flexibilities of BRT for our analysis purposes (De’ath, 2007; Elith et al., 2008). GRIP 2.0 provides the means to generate and propagate uncertainties while decisions regarding how to analyze such outcomes are flexible, with BRT representing one option. The BRT functionality was not coded into GRIP 2.0 but was a supplemental statistical analysis performed once all replicate simulations were created and run using GRIP 2.0 and RAMAS.

We specified a binomial error structure and link function and used untransformed data in the BRT as it does not require data transformations (Elith et al., 2008) . We applied a tree complexity (tc) value of 2, which fits the BRT with up to two-way interactions. Learning rate (lr), which determines the contribution of each tree as it is added to the model, was specified at a value of 0.01, which was optimized to ensure a minimum of 1,000 trees were fit for the model (Elith et al., 2008). Prediction accuracy was measured using the percent deviance explained by the model, which measures the goodness of fit between the prediction and observed values. The relative influence of each predictor variable was assessed by calculating its contribution to reducing the overall model deviance of the BRT model. We identified important modeled interactions by quantifying the strength of pairwise interactions while keeping non-focal variables at their mean values. The BRT model, including relative influence of predictors, and evaluation and visualization of two-way interactions were performed in R (R Core Team, 2015) using the ‘gbm’ package (Ridgeway, 2015) and functions available in Elith et al. (2008).

Influential parameters on whitebark pine conservation status

The baseline model of the whitebark pine in the presence of blister rust predicts a dramatic decline of the metapopulation in Mt. Rainier National Park (Ettl & Cottone, 2004). The two-way interaction BRT model resulted in a high level of explanatory power. The model accounted for 90% of the mean total deviance (1 − mean residual deviance/mean total deviance; 1 − 0.014/0.137). For the 2-way interaction BRT model, the loss function was minimized at 4,250 trees, and the model was optimized with a learning rate of 0.01 and a bag fraction of 0.05. The five most important variables based on the BRT included total habitat amount (40.3%), survival class 4 (13.7%), survival class 3a (12.8%), catastrophe intensity (8.9%), and mean carrying capacity per patch (6.6%) (Table 2).

The partial dependence plots in Fig. 3 illustrate the relationship between the four most important variables and conservation status after accounting for the average effect of all other variables. These plots suggest that a conservation status of ‘Vulnerable’ is a function of low amounts of total habitat in the landscape, low survival rates of trees in classes 3a and 4, and stronger negative effects of catastrophes (i.e., blister rust).

Importance of interactions

The BRT analysis identified a number of strong interactions influencing whitebark pine conservation status. The strongest interaction was between total habitat amount and survival rate of trees in class 4 (Fig. 4A). Specifically, extinction risk increases more rapidly for lower values of survival and lower amounts of habitat in the landscape, but as survival increases above 0.8, less amount of habitat is required for the conservation status to remain ‘not-at-risk’. Based on the second most important interaction, the influence of catastrophe intensity on extinction probability is conditional on the total amount of habitat in the landscape (Fig. 4B). In other words, when the effect of catastrophes is weaker, less amount of habitat is needed for the predicted status to remain not-at-risk whereas when catastrophes have a larger effect more habitat is needed to buffer the risk status rank. Both of these interactions demonstrate the importance of habitat amount in mediating the negative effects of blister rust on whitebark pine extinction risk within this landscape.

Habitat thresholds

Based on qualitative assessments, the 3-dimensional plots revealed thresholds in habitat-based features. For the two strongest interactions (Fig. 4), across all survival rates and catastrophe intensity values, below total habitat amounts of 2 km² the predicted risk of extinction increases more steeply. Despite the strong interaction among these pairs of variables, this suggests the existence of thresholds in total habitat amount for the Mount Rainier whitebark pine metapopulation infected with blister rust. As the number of populations did not interact strongly with other demographic parameters or habitat features, the partial dependence plots may be interpreted directly: thresholds were evident for the number of populations where conserving >20 populations resulted in predicted conservation status of ‘not-at-risk’, based on international listing criteria.

Simulating alternative realizations of landscape and habitat structure

Landscape level experiments involving adequate replication are often difficult to implement due to issues of scale. Consequently, experimentation by simulation using landscape generators presents a viable alternative to develop a better understanding of the relationship between landscape pattern and process (Gardner & Urban, 2007), predict species response to landscape change (Tischendorf, 2001), and more generally to develop and test hypotheses. Because GRIP 2.0 modifies an existing reference landscape, certain landscape elements, such as size and shape of some patches, are retained in replicate simulations. This leads to more realistic simulated landscapes, unlike neutral or multi-fractal models (With & Crist, 1995). Although not applied to the whitebark pine model, GRIP 2.0 allows users to include a landscape mask, enabling habitat creation only in those regions outside of the mask. This aspect increases the functionality GRIP 2.0 and the realism of its outputs; a mask may contain intractable barriers to habitat creation, such as roads, representing landscape constraints that many analysts and managers must contend with.

Replicate landscapes produced by GRIP 2.0 represented a wide range of landscape structural variation allowing an evaluation of the influence of a broad range of structural attributes on extinction probability of whitebark pine, and by the same token, a broad range of scenarios to explore for conservation and management planning. By using Monte Carlo simulations, landscape features modified included the amount of habitat, patch sizes, or the degree and directionality of spatial autocorrelation in replicate landscapes. Such flexibility is important given simulation studies and our results that have identified thresholds related to levels of habitat within a landscape below which species viability rapidly declines (Swift & Hannon, 2010). The GRIP 2.0 code is highly annotated allowing users to understand, scrutinize, and modify the code to reflect a particular species’ biology, sampling distributions, or landscape dynamics.

The original whitebark pine model did not model consequences of climate change on future projections. Thus our implementation of GRIP 2.0 to the whitebark pine model does not address uncertainty in spatio-temporal projections of habitat suitability under future climate change associated with selection of General Circulation Models, for example. However, this is not an inherent limitation of the tool as GRIP 2.0 can be modified to integrate this additional source of uncertainty as in Naujokaitis-Lewis et al. (2013). Additionally, while we did not explicitly address uncertainty in the habitat-demographic relationship, our GSA approach implicitly addresses this source of variation by varying population-specific initial abundances and carrying capacities. Should models incorporate temporal trends in vital rates associated with estimated habitat-demographic relationships, GRIP 2.0 is customizable to reflect this potential source of variation.

The relative influence of parameters and their interactions

By varying multiple parameters simultaneously using probability distribution functions to represent the plausible range of parameter values, GRIP 2.0 performs a global sensitivity analysis, in which parameters are varied concurrently over the plausible range of parameter space (Saltelli et al., 2006). Varying habitat features associated with habitat suitability maps, such as habitat amount, habitat suitability, and number of patches enabled a comparison of the relative influence of habitat-based attributes on whitebark pine meta-population dynamics relative to demographic parameters. Outcomes of our global sensitivity analyses indicated that both demographic and habitat factors influence predictions of whitebark pine persistence. Our GSA identified different influential parameters from the original study, which did not vary habitat-based factors (Ettl & Cottone, 2004). Based on the predicted trajectories of individual whitebark pine subpopulations over time, Ettl & Cottone (2004) concluded that size and distribution of subpopulations influenced the trajectory of individual populations. The authors also concluded that the model was more sensitive to changes in the vital rates of healthy trees than infected trees, and less sensitive to changes in the vital rates of mature trees than of younger trees, and more sensitive to the frequency of invasion than by the effect of blister rust on vital rates. Without applying a GSA, however, the authors were not able to rank the relative influences of these factors on the dynamics of the meta-population as a whole, or identify key interactions among parameters and thresholds that could be used to inform management decisions. By contrast, the results of our GSA indicated that model predictions were more sensitive to older stage classes that were also infected by blister rust disease. We also identified certain habitat variables as highly influential, which broadens the range of information that can be used to develop effective management actions. Our simulations were standardized to 100 years (as per criterion E for IUCN species assessments), but with poor recruitment/low survival, this metapopulation may be at greater risk over the longer term (Ettl & Cottone, 2004) .

The relatively high rank of certain habitat factors, such as total amount of habitat, highlights the importance of spatial habitat-based landscape features as important drivers of whitebark pine persistence as a way to mediate the negative consequences of blister rust disease. Spread of the fungal blister rust pathogens is complex; blister rust has a five-stage life cycle and requires two hosts, whitebark pine (or any other five-needled pines) and any species in the genus Ribes. Predicting the spread of the blister rust through the distribution of whitebark pine individuals and populations is difficult as the pathogen spreads from Ribes hosts to whitebark pine via wind-borne spores, and not from tree to tree (McDonald & Hoff, 2001). Although the original PVA model does not explicitly model transmission through Ribes, the outcomes of our GSA indicate that managing habitat parameters could help offset declines in abundance related to blister rust. Given that proximity to Ribes species contributes to increased rates of blister rust infection (Smith et al., 2011), future research that explicitly explores the spatial context of Ribes’ role in disease transmission dynamics would improve our understanding of the role of multiple interacting parameters that influence whitebark pine persistence.

Thresholds and conservation implications

Theoretical models have been used to generate hypotheses, such as the nonlinear threshold hypothesis, which predicts that species exhibit threshold responses due to the increasing influence of fragmentation below a certain amount of habitat (Andrén, 1994). Although theoretical models predict that below critical threshold points, increased fragmentation results in lowered colonization success and increased extinction probability (With & King, 1999), empirical validation of such responses are uncommon (Fahrig, 2003). Our results corroborate the occurrence of thresholds in response to declines in both habitat amount and quality. Visual assessment of bivariate plots of probability of extinction as a function of landscape variables indicated the presence of thresholds as evidenced by strong negative exponential decay curves (results not shown). Although we did not quantitatively derive threshold points for whitebark pine, these roughly translate to 20 populations, a mean connected distance of approximately 5 km, an average habitat suitability of 2.3 per patch, and a mean patch size of ∼0.35 km². Further insights into the behavior and occurrence of thresholds in real settings may help identify conservation strategies that specify minimum patch size targets, or provide guidance related to landscape level measures, such as habitat amount (Betts, Forbes & Diamond, 2007).