Beyond single-species models: leveraging multispecies forecasts to navigate the dynamics of ecological predictability

Rodent capture data

Our data come from the Portal Project, a long-term monitoring study of a desert rodent community (Brown, 1998; Ernest et al., 2020) that has been undergoing active forecasting since 2016 (White et al., 2019). The Portal Project is based in the Chihuahuan Desert near Portal, Arizona. The sampling design includes 24 experimental plots (50 m × 50 m), each containing a grid of 49 baited traps (Brown, 1998; Ernest et al., 2020). The design uses three experimental treatments. In control plots (N = 10), holes in the fence are large enough to allow free access for all rodents. Full rodent removal plots (N = 6) have fences with no holes. Kangaroo rat exclosures (N = 8) have fences with holes to allow passage of all rodents except kangaroo rats (Dipodomys genus). Investigators close holes during trapping to ensure all captured rodents are residents. Trapping follows the lunar monthly cycle, but weather and other disruptions result in missing observations (~5% on average; Dumandan, Yenni & Ernest, 2023).

The Portal dataset exhibits many of the complexities that confront population forecasting. These include observation errors due to imperfect detection, missing samples due to weather or other issues (e.g., global pandemics), and discrete counts of captured individuals for many species (20 rodent species) that include large numbers of zeros, multi-species dependencies and upper bounds set by the number of traps (Ernest et al., 2020). Environmental drivers, including temperature and measures of primary production, exhibit lagged and nonlinear impacts on rodent breeding, activity rates, and population dynamics (Brown & Ernest, 2002). Moreover, the rodent species at Portal are known to compete for resources in complex ways, and these biotic interactions are postulated to have important consequences for variation in population dynamics. In other words, the Portal Project exhibits many of the complexities that make the ecological forecasting of species populations particularly difficult, making it an ideal real-world test case for exploring whether multi-species models can provide better near-term predictions than single species models.

Open-source software exists to access the Portal Project data (Christensen et al., 2019; Simonis et al., 2022). We used the portalr R package to extract trapping records from the Portal data (version 3.134.0; downloaded October 2022; https://doi.org/10.5281/zenodo.7255488). Our study focused on rodent captures from the long-term control plots for the period December 1996–August 2022. The data has records for 20 rodent species, but some are rarely captured. We excluded species if they were observed in <10% of trapping sessions. We did this to focus inferences on species with the most influence on community dynamics. Each temporal observation was a vector of total captures on long-term control plots for the nine remaining species at a given sampling time (Fig. 1).

Covariate measurements

Rodent populations at Portal, and the associated number of captures recorded during sampling, depend on environmental conditions that reflect resource availability and seasonal breeding signals. We therefore modelled species’ responses to environmental variation using the Normalized Difference Vegetation Index (NDVI, representing a signal of resource availability) and minimum temperature (which strongly reflects variation in seasonal breeding behaviours; Dumandan et al., 2024) as covariates. Hourly air temperature (°C) is recorded by an automated weather station on site, which we used to calculate a daily minimum, while Landsat images are used to calculate NDVI for an area 1,000 m in radius around the geographical centre of the Portal Project site (accessed on a 2-week basis from the US Geological Service Earth Resources Observation and Science Center; https://www.usgs.gov/centers/eros). Measurements for both covariates (daily minimum temperature and bimonthly NDVI) were then converted to monthly averages. We extracted covariate data from one year before the start of captures (from January 1995) so we could calculate lagged and moving average versions. See Ernest et al. (2020) for further details of environmental measurements, including their spatial and temporal resolutions.

Model description

There were several aspects of the data we needed to consider when designing our model. Total rodent captures showed both short- and long-term fluctuations (Fig. S1). Captures for individual species also undulated over multi-annual cycles and were positively autocorrelated at lags up to 20 months (Figs. S2 and S3). To test whether multi-species information improves model performance, we needed to model these dynamics using a multivariate dependence structure. Second, we needed to leverage community information to estimate each species’ time-delayed response to variation in vegetation and temperature. Because species’ responses to environmental change in this system are expected to be delayed and nonlinear (Brown & Ernest, 2002), we used splines to model these responses. Rodent captures were modelled as $Poisson$ observations of a latent state model that was composed of a hierarchical generalized additive model (GAM) component (to capture shared, nonlinear environmental responses) and a multivariate dynamic vector autoregressive component to capture multi-species dependence. This State-Space model was designed to explicitly address several of the key challenges that plague population dynamics forecasting by leveraging (1) a discrete observation model to appropriately deal with count-valued time series with many zeros; (2) a latent real-valued state model with a separate error component to deal with imperfect detection; (3) hierarchical effects and penalized splines to model species’ nonlinear environmental responses and (4) a vector autoregression (VAR) component to capture both lagged and contemporaneous dependencies among species’ estimated latent states. The full description for this model, which we abbreviate to GAM-VAR, is shown in Fig. 2. The GAM component of the model consisted of hierarchical NDVI and minimum temperature effects. The structural forms of these functions were informed by theory and exploration of covariate time series (shown in Figs. S4, S5). We used a 12-month moving average of NDVI ( $NDV{I}_{MA12}$) because we expected rodent populations to respond gradually to vegetation change, whereby productive years (represented as timepoints with relatively high NDVI values) will likely result in delayed population increases as rodents are able to cache more seeds and make use of the higher availability of shelter spaces (Brown, 1998; Brown & Ernest, 2002). Our model assumed linear effects of $NDV{I}_{MA12}$, equivalent to a hierarchical slopes model. The partial pooling properties of this model allowed us to regularize weakly informed slopes toward a community average. Responses to temperature were estimated using a hierarchical distributed lag model in which nonlinear effects of minimum temperature varied smoothly with increasing lag. These effects were constructed as tensor products of four cubic basis functions for lag and three thin plate basis functions for minimum temperature. To allow our model to capitalize on multi-species learning, we included a shared community-level response $f_{global}\left(Mintemp,lag\right)$ and species-level deviation responses $f_{species\left[i\right]}\left(Mintemp,lag\right)$. The sum of these effects allowed each species to show a different temperature response from the wider community, but only if there was information in the data to support such a deviation. We used lags of up to 6 months in the past because this time window is likely sufficient to capture rodent behavioral and breeding responses to seasonal temperature change (Karunarathna, Wells & Clark, 2024).

A VAR of order 1 captured lagged multi-species dependence, where $A$ was a nonsymmetric 9 × 9 matrix of autoregressive coefficients (Fig. 2). Diagonal entries of $A$ described density-dependence, or the effect of a species’ dynamic process (at time $t$) on its own lagged values (at $t-1$). Off-diagonals represented cross-dependencies that could provide useful biological insights into interspecific interactions. For example, the entry in $A$[2,3] described the effect of species $3$’s dynamic state at time $t-1$ on the current state estimate for species $2$ (at time $t$). To encourage stability and prevent forecast variance from increasing indefinitely, we enforced stationarity following methods described in Heaps (2023). Briefly, a multistep process was used to map the constrained $A$ matrix to unconstrained partial autocorrelations $P$, which allows for meaningful prior elicitation about the structure of A while enhancing computational efficiency. Process errors were allowed to be contemporaneously dependent to capture any unmodelled correlations in species’ latent states. Priors for all model components are shown in Fig. 2 and described in detail in the accompanying R code.

Evaluating whether multi-species dependencies improve prediction performance

We formally tested whether learning from multiple species improved our model’s predictions using prediction-based model comparisons. To do so, we estimated a series of benchmark models that acted as natural simplifications of the GAM-VAR by eliminating multi-species components in a stepwise manner. The first benchmark model used the same hierarchical GAM linear predictor as the GAM-VAR but replaced the multi-species VAR(1) dynamics with an AR(1) process. This model (called GAM-AR in subsequent sections) eliminated the covariances and temporal cross-dependencies among species’ latent states, allowing us to ask whether the multivariate dynamic component was supported for improving model fit. Next, we further simplified the GAM-AR by removing the hierarchical environmental response functions from the linear predictor. This forced the model to learn environmental responses for each species without using information from other species in the data (GAM-AR no pooling). The final benchmark, referred to as AR, also used independent AR(1) states but removed the GAM component entirely. Because this model only learned from past observations, comparisons against it helped us understand how covariate responses impacted our models’ predictions and inferences. Each candidate model was trained on all observations (through August 2022, N = 319 timepoints). Models were then compared using Pareto-smoothed importance sampling leave-one-out cross-validation (PSIS-LOO), a method that reweights posterior draws to estimate leave-one-out pointwise prediction accuracy using estimated log predictive density (ELPD) values (Vehtari, Gelman & Gabry, 2017).

To adequately evaluate competing forecast models, it is also necessary to perform out-of-sample validation (Clark et al., 2022; Harris, Taylor & White, 2018; Lewis et al., 2022). This is particularly important because LOO-CV is designed to ask how models would generalize to new observations within the training window. This metric does not evaluate a time series model’s ability to forecast, as information from future timepoints is used to influence predictions for previous time points. To evaluate forecasts in a way that respected the temporal nature of our forecasting exercise, we used exact leave-future-out cross-validation in an iterative expanding window framework. Models were re-trained on the first 273 time points (~22 years), with the subsequent 12 time points (through November 2019; selected to avoid a large sampling gap due to the COVID-19 pandemic) used to evaluate forecasts. This allowed us to gauge how models might perform in a forecast scenario, but it only provided a single comparison. To further scrutinize models, we retrained models on the first 75, 115, 154, 194, and 233 observations, and evaluated the subsequent 12 observations in each cross-validation fold. All forecast comparisons used an evenly weighted combination of two proper multivariate scoring rules. We chose the variogram score, which penalizes distributions that do not adequately capture correlations in test observations, and the energy score, which ignores correlations but penalizes forecasts if they are not well-calibrated (Scheuerer & Hamill, 2015).

Estimation

We estimated posterior distributions with Hamiltonian Monte Carlo in Stan (Carpenter et al., 2017; Stan Development Team, 2022), specifically the cmdstanr interface (Gabry & Češnovar, 2021). Stan’s algorithms provide state-of-the-art diagnostics for probabilistic models (Betancourt, 2017). For example, Hamiltonian Markov chains diverged when attempting to estimate minimum temperature deviations for some species in the GAM-VAR. Our data were not informative enough to learn how, or even if, these species responded to temperature change in ways that differed from the community response. Stan’s diagnostics guided us to a model that could be reliably estimated, which included species’ level distributed lag functions for the four most frequently captured species (D. ordii, D. merriami, Onychomys torridus and C. penicillatus). Posteriors were processed in R 4.3.1 (R Core Team, 2023) with the mvgam R package (Clark & Wells, 2023). Traceplots, rank normalized split-R ̂ and effective sample sizes interrogated convergence of four parallel chains. Each chain was run for 500 warmup and 1,600 sampling iterations. R code to replicate all analyses and produce figures is included in the Supplemental Materials and is permanently archived on Zenodo (https://doi.org/10.5281/zenodo.14607006).

Modeling relationships among species improves prediction performance

Our data included total captures for nine rodent species over 319 time points. All models showed adequate convergence and posterior exploration, and randomized quantile residuals showed no obvious evidence of unmodelled temporal or systematic variation (Figs. S6, S7). However, in-sample performances differed among models, with models that leveraged multi-species information producing higher ELPD scores compared to simpler models (Table 1). This was the case for all stepwise comparisons apart from one: although the GAM-AR, which used partial pooling to learn species’ environmental responses, was favoured over the simpler GAM-AR no pooling, overlapping ELPD standard errors suggested there was still large uncertainty about the magnitude of this difference (Table 1).

We also found that forecast performance differed among models, with more complex multi-species models again tending to score higher for forecast performance than simpler models. Forecasts from the multi-species GAM-VAR were the most accurate when considering all validation points in aggregate and for 4/6 cross-validation folds (Fig. 3; Fig. S8). The GAM-AR and GAM-AR no pooling models gave similar predictions and effectively tied for second in forecast performance, giving the most accurate forecasts in 2/6 cross-validation folds (Fig. 3). The simplest AR model gave the worst forecasts.

The multi-species GAM-VAR model estimated large, positive autoregressive coefficients for most species (diagonal entries in Fig. S9). It also relied strongly on information from multiple species by estimating many non-zero cross-dependence effects (off-diagonal entries in Fig. S9) and process error correlations (Fig. S10), which provided structure that the model leveraged to accurately simulate historical dynamics. The model recovered multiple notable transitions observed in the time-series including a major shift in community composition around the year 2000 following the establishment of Bailey’s pocket mouse C. baileyi, and a second restructuring that happened following a drought in 2008–09 (Fig. S11). It was these multi-species effects that enabled the GAM-VAR to produce more accurate forecasts compared to the benchmarks. For example, Ord’s kangaroo rat (Dipodomys ordii) and silky pocket mouse (Perognathus flavus) had negative cross-dependencies in the GAM-VAR, providing additional information that could be used to make more precise predictions (Fig. 4). The benchmarks, which ignored this structure, produced smoother, less synchronous trends and wider uncertainties (Fig. S12). In the following sections, we use simulations to briefly interpret each of the multi-species effects that allowed the GAM-VAR to outperform simpler models.

Modeling relationships among species provides unique insights into community dynamics

Our cross-validation metrics strongly favoured the GAM-VAR and suggested that the multivariate dynamic component was a particularly important driver of increased performance. Estimates of process error were larger for the benchmarks than the GAM-VAR for nearly all species (Fig. S13), suggesting this model used additional information from multi-species cross-dependencies to produce better predictions. But interpreting this cross-dependence is difficult because VAR effects provide only a marginal view into the complex network of conditional associations. We therefore used impulse response functions (Lütkepohl, 1990) to better understand the model. These functions involve simulating an ‘impulse’ in captures for one species and then evaluating how predicted captures for other species changed over the next 6 months (Fig. 5). Following a simulated impulse of three extra captures for Merriam’s kangaroo rat (D. merriami), the model expected some initial increases (due to the correlated process errors) followed by declines in captures for most of the other species (Fig. 5). The shapes of these declines varied by species. Captures for the two pocket mouse species (C. baileyi and C. penicillatus) showed more immediate declines, while the two grasshopper mouse species (O. leucogaster and O. torridus) declined more gradually (Fig. 5). In contrast, the other kangaroo rat species (D. ordii) was expected to increase following a D. merriami pulse (Fig. 5). Different effects were expected when changing the focal species (Fig. S14).

Positive NDVI associations and hierarchical minimum temperature effects

We found broad support for positive associations with $NDV{I}_{MA12}$ (Fig. 6). Conditional simulations, which asked how rodents might respond if moved from a relatively dry/brown vegetation state to a relatively moist/green vegetation state, gave higher probability to increased captures in the moist/green scenario for all species. But uncertainties about this effect varied among species. Greatest increases were expected for Ord’s kangaroo rat (D. ordii), Western harvest mouse (R. megalotis) and cactus mouse (Peromyscus eremicus). The model was less confident about the direction of effect for Northern grasshopper mouse (O. leucogaster) and for one of the most dominant species in the study, Meriam’s kangaroo rat (D. merriami). For these species, the model showed a ~70% chance of increasing abundance in the higher $NDV{I}_{MA12}$ scenario (Fig. 6). While primary conclusions were generally similar when using the GAM-AR no pooling model, which did not leverage multi-species learning, the estimates of these contrasts were much more variable.

Interpreting minimum temperature distributed lag effects also required simulation. We visualized 1,000 simulated functions for each species using temperatures from the year 1997 (Fig. S15). There was large uncertainty in function shapes for all species except the desert pocket mouse (C. penicillatus). Captures for this species were expected to increase from May to October and decrease sharply in winter. For seven of the other eight species, the model generally expected more captures in spring (March–May) and fewer in late summer/autumn (July–October). But the shapes of these responses varied. The two kangaroo rats (D. merriami and D. ordii) had smoother shapes that decreased gradually from mid-summer to mid-winter. But the model expected D. ordii captures to peak slightly later (May as opposed to March for D. merriami). The Southern grasshopper mouse (O. torridus) was the only species that was expected to show higher captures in late autumn/early winter (Fig. S16). The five species that relied solely on the global distributed lag minimum temperature function (O. leucogaster, C. baileyi, P. eremicus, P. flavus and R. megalotis) were expected to show tighter spring peaks (higher captures in April and May) and autumn troughs (fewer captures in August and September). When simulating from the GAM-AR no pooling model, the lack of multi-species learning was immediately obvious. There was not enough information to learn nonlinear distributed lag functions for these five species, with the model instead estimating flat functions centred on zero for all five species (Fig. S16).

Leveraging relationships between species for ecological forecasting

Interactions and dependencies among multiple species are hypothesized to play pivotal roles in the assembly of ecological communities and for broader ecosystem functions (Dobzhansky, 1950; Fecchio et al., 2019; Mayfield & Stouffer, 2017; Mutshinda, O’Hara & Woiwod, 2009). This study shows why models that target multi-species effects in both their environmental responses and their biotic dependencies should also be strongly considered when studying community dynamics. Our approach to constructing hierarchical dynamic GAMs and evaluating forecasts using multivariate proper scoring rules offers a way to quantitatively assess multi-species forecasts and scrutinize their value in real-world ecological forecasting applications. We also demonstrate how inferences from these models provide deeper insights into why they may or may not perform better. For example, the GAM-VAR’s process error estimates were smaller than those from the benchmarks because it used more information from the data. By learning about the relationships between species the model could better capture both shared responses to environmental factors (e.g., a wet year in the desert is good for most species) and direct temporal effects (e.g., dynamic competition for seeds). These relationships between species can allow forecasts for less common species to borrow strength from more common species, yielding better hindcasts and forecasts compared to simpler single-species models. But like other multivariate autoregressive models (Hannaford et al., 2023; Holmes, Ward & Scheuerell, 2014; Ives et al., 2003) the VAR parameters of the GAM-VAR should not be interpreted as a species interaction matrix, because these relationships can result from multiple sources (i.e., shared environmental responses, shared biotic responses and/or direct interactions). While the parameters are not interpretable as direct interactions, this approach does make it possible to gain a more detailed understanding of population dynamics. Conducting simulations from this type of model allows exploring which species have the strongest cascading effects, what changes might we expect if management increases or decreases abundance for target species, and how these effects relate to regime transitions.

Our hierarchical modelling approach also makes it possible to partition variance into contributions from observation error, environmental responses, and multi-species dependence to guide future efforts to improve ecological forecasting. In our study, forecasts were dominated by uncertainty in the dynamic process model, but using a Vector Autoregressive process allowed us to dissect this uncertainty in meaningful ways (Ives et al., 2003; Lütkepohl, 1990). Simulated responses to sudden impulses in captures were often delayed and nonlinear. Despite the restriction to a VAR of lag of 1 month, these responses resulted in cascading changes that lasted up to 6 months. Our model’s ability to simulate and dissect community change in this way offers a useful avenue for ecologists to better understand, and expand on, theoretical predictions from both classical and more recent empirical studies that have described strong interspecific interactions in ecological systems (Dumandan et al., 2024; Ebersole, 1977; Mayfield & Stouffer, 2017; Volterra, 1928).

Learning hierarchical nonlinear effects from community data

Our model captured linear, nonlinear, and lagged responses to environmental and climatic covariates that were informed by data from all species at once. We found positive linear associations between capture rates and a 12-month moving average of NDVI. This response was expected because the rodents at Portal depend on plants for food and other resources (Brown & Ernest, 2002; Ernest et al., 2020) and NDVI measures vegetation greenness in the landscape. But within this overarching community pattern there were interesting patterns of variation in these responses among species. The strongest positive association was shown by Ord’s kangaroo rat (D. ordii), a species that field evidence suggests consumes and harvests grasses (Kerley, Whitford & Kay, 1997). In contrast, Merriam’s kangaroo rat (D. merriami) showed weaker associations with NDVI. This species has been predicted to increase in relative abundance in the study region with more severe droughts, in part due to a preference for more open foraging habitat with less vegetation (Cárdenas et al., 2021).

Distributed lag functions determine the best combination of lags for environmental covariates but are not commonly used in ecology (but see Karunarathna, Wells & Clark, 2024; Ogle et al., 2015; Wells et al., 2016). Our study shows how these effects can be learned hierarchically and provides useful insights into delayed responses to temperature change for rodent species at Portal. Most species showed higher captures when minimum temperatures were low 3–4 months prior, suggesting increases begin during mid to late spring when resources such as seeds become available. But others, such as Merriam’s kangaroo rat and Southern grasshopper mouse, showed increases during cooler months in autumn and winter. Asynchronous phenology, where species show different reproductive timing, is sometimes expected in competitive communities (Carter & Rudolf, 2022). Analysis of individual reproductive status in different biotic contexts suggests that some species shift their reproductive timing in the presence of strong competitors in the Portal system (Dumandan, Yenni & Ernest, 2023). Do these competitive forces play a role in seasonal capture variation over the long-term? Comparing temperature responses on control vs. experimental plots would be one interesting way to tackle this question.

Interestingly, despite the relatively large number of observations our data contained for each species, estimates of environmental responses were still more precise and informative when using hierarchical models (which use partial pooling) as opposed to a no-pooling model that only considers species’ effects in isolation. While hierarchical intercepts and slopes are commonly used in ecological models, there has been less emphasis on hierarchical nonlinear functions (but see Pedersen et al., 2019). Open access to new software that makes it easy to construct and estimate these types of functions, such as the mvgam R package that we used here (Clark & Wells, 2023), should improve their uptake in ecological forecasting exercises.

But despite the power of hierarchical environmental effects to improve predictions, we cannot interpret environmental response estimates as directly causal for several reasons. First, we know NDVI is not a perfect measure of changes in seed production. Second, it is likely that changes to NDVI and minimum temperature are both related to other unmeasured environmental drivers that may also influence rodent abundance. Major storms, the El Niño Southern Oscillation and other factors that influence moisture levels can all influence temperature and vegetation change (Sun & Kafatos, 2007). These other drivers could act as unmeasured confounds, biasing estimates in a causal inference framework (McElreath, 2020).

Future directions

Two additional steps would be useful to fully assess the value of multi-species models for ecological forecasting, both in this system and more broadly as an ecological application. First, a more diverse suite of candidate models could be estimated to determine how forecasts could be combined into an ensemble to provide the best predictions in situations where prediction accuracy is the primary goal (Clark et al., 2022; Powell-Romero et al., 2023). This could be especially useful for detecting changes in the system. For example, GAM-VAR gave better forecasts in most cross-validation tests, but its performance was slightly worse than the simpler GAM-AR when the training window stopped just prior to a major restructuring of rodent abundances that was taking place in response to a drought. Second, determining which models are best for true forecasting requires evaluating forecasts in the presence of uncertainty in future covariate values. In this study we were hindcasting and therefore used the actual observed environmental measurements for the period reserved for model evaluation. Fortunately, the system is undergoing active forecasting involving a suite of simpler models and leveraging actual forecasts for environmental covariates (Simonis et al., 2022; White et al., 2019). A natural next step for this work is to compare the performance of the GAM-VAR model to simpler models both using hindcasting with observed covariates and when making true forecasts relying on predictions instead of observations for NDVI and minimum temperature.

The Portal Project also provides a unique opportunity to disentangle the combined influence of shared environmental responses and direct species interactions in driving observed relationships between species. The site includes a long-term experimental manipulation where kangaroo rats (Dipodomys species) are excluded from some plots. Recent work shows that single species forecasting models for C. baileyi do not transfer well between the control plots and this experimental manipulation, likely due to the different competitive environment experienced in the absence of the behaviorially dominant kangaroo rats (Dumandan, Yenni & Ernest, 2023). Multi-species models like the GAM-VAR have the potential to transfer better in situations where one or more species are removed from the system by accurately predicting the response of the other species to this removal. Therefore, a key next step in evaluating the potential strengths of our models is to determine if they can more effectively transfer to make accurate predictions on the plots with the experimentally manipulated species composition. More broadly, we hope that the ability to estimate multi-species dependence and species-level variation in nonlinear environmental responses result in more accurate forecasts, inspire new questions, and lead to an improved understanding of the factors that govern ecological community dynamics.