Estimation methods for nonlinear state-space models in ecology

The use of nonlinear state-space models for analyzing ecological systems is increasing. A wide range of estimation methods for such models are available to ecologists, however it is not always clear, which is the appropriate method to choose. To this end, three approaches to estimation in the theta logistic model for population dynamics were benchmarked by Wang (2007). Similarly, we examine and compare the estimation performance of three alternative methods using simulated data. The first approach is to partition the state-space into a finite number of states and formulate the problem as a hidden Markov model (HMM). The second method uses the mixed effects modeling and fast numerical integration framework of the AD Model Builder (ADMB) open-source software. The third alternative is to use the popular Bayesian framework of BUGS. The study showed that state and parameter estimation performance for all three methods was largely identical, however with BUGS providing overall wider credible intervals for parameters than HMM and ADMB confidence intervals.     

The Extended Kalman Filter (EKF) as a tool for the assimilation of high frequency water quality data

The Extended Kalman Filter (EKF) was applied to the analysis of high frequency field measurements of dissolved oxygen (DO), water temperature, salinity, collected by multiparametric sensors in the lagoon of Venice. This paper focuses on the practical aspects of the implementation of the EKF as a data assimilation technique and does not deal with the problems associated with the identification of the model. In this regard, the EKF has proved to be a useful tool for the updating of the estimates of the parameters of a simple DO-chlorophyll model, which can be used for linking the high frequency data to meteorological forcings, such as solar radiation and wind, and to other low frequency measurements of water quality parameters, such as the concentrations of Chlorophyll a and nutrients. The model can subsequently be used as a tool for checking the consistency of all this data, and may also be employed for controlling the quality of the data collected by the multiparametric sensors.     

Estimation of population parameters from aerial counts of North American mallards: a cautionary tale.

The accuracy of population estimates strongly interferes with our ability to obtain unbiased estimates of population parameters based on analyses of time series of population fluctuations. Here we use long-term data on fluctuations in the size of Mallard populations collected as part of the May Breeding Waterfowl Survey covering a large section of North America. We assume a log-linear model of density dependence and use a hierarchical Bayesian state-space approach in which all parameters are assumed to be realizations from a common underlying distribution. Thus, parameters for different populations are not allowed to vary independently of each other. We then simulated independent time series of aerial counts, using the estimated parameters and adding various levels of observation error. These simulations showed that the estimates of stochastic population growth rate and strength of density dependence were biased even when moderate sampling errors were present. In contrast, the estimates of the environmental stochasticity and the carrying capacity were unbiased even for short time series and large observation error. Our results underline the importance of reducing the magnitude of sampling error in the design of large-scale monitoring programs of population fluctuations.     

Field-scale roles of density, temperature, nitrogen, and predation on aphid population dynamics

Robust analyses of noisy, stage-structured, irregularly spaced, field-scale data incorporating multiple sources of variability and nonlinear dynamics remain very limited, hindering understanding of how small-scale studies relate to large-scale population dynamics. We used a novel, complementary Bayesian and frequentist state-space model analysis to ask how density, temperature, plant nitrogen, and predators affect cotton aphid (Aphis gossypii) population dynamics in weekly data from 18 field-years and whether estimated effects are consistent with small-scale studies. We found clear roles of density and temperature but not of plant nitrogen or predators, for which Bayesian and frequentist evidence differed. However, overall predictability of field-scale dynamics remained low. This study demonstrates stage-structured state-space model analysis incorporating bottom-up, top-down, and density-dependent effects for within-season (nearly continuous time), nonlinear population dynamics. The analysis combines Bayesian posterior evidence with maximum-likelihood estimation and frequentist hypothesis testing using average one-step-ahead residuals.     

Fitting complex population models by combining particle filters with Markov chain Monte Carlo

We show how a recent framework combining Markov chain Monte Carlo (MCMC) with particle filters (PFMCMC) may be used to estimate population state-space models. With the purpose of utilizing the strengths of each method, PFMCMC explores hidden states by particle filters, while process and observation parameters are estimated using an MCMC algorithm. PFMCMC is exemplified by analyzing time series data on a red kangaroo (Macropus rufus) population in New South Wales, Australia, using MCMC over model parameters based on an adaptive Metropolis-Hastings algorithm. We fit three population models to these data; a density-dependent logistic diffusion model with environmental variance, an unregulated stochastic exponential growth model, and a random-walk model. Bayes factors and posterior model probabilities show that there is little support for density dependence and that the random-walk model is the most parsimonious model. The particle filter Metropolis-Hastings algorithm is a brute-force method that may be used to fit a range of complex population models. Implementation is straightforward and less involved than standard MCMC for many models, and marginal densities for model selection can be obtained with little additional effort. The cost is mainly computational, resulting in long running times that may be improved by parallelizing the algorithm.     

Population growth in snow geese: a modeling approach integrating demographic and survey information

There are few analytic tools available to formally integrate information coming from population surveys and demographic studies. The Kalman filter is a procedure that facilitates such integration. Based on a state-space model, we can obtain a likelihood function for the survey data using a Kalman filter, which we may then combine with a likelihood for the demographic data. In this paper, we used this combined approach to analyze the population dynamics of a hunted species, the Greater Snow Goose (Chen caerulescens atlantica), and to examine the extent to which it can improve previous demographic population models. The state equation of the state-space model was a matrix population model with fecundity and regression parameters relating adult survival and harvest rate estimated in a previous capture-recapture study. The observation equation combined the output from this model with estimates from an annual spring photographic survey of the population. The maximum likelihood estimates of the regression parameters from the combined analysis differed little from the values of the original capture-recapture analysis, though their precision improved. The model output was found to be insensitive to a wide range of coefficient of variation (CV) in fecundity parameters. We found a close match between the surveyed and smoothed population size estimates generated by the Kalman filter over an 18-year period, and the estimated CV of the survey (0.078-0.150) was quite compatible with its assumed value (approximately 0.10). When we used the updated parameter values to predict future population size, the model underestimated the surveyed population size by 18% over a three-year period. However, this could be explained by a concurrent change in the survey method. We conclude that the Kalman filter is a promising approach to forecast population change because it incorporates survey information in a formal way compared with ad hoc approaches that either neglect this information or require some parameter or model tuning.     

Nonlinear relationships between vital rates and state variables in demographic models

To accurately estimate population dynamics and viability, structured population models account for among-individual differences in demographic parameters that are related to individual state. In the widely used matrix models, such differences are incorporated in terms of discrete state categories, whereas integral projection models (IPMs) use continuous state variables to avoid artificial classes. In IPMs, and sometimes also in matrix models, parameterization is based on regressions that do not always model nonlinear relationships between demographic parameters and state variables. We stress the importance of testing for nonlinearity and propose using restricted cubic splines in order to allow for a wide variety of relationships in regressions and demographic models. For the plant Borderea pyrenaica, we found that vital rate relationships with size and age were nonlinear and that the parameterization method had large effects on predicted population growth rates, X (linear IPM, 0.95; nonlinear IPMs, 1.00; matrix model, 0.96). Our results suggest that restricted cubic spline models are more reliable than linear or polynomial models. Because even weak nonlinearity in relationships between vital rates and state variables can have large effects on model predictions, we suggest that restricted cubic regression splines should be considered for parameterizing models of population dynamics whenever linearity cannot be assumed.     

Performance of a Bayesian state-space model of semelparous species for stock-recruitment data subject to measurement error

Measurement errors in spawner abundance create problems for fish stock assessment scientists. To deal with measurement error, we develop a Bayesian state-space model for stock-recruitment data that contain measurement error in spawner abundance, process error in recruitment, and time series bias. Through extensive simulations across numerous scenarios, we compare the statistical performance of the Bayesian state-space model with that of standard regression for a traditional stock-recruitment model that only considers process error. Performance varies depending on the information content in data, as determined by stock productivity, types of harvest situations, and amount of measurement error. Overall, in terms of estimating optimal spawner abundance S_MSY, the Ricker density-dependence parameter β, and optimal harvest rate h_MSY, the Bayesian state-space model works best for informative data from low and variable harvest rate situations for high-productivity salmon stocks. The traditional stock-recruitment model (TSR) may be used for estimating α and h_MSY for low-productivity stocks from variable and high harvest rate situations. However, TSR can severely overestimate S_MSY when spawner abundance is measured with large error in low and variable harvest rate situations. We also found that there is substantial merit in using h_MSY (or benchmarks derived from it) instead of S_MSY as a management target.     

Estimating age from recapture data: integrating incremental growth measures with ancillary data to infer age-at-length

Estimating the age of individuals in wild populations can be of fundamental importance for answering ecological questions, modeling population demographics, and managing exploited or threatened species. Significant effort has been devoted to determining age through the use of growth annuli, secondary physical characteristics related to age, and growth models. Many species, however, either do not exhibit physical characteristics useful for independent age validation or are too rare to justify sacrificing a large number of individuals to establish the relationship between size and age. Length-at-age models are well represented in the fisheries and other wildlife management literature. Many of these models overlook variation in growth rates of individuals and consider growth parameters as population parameters. More recent models have taken advantage of hierarchical structuring of parameters and Bayesian inference methods to allow for variation among individuals as functions of environmental covariates or individual-specific random effects. Here, we describe hierarchical models in which growth curves vary as individual-specific stochastic processes, and we show how these models can be fit using capture-recapture data for animals of unknown age along with data for animals of known age. We combine these independent data sources in a Bayesian analysis, distinguishing natural variation (among and within individuals) from measurement error. We illustrate using data for African dwarf crocodiles, comparing von Bertalanffy and logistic growth models. The analysis provides the means of predicting crocodile age, given a single measurement of head length. The von Bertalanffy was much better supported than the logistic growth model and predicted that dwarf crocodiles grow from 19.4 cm total length at birth to 32.9 cm in the first year and 45.3 cm by the end of their second year. Based on the minimum size of females observed with hatchlings, reproductive maturity was estimated to be at nine years. These size benchmarks are believed to represent thresholds for important demographic parameters; improved estimates of age, therefore, will increase the precision of population projection models. The modeling approach that we present can be applied to other species and offers significant advantages when multiple sources of data are available and traditional aging techniques are not practical.     

A Bayesian state-space formulation of dynamic occupancy models

Species occurrence and its dynamic components, extinction and colonization probabilities, are focal quantities in biogeography and metapopulation biology, and for species conservation assessments. It has been increasingly appreciated that these parameters must be estimated separately from detection probability to avoid the biases induced by non-detection error. Hence, there is now considerable theoretical and practical interest in dynamic occupancy models that contain explicit representations of metapopulation dynamics such as extinction, colonization, and turnover as well as growth rates. We describe a hierarchical parameterization of these models that is analogous to the state-space formulation of models in time series, where the model is represented by two components, one for the partially observable occupancy process and another for the observations conditional on that process. This parameterization naturally allows estimation of all parameters of the conventional approach to occupancy models, but in addition, yields great flexibility and extensibility, e.g., to modeling heterogeneity or latent structure in model parameters. We also highlight the important distinction between population and finite sample inference; the latter yields much more precise estimates for the particular sample at hand. Finite sample estimates can easily be obtained using the state-space representation of the model but are difficult to obtain under the conventional approach of likelihood-based estimation. We use R and WinBUGS to apply the model to two examples. In a standard analysis for the European Crossbill in a large Swiss monitoring program, we fit a model with year-specific parameters. Estimates of the dynamic parameters varied greatly among years, highlighting the irruptive population dynamics of that species. In the second example, we analyze route occupancy of Cerulean Warblers in the North American Breeding Bird Survey (BBS) using a model allowing for site-specific heterogeneity in model parameters. The results indicate relatively low turnover and a stable distribution of Cerulean Warblers which is in contrast to analyses of counts of individuals from the same survey that indicate important declines. This discrepancy illustrates the inertia in occupancy relative to actual abundance. Furthermore, the model reveals a declining patch survival probability, and increasing turnover, toward the edge of the range of the species, which is consistent with metapopulation perspectives on the genesis of range edges. Given detection/non-detection data, dynamic occupancy models as described here have considerable potential for the study of distributions and range dynamics.     

Poor-data and data-poor species stock assessment using a Bayesian hierarchical approach

Appropriate inference for stocks or species with low-quality data (poor data) or limited data (data poor) is extremely important. Hierarchical Bayesian methods are especially applicable to small-area, small-sample-size estimation problems because they allow poor-data species to borrow strength from species with good-quality data. We used a hammerhead shark complex as an example to investigate the advantages of using hierarchical Bayesian models in assessing the status of poor-data and data-poor exploited species. The hammerhead shark complex (Sphyrna spp.) along the Atlantic and Gulf of Mexico coasts of the United States is composed of three species: the scalloped hammerhead (S. lewini), the great hammerhead (S. mokarran), and the smooth hammerhead (S. zygaena) sharks. The scalloped hammerhead comprises 70-80% of the catch and has catch and relative abundance data of good quality, whereas great and smooth hammerheads have relative abundance indices that are both limited and of low quality presumably because of low stock density and limited sampling. Four hierarchical Bayesian state-space surplus production models were developed to simulate variability in population growth rates, carrying capacity, and catchability of the species. The results from the hierarchical Bayesian models were considerably more robust than those of the nonhierarchical models. The hierarchical Bayesian approach represents an intermediate strategy between traditional models that assume different population parameters for each species and those that assume all species share identical parameters. Use of the hierarchical Bayesian approach is suggested for future hammerhead shark stock assessments and for modeling fish complexes with species-specific data, because the poor-data species can borrow strength from the species with good data, making the estimation more stable and robust.     

A Bayesian hierarchical framework for calibrating aquatic biogeochemical models

Model practitioners increasingly place emphasis on rigorous quantitative error analysis in aquatic biogeochemical models and the existing initiatives range from the development of alternative metrics for goodness of fit, to data assimilation into operational models, to parameter estimation techniques. However, the treatment of error in many of these efforts is arguably selective and/or ad hoc. A Bayesian hierarchical framework enables the development of robust probabilistic analysis of error and uncertainty in model predictions by explicitly accommodating measurement error, parameter uncertainty, and model structure imperfection. This paper presents a Bayesian hierarchical formulation for simultaneously calibrating aquatic biogeochemical models at multiple systems (or sites of the same system) with differences in their trophic conditions, prior precisions of model parameters, available information, measurement error or inter-annual variability. Our statistical formulation also explicitly considers the uncertainty in model inputs (model parameters, initial conditions), the analytical/sampling error associated with the field data, and the discrepancy between model structure and the natural system dynamics (e.g., missing key ecological processes, erroneous formulations, misspecified forcing functions). The comparison between observations and posterior predictive monthly distributions indicates that the plankton models calibrated under the Bayesian hierarchical scheme provided accurate system representations for all the scenarios examined. Our results also suggest that the Bayesian hierarchical approach allows overcoming problems of insufficient local data by “borrowing strength” from well-studied sites and this feature will be highly relevant to conservation practices of regions with a high number of freshwater resources for which complete data could never be practically collected. Finally, we discuss the prospect of extending this framework to spatially explicit biogeochemical models (e.g., more effectively connect inshore with offshore areas) along with the benefits for environmental management, such as the optimization of the sampling design of monitoring programs and the alignment with the policy practice of adaptive management.     

Lifting a veil on diversity: a Bayesian approach to fitting relative-abundance models.

Bayesian methods incorporate prior knowledge into a statistical analysis. This prior knowledge is usually restricted to assumptions regarding the form of probability distributions of the parameters of interest, leaving their values to be determined mainly through the data. Here we show how a Bayesian approach can be applied to the problem of drawing inference regarding species abundance distributions and comparing diversity indices between sites. The classic log series and the lognormal models of relative- abundance distribution are apparently quite different in form. The first is a sampling distribution while the other is a model of abundance of the underlying population. Bayesian methods help unite these two models in a common framework. Markov chain Monte Carlo simulation can be used to fit both distributions as small hierarchical models with shared common assumptions. Sampling error can be assumed to follow a Poisson distribution. Species not found in a sample, but suspected to be present in the region or community of interest, can be given zero abundance. This not only simplifies the process of model fitting, but also provides a convenient way of calculating confidence intervals for diversity indices. The method is especially useful when a comparison of species diversity between sites with different sample sizes is the key motivation behind the research. We illustrate the potential of the approach using data on fruit-feeding butterflies in southern Mexico. We conclude that, once all assumptions have been made transparent, a single data set may provide support for the belief that diversity is negatively affected by anthropogenic forest disturbance. Bayesian methods help to apply theory regarding the distribution of abundance in ecological communities to applied conservation.     

Estimating wildfire growth from noisy and incomplete incident data using a state space model

Wildfire behaviors are complex and are of interest to fire managers and scientists for a variety of reasons. Many of these important behaviors are directly measured from the cumulative burn area time series of individual wildfires; however, estimating cumulative burn area time series is challenging due to the magnitude of measurement errors and missing entries. To resolve this, we introduce two state space models for reconstructing wildfire burn area using repeated observations from multiple data sources that include different levels of measurement error and temporal gaps. The constant growth parameter model uses a few parameters and assumes a burn area time series that follows a logistic growth curve. The non-constant growth parameter model uses a time-varying logistic growth curve to produce detailed estimates of the burn area time series that permit sudden pauses and pulses of growth. We apply both reconstruction models to burn area data from 13 large wildfire incidents to compare the quality of the burn area time series reconstructions and computational requirements. The constant growth parameter model reconstructs burn area time series with minimal computational requirements, but inadequately fits observed data in most cases. The non-constant growth parameter model better describes burn area time series, but can also be highly computationally demanding. Sensitivity analyses suggest that in a typical application, the reconstructed cumulative burn area time series is fairly robust to minor changes in the prior distributions.     

Impact of bias in predicted height on tree volume estimation: A case-study of intrinsic nonlinearity

Bias originating from intrinsic nonlinearity in nonlinear models is caused by excess curvature in the solution locus of parameter estimates derived from least squares procedures. Bias due to intrinsic nonlinearity varies according to sample size as well as model specification. This paper analyses consequences of fractionising data into smaller sub-samples. Based on measurements of stem diameter and total tree height from the first Danish national forest inventory, it is demonstrated how data splitting at random may cause the intrinsic nonlinear curvature to exceed the critical F-value. Application of a Taylor-series expansion shows that, for all practical purposes, the bias in predictions of individual tree volume (based on stem diameter and tree height) is negligible. To minimize residual variance, intrinsic curvature and, in turn, prediction bias, it is recommended that data be stratified according to site conditions, stand characteristics or other relevant criteria. Finally, the preferred model should exhibit close-to-linear behaviour.     

Optimization of ecosystem model parameters using spatio-temporal soil moisture information

Parameters in process-based terrestrial ecosystem models are often nonlinearly related to the water flux to the atmosphere, and they also change temporally and spatially. Therefore, for estimating soil moisture, process-based terrestrial ecosystem models inevitably need to specify spatially and temporally variant model parameters. This study presents a two-stage data assimilation scheme (TSDA) to spatially and temporally optimize some key parameters of an ecosystem model which are closely related to soil moisture. At the first stage, a simplified ecosystem model, namely the Boreal Ecosystem Productivity Simulator (BEPS), is used to obtain the prior estimation of daily soil moisture. After the spatial distribution of 0–10 cm surface soil moisture is derived from remote sensing, an Ensemble Kalman Filter is used to minimize the difference between the remote sensing model results, through optimizing some model parameters spatially. At the second stage, BEPS is reinitialized using the optimized parameters to provide the updated model predictions of daily soil moisture. TSDA has been applied to an arid and semi-arid area of northwest China, and the performance of the model for estimating daily 0–10 cm soil moisture after parameter optimization was validated using field measurements. Results indicate that the TSDA developed in this study is robust and efficient in both temporal and spatial model parameter optimization. After performing the optimization, the correlation (r²) between model-predicted 0–10 cm soil moisture and field measurement increased from 0.66 to 0.75. It is demonstrated that spatial and temporal optimization of ecosystem model parameters can not only improve the model prediction of daily soil moisture but also help to understand the spatial and temporal variation of some key parameters in an ecosystem model and the corresponding ecological mechanisms controlling the variation.     

A two-stage ensemble Kalman filter for smooth data assimilation

The ensemble Kalman Filter (EnKF) applied to a simple fire propagation model by a nonlinear convection-diffusion-reaction partial differential equation breaks down because the EnKF creates nonphysical ensemble members with large gradients. A modification of the EnKF is proposed by adding a regularization term that penalizes large gradients. The method is implemented by applying the EnKF formulas twice, with the regularization term as another observation. The regularization step is also interpreted as a shrinkage of the prior distribution. Numerical results are given to illustrate success of the new method.     

Sampling variability and estimates of density dependence: a composite-likelihood approach

It is well known that sampling variability, if not properly taken into account, affects various ecologically important analyses. Statistical inference for stochastic population dynamics models is difficult when, in addition to the process error, there is also sampling error. The standard maximum-likelihood approach suffers from large computational burden. In this paper, I discuss an application of the composite-likelihood method for estimation of the parameters of the Gompertz model in the presence of sampling variability. The main advantage of the method of composite likelihood is that it reduces the computational burden substantially with little loss of statistical efficiency. Missing observations are a common problem with many ecological time series. The method of composite likelihood can accommodate missing observations in a straightforward fashion. Environmental conditions also affect the parameters of stochastic population dynamics models. This method is shown to handle such nonstationary population dynamics processes as well. Many ecological time series are short, and statistical inferences based on such short time series tend to be less precise. However, spatial replications of short time series provide an opportunity to increase the effective sample size. Application of likelihood-based methods for spatial time-series data for population dynamics models is computationally prohibitive. The method of composite likelihood is shown to have significantly less computational burden, making it possible to analyze large spatial time-series data. After discussing the methodology in general terms, I illustrate its use by analyzing a time series of counts of American Redstart (Setophaga ruticilla) from the Breeding Bird Survey data, San Joaquin kit fox (Vulpes macrotis mutica) population abundance data, and spatial time series of Bull trout (Salvelinus confluentus) redds count data.     

Nonlinearities in mathematical ecology: Phenomena and models: Would we live in Volterra's world?

Presented is a critical survey of canonical nonlinear models in theoretical population ecology, namely single-species population, prey–predator, competition, migration within a metapopulation, and trophic chains. Various nonlinear effects, like hysteresis, structural instability, dissipative structures, dynamic chaos, etc., do exist in these models, but the problem how to detect these phenomena in real ecosystems is not yet solved. In the mathematics of nonlinear models, the central question is whether the simplest, i.e., Volterra-type, nonlinearity is sufficient to reproduce a variety of nonlinear phenomena in a given model or we need a more sophisticated formalism. Examples are considered where the Volterra models fail. Although fundamental physical principles, like, e.g., the mass conservation law, should work in ecology too, the ecological origin of the models often causes mathematical effects which are distinct from those in theoretical physics. For example, the trophic-chain model does reveal a kind of chaotic behaviour, but the “ecological strange attractor” occupies an intermediate position between Lorenz's and Feigenbaum's attractors; moreover, the phase volume of our system contracts, hence the system is dissipative (like a Lorenz's one) in spite of its matter conservation property. Nevertheless, when applied properly, physical concepts, like, e.g., the thermodynamic notion of exergy, give better insight both to the patterns of nonlinear ecosystem behaviour and to comparison of the patterns.     

Hierarchical Bayesian Spatial Models for Multispecies Conservation Planning and Monitoring

Abstract: Biologists who develop and apply habitat models are often familiar with the statistical challenges posed by their data's spatial structure but are unsure of whether the use of complex spatial models will increase the utility of model results in planning. We compared the relative performance of nonspatial and hierarchical Bayesian spatial models for three vertebrate and invertebrate taxa of conservation concern (Church's sideband snails [Monadenia churchi], red tree voles [Arborimus longicaudus], and Pacific fishers [Martes pennanti pacifica]) that provide examples of a range of distributional extents and dispersal abilities. We used presence–absence data derived from regional monitoring programs to develop models with both landscape and site‐level environmental covariates. We used Markov chain Monte Carlo algorithms and a conditional autoregressive or intrinsic conditional autoregressive model framework to fit spatial models. The fit of Bayesian spatial models was between 35 and 55% better than the fit of nonspatial analogue models. Bayesian spatial models outperformed analogous models developed with maximum entropy (Maxent) methods. Although the best spatial and nonspatial models included similar environmental variables, spatial models provided estimates of residual spatial effects that suggested how ecological processes might structure distribution patterns. Spatial models built from presence–absence data improved fit most for localized endemic species with ranges constrained by poorly known biogeographic factors and for widely distributed species suspected to be strongly affected by unmeasured environmental variables or population processes. By treating spatial effects as a variable of interest rather than a nuisance, hierarchical Bayesian spatial models, especially when they are based on a common broad‐scale spatial lattice (here the national Forest Inventory and Analysis grid of 24 km² hexagons), can increase the relevance of habitat models to multispecies conservation planning.