Matching the forecast horizon with the relevant ecological processes

Most models used to generate ecological forecasts take either a time-series approach, based on long-term data from one location, or a space-for-time approach, based on data describing spatial patterns across environmental gradients. Here we consider how the forecast horizon determines whether the most accurate predictions come from the time-series approach, the space-for-time approach, or a combination of the two. We use two simulation case studies to show that forecasts for short and long-time scales need to focus on different ecological processes, which are reflected in different kinds of data. In the short-term, dynamics reflect initial conditions and fast processes such as birth and death, and the phenomenological time-series approach makes the best predictions. In the long-term, dynamics reflect the additional influence of slower processes such as evolutionary and ecological selection, colonization and extinction, which the space-for-time approach can effectively capture. At intermediate time-scales, a weighted average of the two approaches shows promise. However, making this weighted model operational will require new research to predict the rate at which slow processes begin to influence dynamics.


Introduction
1. For short-term forecasts, phenomenological time-series approaches are hard to beat, whereas 84 longer-term forecasts require accounting for the influence of slow processes such as evolu-85 tionary and ecological selection as well as dispersal. Fast processes, such as births, deaths, and individual growth, operate at all time scales, but are the exclusive drivers of the short-term dynamics captured in most time series datasets. Slower processes, such as evolutionary selection on genotype frequencies, ecological selection on species abundances, and colonization and extinction, interact with fast processes to drive dynamics over the long-term. The influence of these slow processes is seen in very long time series, or in spatial gradients. Understanding dynamics at intermediate time scales requires integrating information from spatial and temporal data sources. We propose a model weighting approach; mechanistic spatiotemporal modeling is another alternative. The time scales shown here were chosen with vascular plants in mind, but the same concepts would apply for much shorter-lived organisms but at shorter time scales.

3.
A key challenge for future research is determining the rate at which slow processes begin to 92 influence dynamics.

93
Modeling approach 94 In two case studies, we simulated the effects of an increase in temperature on simple systems 95 with known dynamics. The "truth" is represented by a model that is mechanistic for at least one 96 important process, but we treat the model as unknown when analyzing the data and we assumed 97 that perfectly recovering this model would not be possible in practice. We began each simulation 98 under stationary temperature, allowing the system to equilibrate; we call this the baseline phase. 99 We then increased temperature progressively over a period of time, followed by a second period 100 of stationary, now elevated, temperature. The objective was to forecast the response of the system 101 to the temperature increase based on data gathered during the baseline period. 102 correlated interannual variation in an ecological response with interannual variation in the envi-105 ronment at just one site. The other model represents the space-for-time substitution approach, 106 which we call the "spatial approach" for brevity. We correlated the mean environment with the 107 mean of an ecological state or rate across many sites. We compared forecasts from both models 108 to the simulated dynamics to determine how well the two approaches performed at different 109 forecast horizons. We also assessed the potential for combining the information available in tem-110 poral and spatial patterns by using a weighted average of the forecasts from the temporal and 111 spatial approaches optimized to best match the (simulated) observations. We then studied how 112 the optimal model weights changed over time. We expected the temporal model to best predict 113 short-term dynamics, the spatial model to best predict long-term dynamics, while the weighted 114 model would show potential to provide the best forecasts at transitional, intermediate time scales.

115
The three statistical models are described in Supporting Information (Appendix A) and all code 116 for both case studies is available at Github (https://github.com/pbadler/space-time-forecast). were produced and reaches all other sites with equal probability. We provide a more detailed 134 description of the model in SI Appendix B.

135
We first simulated a baseline period with variable but stationary temperature, followed by 136 a period of rapid temperature increase, and then a final period of stationary temperature. In-137 terannual variation in temperature is the same at all sites, but mean temperature varies among 138 sites. All sites experienced the same absolute increase in mean temperature. We focused on the 139 biomass dynamics of one focal species that dominated the central site during the baseline period. During the baseline period there were strong spatial patterns across the mean temperature 141 gradient. Individual species, including our focal species, showed classic, unimodal "Whittaker" 142 patterns of abundances across the gradient ( Fig. 2A). These spatial patterns are the basis for 143 our "spatial model" of the temperature-biomass relationship for our focal species ( Fig. 2A).

144
In contrast to the strong spatial patterns, population and community responses to interannual 145 variation in temperature within sites were weak. At our focal site in the center of the gradient, the 146 biomass of the focal species was quite insensitive to interannnual variation in temperature, but 147 showed strong temporal autocorrelation (Fig. 2B). Our "temporal model" estimates this weak, 148 linear temperature effect, along with the strong lag effect of biomass in the previous year. 149 We fit both a temporal and a spatial statistical model to forecast the effect of a temperature 150 increase (Fig. 3A) on the focal species' biomass at one location in the center of the tempera-151 ture gradient. The predictions from the spatial and temporal models contrasted markedly, with 152 the temporal model predicting a large increase in biomass and the spatial model predicting a 153 decrease. Initially, the simulated abundances followed the increase predicted by the temporal 154 model, but as faster-growing species colonized and increased in abundance at the focal site, the 155 biomass of the focal species decreased, eventually falling below its baseline level (Fig. 3B).

156
To combine the temporal and spatial model into a single forecast, we fit a weighting parame-on K (see SI Appendix A for a full description of the approach). The weighted model accurately predicts the simulated dynamics across the full forecast horizon (Fig. 3B). It also shows that the 162 most rapid shifts in the model weights occurred during the period when warm-adapted, faster 163 growing species were increasing most rapidly in abundance (Fig. 3C). However, the reason the 164 weighted models works so well is that the weights were determined by fitting directly to the 165 data. Unlike our spatial and temporal model forecasts, we did not generate out-of-sample pre-166 dictions from the weighted model; it merely provides a convenient way to quantify how rapidly 167 dynamics shift from being dominated by interannual variation captured in the temporal model  The compositional turnover affecting our focal species also influences total biomass, linking 172 community and ecosystem dynamics. We repeated our focal species analysis for total commu-173 nity biomass, and the results were similar: the temporal model initially made the best forecasts 174 immediately following the onset of the temperature increase, but as the identity and abundances 175 of species at the study site changed, the model weights rapidly shifted to the spatial model (SI   To create "spatial data", we simulated the equilibrium density of each genotype under differ-191 ent mean temperatures. The pattern of equilibrium densities across a gradient in mean annual 192 temperature defines our spatial forecast: cold sites will be dominated by the cold-adapted ho-193 mozygous genotype, warm sites will be dominated by the heat-adapted homozygous genotype, 194 and intermediate sites will be dominated by the heterozygous genotype (Fig. 4B). The full model 195 description is provided in SI Appendix C.  The spatial pattern of individual genotypes (colors) and total population abundance (black) at sites arrayed across a gradient of mean annual temperature. The dashed line shows predictions from an empirical "spatial model," a linear regression that describes mean population size as a function of mean temperature. (C) The relationship between annual temperature and per capita growth rate at a location with a mean temperature that favors the cold-adapted genotype. Colors show population size (the green to brown gradient depicting low to high population density), which influences the population growth rate through density dependence.
The spatial pattern shown in Fig. 4B is the outcome of steady-state conditions. But at any one 197 site, the population's short-term response to temperature will be determined by the dominant 198 genotype's reaction norm (Fig. 4A). For example, at a cold site dominated by the cold-adapted homozygous genotype, a warmer than average year would cause a decrease in population size 200 due to decreases in fecundity (blue line in Fig. 4A), even though the heat-adapted homozygote 201 might perform optimally at that temperature. However, if warmer than normal conditions persist 202 for many years, then genotype frequencies should shift, and the heat-adapted homozygote will 203 compensate for the decreases of the cold-adapted genotype.

204
To demonstrate these dynamics, we simulated a diploid annual plant population at a colder 205 than average site. During the baseline period, the population is dominated by the cold-adapted 206 genotype. We used the simulated data from this baseline period to fit an empirical model that 207 assumes no knowledge of the underlying eco-evolutionary process. This empirical temporal 208 model (Appendix A) predicts population growth rate as a function of annual temperature and 209 population size (Fig. 4C). We then imposed a period of warming, followed by a final period of

218
As in the community turnover example, we also fit a weighted average of the spatial and 219 temporal model, with the weights changing over time. This weighted model initially reflected 220 the temporal model (decrease from t = 500 to t = 600), but then rapidly transitioned to reflect the 221 spatial model (t ≥ 700). The rapid transition in the weighting term, ω, occurred during the period 222 of most rapid change in genotype frequencies ( Fig. S-3). The weighted model's predictions look 223 impressively accurate, but, as in the community turnover example, that is because we used the 224 full, simulated time series to fit the weighting term. A true forecast would require an independent 225 method to predict how the model weights shift over time.

227
Ecological forecasts are typically made using either a space-for-time substitution approach based 228 on models fit to spatial data or using dynamic models fit to time-series data. Our results demon-229 strate that these two approaches can make very different predictions about the future state of eco-230 logical systems. Which approach provides the most accurate forecasts depends on the forecast-231 horizon. In our simulations, time-series approaches performed best for short-term forecasts, 232 whereas models based on spatial data made more accurate long-term forecasts. In addition, our 233 simulations demonstrate extended transitional periods during which neither the time-series or 234 the spatial approach is effective on its own. The challenge is determining what is "short-term," 235 what is "long-term," and how to handle the many forecasts we need in ecology which fall in 236 Figure 5: (Top) Simulated annual temperatures (grey) and expected temperature (black), which was used to make forecasts. (Bottom) Simulated population size and forecasts from the spatial, temporal and weighted models.
between. We have proposed that a weighted combination of the time-series and space-for-time 237 approaches may produce better forecasts at these intermediate forecast horizons. 238 We designed our simulation studies to illustrate how the change in model performance with 239 increasing forecast horizon reflects differences in the types and scales of processes captured by 240 spatial and temporal data sets. How could these hypotheses be tested with empirical data?

241
The hypothesis that time-series models will be most effective for near-term forecasts already has

262
The greatest empirical challenge will be testing our hypothesis that a weighted average of 263 spatial and temporal models will make the best forecasts at intermediate time scales. There are 264 two problems: finding appropriate data and determining the model weights a priori. Many data 265 sets have both a longitudinal and spatial dimension, but we could not think of one which also 266 featured a clear ecological response to significant environmental change. Surely such datasets 267 exist, and we hope researchers who work with them will test our proposed weighted model. De-268 termining model weights may be more difficult. In our simulations, we fit the weights directly to 269 the simulated data, which is impossible to do for actual forecasting when the future is unknown. 270 We need new theory or empirical case studies in order to assign these weights a priori. known, estimating the high number of parameters and quantifying how they vary across ecosys-302 tems typically requires more data than is currently available even for well studied systems. As 303 a result, models used for ecological forecasting will include at least some phenomenological 304 components. But that does not mean that phenomenological forecast models do not benefit 305 from process-based understanding. The message from our simulations is that different processes 306 should be considered for different forecast time-scales, and this can be done by fitting models 307 to different kinds of datasets. Even when process-level understanding does not enable a fully 308 mechanistic model, it can improve the specification of phenomenological models.

309
While fully process-based models may not be practical, there are more mechanistic alterna-310 tives to our phenomenological, weighted model for integrating spatial and temporal information.

311
cesses of interest to ecological forecasters, such the spread of an invasive species or population 313 status of a threatened species (Wikle, 2003;Williams et al., 2017;Schliep et al., 2018). Because 314 these models include both fast processes, such as births and deaths, and slower processes, such as 315 colonization and extinction dynamics, they have the potential to make better predictions at inter-316 mediate forecast horizons than purely spatial or temporal models. However, these spatiotempo-317 ral models have rarely been used in a forecasting context, due to a combination of data limitation 318 and computational challenges. Many data sources contain either spatial or temporal variation, 319 but not both, and when spatiotemporal datasets are available they often involve irregular sam- 347 Appendices A Spatial, temporal and spatial-temporal-weighted models 441 The two simulation models in the main text describe how population size, N(x, t), at location x 442 changes over time (t). We assume that the temperature, K(x, t), at each location can vary in time 443 and space. To forecast the dynamics generated by these simulations models, we fit a series of 444 statistical models.

445
The spatial model, which we refer to as S, is a quadratic regression of the mean long-term 446 population density at a location (N(x)) against the mean temperature at that location (K(x)).

447
The quadratic term describes the unimodal relationship betweenN andK. The spatial statistical 448 model is The temporal model, which we call T, starts with a time-series of "observed" population 450 sizes, or total biomasses, at one location, N(t), for t = 1...n (the spatial index is suppressed 451 because we only focus on one location at a time). In the community turnover example, we fit the 452 following regression, which predicts biomass at time t + 1 as a function of biomass (N(t)) and 453 annual temperature (K(t)) at time t, In the eco-evolutionary example, the response variable is the log of the population growth rate.

455
The regression is This version of the temporal model returns a per capita growth rate on the log scale. To predict 457 population size at the next time step, we exponentiate the growth rate and multiply it by the 458 current population size: exp(T(N(t), K(t)))N(t).

459
The weighted model is a weighted average of predictions from the spatial and temporal 460 models, with the weights changing as a function of time, here expressed as the forecast horizon.

461
The weights change as a function of the square root of the forecast horizon, to allow rapid shifts 462 in the model weights.
For the community turnover example, the predicted biomass from the weighted model is: one location. For the eco-evolutionary example, the predicted population size from the weighted 466 model is: We used the optim function to estimate the β W s that minimize the sum of squared errors,

469
In the main text, we show the point forecasts but not the uncertainty around the forecasts.

470
After exploring that uncertainty, we decided that presenting it would be misleading. For the spa-471 tial and, especially, the temporal statistical models, the uncertainty is unrealistically low, because cesses: temperature-dependent growth, competition, and dispersal. Here we adapt their notation 483 to be consistent our own.

484
The population size of species i in cell x at time t + 1, N i (x, t + 1), is computed in two 485 steps. The first step accounts for changes in local population sizes due to dispersal. In each 486 local community, all species export a fraction (d) of their local population to the two adjacent 487 communities in the 1-dimensional landscape: Here N distinguishes the post-dispersal population size from the pre-dispersal population size.

489
The second step computes population growth, taking into account competition: In the absence of competition, the growth rate (g i ) is determined by the difference between the 491 temperature at site x (K(x)) and the focal species' minimum temperature tolerance, Kmin i , the ]N 1,t + λ 1 g 1 (K(t))N 1,t 1 + α 11 g 1 (K(t))N 1,t + α 12 g 2 (K(t))N 2,t N 2,t+1 = s 2 [1 − g 2 (K(t))]N 2,t + λ 2 g 2 (K(t))N 2,t 1 + α 21 g 1 (K(t))N 1,t + α 22 g 2 (K(t))N 2,t where g i (K(t)) is the probability of germination, K(t) is the temperature at time t, s i is the seed 502 survival probability for species i, and λ i is the seed production rate per plant. Below we refer to 503 the α ij as intra-and inter-genotype competition coefficients. and aa. The number of each genotypes at time t is N AA (t), N Aa (t), and N aa (t). The germination 506 rates for each genotype are g AA (K(t)), g Aa (K(t)), and g aa (K(t)). The seed survival probability 507 and seed production rate for genotype AA are s AA and λ AA , respectively. The analogous param-508 eters for the other genotypes are similarly denoted. The competition coefficients are denoted by 509 α i,j , e.g., α AA,AA or α AA,Aa . Throughout we assume that gametes mix randomly in the population.

510
First consider the case where the competition coefficients are zero (α i,j = 0). Let T denote the 511 total number of gamete-pairs produced in a given year, 512 T = λ AA N AA (t)g AA (K(t)) + λ Aa N Aa (t)g Aa (K(t)) + λ aa N aa (t)g aa (K(t)).
The first term is the number of gamete-pairs produced by AA individuals. The second and third terms are the numbers of gamete-pairs produced by Aa and aa individuals, respectively. The proportion of A gametes (φ A ) and the proportion of a gametes (φ a ) are given by φ A = λ AA N AA (t)g AA (K(t)) + 1 2 λ Aa N Aa (t)g Aa (K(t)) T and φ a = 1 − φ A .
Note that the T in the denominator of φ A shows up because we are computing proportions.
The first line is the number of gamete-pairs produced by AA individuals after accounting for the effects of competition. The second and third lines are the numbers of gamete-pairs produced by Aa and aa individuals, respectively. The proportions of A gametes and a gametes are φ A = 1 T λ AA N AA (t)g AA (K(t)) 1 + α AA,AA g AA (K(t))N AA (t) + α AA,Aa g Aa (K(t))N Aa (t) + α AA,aa g aa (K(t))N aa (t) λ Aa N Aa (t)g Aa (K(t)) 1 + α Aa,AA g AA (K(t))N AA (t) + α Aa,Aa g Aa (K(t))N Aa (t) + α Aa,aa g aa (K(t))N aa (t) Combining all of this results in the same model as above, N AA (t + 1) = s AA [1 − g AA (K(t))]N AA (t) + φ 2 A T N Aa (t + 1) = s Aa [1 − g Aa (K(t))]N Aa (t) + 2φ A φ a T N aa (t + 1) = s aa [1 − g aa (K(t))]N aa (t) + φ 2 a T, but the definitions of T, φ A , and φ a are given by equations (13) and (14) .