Matching the forecast horizon with the relevant spatial and temporal processes and data sources

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


38
Forecasting is increasingly recognized as important to the application and advancement of eco-39 logical research. Forecasts are necessary to guide environmental policy and management deci-40 sions about mitigation and adaption to global change (Clark et al., 2001;Mouquet et al., 2015;41 Dietze et al., 2018). But forecasts can also advance understanding of the processes governing example, meteorological models for short-term weather forecasts differ substantially in spatial 97 and temporal resolution and extent from the global circulation models used to predict long-term 98 changes in climate. 99 Why not simply use process-based models to avoid the difficulties posed by phenomenolog-100 ical time-series and space-for-time modeling approaches? If we could accurately characterize all 101 of the processes governing a system, then a model based on that understanding should make 102 accurate predictions at all forecast horizons. Process-based models should also be more robust 103 for making predictions outside of historically observed conditions and even beyond the con-104 ditions observed across spatial gradients, which will be especially important in a future with 105 increasingly novel combinations of environment and species interactions (Williams and Jackson, 106 2007). Unfortunately, in most cases this approach is not currently feasible because we lack a 107 detailed knowledge of all the complex and interacting processes influencing the dynamics of real 108 ecological systems. Even if the general form of the models were known, estimating the high 109 number of parameters and quantifying how they vary across ecosystems typically requires more 110 data than is currently available even for well-studied systems. Furthermore, the high complexity 111 and corresponding parameter uncertainty of such models can increase predictive errors; simpler 112 time-series models may actually perform better (Ward et al., 2014), though spatial replication can 113 reduce the cost of complexity (Chevalier and Knape, 2020). As a result, models used for eco-114 logical forecasting will include at least some phenomenological components. But that does not 115 mean that phenomenological forecast models cannot benefit from process-based understanding. 116 Even if process-level understanding does not enable a fully mechanistic model, it can improve 117 the specification of phenomenological models. Our hypothesis is that different processes may be 118 relevant for different forecast horizons, and that we can act on this knowledge by fitting models 119 to different kinds of datasets. 120 Here we use two simulated case studies to 1) demonstrate that the best model-building ap-121 proaches for ecological forecasting depend on the time horizon of the forecast, and 2) explore  2. Different kinds of data reflect the operation of different processes: longitudinal data cap-130 ture autocorrelation and fast responses of current assemblages to interannual environmental 131 variation, while data spanning spatial gradients capture the long-term outcome of interac-132 tions between fast and slow processes. Whether predictive models should be trained using longitudinal or spatial data sets, or both, depends on the time-scale of the desired forecast. 134 3. A key challenge for future research is determining the rate at which slow processes begin to influence dynamics.

137
In each case study, we simulated the effects of an increase in temperature on simple systems with 138 known dynamics. The truth was represented by a simulation model that was mechanistic for at 139 least one important process, but we treated the data-generating model as unknown when ana-140 lyzing the data and we assumed that perfectly recovering the mechanisms it contains would not 141 be possible in practice. We began each simulation under a stationary distribution of annual tem-142 peratures, allowing the system to equilibrate; we call this the baseline phase. We then increased 143 temperature progressively over a period of time, followed by a second period of stationary, now 144 elevated, temperature. The objective was to forecast the response of the system to the tempera-145 ture increase based on spatial and/or temporal data "sampled" from the simulation during the 146 baseline period. 147 We made forecasts based on two phenomenological statistical models, each representing pro-148 cesses operating at different time scales. One statistical model represents the time-series or "tem-149 poral approach." We correlated interannual variation in an ecological response with interannual 150 variation in temperature at just one site. The other statistical model relies on a space-for-time 151 substitution, which we call the "spatial approach" for brevity. We correlated the mean tempera-152 ture with the mean of an ecological state or rate across many sites. We compared forecasts from 153 both statistical models to the simulated dynamics to determine how well the two approaches 154 performed at different forecast horizons. We also assessed the potential for combining the infor-155 mation available in temporal and spatial patterns by using a weighted average of the forecasts 156 from the temporal and spatial approaches optimized to best match the (simulated) observations. 157 We then studied how the optimal model weights changed over time. We expected the tem-158 poral approach to best predict short-term dynamics, the spatial approach to best predict long-159 term dynamics, while the weighted model would show potential to provide the best forecasts at  We first simulated a baseline period with variable but stationary temperature, followed by 183 a period of rapid temperature increase, and then a final period of stationary temperature. In-184 terannual variation in temperature is the same at all sites, but mean temperature varies among 185 sites. All sites experienced the same absolute increase in mean temperature. We focused on the 186 biomass dynamics of one focal species that dominated the central site during the baseline period.

187
During the baseline period there were strong spatial patterns across the mean temperature 188 gradient. Individual species, including our focal species, showed classic, unimodal "Whittaker" 189 patterns of abundances across the gradient ( Fig. 2A). These spatial patterns are the basis for our 190 spatial statistical model of the temperature-biomass relationship for our focal species ( Fig. 2A).

191
In contrast to the strong spatial patterns, population and community responses to interannual year. 197 We used both the temporal and spatial statistical models to forecast the effect of a temperature 198 increase (Fig. 3A) on the focal species' biomass at one location in the center of the temperature 199 gradient. The predictions from these two models contrasted markedly, with the temporal sta-200 tistical model predicting a large increase in biomass and the spatial statistical model predicting 201 a decrease. Initially, the simulated abundances followed the increase predicted by the temporal 202 model, but as faster-growing species colonized and increased in abundance at the focal site, the 203 biomass of the focal species decreased, eventually falling below its baseline level (Fig. 3B).

204
To combine information from the temporal and spatial statistical models into a single pre-205 diction, we fit a weighting parameter, ω, which varies over time and is bounded between 0 and 206 1. At any time point, t, this weighted forecast is ω · T(N t−1 , K t ) + (1 − ω) · S(K t ) where T is the 207 temporal statistical model, which depends on population size, N, and expected temperature, K, 208 and S is the spatial statistical model, which depends only on K (see SI Appendix A for a full 209 description of the approach). The weighted model accurately predicts the simulated dynamics 210 across the full forecast horizon (Fig. 3B). It also shows that the most rapid shifts in the model 211 weights occurred during the period when warm-adapted, faster growing species were increasing 212 most rapidly in abundance (Fig. 3C). However, the reason the weighted models works so well 213 is that the weights were determined by fitting directly to the data. Unlike the forecasts from the 214 spatial and temporal statistical models, we did not generate out-of-sample predictions from the 215 weighted model; it merely provides a convenient way to quantify how rapidly dynamics shift  The compositional turnover affecting our focal species also influences total biomass, linking 221 community and ecosystem dynamics. We repeated our focal species analysis for total community 222 biomass, and the results were similar: the temporal statistical model initially made the best 223 forecasts immediately following the onset of the temperature increase, but as the identity and 224 abundances of species at the study site changed, the model weights rapidly shifted to the spatial  rium densities across the mean annual temperature gradient is the basis for our spatial statistical 243 model: cold sites will be dominated by the cold-adapted homozygous genotype, warm sites 244 will be dominated by the heat-adapted homozygous genotype, and intermediate sites will be 245 dominated by the heterozygous genotype (Fig. 4B). The full description of the eco-evolutionary 246 simulation model is provided in SI Appendix C.

247
The spatial pattern shown in Fig. 4B is the outcome of steady-state conditions. But at any one 248 site, the population's short-term response to temperature will be determined by the dominant 249 genotype's reaction norm (Fig. 4A). For example, at a cold site dominated by the cold-adapted 250 homozygous genotype, a warmer than average year would cause a decrease in population size 251 due to decreases in fecundity (blue line in Fig. 4A), even though the heat-adapted homozygote 252 might perform optimally at that temperature. However, if warmer than normal conditions persist 253 for many years, then genotype frequencies should shift, and the heat-adapted homozygote will 254 compensate for the decreases of the cold-adapted genotype.

255
To demonstrate these dynamics, we simulated a diploid annual plant population at a colder 256 than average site. During the baseline period, the population is dominated by the cold-adapted 257 genotype. We used the simulated data from this baseline period to fit a temporal statistical model

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

278
Ecological forecasts are typically made using either a space-for-time substitution approach based 279 on models fit to spatial data or using dynamic models fit to time-series data. Our results demon-280 strate that these two approaches can make very different predictions about the future state 281 of ecological systems. Which approach provides the most accurate forecasts depends on the 282 forecast-horizon. In our simulations, time-series approaches performed best for short forecast 283 horizons, whereas models based on spatial data made more accurate forecasts at long horizons.

284
In addition, our simulations demonstrate extended transitional periods during which neither the 285 time-series or the spatial approach was effective on its own. The challenge is determining what 286 is "short-term," what is "long-term," and how to handle the many forecasts we need in ecology 287 which fall in between. We have proposed that a weighted combination of the time-series and 288 space-for-time approaches may produce better forecasts at these intermediate forecast horizons. 289 We designed our simulation studies to illustrate how the change in statistical model perfor-290 mance with increasing forecast horizon reflects differences in the types and scales of processes 291 captured by spatial and temporal data sets. How could these hypotheses be tested with empir-292 ical data? The hypothesis that time-series models will be most effective for near-term forecasts   The greatest empirical challenge will be testing our hypothesis that a weighted average of 314 spatial and temporal statistical models will make the best forecasts at intermediate time scales.

315
There are two problems: finding appropriate data and determining the model weights a priori.

316
Many data sets have both a longitudinal and spatial dimension, but we could not think of one 317 which also featured a clear ecological response to significant environmental change. Surely such 318 datasets exist, and we hope researchers who work with them will test our proposed weighted 319 model. Determining model weights may be more difficult. In our simulations, we fit the weights 320 directly to the simulated data, which is impossible to do for actual forecasting when the future is 321 unknown. We need new theory or empirical case studies in order to assign these weights a priori.

322
Theory could explore the influence of different parameters on the rate at which slow processes 323 begin to influence dynamics. The effects of some parameters are intuitive: in the community 324 turnover example, increasing the fraction of dispersing individuals caused a more rapid shift in 325 species composition and in model weights (Fig. 6A). Other parameters have less intuitive effects: such as births and deaths, and slower processes, such as colonization and extinction dynamics, 347 they have the potential to make better predictions at intermediate forecast horizons than purely 348 spatial or temporal models. However, these spatiotemporal models have rarely been used in a 349 forecasting context, due to a combination of data limitation and computational challenges. Many 350 data sources contain either spatial or temporal variation, but not both, and when spatiotempo-351 ral datasets are available they often involve irregular sampling, creating challenges for modeling.

352
Fitting and generating predictions from spatiotemporal models is also computationally intensive,  Weather forecasts based on regional-scale meteorological models are very effective for forecast-369 ing a week to ten days in advance, but then become largely uninformative. Forecasting these 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.   The spatial pattern of individual genotypes (colors) and total population abundance (black) at sites arrayed across a gradient of mean annual temperature. The dashed black line (almost entirely hidden by the slid black line) shows predictions from an empirical, spatial statistical 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 coldadapted 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. 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 statistical models.

Appendices
A Spatial, temporal and spatial-temporal-weighted models 455 The two simulation models in the main text describe how population size, N(x, t), at location x 456 changes over time (t). We assume that the temperature, K(x, t), at each location can vary in time 457 and space. To forecast the dynamics generated by these simulations models, we fit a series of 458 statistical models.

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

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

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

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

475
The weights change as a function of the square root of the forecast horizon, to allow rapid shifts 476 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 480 model is: We used the optim function to estimate the β W s that minimize the sum of squared errors,

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

484
After exploring that uncertainty, we decided that presenting it would be misleading. For the spa-485 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 497 to be consistent our own.

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

503
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 505 temperature at site x (K(x)) and the focal species' minimum temperature tolerance, Kmin i , the 506 lowest temperature at which a species can maintain a positive growth rate. Growth is further 507 reduced by intraspecific and interspecific competition, parameterized by c i and l i . All species are 508 assigned the same value of c i , which represents an additional effect of intraspecific competition 509 on top of interspecific competition. This stabilizes coexistence, since every species will exert 510 stronger intra-than interspecific competition. However, values of l vary among species to create 511 a trade-off between growth rates and competitive ability versus low temperature tolerance: fast-512 growing species (high g i ) are more tolerant of interspecific competition (low l i ) but are more 513 limited by temperature (high Kmin i ).

518
Diploid Model: Consider a one-species diploid model. The genotypes are denoted by AA, Aa, 519 and aa. The number of each genotypes at time t is N AA (t), N Aa (t), and N aa (t). The germination 520 rates for each genotype are g AA (K(t)), g Aa (K(t)), and g aa (K(t)). The seed survival probability 521 and seed production rate for genotype AA are s AA and λ AA , respectively. The analogous param-522 eters for the other genotypes are similarly denoted. The competition coefficients are denoted by 523 α i,j , e.g., α AA,AA or α AA,Aa . Throughout we assume that gametes mix randomly in the population.

524
First consider the case where the competition coefficients are zero (α i,j = 0). Let T denote the 525 total number of gamete-pairs produced in a given year, 526 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) .