Seasonal difference in temporal transferability of an ecological model: near-term predictions of lemming outbreak abundances

Ecological models have been criticized for a lack of validation of their temporal transferability. Here we answer this call by investigating the temporal transferability of a dynamic state-space model developed to estimate season-dependent biotic and climatic predictors of spatial variability in outbreak abundance of the Norwegian lemming. Modelled summer and winter dynamics parametrized by spatial trapping data from one cyclic outbreak were validated with data from a subsequent outbreak. There was a distinct difference in model transferability between seasons. Summer dynamics had good temporal transferability, displaying ecological models’ potential to be temporally transferable. However, the winter dynamics transferred poorly. This discrepancy is likely due to a temporal inconsistency in the ability of the climate predictor (i.e. elevation) to reflect the winter conditions affecting lemmings both directly and indirectly. We conclude that there is an urgent need for data and models that yield better predictions of winter processes, in particular in face of the expected rapid climate change in the Arctic.

lemming (Lemmus lemmus) -either across space within a single outbreak 33 or across several outbreaks for a single locality 34 . Here we investigate to what extent a model based on abiotic (elevation as a proxy for local climate) and biotic predictors (intra -and interspecific density dependence) of lemming outbreak abundances in space (i.e. Ims et al. 33 ) has temporal transferability. Specifically, we do this by assessing the ability of a season-specific state-space model, parameterized by spatial data from one outbreak cycle to predict outbreak amplitudes in the subsequent cycle. By doing so we apply the framework of near-term forecasting of Dietze et al. 13 to advance our knowledge about what causes variability in lemming outbreak abundance.

Results
Overall characteristics of the abundance dynamics. Comparing the two consecutive cyclic peaks, there were some differences in estimated abundances (Table 1). While estimated lemming abundances were similar in the two peak springs, the spatial variability in spring densities (cf. CV values in Table 1) was substantially larger in 2007 than in 2011. Moreover, the growth of the lemming population over the peak summer was higher in 2007 than in 2011 (Table 1), leading to 40% higher estimated autumn abundance in the first peak compared to the second peak. For the grey-sided vole that co-occur with the lemming in our study region, the estimated pre-peak autumn abundances were similar between the two peaks, while the estimated spring abundance was somewhat higher in the second compared to the first peak ( Table 1).
Predictors of lemming peak abundances. The estimated coefficients of the fitted models for each of the two lemming peaks are given in Table 2. Lemming abundances were generally negatively affected by intra-specific density dependence. The elevation effect was consistently positive for the autumn abundances (i.e. increasing abundance with elevation). A positive effect of elevation was also found in the spring of the first peak, while there was no such over-winter effect during the second peak. Considering the interaction between the rodent species (inter-specific density dependence), there was a tendency for a positive effect of grey-sided vole autumn abundance on lemming spring abundance in the first peak, while we found no such effect in the second peak (credible intervals encompassed 0 with good margins; Table 2). The estimates for a grey-sided vole effects on autumn lemming abundance were smaller and less certain (credible intervals encompassed 0 with good margins; Table 2) in both peaks.
Transferability. The temporally consistent parameter estimates of the autumn dynamics in both peaks contributed to a relatively good ability of the estimated autumn dynamics in 2007 to predict the site-specific autumn lemming abundances for 2011 (Fig. 2). The mean absolute error (MAE) was 0.55 individual per site for this model and there was no apparent bias (Fig. 2). In contrast, the spring part of the model parameterized by data from the first peak (2007) performed poorly in terms of its temporal transferability to the spring of the next peak (2011); as could be expected from the inconsistent estimates of the spring dynamics in the two peaks ( Table 2). The MAE for the spring predictions was 0.70 individual per site and a strong bias was evident (Fig. 2). The MAE values have to be seen relative to the mean season specific abundance, which is more than 3 times higher in autumn than in spring. Therefore, the calculated MAE values show a clear difference in model transferability.

Discussion
We investigated the temporal transferability of a dynamical state-space model that was developed to identify season-specific biotic and abiotic predictors of cyclic lemming outbreaks. Based on spatial data from one lemming outbreak, Ims et al. 33 found that a relatively simple model (i.e. with intra -and interspecific-density dependence and elevation as predictors) explained well the spatial variation in outbreak abundances. However, our results show that the temporal transferability of the model with respect to the subsequent cyclic lemming outbreak was only partial. That is, the model part projecting autumn abundances (i.e. reflecting population change over the summer) exhibited good transferability, whereas the model part predicting spring abundances (i.e. reflecting population change over the winter) performed poorly.
Previous studies have claimed that highly detailed knowledge about a modest number of interactions would be most beneficial regarding model transferability 35 . Indeed, our state-space model included few biotic interactions -i.e. only direct inter-and intra-specific density dependence -yet it appeared to be temporally transferable with respect to predicting population changes over lemming outbreak summers. This means that both the data (site-specific rodent abundance data and elevation) and the model appears adequate for near-term forecasting of lemming outbreak abundances in the autumn. However, the forecast of lemming spring abundances performed poorly, meaning that the data/model used for this purpose did not prove to be adequate/transferable. This seasonal difference in the degree of model transferability is interesting. Arctic ecosystems are known to experience high temporal environmental variability both within (i.e. seasonality) and between years 36,37 . However, climatic variation between summers are known to be lower than between winters, in particular for the Atlantic sector of the Arctic 36 . Especially, large inter-annual differences in qualitative snow characteristics, towards which lemmings exhibit high sensitivity 34 , adds significantly to the winter variability. In particular, mild winters increase the hardness and humidity of the snow that impact lemming survival negatively 34 . Thus it appears that detailed temporal data and understanding of the impact of winter must be incorporated in models to provide temporally consistent predictions. Kausrud et al. 34 did this with good result when they projected lemming dynamics for a single site. In the present study we attempted to make projection across a large number of sites within an area of approximately 10 000 km 2 with elevation as a proxy of spatial climatic variation. Elevation has previously been used as a proxy for spatial variation in climate in many ecological studies 38 including studies of outbreak amplitude of cyclic herbivore populations 39,40 . However, while elevation gradients may reflect spatial differences in snow conditions in some winters (i.e. winter 2006/2007), it may not in other winters (i.e. winter 2010/2011). Hansen et al. 41 recently demonstrated that climatic extreme events during the winter in the high arctic could disconnect the association between snow quality and elevation. Generally, transferability of ecological models based on spatial data has been found to decrease whenever the magnitude and nature of the spatial and environmental variation differs between temporal domains 26,42 .
This study based on data from a relatively new program to monitor population dynamics of tundra rodents, should be seen as an initial loop of the iterative near-term forecasting cycle of Dietze et al. 13 . Thereby we have learned that the summer dynamics of outbreaking Norwegian lemming populations is near-term predictable based on the trapping data and elevation, whereas clearly more information is needed to be able to predict the pre-outbreak winter dynamics. We consider this lesson to be particular important in face of the ongoing rapid change in winter climate in the Arctic 43 .  System and sampling. The study was conducted within a tundra area of about 10 000 km 2 at the north-easternmost tip of the Scandinavian Peninsula (70°N to 71°N). Rodent cycles with periodicity of 4-5 years prevail in the focal tundra ecosystem, with grey-sided voles (Myodes rufocanus) and the Norwegian lemming as the most abundant species 44 .

Methods
Since spring 2004, small rodent snap trapping has been performed on 98-109 sites according to the small quadrat method 45 , with one quadrat (i.e. 12 traps) per site (see Ims et al. 33 for details). In order to include spatial variability in environmental conditions, the design contains trapping sites that span a range of 30 to 346 m.a.s.l. (mean of ~ 200 m.a.s.l.). The orographic effect of elevation amounts to a decrease of approximately 0.6 °C per 100 m 46 , making elevation a proxy for spatial variation in temperature. Trapping was conducted twice annually; 2 days in late June (spring) and two days in early September (autumn) before the onset of winter.
Statistical modelling: analyses and validation. Following Ims et al. 33 , the trapping data were analyzed at the site level (i) including data for the two cyclic lemming peaks (t = 2 peaks) contained in the time series (see Typically, Norwegian lemmings are mostly absent in trapping data between peaks (e.g. Turchin et al. 47 , Fig. 1) thus we included only data from the pre-peak autumn (k-2) together with spring (k-1) and autumn (k) in the lemming peak years (k = 3 trapping seasons). The data for the two peaks were analyzed separately with the purpose of assessing the transferability of predictors of lemming outbreak abundance across different cycles. The predictors investigated were the same of those identified by Ims et al. 33 ; elevation as a proxy for spatial climate variation and lemming and grey-sided vole abundance to model intra -and interspecific density dependence. The interspecific density dependence is most likely due to the influence of shared predators 33,48 . Ims et al. 33 found that there was no residual spatial autocorrelation in the lemming abundance data so we did not include any extra spatial terms in the models. Moreover, previous time series analyses (e.g. Stenseth et al. 49 ) have shown that there is no time-lags >2 years in small rodent population dynamics, meaning that consecutive cycles can be considered independent.
Small rodent trapping data includes stochastic sampling variability, therefore we analyzed the data using a state space model. We modelled the sampling variance in the number of trapped lemmings y ( ) i k t , , and grey-sided voles (x ) i,k,t per site (i), season (k) and peak (t) as a Poisson process (λ) 33,49 . We used the mean absolute predictive error (MAE) 50 to evaluate model fit (Appendix S1). We also plotted the estimated counts against observed counts to investigate whether there were some systematic differences between raw counts and estimated abundance.
With some small modifications from the model of Ims et al. 33 (see Appendix S3), we then applied the following state-space model to estimate the season-(spring (k = 2) and fall (k = 3)) and peak outbreak-specific (year 2007 (t = 1) and 2011 (t = 2)) effects of elevation β ( elev ), inter-specific (β dvole ) and intra-specific (β dlem ) density dependence on lemming abundance: where σ t is the standard deviation. For lemming abundance in the initial season (k = 1) and for grey sided vole abundance in all seasons (k = 1:3), the trapping data is also assumed to follow a Poisson process with mean λ i,t . However, since delayed effects cannot be included in the initial season, log (λ i,t ) is modelled as Norm (µ i,t , σ t ) where µ i,t is a site-specific intercept. We checked that this difference between the autumn and spring models did not affect our conclusions regarding transferability, by fitting also a model with only the spring densities as a predictor of the fall densities (see Appendix S5). Finally, to evaluate the temporal transferability, the model parameter estimates obtained based on data from the first peak (t = 1) was applied to the predictor data for the second peak (t = 2) to derive predicted lemming abundance ( λ p y i k , ): The predicted lemming abundance ( λ p y i k , ) was validated against the estimated posterior means for lemming abundance for the second peak (λ = y i k t , , 2 ) by means of the mean absolute error 50 : with P being the predicted abundance ( λ p y i k , ) and O the abundance estimated with the Poisson state-space model for the second peak (λ = y i k t , , 2 ). The mean (P and O) is subtracted to account for seasonal differences in abundance.
The state space models were specified in a Bayesian framework and priors were kept uninformative 51 . Posterior distributions were obtained using Markov Chain Monte Carlo (MCMC) techniques computed through Jags run from R (R Core Team 2015) using the jagsUI package. We used 4 chains, each of 50 000 iterations, with a burn-in of 15 000 (see Appendix S2 for details). To assess convergence of the chains, trace plots for all parameters where investigated graphically as well as from the Gelman-Rubin statistics (where R <1.1 indicates convergence) 52 .

Data Availability
All data analyzed in this paper is available on DRYAD, doi:10.5061.