Population response to climate change: linear vs. non-linear modeling approaches.

BACKGROUND: Research on the ecological consequences of global climate change has elicited a growing interest in the use of time series analysis to investigate population dynamics in a changing climate. Here, we compare linear and non-linear models describing the contribution of climate to the density fluctuations of the population of wolves on Isle Royale, Michigan from 1959 to 1999.
RESULTS: The non-linear self excitatory threshold autoregressive (SETAR) model revealed that, due to differences in the strength and nature of density dependence, relatively small and large populations may be differentially affected by future changes in climate. Both linear and non-linear models predict a decrease in the population of wolves with predicted changes in climate.
CONCLUSIONS: Because specific predictions differed between linear and non-linear models, our study highlights the importance of using non-linear methods that allow the detection of non-linearity in the strength and nature of density dependence. Failure to adopt a non-linear approach to modelling population response to climate change, either exclusively or in addition to linear approaches, may compromise efforts to quantify ecological consequences of future warming.

Background

Several recent studies have suggested that species from many taxa and environments are already responding to regional climate change with, for example, poleward and elevational range shifts, changes in abundance, changes in body size, and changes in the timing of life history events [1-3]. Because these changes may influence important stabilizing mechanisms in communities and ecosystems, and, thus, the goods and services they provide to our society, and because the rates of such changes may continue to increase [3], they have moved to the forefront of political and environmental concern [1]. Nonetheless, there is still considerable uncertainty regarding the magnitude of future changes in climate on both local and regional scales and the impact they will have on populations, communities, and ecosystems. Clearly, studying the consequences of climate change for the dynamics and stability of populations will improve our ability to predict these effects and respond accordingly.

One growing body of research investigating this issue has used time series analysis of long-term datasets to quantify the contribution of large-scale climatic fluctuation to population dynamics (reviewed in [4] and [5]). The approach commonly used in such studies is log-linear autoregressive modeling. The dynamical properties of this approach make it a logical, tractable, and appropriate candidate for examining and quantitatively predicting the influence of climate change on mean population density [6].

While appropriate in some situations, linear models may obscure potentially important properties of the population, which may ultimately influence the accuracy of predictions drawn from them. For example, linear models cannot detect non-linear density dependence and the differential influence of abiotic factors above and below certain density thresholds (see, e.g., [7]). While several studies have indicated the potential for non-linear density-dependence in natural populations, the consequences of incorrectly assuming linear density-dependence for the accuracy of model predictions is not well understood.

Here, we present a comparative analysis of linear and non-linear models describing the contribution of climate to the interannual density fluctuations of the population of wolves (Canis lupus) on Isle Royale, Michigan, USA, from 1959 to 1999. Our analysis builds on previous studies indicating that the dynamics of the focal wolf population may be influenced by large-scale climatic variation [8,9], and that the population may exhibit phase-dependence [10-12] that could generate non-linearity in the dynamics [13]. Thus, there are indications that climate change may have important implications for the dynamics of wolves on Isle Royale, and, in addition, there are convincing arguments suggesting that non-linear modeling may be appropriate for investigating this issue. We proceed by first developing a non-linear self-excitatory threshold autoregressive (SETAR) model for the focal population. We then develop a linear autoregressive model using similar methods and compare the stability properties, statistical fit, and predictions obtained by both models. Our analysis indicates that linear and non-linear models may, in fact, lead to significantly different predictions of population response to climate change.

Results

The fit of the models

The most parsimonious non-linear model identified here includes the number of packs in the previous year and both temperature and snow in the current year as covariates in the low-density regime; and pack size in the previous year, the number of packs in the previous year, and snow of the current year as covariates in the high-density regime (Table 1). Piecewise non-linear modeling indicated that the optimal density threshold for the SETAR model was 3.4, which on a raw scale is a density of approximately 30 wolves (Table 1). The non-linear SETAR model is discontinuous at this threshold. In the lower regime, only the parameters for the delayed density and temperature covariates are significant at α = 0.1 Table 1). In contrast, only the parameters for delayed density and delayed number of packs covariates are significant in the upper regime (Table 1).

Table 1

The most parsimonious SETAR models of the population dynamics of wolves on Isle Royale, Michigan, 1959–99. Covariates included winter snow accumulation (SW), Northern Hemisphere temperature anomalies (T), number of packs in the previous year (PA), and mean pack size in the previous year (PS). Xis loge-transformed density, and aare statistical parameters (i = 1 and 2 corresponds to lower and upper regimes, respectively, j = 0, 1, 2, 3 correspond to the constant, lag-1 density coefficients, number of packs coefficients, and climatic coefficients, respectively). n indicates the number of data points in each regime. Parameter estimates were obtained by the method of conditional least squares. X* is the equilibrium point on the log-scale.

	Full Model	θ	Coefficients	SE	p	n	AIC	R2	Equilibrium
	Xt = a1,0 + a1,1Xt-1 + a1,2PAt-1 + a1,3Tt + a1,4SWt	Xt-1 < 3.40	a1,0 = 1.89	0.43	0.0002	32	14.18	0.58	X* = 2.40
Non-Linear Model			a1,1 = 0.50	0.14	0.001				Stable
			a1,2 = -0.06	0.04	0.17
			a1,3 = -0.47	0.16	0.009
			a1,4 = -.0008	0.0005	0.15
	Xt = a2,0 + a2,1Xt-1 + a2,2PSt-1 + a2,3PAt-1 + a2,4SWt	Xt-1 ≥ 3.40	a2,0 = -5.82	1.43	0.03	8		0.96	X* = 3.58
			a2,1 = 2.63	0.46	0.01				Unstable
			a2,2 = 0.06	0.03	0.11
			a2,3 = -0.24	0.05	0.02
			a2,4 = 0.002	0.0009	0.12
Linear Model	Xt = a0+ a1Xt-1 + a2PAt-1 + a3Tt		a0 = 0.86	0.33	0.01	40	27.98	0.69	X* = 2.97
			a1 = 0.81	0.12	0.00				Stable
			a2 = -0.07	0.03	0.06
			a3 = -0.39	0.18	0.04

The most parsimonious linear model includes the number of packs in the previous year and temperature of the current year as covariates and the parameters for these covariate are significant at α = 0.1 (Table 1). Given only the first observed value of the density time series and the observed values of snow, temperature, the number of packs, and pack size, both models simulate dynamics reasonably well (Fig. 2). However, AIC scores indicate a better fit for the non-linear model (Table 1).

Figure 2

Dynamics predicted by non-linear SETAR model and linear model. Predicted dynamics were calculated from the full models in Table 1 using one seed density estimates from 1960 and the observed climate and pack time series (one step ahead simulations). (a). Dynamics predicted by SETAR model with winter snow accumulation, SW, as a covariate. (b). Dynamics predicted by linear model with the number of packs, PA, and the Northern Hemisphere winter temperature anomaly, T, as covariates. (c). Predicted versus observed values for the non-linear model. (d). Predicted versus observed values for the linear model.

Stability properties

The stability of the dynamics within the each regime of the non-linear model and for the linear model were determined by the values of the coefficients of density dependence (a1,1, a2,1, and a1 in Table 1). Ninety-five percent confidence intervals for the coefficient of direct density dependence in the lower regime (0.05 ± 0.28) and the upper regime (2.63 ± 0.13) do not overlap suggesting that the strength of density dependence differs between the two regimes. The equilibrium point in the lower regime is stable, while that for the upper regime is not (Table 1). The equilibrium point in the linear model is stable (Table 1). The location of the equilibrium points in both regimes of the non-linear model and in the linear model, however, depends on climate (Table 1). If the average values of pack size (8.24), number of packs (3.7), temperature (0.10), and snow (190.4) from the observed time series are used, then the equilibrium points of the lower and upper regime of the non-linear model and of the linear model are approximately 11, 36, and 20 wolves on the raw scale, respectively (Table 1).

Discussion

Non-linear SETAR model

Previous work has revealed that density dependence in the wolf population on Isle Royale is phase-dependent [12], which may contribute to non-linearity in the dynamics of this population. Importantly, our results also indicate that there are differences in the nature of density dependence between low- and high-density regimes (Table 1, Fig. 1). Analysis of the model suggests that populations below the density threshold experience negative density dependence and exhibit monotonic damping, while populations above the threshold display positive density dependence.

Figure 1

Population growth rate (R= X- X) of wolves versus ln-transformed density in the previous year (X). Solid lines indicate non-linear spline functions, and dashed lines indicate 95% confidence bands, estimated with generalized additive models. These plots suggest a threshold at ~3.25–3.50 wolves on a loge scale. The linear relationship below this threshold suggests weak density dependence while the relationship above this threshold suggests that density dependence is stronger at higher densities.

The social behavior of wolves provides a convincing biological explanation for the presence of a density dependence threshold and positive density dependence above the threshold. Because they live on an island and, thus, have a finite amount of space, growth in the population must be accompanied by an increase in the number of wolves per pack and/or the number of packs on the island. A closer examination of the data demonstrates that at low densities, population growth on Isle Royale is accompanied by an increase in average pack size (Fig. 3). The number of packs on the island remains relatively constant (3 to 5 packs). When the total population size reaches 30 wolves, however, the number of packs begins to increase (Fig. 3). The average pack size remains relatively constant at 10–15 wolves per pack. These data suggest that there may be an upper limit to the size of a pack, possibly due to difficulties obtaining enough to food to feed all pack members or due to an increase in aggressive behaviors. This upper limit in pack size may generate a density threshold. In addition, this shift in behavior may lead to positive density dependence at high densities. An increase in the number of packs corresponds to an increase in the number of breeding pairs on the island. As density continues to increase, more packs are formed, leading to higher population growth. However, and importantly, this positive feedback loop is inherently unstable. As the number of packs on the island increases, the amount of territory per pack decreases and the frequency of aggressive encounters among packs must increase. Together, these processes should lead to increased mortality. In fact, the observed dynamics of the population support this theory. The first time the population reached 30 wolves in 1965, it declined over the next four years, and the two times it has crossed this threshold, between 1978 and 1980, it has crashed precipitously [12].

Figure 3

(a) Number of packs (PA) versus Number of wolves (N). (b) Average pack size (PS) versus Number of wolves (N). Threshold for non-linear model is ~30 wolves on the raw scale.

Population response to climate change

Differences in the nature of density dependence above and below the threshold suggest that the consequences of climate change for the dynamics of Isle Royale wolves differ depending on if the population is relatively small or large. At low densities dynamics are largely determined by within pack dynamics, specifically their ability to capture prey and feed all pack members. Because changes in climate influence the ability of wolves to capture prey, climate change should have significant effects when the population is small. Because dynamics maybe determined more by interpack dynamics and less by the ability of individual packs to capture prey at high densities, climate and climate change should be less important in determining population dynamics when the population is large. Climate may influence the dynamics of large populations, but the signature of climatic effects may be diluted by the effects of interpack aggression. These predictions are supported by the fact that though the most parsimonious non-linear model contains climate as covariates in both regimes, the only significant coefficient of a climatic variable occurs in the low-density regime (Table 1). Thus, our results suggest that the dynamics of this population, when its density is low, are stable, determined mainly by direct density dependence, and affected by changes in climate. At higher densities, the dynamics are intrinsically unstable and the population is less influenced by climate.

In addition to predicting when climate and climate change is important, the non-linear SETAR model suggests that predicted changes in climate will cause the population of wolves on Isle Royale to decline, possibly to extinction. While wolf density is negatively related to both temperature and snow in the low density regime, the coefficient for snow is not significant (Table 1). This suggests that below the density threshold, changes in snow depth will have a lesser effect on the population than changes in temperature. Hence, predicted increases in temperature may lead to a decrease in the population of wolves on Isle Royale, possibly contributing to extinction risk. At densities above the threshold, the population is unlikely to be systematically affected directly by changes in climate. However, because dynamics are unstable at high densities, large populations are likely to crash, falling below the threshold where predicted changes in climate will cause the population to decline.

While increases in mean temperature may cause the population to decline, this warming trend may also be accompanied by increases in variability that could have additional effects on the population. Extreme variations in climate may cause the population of wolves on Isle Royale to fluctuate above and below the density threshold. This may cause extreme fluctuations in the size of the population greatly increasing the chance of extinction.

Comparison to linear model

Qualitatively, the linear model obtained here predicts similar responses to climate change in comparison to the non-linear model: increasing temperatures will lead to declines in population density. However, the linear model indicates negative density dependence for the whole range of observed densities and suggests that the number of packs is important in determining the dynamics of the population. Importantly, the absence of a density threshold and positive density dependence in this model suggests that the population should experience relatively small fluctuations in density in constant and variable environments. Negative density dependence should regulate the population at all densities and buffer the population from changes in climate. Thus, in contrast to the non-linear SETAR model, the linear model predicts similar effects of climate and climate change on relatively small and large populations, predicts that the population is relatively stable, and suggests that population fluctuations will be relatively small.

In terms of statistical fit to the data, the non-linear SETAR model performs better than the linear model, indicated by the fact that the AIC of the linear model is almost double that of the non-linear model (Table 1). Because these are phenomenological models, this criterion does not necessarily mean that the non-linear model more accurately portrays the biology of the system or that it is better for the purposes of predicting response to climate change. In fact, we urge caution in interpretations of both models. For the non-linear model in particular, there are only eight data points, but relatively many parameters in the upper regime. This decreases certainty surrounding parameter estimates and thus, our conclusions regarding stability and response to climate change for large populations. In addition, despite differences in AIC scores, the linear model still provides a reasonably good fit to the data. It may, therefore, give just as reliable predictions as the non-linear model. The main contribution of our analysis is its demonstration of how linear versus non-linear time series modeling approaches may influence conclusions regarding population response to climate change.

The importance of non-linear time-series modelling

Using time series data for examining population dynamics has both disadvantages and advantages, regardless of whether linear or non-linear approaches are applied. Our ability to model the dynamics of a population and the reliability of the predictions that can be drawn from it depend on the data available and the length of the time series. As mentioned above, trophic interactions were not directly included in either model but could be very important in determining not only the dynamics of the population of wolves on Isle Royale, but also its response to climate change. In addition, though this data set represents one of the longest time series of a three-trophic level system, it spans only 40 years and teasing out the relative influences of intrinsic versus extrinsic processes on the dynamics of long-lived species may require longer time series [14]. In addition, it is possible that the factors most important in determining dynamics vary through time. For example, it is possible that the population was exposed to canine parvovirus in the 1980s [11]. The dynamics of the population during this time period may differ in important ways that cannot be captured in a simple model of the full time series lacking this information.

Despite the potential difficulties, the analysis of time series data offers several advantages. First, these methods provide one of the only means for examining population dynamics when experimentation and detailed demographic analyses are impractical or impossible, as they are on Isle Royale. Second, these methods allow analyses of population dynamics when only limited data are available. Finally, time series analysis is useful when examining processes that occur over large spatio-temporal scales such as climate change.

Conclusions

While time series methods may allow us to predict the effects of climate change on population dynamics, the reliability of those predictions and our mechanistic understanding of population dynamics is influenced by whether linear or non-linear approaches are taken. Linear models may be more appropriate for determining response of the whole system to changes in climate, however, linear models assume that the strength of density dependence and, thus, the susceptibility of the population to climate change, are independent of density (sensu [6]). The non-linear model presented here, however, suggests that the nature and strength of density dependence may be non-linear (i.e. differs between high- and low-density regimes) (Table 1) and, because of this, the manner in which the population responds to climate change may differ depending on if the population is at relatively high or relatively low densities. By specifically examining the possibility of non-linearity, we were able to suggest mechanisms through which changes in wolf behavior may influence population dynamics and their response to changes in climate. Thus, we suggest that because density dependence determines, in part, population response to climate change, it is important to investigate methods that allow the detection of non-linearity in the nature and strength of this force and other potentially important processes.

Methods

Developing the model

Wolves on Isle Royale are the top predator in a straight-chain, three trophic-level system including moose (Alces alces) and balsam fir (Abies balsamea) [16]. Therefore, the dynamics of the wolves in this system are most broadly influenced by trophic interactions with moose and intrinsic density dependence. Because wolves prey primarily on young and old moose, their population dynamics may depend more on the density of these age groups rather than the total density of moose. Unfortunately, there is very little information regarding the age structure of the moose population. In addition, estimates of moose density are based upon cohort reconstruction and, therefore, may be less accurate than those of wolf density, which are based on direct counts. For these reasons, in addition to previous studies suggesting that a uni-variate equation in delay coordinates may indirectly incorporate trophic interactions [17], we begin with a simple single-species linear autogressive model:

X= α0 + α1X+ α2X+ .... + αX+ ε    (1)

where Xis loge-transformed wolf density at time i and εis white noise. Although this simplification compromises detail in terms of the ecology behind wolf dynamics, it may facilitate further manipulations while, at some level, incorporating trophic interactions between wolves and moose (see, e.g. [13,18]).

Partial autocorrelation function (PACF) plots and recent first- and second-order linear and piecewise non-linear modeling of the wolf time series indicate that a first-order model provides a better fit in terms of statistical parsimony than second-order or higher models [9,12]; A. Ellis, unpublished results) giving:

X= α0 + α1X+ ε    (2)

To examine the structure of density dependence we set X- X= R, simplify coefficients, and express equation (2) in terms of the population growth rate of wolves:

R= β0 + β1X+ ε    (3)

The scatter plot of Rversus Xfit with a spline function (Fig. 1) indicates that the relationship between population growth rate in the current year and density in the previous year may be non-linear. While this inference is supported by the reduction in residual error achieved by applying a generalized additive model [19] with a non-linear spline function (deviance = 1.44) compared to a linear regression (deviance = 2.19) between Rand lagged density, more data points are needed to verify the robustness of this claim. For the purposes of this study, this plot suggests that the strength of density dependence may differ between relatively small and relatively large populations around some threshold (in X) and that this non-linearity may appropriately be expressed with simple linear piecewise models.

A first-order self-excitatory threshold autoregressive (SETAR) model with a low- and high-density regime and a threshold in Xhas the general form:

X= a1,0 + a1,1X+ ε1,    X< θ     (4)

X= a2,0 + a2,1X+ ε2,    X≥ θ

where, Xand Xare loge-transformed densities of the current and previous years, respectively, aare statistical parameters (i = 1 and 2 corresponds to lower and upper density regimes, respectively, j = 0 and 1 correspond to the constant and lag-1 coefficients, respectively), εis white noise (i = 1 and 2 corresponds to lower and upper density regimes, respectively), and θ is the threshold set in terms of density [7,20].

Covariates

Clearly, wolf dynamics are not determined exclusively by the direct influence of density dependence and the indirect influence of moose density. Wolf dynamics may also be influenced by social structure [9,12,21], local climate [22], and large-scale climate [9]. There is, for example, typically one breeding pair per pack [23], so that the birth rate of the population depends on the number of packs in the population. A closely related aspect of social structure, pack size, may also influence dynamics because pack size may affect ability to capture prey [9,12]. Climate, specifically snow depth and the North Atlantic Oscillation (NAO) [24], may also influence wolf dynamics through its influence on prey vulnerability and kill rate [8,22]. Because snow depth itself is influenced by winter temperature (R2 = 0.14 for the period of 1959 to 1999), and because the NAO may represent the vehicle through which large-scale warming occurs [25], it is reasonable to suggest that winter temperature may also influence wolf dynamics on Isle Royale though, to our knowledge, this influence has not been examined.

Hence, we examined all possible non-linear SETAR models of the form specified in equation (4) with covariates of most-recent lags including pack size of the previous year, number of packs of the previous year, cumulative winter snow depth of the current year (Ontario Climate Center), the Northern Hemisphere winter temperature anomaly of the current year, and the North Atlantic Oscillation of the current year [26]. Because estimates of cumulative snowfall and temperature on Isle Royale are not available, we used snowfall measurements from nearby Thunder Bay, Ontario and the Northern Hemisphere winter temperature anomaly.

Threshold selection

Before obtaining parameter estimates and comparing models, we determined the optimal threshold for each possible model. The optimal threshold for a particular model based upon the skeleton (Eq. 4) is the lag-one density term (X) from the observed unique values of Xthat minimizes the negative log-likelihood value for the entire model. This was determined by analyzing each model with θ in every unique value of X, except for the lowest and highest values to ensure minimum degrees of freedom to perform a linear regression. The threshold that minimized the negative log-likelihood was the optimal threshold for that particular model.

Model selection

All possible SETAR models with the form of equation (4) and the covariates mentioned above were analyzed to find the most statistically parsimonious model. For comparison, all possible linear models of the form of equation (2) and the covariates mentioned above were also analyzed and the most parsimonious model selected. We examined models with the separate effects of snow, temperature, or the NAO and models that included two or three climate variables to investigate the interaction between them. The most parsimonious model was that which minimized the corrected Akaike's Information Criterion (AIC) score [27]. Selected models were then tested for multicollinearity. Pairs of covariates resulting in variance inflation factors greater than 10 were not simultaneously included in the model, excluding the covariate with lower explanatory power when necessary.

Stability

Stability theory predicts that the greater the strength of density dependence, the more likely a population is to remain within a range of densities despite environmental change [28]. We inferred the stability of the system on the basis of the value of the coefficient of direct density dependence (a1,1, a2,1 in Eq. (4) and α1 in Eq. (2)) in the regime where the equilibrium point lies [18,29,30].

Authors' contributions

AE designed the study, conducted the analyses, and did most of the writing. EP conceived of the study, supervised its design, and contributed to the writing. Both authors read and approved of the final manuscript.