Predicting the process of extinction in experimental microcosms and accounting for interspecific interactions in single-species time series

Predicting population extinction risk is a fundamental application of ecological theory to the practice of conservation biology. Here, we compared the prediction performance of a wide array of stochastic, population dynamics models against direct observations of the extinction process from an extensive experimental data set. By varying a series of biological and statistical assumptions in the proposed models, we were able to identify the assumptions that affected predictions about population extinction. We also show how certain autocorrelation structures can emerge due to interspecific interactions, and that accounting for the stochastic effect of these interactions can improve predictions of the extinction process. We conclude that it is possible to account for the stochastic effects of community interactions on extinction when using single-species time series.

The density dependence can also be placed in the fecundity term, with survival as a constant, leading 1 to a similar form for the mean and variance terms. It is also important to note that the terms λ 2 and p 0 are not identifiable under the Poisson assumption for demographic stochasticity, φ(λ) 2 = λ.

3
Under the Poisson assumption V [N t | N t−1 = n t−1 ] = λn t−1 p t−1 = λp 0 n t−1 p(n t−1 ). Therefore, 4 we explored whether the parameter p 0 was identifiable in our data when the Poisson assumption 5 is relaxed by fitting the joint-likelihood profile of p 0 and λ for the first microcosm experimental 6 dataset (the likelihoods are fully described in Appendix C). As in the main manuscript we let 7 the demographic stochasticity in reproduction be a constant, φ(λ) 2 = φ 2 . We found that λ and 8 p 0 were collinear ( Figure A1) showing that it is not possible to jointly estimate both parameters 9 from abundance data. Therefore we follow Melbourne & Hastings (2008) and treat λp 0 as a single 10 parameter r. This leads to the following approximation for the demographic variance which ignores 11 the density independent contribution in the binomial variance term: This approximation can lead to estimates that predict higher or lower variance than the true value The lack of a well defined peak indicates that the parameters r and p 0 are collinear and are not independently identifiable.
the mean and variance conditional on the state of the environment, W t , are: Finally, averaging over the environmental state W t gives a model with both demographic and envi-2 ronmental stochasticity. The mean is given by = rn t−1 p(n t−1 ), and the variance by where we used the approximation for the demographic term described above.

2
For estimation purposes in the manuscript, we assumed that the demographic variance in 3 reproduction is proportional to the yearly realized reproductive value, in equation A4, we can then plugin πλ . The variance is then given by, ≈ n t−1 rp(n t−1 )(1 − p(n t−1 )) + n t−1 πrp(n t−1 ) 2 + n 2 t−1 p(n t−1 ) 2 ϕ 2 , where we have used the approximation due to the nonidentifiability of the density independent 6 survival that was used above. Additionally, π and r in the demographic stochasticity are non-7 identifiable. Therefore we define the term, φ 2 ≡ rπ. We then get the form of the variance used in The dynamics of a single species in a linear interacting system can be re-expressed through a combi-2 nation of lagged observations of the population abundance and lagged error terms (Royama, 1981).

3
Here we show how to obtain these dynamics in a two-species system by substituting the dynam-4 ics of one species into the other. Beginning with the stochastic discrete-time Gompertz model, , we can re-express the dynamics as a linear function with  This gives strength of density dependence is then given by the parameter a which for a stationary process can 9 go between −1 and 1, with a = 0 corresponding to a random walk. The range for the parameter b 10 is then given by 0 ≤ b ≤ 2. Consider a linear interacting stochastic system with this formulation, 11 written as: The error terms are independent nonidentical-normal distributions such that: ε n (t) ∼ 13 Norm(0, σ 2 n (t)). We remove the constants c 1 and c 2 and reintroduce them later in order to sim- dynamics of X 1 (t) or X 2 (t) through a substitution process. The substitution for X 1 (t) is as follows, 16 first solve equation (B1) for X 2 (t − 1) in terms of the X 1 (t)'s and plug this into equation (B2): Note that this procedure implies that b 12 = 0 or else there is no solution. Now, plugging 18 this result back into equation (B1), then collecting terms gives Jenkins (1970)), the process can be rewritten as: The right hand side of equation B4 is a new moving average process, ξ, with unknown moving 10 average parameter, θ(t) and variance, σ 2 ξ (t − 1). The variance and autocovariance of left-hand side 11 of equation B4 are known and given by 12 while the variance and autocovariance of the right hand side is given by We can then equate the autocovariances to give Plugging this in for σ 2 ξ (t) in the variances, We can then solve for θ(t) using the quadratic formula, where B(t) = (1 + b 2 22 )σ 2 1 (t) + b 12 σ 2 2 (t). For estimation purposes we assume θ(t) = θ. This is likely to 3 be a reasonable assumption when demographic stochasticity is low but becomes worse as it increases.

4
A time-dependency also enters the variance of the process σ 2 ξ (t) which is also dependent on σ 2 1 (t) and assumptions we apply this model as a rough, but potentially useful, approximation to the complex 9 stochastic dynamics that can be induced through predator-prey interactions. 1 We fit gamma, log-normal, and negative binomial distributions to the data by setting mean and 2 variance for the distributions to expressions derived in Appendix A. Here we use some shortcut 3 notation for convenience, We fit all models by matching the mean and variance terms in (C6) to the distribution mean 5 and variance. For the lognormal distribution this leads to: Similarly we fit the gamma distribution to the abundance transitions with shape (k) and scale (θ) 7 parameters given by The negative binomial distribution was also fit to the abundance data. This model was parameterized 9 by the mean (µ NB ) and size (k) parameters, In the negative binomial the demographic variance is assumed to be equal to the growth rate, 11 φ(λ) 2 = r, while the size parameter, k, is related the environmental variance.

12
We used the one-step formulations of the AR, MA, and ARMA processes (Shumway & 13 Stoffer, 2006). An exploratory analysis suggested that the autocorrelation processes operated on 14 the scale of the population growth rate, ln Nt nt−1 , therefore we included the autocorrelation as multiplicative processes such that they were additive to the density dependence on the log-scale.

1
This analyses also suggested that log-population abundances performed better than abundances so 2 all autocorrelation models are formulated with this transformation. All models contained an AR(1) 3 contribution through population growth and regulation processes contained in rn t−1 p(n t−1 ). The 4 AR(2) process was modeled as a contribution to the population growth rate given by ψ ln(n t−2 ).

5
The MA(1) process was also modeled as a contribution to the population growth rate, this term was 6 modeled as, θ (ln(n t−1 ) − ln(µ t−1 )). The mean of the process given in (C6) with the full ARMA 7 process therefore becomes The ARMA processes are additive combinations of the AR and MA components so we can remove 9 the AR(2) term by setting ψ = 0 and remove the MA(1) term by setting θ = 0 to examine models 10 nested to the ARMA(2,1). Although even higher order ARMA processes can be fit to the data using 11 similar methods, we limited ourselves to ARMA(2,1) processes due to the patterns observed in the 12 empirical autocorrelation function of the population growth rate ( Figure C2).

13
The likelihoods of the abundances for these data were calculated using the transition prob- . We also note that the logistic 1 form of density dependence, p(n t−1 ) = 1 − bn t−1 , was constrained in the optimization by applying 2 a constraint to the likelihood such that 0 < p(n t−1 ) < 1. The ∆AIC values for the abundances are 3 presented in Table C1.