Inference for nonlinear epidemiological models using genealogies and time series

Phylodynamics - the field aiming to quantitatively integrate the ecological and evolutionary dynamics of rapidly evolving populations like those of RNA viruses - increasingly relies upon coalescent approaches to infer past population dynamics from reconstructed genealogies. As sequence data have become more abundant, these approaches are beginning to be used on populations undergoing rapid and rather complex dynamics. In such cases, the simple demographic models that current phylodynamic methods employ can be limiting. First, these models are not ideal for yielding biological insight into the processes that drive the dynamics of the populations of interest. Second, these models differ in form from mechanistic and often stochastic population dynamic models that are currently widely used when fitting models to time series data. As such, their use does not allow for both genealogical data and time series data to be considered in tandem when conducting inference. Here, we present a flexible statistical framework for phylodynamic inference that goes beyond these current limitations. The framework we present employs a recently developed method known as particle MCMC to fit stochastic, nonlinear mechanistic models for complex population dynamics to gene genealogies and time series data in a Bayesian framework. We demonstrate our approach using a nonlinear Susceptible-Infected-Recovered (SIR) model for the transmission dynamics of an infectious disease and show through simulations that it provides accurate estimates of past disease dynamics and key epidemiological parameters from genealogies with or without accompanying time series data.     

Estimation of growth regulation in natural populations by extended family of growth curve models with fractional order derivative: Case studies from the global population dynamics database

Estimating the trend in population time series data using growth curve models is a central idea in population ecology. Several models, mainly governed by differential or difference equations, have been applied to real data sets to identify general growth pattern and make predictions. In this article, we analyze ecological time series data by fitting mathematical models governed by fractional differential equations (FDE). The order of the FDE (α) is used to quantify the evidence of memory in the population processes. The application of FDE is exemplified by analyzing time series data on two bird species Phalacrocorax carbo (Great cormorant) and Parus bicolor (Tufted titmouse) and two mammal species Castor canadensis (Beaver) and Ursus americanus (American black bear) extracted from the global population dynamics database. Five different population growth models were fitted to these data; density-independent exponential, negative density-dependent logistic and θ-logistic model, positive density-dependent exponential Allee and strong Allee model. Both ordinary and fractional derivative representations of these models were fitted to the time series data. Markov chain Monte Carlo (MCMC) method was used to estimate the model parameters and Akaike information criterion was used to select the best model. By estimating the return rate for each of the time series, we have shown that populations governed by FDE with a small value of α (high level of memory) return to the stable equilibrium faster. This demonstrates a synergistic interplay between memory and stability in natural populations.     

Analysis of cyclic fluctuations in larch bud moth populations by means of discrete-time dynamic models

Analysed are the data of larch bud moth (Zeiraphera diniana Gn.) fluctuations in Swiss Alps. The analysis applies simplest mathematical models of isolated population dynamics (in particular, Kostitzin model, Skellam model, the discrete logistic model, and some other ones), which include the minimal number of unknown parameters. The parameters have been estimated, for all the models in hand, by the least-squares method, to fit certain data from the Global Population Dynamics Database (N 1407 and N 6195), the sequences of the data deviations from the model trajectories being treated as well. The best approximations are shown to be achieved with Moran-Ricker model and the discrete logistic model. Statistical criteria (Kolmogorov-Smirnov and Shapiro-Wilk tests) reveal that the hypotheses of normal distribution of residuals must be rejected for one of the time series (N 1407); some models demonstrate serial correlations in the sequence of residuals (according to Durbin-Watson test). This leads to the conclusion that periodic fluctuations in the larch bud moth population (N 1407) can hardly be explained by self-regulation mechanisms alone. For another time series (N 6195), the modified discrete logistic model has appeared to be acceptable as a mode of fluctuations.     

A statistical approach to quasi-extinction forecasting

Forecasting population decline to a certain critical threshold (the quasi-extinction risk) is one of the central objectives of population viability analysis (PVA), and such predictions figure prominently in the decisions of major conservation organizations. In this paper, we argue that accurate forecasting of a population's quasi-extinction risk does not necessarily require knowledge of the underlying biological mechanisms. Because of the stochastic and multiplicative nature of population growth, the ensemble behaviour of population trajectories converges to common statistical forms across a wide variety of stochastic population processes. This paper provides a theoretical basis for this argument. We show that the quasi-extinction surfaces of a variety of complex stochastic population processes (including age-structured, density-dependent and spatially structured populations) can be modelled by a simple stochastic approximation: the stochastic exponential growth process overlaid with Gaussian errors. Using simulated and real data, we show that this model can be estimated with 20-30 years of data and can provide relatively unbiased quasi-extinction risk with confidence intervals considerably smaller than (0,1). This was found to be true even for simulated data derived from some of the noisiest population processes (density-dependent feedback, species interactions and strong age-structure cycling). A key advantage of statistical models is that their parameters and the uncertainty of those parameters can be estimated from time series data using standard statistical methods. In contrast for most species of conservation concern, biologically realistic models must often be specified rather than estimated because of the limited data available for all the various parameters. Biologically realistic models will always have a prominent place in PVA for evaluating specific management options which affect a single segment of a population, a single demographic rate, or different geographic areas. However, for forecasting quasi-extinction risk, statistical models that are based on the convergent statistical properties of population processes offer many advantages over biologically realistic models.     

Complexity is costly: a meta‐analysis of parametric and non‐parametric methods for short‐term population forecasting

Short‐term forecasts based on time series of counts or survey data are widely used in population biology to provide advice concerning the management, harvest and conservation of natural populations. A common approach to produce these forecasts uses time‐series models, of different types, fit to time series of counts. Similar time‐series models are used in many other disciplines, however relative to the data available in these other disciplines, population data are often unusually short and noisy and models that perform well for data from other disciplines may not be appropriate for population data. In order to study the performance of time‐series forecasting models for natural animal population data, we assembled 2379 time series of vertebrate population indices from actual surveys. Our data were comprised of three vastly different types: highly variable (marine fish productivity), strongly cyclic (adult salmon counts), and small variance but long‐memory (bird and mammal counts). We tested the predictive performance of 49 different forecasting models grouped into three broad classes: autoregressive time‐series models, non‐linear regression‐type models and non‐parametric time‐series models. Low‐dimensional parametric autoregressive models gave the most accurate forecasts across a wide range of taxa; the most accurate model was one that simply treated the most recent observation as the forecast. More complex parametric and non‐parametric models performed worse, except when applied to highly cyclic species. Across taxa, certain life history characteristics were correlated with lower forecast error; specifically, we found that better forecasts were correlated with attributes of slow growing species: large maximum age and size for fishes and high trophic level for birds. Synthesis Evaluating the data support for multiple plausible models has been an integral focus of many ecological analyses. However, the most commonly used tools to quantify support have weighted models’ hindcasting and forecasting abilities. For many applications, predicting the past may be of little interest. Concentrating only on the future predictive performance of time series models, we performed a forecasting competition among many different kinds of statistical models, applying each to many different kinds of vertebrate time series of population abundance. Low‐dimensional (simple) models performed well overall, but more complex models did slightly better when applied to time series of cyclic species (e.g. salmon).     

A comparative evaluation of models of cinnabar moth dynamics

Summary 1. A complex model of cinnabar moth dynamics proposed by Dempster and Lakhani (1979) with 23 parameters is reduced to a single equation with five parameters, and the behaviour of the reduced model shown to explain most features of the full model. 2. The efficiency of the full model is compared with the reduced model and with two even simpler models (the two parameter discrete logistic and a four parameter model based on a step-function for mortality) in their abilities to describe time series data of cinnabar moth population densities from Weeting Heath. Models with more parameters were not significantly better than few-parameter models in describing population trajectories. 3. Models that included a driving variable (in this case observed rainfall data) were no better at describing the data than simpler models without driving variables. It appears, therefore, that the routine inclusion of driving variables may be counterproductive, unless there is compelling empirical or theoretical evidence of their importance and the mode of action of the driving variables can be modelled mechanistically. For example, the regression model used to describe the relationship between rainfall and plant biomass in Dempster and Lakhani (1979), breaks down if rainfall is assumed to be constant, because there is no explicit model for the regulation of plant biomass. 4. The parameter values of the cinnabar-ragwort interaction suggest that cinnabar moth dynamics may be chaotic. Whether or not field data exhibit chaos or environmental stochasticity (or a mixture of both) is impossible to determine from inspection of time series data on population density. There is an urgent need for experimental and theoretical protocols to disentangle these two sources of population fluctuation.     

Sturdy cycles in the chaotic Tribolium castaneum data series

Because of the inherent discreteness of individuals, population dynamical models must be discrete variable systems. In case of strong nonlinearity, such systems interacting with noise can generate a great variety of patterns from nearly periodic behavior through complex combination of nearly periodic and chaotic patterns to noisy chaotic time series. The interaction of a population consisting of discrete individuals and demographic noise has been analyzed in laboratory population data Henson et al. (Science 294 (2001) 602; Proc. Roy. Soc. Ser. B 270 (2003) 1549). In this paper we point out that some of the cycles are fragile, i.e. they are sensitive to the discretization algorithm and to small variation of the model parameters, while others remain "sturdy" against the perturbations. We introduce a statistical algorithm to detect disjoint, nearly-periodic patterns in data series. We show that only the sturdy cycles of the discrete variable models appear in the data series significantly. Our analysis identified the quasiperiodic 11-cycle (emerging in the continuous model) to be present significantly only in one of the three experimental data series. Numerical simulations confirm that cycles can be detected only if noise is smaller than a certain critical level and population dynamics display the largest variety of nearly-periodic patterns if they are on the border of "grey" and "noisy" regions, defined in Domokos and Scheuring (J. Theor. Biol. 227 (2004) 535).     

Likelihood based population viability analysis in the presence of observation error

Population viability analysis (PVA) entails calculation of extinction risk, as defined by various extinction metrics, for a study population. These calculations strongly depend on the form of the population growth model and inclusion of demographic and/or environmental stochasticity. Form of the model and its parameters are determined based on observed population time series data. A typical population time series, consisting of estimated population sizes, inevitably has some observation error and likely has missing observations. In this paper, we present a likelihood based PVA in the presence of observation error and missing data. We illustrate the importance of incorporation of observation error in PVA by reanalyzing the population time series of song sparrow Melospiza melodia on Mandarte Island, British Columbia, Canada from 1975–1998. Using Akaike information criterion we show that model with observation error fits the data better than the one without observation error. The extinction risks predicted by with and without observation error models are quite different. Further analysis of possible causes for observation error revealed that some component of the observation error might be due to unreported dispersal. A complete analysis of such data, thus, would require explicit spatial models and data on dispersal along with observation error. Our conclusions are, therefore, two‐fold: 1) observation errors in PVA matter and 2) integrating these errors in PVA is not always enough and can still lead to important biases in parameter estimates if other processes such as dispersal are ignored.     

A cure for the plague of parameters: constraining models of complex population dynamics with allometries

A major goal of ecology is to discover how dynamics and structure of multi-trophic ecological communities are related. This is difficult, because whole-community data are limited and typically comprise only a snapshot of a community instead of a time series of dynamics, and mathematical models of complex system dynamics have a large number of unmeasured parameters and therefore have been only tenuously related to real systems. These are related problems, because long time-series, if they were commonly available, would enable inference of parameters. The resulting ‘plague of parameters’ means most studies of multi-species population dynamics have been very theoretical. Dynamical models parametrized using physiological allometries may offer a partial cure for the plague of parameters, and these models are increasingly used in theoretical studies. However, physiological allometries cannot determine all parameters, and the models have also rarely been directly tested against data. We confronted a model of community dynamics with data from a lake community. Many important empirical patterns were reproducible as outcomes of dynamics, and were not reproducible when parameters did not follow physiological allometries. Results validate the usefulness, when parameters follow physiological allometries, of classic differential-equation models for understanding whole-community dynamics and the structure–dynamics relationship.     

Bayesian meta-analysis of demographic parameters in three small, temperate passerines

Accurate estimates of population parameters are vital for estimating extinction risk. Such parameters, however, are typically not available for threatened populations. We used a recently developed software tool based on Markov Chain Monte Carlo methods for carrying out Bayesian inference (the BUGS package) to estimate four demographic parameters; the intrinsic growth rate, the strength of density dependence, and the demographic and environmental variance, in three species of small temperate passerines from two sets of time series data taken from a dipper and a song sparrow population, and from previously obtained frequentist estimates of the same parameters in the great tit. By simultaneously modeling variation in these demographic parameters across species and using the resulting distributions as priors in the estimation for individual species, we improve the estimates for each individual species. This framework also allows us to make probabilistic statements about plausible parameter values for small passerines temperate birds in general which is often critically needed in management of species for which little or no data are available. We also discuss how our work relates to recently developed theory on dynamic stochastic population models, and finally note some important differences between frequentist and Bayesian methods.     

Extinction risk of a meta-population: aggregation approach

Aggregation of variables of a complex mathematical model with realistic structure gives a simplified model which is more suitable than the original one when the amount of data for parameter estimation is limited. Here we explore use of a formula derived for a single unstructured population (canonical model) in predicting the extinction time for a population living in multiple habitats. In particular we focus multiple populations each following logistic growth with demographic and environmental stochasticities, and examine how the mean extinction time depends on the migration and environmental correlation. When migration rate and/or environmental correlation are very large or very small, we may express the mean extinction time exactly using the formula with properly modified parameters. When parameters are of intermediate magnitude, we generate a Monte Carlo time series of the population size for the realistic structured model, estimate the "effective parameters" by fitting the time series to the canonical model, and then calculate the mean extinction time using the formula for a single population. The mean extinction time predicted by the formula was close to those obtained from direct computer simulation of structured models. We conclude that the formula for an unstructured single-population model has good approximation capability and can be applicable in estimating the extinction risk of the structured meta-population model for a limited data set.     

Stochastic analogues of deterministic single-species population models

Although single-species deterministic difference equations have long been used in modeling the dynamics of animal populations, little attention has been paid to how stochasticity should be incorporated into these models. By deriving stochastic analogues to difference equations from first principles, we show that the form of these models depends on whether noise in the population process is demographic or environmental. When noise is demographic, we argue that variance around the expectation is proportional to the expectation. When noise is environmental the variance depends in a non-trivial way on how variation enters into model parameters, but we argue that if the environment affects the population multiplicatively then variance is proportional to the square of the expectation. We compare various stochastic analogues of the Ricker map model by fitting them, using maximum likelihood estimation, to data generated from an individual-based model and the weevil data of Utida. Our demographic models are significantly better than our environmental models at fitting noise generated by population processes where noise is mainly demographic. However, the traditionally chosen stochastic analogues to deterministic models--additive normally distributed noise and multiplicative lognormally distributed noise--generally fit all data sets well. Thus, the form of the variance does play a role in the fitting of models to ecological time series, but may not be important in practice as first supposed.     

Incorporating movement into models of grey seal population dynamics

1. One of the most difficult problems in developing spatially explicit models of population dynamics is the validation and parameterization of the movement process. We show how movement models derived from capture-recapture analysis can be improved by incorporating them into a spatially explicit metapopulation model that is fitted to a time series of abundance data. 2. We applied multisite capture-recapture analysis techniques to photo-identification data collected from female grey seals at the four main breeding colonies in the North Sea between 1999 and 2001. The best-fitting movement models were then incorporated into state-space metapopulation models that explicitly accounted for demographic and observational stochasticity. 3. These metapopulation models were fitted to a 20-year time series of pup production data for each colony using a Bayesian approach. The best-fitting model, based on the Akaike Information Criterion (AIC), had only a single movement parameter, whose confidence interval was 82% less than that obtained from the capture-recapture study, but there was some support for a model that included an effect of distance between colonies. 4. The state-space modelling provided improved estimates of other demographic parameters. 5. The incorporation of movement, and the way in which it was modelled, affected both local and regional dynamics. These differences were most evident as colonies approached their carrying capacities, suggesting that our ability to discriminate between models should improve as the length of the grey seal time series increases.     

A general framework for modeling population abundance data

Time-series data resulting from surveying wild animals are often described using state-space population dynamics models, in particular with Gompertz, Beverton-Holt, or Moran-Ricker latent processes. We show how hidden Markov model methodology provides a flexible framework for fitting a wide range of models to such data. This general approach makes it possible to model abundance on the natural or log scale, include multiple observations at each sampling occasion and compare alternative models using information criteria. It also easily accommodates unequal sampling time intervals, should that possibility occur, and allows testing for density dependence using the bootstrap. The paper is illustrated by replicated time series of red kangaroo abundances, and a univariate time series of ibex counts which are an order of magnitude larger. In the analyses carried out, we fit different latent process and observation models using the hidden Markov framework. Results are robust with regard to the necessary discretization of the state variable. We find no effective difference between the three latent models of the paper in terms of maximized likelihood value for the two applications presented, and also others analyzed. Simulations suggest that ecological time series are not sufficiently informative to distinguish between alternative latent processes for modeling population survey data when data do not indicate strong density dependence.     

Data estimation and the colour of time series.

There has been a long debate on the source of temporal fluctuations in natural population densities. The difficulty is that unpredictable irregularities might be attributed either to external environmental factors or to chaotic dynamics of populations, or even to the interaction of these two factors. Some years ago Cohen (1995) pointed out that real time series follow redshifted Fourier power spectra, while the simplest chaotic population dynamical models are mostly blueshifted. Since then, the controversy has focused on comparisons of Fourier spectra originating from different models and data. Here, we show experimentally that estimation process by human observers shifts power spectra to the red. This result implies that because of estimation distortion, real population data must be less redshifted than many recorded time series suggest.     

Extinction risk and the 1/f family of noise models

In order to predict extinction risk in the presence of reddened, or correlated, environmental variability, fluctuating parameters may be represented by the family of 1/f noises, a series of stochastic models with different levels of variation acting on different timescales. We compare the process of parameter estimation for three 1/f models (white, pink and brown noise) with each other, and with autoregressive noise models (which are not 1/f noises), using data from a model time-series (length, T) of population. We then calculate the expected increase in variance and the expected extinction risk for each model, and we use these to explore the implication of assuming an incorrect noise model. When parameterising these models, it is necessary to do so in terms of the measured ("sample") parameters rather than fundamental ("population") parameters. This is because these models are non-stationary: their parameters need not stabilize on measurement over long periods of time and are uniquely defined only over a specified "window" of timescales defined by a measurement process. We find that extinction forecasts can differ greatly between models, depending on the length, T, and the coefficient of variability, CV, of the time series used to parameterise the models, and on the length of time into the future which is to be projected. For the simplest possible models, ones with population itself the 1/f noise process, it is possible to predict the extinction risk based on CV of the observed time series. Our predictions, based on explicit formulae and on simulations, indicate that (a) for very short projection times relative to T, brown and pink noise models are usually optimistic relative to equivalent white noise model; (b) for projection timescales equal to and substantially greater than T, an equivalent brown or pink noise model usually predicts a greater extinction risk, unless CV is very large; and (c) except for very small values of CV, for timescales very much greater than T, the brown and pink models present a more optimistic picture than the white noise model. In most cases, a pink noise is intermediate between white and brown models. Thus, while reddening of environmental noise may increase the long-term extinction probability for stationary processes, this is not generally true for non-stationary processes, such as pink or brown noises.     

Hyperparameter estimation using stochastic approximation with application to population pharmacokinetics

A stochastic approximation algorithm is proposed for recursive estimation of the hyperparameters characterizing, in a population, the probability density function of the parameters of a statistical model. For a given population model defined by a parametric model of a biological process, an error model, and a class of densities on the set of the individual parameters, this algorithm provides a sequence of estimates from a sequence of individuals' observation vectors. Convergence conditions are verified for a class of population models including usual pharmacokinetic applications. This method is implemented for estimation of pharmacokinetic population parameters from drug multiple-dosing data. Its estimation capabilities are evaluated and compared to a classical method in population pharmacokinetics, the first-order method (NONMEM), on simulated data.     

Assessing the spatial variability of density dependence in waterfowl populations

Population models commonly assume that the demographic parameters are spatially invariant, but there is considerable evidence that population growth rate (r) and the strength of density dependence (β) can vary over a species' range. To address this issue we developed a spatially explicit Gompertz population model based on the spatially varying coefficients approach to assess the spatial variation in population drivers. The model was fit to spatially stratified time series population estimates of the mallard Anas platyrhynchos in western North America. We included precipitation during the previous year and spring maximum temperature in the current year as environmental factors in the density dependent population model. Because density dependent models can give biased estimates for time series of abundance data, we fit a naïve model without informative priors and a model where we constrained the mean and variance of r to biologically realistic values that were derived via a comparative demography approach. In the naïve model, r and β were not separately identifiable and their values were overestimated, leading to unrealistic population growth. The naïve model also implied spatial variation in population r and the return time to equilibrium [1?(– β)] across the survey area. In contrast, in the informative model, r and the return time to equilibrium did not vary markedly among populations and were generally equal across populations. The effects of the climatic factors were similar across models. Population growth rates in the Prairie‐pothole region were positively correlated with precipitation, while in Alaska rates were positively correlated with spring temperature. Although it has been argued in the past that adding ecological realism could help avoid the pitfalls associated with density dependent models, our results demonstrate that imposing constraints on the population parameters is still the best course of action.     

Time series analysis of the epidemiological transition in Minorca, 1634-1997

Autoregressive integrated moving average (ARIMA) models provide a powerful tool for detecting seasonal patterns in mortality statistics. The strength of ARIMA models lies in their ability to reveal complex structures of temporal interdependence in time series. Moreover, changes in model parameters provide an empirical basis for detecting secular trends and death seasonality patterns. This approach is illustrated by our analysis of changes in the mortality patterns of the population of the town of Es Mercadal on the island of Minorca between 1634 and 1997. These data reveal a transition from an early mortality pattern requiring a complex ARIMA model that accounts for a strong seasonal death pattern and periodic epidemic-related mortality crises to a much simpler 20th-century pattern that can be described by a simple single-parameter ARIMA model. These same data were also analyzed using standard seasonality tests. The results show that the reduction in the number of parameters required to fit the Es Mercadal mortality data coincides with the epidemiological transition in which the predominant causes of morbidly and mortality shift from infectious to degenerative causes.     

Semiparametric models and inference for biomedical time series with extra-variation

Biomedical trials often give rise to data having the form of time series of a common process on separate individuals. One model which has been proposed to explain variations in such series across individuals is a random effects model based on sample periodograms. The use of spectral coefficients enables models for individual series to be constructed on the basis of standard asymptotic theory, whilst variations between individuals are handled by permitting a random effect perturbation of model coefficients. This paper extends such methodology in two ways: first, by enabling a nonparametric specification of underlying spectral behaviour; second, by addressing some of the tricky computational issues which are encountered when working with this class of random effect models. This leads to a model in which a population spectrum is specified nonparametrically through a dynamic system, and the processes measured on individuals within the population are assumed to have a spectrum which has a random effect perturbation from the population norm. Simulation studies show that standard MCMC algorithms give effective inferences for this model, and applications to biomedical data suggest that the model itself is capable of revealing scientifically important structure in temporal characteristics both within and between individual processes.