On the latent state estimation of nonlinear population dynamics using Bayesian and non-Bayesian state-space models

Nonlinear state-space models have been increasingly applied to study population dynamics and data assimilation in environmental sciences. State-space models can account for process error and measurement error simultaneously to correct for the bias in the estimates of system state and model parameters. However, few studies have compared the performance of different nonlinear state-space models for reconstructing the state of population dynamics from noisy time series. This study compared the performance of the extended Kalman filter (EKF), unscented Kalman filter (UKF) and Bayesian nonlinear state-space models (BNSSM) through simulations. Synthetic population time series were generated using the theta logistic model with known parameters, and normally distributed process and measurement errors were introduced using the Monte Carlo simulations. At higher levels of nonlinearity, the UKF and BNSSM had lower root mean square error (RMSE) than the EKF. The BNSSM performed reliably across all levels of nonlinearity, whereas increased levels of nonlinearity resulted in higher RMSE of the EKF. The Metropolis–Hastings algorithm within the Gibbs algorithm was used to fit the theta logistic model to synthetic time series to estimate model parameters. The estimated posterior distribution of the parameter θ indicated that the 95% credible intervals included the true values of θ (=0.5 and 1.5), but did not include 1.0 and 0.0. Future studies need to incorporate the adaptive Metropolis algorithm to estimate unknown model parameters for broad applications of Bayesian nonlinear state-space models in ecological studies.     

Hidden process models for animal population dynamics.

Hidden process models are a conceptually useful and practical way to simultaneously account for process variation in animal population dynamics and measurement errors in observations and estimates made on the population. Process variation, which can be both demographic and environmental, is modeled by linking a series of stochastic and deterministic subprocesses that characterize processes such as birth, survival, maturation, and movement. Observations of the population can be modeled as functions of true abundance with realistic probability distributions to describe observation or estimation error. Computer-intensive procedures, such as sequential Monte Carlo methods or Markov chain Monte Carlo, condition on the observed data to yield estimates of both the underlying true population abundances and the unknown population dynamics parameters. Formulation and fitting of a hidden process model are demonstrated for Sacramento River winter-run chinook salmon (Oncorhynchus tshawytsha).     

Uncertainty analyses for calibrating a soil carbon balance model to agricultural field trial data in Sweden and Kenya

How do additional data of the same and/or different type contribute to reducing model parameter and predictive uncertainties? Most modeling applications of soil organic carbon (SOC) time series in agricultural field trial datasets have been conducted without accounting for model parameter uncertainty. There have been recent advances with Monte Carlo-based uncertainty analyses in the field of hydrological modeling that are applicable, relevant and potentially valuable in modeling the dynamics of SOC. Here we employed a Monte Carlo method with threshold screening known as Generalized Likelihood Uncertainty Estimation (GLUE) to calibrate the Introductory Carbon Balance Model (ICBM) to long-term field trail data from Ultuna, Sweden and Machang’a, Kenya. Calibration results are presented in terms of parameter distributions and credibility bands on time series simulations for a number of case studies. Using these methods, we demonstrate that widely uncertain model parameters, as well as strong covariance between inert pool size and rate constant parameters, exist when root mean square simulation errors were within uncertainties in input estimations and data observations. We show that even rough estimates of the inert pool (perhaps from chemical analysis) can be quite valuable to reduce uncertainties in model parameters. In fact, such estimates were more effective at reducing parameter and predictive uncertainty than an additional 16 years time series data at Ultuna. We also demonstrate an effective method to jointly, simultaneously and in principle more robustly calibrate model parameters to multiple datasets across different climatic regions within an uncertainty framework. These methods and approaches should have benefits for use with other SOC models and datasets as well.     

Sensitivity analysis of North American bird population estimates

The Partners in Flight North American Landbird Conservation Plan provided estimates of population sizes for 448 landbird species using a multiplicative model. Input parameters in this calculation included the area of state × Bird Conservation Region polygons, area-specific mean Breeding Bird Survey counts circa 1995, and adjustment factors for the distance over which species may presumably be correctly counted, the assumed pairing of singing males with non-singing females, and variability in the propensity of birds to sing over the course of the survey day. I assessed the sensitivity of this population calculation to changes in the input parameters. I assessed both local and global sensitivity of the model to changes in the parameters with Monte Carlo one-at-a-time simulations and the Fourier amplitude sensitivity test (FAST). Monte Carlo simulations were an estimate of local model sensitivity whereas FAST estimated global model sensitivity, accommodating the potential shared variance between model parameters. Monte Carlo simulations suggested population estimates were 39% more sensitive to changes in the detection distance adjustment than to the other parameters; the other parameters were nearly equal in their contribution to model sensitivity. Conversely, FAST analysis determined that each of the input variables aside from the pair adjustment provided roughly equal contributions to variability in population estimates. The most efficient means for improving continental population estimates for birds surveyed by the Breeding Bird Survey will be through increased scrutiny of the species-specific distance detection and time-of-day adjustments and improved understanding in the spatial and temporal variability in the mean Breeding Bird Survey count.     

Determining patterns of variability in ecological communities: time lag analysis revisited

All ecological communities experience change over time. One method to quantify temporal variation in the patterns of relative abundance of communities is time lag analysis (TLA). It uses a distance-based approach to study temporal community dynamics by regressing community dissimilarity over increasing time lags (one-unit lags, two-unit lags, three-unit lags). Here, we suggest some modifications to the method and revaluate its potential for detecting patterns of community change. We apply Hellinger distance based TLA to artificial data simulating communities with different levels of directional and stochastic dynamics and analyse their effects on the slope and its statistical significance. We conclude that statistical significance of the TLA slope (obtained by a Monte Carlo permutation procedure) is a valid criterion to discriminate between (i) communities with directional change in species composition, regardless whether it is caused by directional abundance change of the species or by stochastic change according to a Markov process, and (ii) communities that are composed of species with population sizes oscillating around a constant mean or communities whose species abundances are governed by a white noise process. TLA slopes range between 0.02 and 0.25, depending on the proportions of species with different dynamics; higher proportions of species with constant means imply shallower slopes; and higher proportions of species with stochastic dynamics or directional change imply steeper slopes. These values are broadly in line with TLA slopes from real world data. Caution must be exercised when TLA is used for the comparison of community time series with different lengths since the slope depends on time series length and tends to decrease non-linearly with it.     

A comparison of sensitivity analysis and error analysis based on a stream ecosystem model

When individual model parameters must be measured in field or laboratory experiments, the provision of feedback information for allocation of research efforts is an important function of modeling. Both sensitivity analysis and Monte Carlo error analysis can be used to determine which parameters require intensified measurement effort. When both methods are applied to a stream ecosystem model, the assumptions of sensitivity analysis are violated if reasonable estimates of measurement errors on parameters are used. Sensitivity analysis estimates a linear relationship between a state variable and a parameter and largely ignores higher order effects.In the model investigated in this study, higher-order effects dominate prediction error, and the results of sensitivity analysis are misleading. It is suggested that the simple correlation coefficient derived from analysis of Monte Carlo simulations is a more reasonable way to rank model parameters according to their contribution to prediction uncertainty. For the stream model used in this study, halving variance on the four parameters, indicated as most important by sensitivity analysis, reduces prediction errors by only 2–6%. Halving variance on the four, completely different, parameters with the largest simple correlation coefficients reduces prediction errors by 17–31%.     

The Importance of Being Spatial (and Reserved): Assessing Northern Spotted Owl Habitat Relationships with Hierarchical Bayesian Models

Abstract:  Regional conservation planning increasingly draws on habitat suitability models to support decisions regarding land allocation and management. Nevertheless, statistical techniques commonly used for developing such models may give misleading results because they fail to account for 3 factors common in data sets of species distribution: spatial autocorrelation, the large number of sites where the species is absent (zero inflation), and uneven survey effort. We used spatial autoregressive models fit with Bayesian Markov Chain Monte Carlo techniques to assess the relationship between older coniferous forest and the abundance of Northern Spotted Owl nest and activity sites throughout the species' range. The spatial random-effect term incorporated in the autoregressive models successfully accounted for zero inflation and reduced the effect of survey bias on estimates of species–habitat associations. Our results support the hypothesis that the relationship between owl distribution and older forest varies with latitude. A quadratic relationship between owl abundance and older forest was evident in the southern portion of the range, and a pseudothreshold relationship was evident in the northern portion of the range. Our results suggest that proposed changes to the network of owl habitat reserves would reduce the proportion of the population protected by up to one-third, and that proposed guidelines for forest management within reserves underestimate the proportion of older forest associated with maximum owl abundance and inappropriately generalize threshold relationships among subregions. Bayesian spatial models can greatly enhance the utility of habitat analysis for conservation planning because they add the statistical flexibility necessary for analyzing regional survey data while retaining the interpretability of simpler models.     

Assessing plausible rates of population growth in humpback whales from life-history data

The rate of growth of any population is a quantity of interest in conservation and management and is constrained by biological factors. In this study, recent data on life-history parameters influencing rates of population growth in humpback whales, including survival, age at first parturition and calving rate are reviewed. Monte Carlo simulations are used to compute a distribution of rates of increase (ROIs) taking into account uncertainty in biological parameter estimates. Two approaches for computing juvenile survival are proposed, which taken into account along with other life-history data, resulted in the following estimates of the rate of population growth: Approach A: mean of 7.3%/year (95% CI = 3.5–10.5%/year) and Approach B: mean of 8.6%/year (95% CI = 5.0–11.4%/year). It is proposed that the upper 99% quantile of the resulting distribution of the ROI for Approach B (11.8%/year) be established as the maximum plausible ROI for humpback whales and be used in population assessment of the species. Possible sources of positive and negative biases in the present estimates are presented and include measurement error in estimation of life-history parameters, changes in the environment within the period these quantities are measured, density dependence or other natural factors. However, it is difficult to evaluate potential biases without additional data. The methods presented in this study can be applied to other species for which life-history parameters are available and are useful in assessing plausibility in the estimation of population growth rates from time series of abundance estimates.     

Spatial statistics for modeling of abundance and distribution of wildlife species in the Masai Mara ecosystem, Kenya

This study illustrates the use of modern statistical procedures for better wildlife management by addressing three key issues: determination of abundance, modeling of animal distributions and variability of diversity in space and time. Prior information in Markov Chain Monte Carlo (MCMC) methods is used to improve estimates of abundance. Measures of autocorrelation are included when modeling distributions of animal counts, and a diversity index to indicate species abundance and richness for large herbivores is developed. Data from the Masai Mara ecosystem in Kenya are used to develop and demonstrate these procedures. The new abundance estimates are up to 35% more accurate than those obtained by existing methods. Significant temporal changes in spatial patterns are found from a space-time analysis of elephant counts over a 20-year period, with strong interactions over 5 km and 6 months space and time separations, respectively. The new diversity index is sensitive to both high abundance and species richness and is also able to capture year to year variation. It indicates an overall marginal decrease in diversity for large herbivores in the Mara ecosystem. The space-time analyses and diversity index can easily be computed thereby providing tools for rapid decision making.     

Effects of site-selection bias on estimates of biodiversity change

Estimates of biodiversity change are essential for the management and conservation of ecosystems. Accurate estimates rely on selecting representative sites, but monitoring often focuses on sites of special interest. How such site-selection biases influence estimates of biodiversity change is largely unknown. Site-selection bias potentially occurs across four major sources of biodiversity data, decreasing in likelihood from citizen science, museums, national park monitoring, and academic research. We defined site-selection bias as a preference for sites that are either densely populated (i.e., abundance bias) or species rich (i.e., richness bias). We simulated biodiversity change in a virtual landscape and tracked the observed biodiversity at a sampled site. The site was selected either randomly or with a site-selection bias. We used a simple spatially resolved, individual-based model to predict the movement or dispersal of individuals in and out of the chosen sampling site. Site-selection bias exaggerated estimates of biodiversity loss in sites selected with a bias by on average 300–400% compared with randomly selected sites. Based on our simulations, site-selection bias resulted in positive trends being estimated as negative trends: richness increase was estimated as 0.1 in randomly selected sites, whereas sites selected with a bias showed a richness change of −0.1 to −0.2 on average. Thus, site-selection bias may falsely indicate decreases in biodiversity. We varied sampling design and characteristics of the species and found that site-selection biases were strongest in short time series, for small grains, organisms with low dispersal ability, large regional species pools, and strong spatial aggregation. Based on these findings, to minimize site-selection bias, we recommend use of systematic site-selection schemes; maximizing sampling area; calculating biodiversity measures cumulatively across plots; and use of biodiversity measures that are less sensitive to rare species, such as the effective number of species. Awareness of the potential impact of site-selection bias is needed for biodiversity monitoring, the design of new studies on biodiversity change, and the interpretation of existing data.     

Statistical inference on patch-specific survival and movement rates from marked animals

Empirical estimates of patch-specific survival and movement rates are needed to parametrize spatially explicit population models, and for inference on the effects of habitat quality and fragmentation on populations. Data from radio-marked animals, in which both the fates and habitat locations of animals are known over time, can be used in conjunction with continuous-time proportional hazards models to obtain inferences on survival rates. Discrete-time conditional logistic models may provide inference on both survival and movement rates. We use Monte Carlo simulation to investigate accuracy of estimates of survival from both approaches, and movement rates from conditional logistic regression, for two habitats. Bias was low (relative bias < 0.04) and interval coverage accurate (close to the nominal 0.95) for estimates of habitat effect on survival based on proportional hazards. Bias was high ( relative bias 0.60) and interval coverage poor ( = 0.26 vs. nominal 0.95) for estimates of habitat effect based on conditional logistic regression; bias was especially influenced by heterogeneity in survival and the shape of the hazard function, whereas both bias and coverage were affected by ‘memory’ effects in movement patterns. Bias estimates of movement rate was low ( relative bias < 0.05), but interval coverage was poor ( = 0.48–0.80), possibly as a result of poor performance of a Taylor series estimate of variance. An example is provided from a radio-telemetry study of 47 wintering American woodcock (Scolopax minor), illustrating practical difficulties in field studies to parametrize these models. We also discuss extensions of continuous-time models to explicitly include a movement process, and further examine tradeoffs between continuous and discrete models.     

Comparing model averaging with other model selection strategies for benchmark dose estimation

Model averaging (MA) has been proposed as a method of accommodating model uncertainty when estimating risk. Although the use of MA is inherently appealing, little is known about its performance using general modeling conditions. We investigate the use of MA for estimating excess risk using a Monte Carlo simulation. Dichotomous response data are simulated under various assumed underlying dose–response curves, and nine dose–response models (from the USEPA Benchmark dose model suite) are fit to obtain both model specific and MA risk estimates. The benchmark dose estimates (BMDs) from the MA method, as well as estimates from other commonly selected models, e.g., best fitting model or the model resulting in the smallest BMD, are compared to the true benchmark dose value to better understand both bias and coverage behavior in the estimation procedure. The MA method has a small bias when estimating the BMD that is similar to the bias of BMD estimates derived from the assumed model. Further, when a broader range of models are included in the family of models considered in the MA process, the lower bound estimate provided coverage close to the nominal level, which is superior to the other strategies considered. This approach provides an alternative method for risk managers to estimate risk while incorporating model uncertainty.

Model correlation in stochastic forest simulators—A case of multilevel multivariate model for seedling establishment

Forest development can be predicted by the use of forest simulators based on various statistical models describing the forest and its dynamics. One potential approach to study the reliability of the simulators is to utilise Monte Carlo simulation techniques to generate a predictive distribution of a forest characteristic. One problem in examining the effect of model uncertainty in forestry decision making, however, is correlation between the models. If this is not taken into account, predictions of the model systems may become biased, and the effect of errors on decision making may be underestimated. In reality, the models often are interdependent, but the correlations usually are not known because the models have been estimated in separate studies. The aim of this paper is to study the impacts of between-model dependencies on the predictive distribution of forest characteristics by Monte Carlo simulation techniques. We utilise a case of predicting seedling establishment of planted Norway spruce (Picea abies (L.) Karst.) stands as an example with multivariate multilevel model structures. Regardless of low cross-correlations between the models, ignoring them led to significant underestimation of the amount of competing broadleaves to be removed in pre-commercial thinning. Therefore, we recommend that between-model dependencies are clarified and considered in stochastic simulations. In our case, between-model interdependencies can be reliably estimated with a limited dataset. In addition, estimating the models separately and using the model residuals to estimate interdependencies between models were also sufficient to take the between-model dependencies into account when producing stochastic predictions for silvicultural decision making.     

Rarefaction and nonrandom spatial dispersion patterns

Rarefaction estimates how many species are expected in a random sample of individuals from a larger collection and allows meaningful comparisons among collections of different sizes. It assumes random spatial dispersion. However, two common dispersion patterns, within-species clumping and segregation among species, can cause rarefaction to overestimate the species richness of a smaller continuous area. We use field studies and computer simulations to determine (1) how robust rarefaction is to nonrandom spatial dispersion and (2) whether simple measures of spatial autocorrelation can predict the bias in rarefaction estimates. Rarefaction does not estimate species richness accurately for many communities, especially at small sample sizes. Measures of spatial autocorrelation of the more abundant species do not reliably predict amount of bias. Survey sites should be standardized to equal-sized areas before sampling. When sites are of equal area but differ in number of individuals sampled, rarefaction can standardize collections. When communities are sampled from different-sized areas, the mean and confidence intervals of species accumulation curves allow more meaningful comparisons among sites. Electronic supplementary material  The online version of this article (doi:) contains supplementary material, which is available to authorized users.

Coregionalization analysis with a drift for multi-scale assessment of spatial relationships between ecological variables 1. Estimation of drift and random components

In two articles, we present ‘coregionalization analysis with a drift’ (CRAD), a method to assess the multi-scale variability of and relationships between ecological variables from a multivariate spatial data set. In phase I of CRAD (the first article), a deterministic drift component representing the large-scale pattern and a random component modeled as a second-order stationary process are estimated for each variable separately. In phase II (this article), a linear model of coregionalization (LMC) is fitted by estimated generalized least squares to the direct and cross experimental variograms of residuals (i.e., after the removal of estimated drifts). Structural correlations and coefficients of determination at smaller scales are then computed from the estimated coregionalization matrices, while the estimated drifts are used to calculate pseudo coefficients at large scale. The performance of five procedures in estimating correlations and coefficients of determination was compared using a Monte Carlo study. In four CRAD procedures, drift estimation was based on local polynomials of order 0, 1, 2 (L₀, L₁, L₂) or a global polynomial with forward selection of the basis functions; the fifth procedure was coregionalization analysis (CRA), in which large-scale patterns were modeled as a supplemental component in the LMC. In bivariate and multivariate analyses, the uncertainty in the estimation of correlations and coefficients of determination could be related to the interference between spatial components within a bounded sampling domain. In the bivariate case, most procedures provided acceptable estimates of correlations. In regionalized redundancy analysis, uncertainty was highest for CRA, while L₁ provided the best results overall. In a forest ecology example, the identification of scale-specific correlations between plant species diversity and soil and topographical variables illustrated the potential of CRAD to provide unique insight into the functioning of complex ecosystems.     

Correlations between forest fires in British Columbia, Canada, and sea surface temperature of the Pacific Ocean

Correlations and cross-correlations between forest fires in the province of British Columbia, Canada, and sea surface temperatures in the Pacific Ocean were evaluated. British Columbia has a long Pacific Ocean coastline; given that there may be teleconnections between the province's forest fires and climate variability over the ocean, significant correlations may exist between forest fires and the sea surface temperature of the Pacific Ocean. Fire occurrences and areas burned through lightning-caused and human-caused fires were analyzed against individual 1° × 1° grid cells of anomalies in the sea surface temperature to determine correlations for the period 1950-2006. Significant correlations (p < 0.05) for vast areas of the ocean were found between occurrences of lightning-caused fires and sea surface temperature anomalies for time lags of 1 and 2 years, whereas significant correlations between occurrences of human-caused fires and sea surface temperature anomalies occurred extensively for many time lags. To support the results of this approach, correlations between fire data and the Niño 3.4, Pacific Decadal Oscillation, and Arctic Oscillation indices were tested for the same period. Significant correlations were found between fire occurrences and these indices at certain time lags. Overall, fire occurrence appeared to be more extensively correlated with sea surface temperature anomalies than was area burned. These results support the hypothesis that teleconnections exist between fire activity in British Columbia and sea surface temperatures in the Pacific Ocean, and the correlations suggest that linear regression models or other regression techniques may be appropriate for predicting fire severity from the sea surface temperatures of one or more previous years.     

Analysis of binary data via a centered spatial-temporal autologistic regression model

A centered spatial-temporal autologistic model is developed for analyzing spatial-temporal binary data observed on a lattice over time. We propose expectation-maximization pseudolikelihood and Monte Carlo expectation-maximization likelihood as well as consider Bayesian inference to obtain the estimates of model parameters. Further, we compare the statistical efficiency of the three approaches for various sizes of sampling lattices and numbers of sampling time points. Regarding prediction, we use Monte Carlo to obtain predictive distributions at future time points and compare the performance of the model with the uncentered spatial-temporal autologistic regression model. The methodology is demonstrated via simulation studies and a real data example concerning southern pine beetle outbreak in North Carolina.     

Kernel density estimators of home range: smoothing and the autocorrelation red herring

Two oft-cited drawbacks of kernel density estimators (KDEs) of home range are their sensitivity to the choice of smoothing parameter(s) and their need for independent data. Several simulation studies have been conducted to compare the performance of objective, data-based methods of choosing optimal smoothing parameters in the context of home range and utilization distribution (UD) estimation. Lost in this discussion of choice of smoothing parameters is the general role of smoothing in data analysis, namely, that smoothing serves to increase precision at the cost of increased bias. A primary goal of this paper is to illustrate this bias-variance trade-off by applying KDEs to sampled locations from simulated movement paths. These simulations will also be used to explore the role of autocorrelation in estimating UDs. Autocorrelation can be reduced (1) by increasing study duration (for a fixed sample size) or (2) by decreasing the sampling rate. While the first option will often be reasonable, for a fixed study duration higher sampling rates should always result in improved estimates of space use. Further, KDEs with typical data-based methods of choosing smoothing parameters should provide competitive estimates of space use for fixed study periods unless autocorrelation substantially alters the optimal level of smoothing.     

Early warnings of unknown nonlinear shifts: a nonparametric approach

Early warning signals (EWS) of regime shifts are challenging in cases where the true natural data-generating process is uncertain. Nonparametric drift-diffusion-jump models address this problem by fitting a general model that can approximate a wide range of data-generating processes. Drift measures the local rate of change. Diffusion measures relatively small shocks that occur at each time step. Jumps are large intermittent shocks. Total variance combines the contributions of diffusion and jumps. Nonparametric methods are well suited to emerging technology for automated, high-frequency sensors. Total variance is the most precisely measured indicator. Jump intensity appears to be a useful EWS. Estimates of the drift are highly uncertain unless long time series with many regime shifts are available. EWS computed from drift estimates (such as autocorrelation coefficients or return rates) have low precision and should be used with caution. Nonetheless, in the current state of knowledge, it is premature to disregard any potential EWS.     

The interpretation of capitalized hedonic prices in a dynamic environment

This paper considers the biases in hedonic price estimates of environmental variables. It is shown that using the current values of such variables can seriously bias the measure of the true effect of that variable when the future is changing. Although the bias can be in either direction most plausible examples considered show that we tend to underestimate such prices. A discussion of the relation between rental and capital value of hedonic prices follows from the distinction between present and future levels of environmental effects.