A generalized estimating equations approach for capture–recapture closed population models

The estimation of population density animal population parameters, such as capture probability, population size, or population density, is an important issue in many ecological applications. Capture–recapture data may be considered as repeated observations that are often correlated over time. If these correlations are not taken into account then parameter estimates may be biased, possibly producing misleading results. We propose a generalized estimating equations (GEE) approach to account for correlation over time instead of assuming independence as in the traditional closed population capture–recapture studies. We also account for heterogeneity among observed individuals and over-dispersion, modelling capture probabilities as a function of covariates. The GEE versions of all closed population capture–recapture models and their corresponding estimating equations are proposed. We evaluate the effect of accounting for correlation structures on capture–recapture model selection based on the quasi-likelihood information criterion (QIC). An example is used for an illustrative application and for comparison to currently used methodology. A Horvitz–Thompson-like estimator is used to obtain estimates of population size based on conditional arguments. A simulation study is conducted to evaluate the performance of the GEE approach in capture-recapture studies. The GEE approach performs well for estimating population parameters, particularly when capture probabilities are high. The simulation results also reveal that estimated population size varies on the nature of the existing correlation among capture occasions.     

Analysis of non-separable spatiotemporal data with application to enterovirus study in Taiwan

We propose an estimation procedure to incorporate non-separable spatiotemporal correlation into a generalized linear mixed model. The motivation of this paper is from a study of enterovirus infection with spatial-temporal correlation. The proposed method underlying a working estimating equation comes from a generalization of weighted least squares approaches. With an iterative two-stage estimation procedure, we may address the non-identifiability problem caused by latent random effects. Under certain regularity conditions, we show that the proposed estimate has consistency and asymptotic normality for spatiotemporal data. We also conduct a model-based simulation and apply the method to the enterovirus data in Taiwan.     

Compositional analysis of overdispersed counts using generalized estimating equations

Multivariate abundance data are commonly collected in ecology, and used to explore questions of “community composition”—how relative abundance of different taxa changes with environmental conditions. In this paper, we propose a log-linear marginal modeling approach for analyzing such compositional count data, via generalized estimating equations. This method exploits the multiplicative nature of log-linear models for counts, by reparameterizing models that describe marginal effects on mean abundance. This allows partitioning into “main effects” and compositional effects, which is appealing for interpretation. We apply the proposed approach to reanalyze compositional counts of benthic invertebrates from Delaware Bay, and data of invertebrate communities inhabiting Acacia plants in eastern Australia. In both cases we resort to a resampling approach to make inferences about regression parameters, because the number of clusters was not large compared to cluster size.     

Methods for estimating peak physiological performance and correlating performance measures

Estimates of animal performance often use the maximum of a small number of laboratory trials, a method which has several statistical disadvantages. Sample maxima always underestimate the true maximum performance, and the degree of the bias depends on sample size. Here, we suggest an alternative approach that involves estimating a specific performance quantile (e.g., the 0.90 quantile). We use the information on within-individual variation in performance to obtain a sampling distribution for the residual performance measures; we use this distribution to estimate a desired performance quantile for each individual. We illustrate our approach using simulations and with data on sprint speed in lizards. The quantile method has several advantages over the sample maximum: it reduces or eliminates bias, it uses all of the data from each individual, and its accuracy is independent of sample size. Additionally, we address the estimation of correlations between two different performance measures, such as sample maxima, quantiles, or means. In particular, because of sampling variability, we propose that the correlation of sample means does a better job estimating the correlation of population maxima than the estimator which is the correlation of sample maxima.     

Efficiency of composite sampling for estimating a lognormal distribution

In many environmental studies measuring the amount of a contaminant in a sampling unit is expensive. In such cases, composite sampling is often used to reduce data collection cost. However, composite sampling is known to be beneficial for estimating the mean of a population, but not necessarily for estimating the variance or other parameters. As some applications, for example, Monte Carlo risk assessment, require an estimate of the entire distribution, and as the lognormal model is commonly used in environmental risk assessment, in this paper we investigate efficiency of composite sampling for estimating a lognormal distribution. In particular, we examine the magnitude of savings in the number of measurements over simple random sampling, and the nature of its dependence on composite size and the parameters of the distribution utilizing simulation and asymptotic calculations.     

Robust benchmark dose determination based on profile score methods

We investigate several methods commonly used to obtain a benchmark dose and show that those based on full likelihood or profile likelihood methods might have severe shortcomings. We propose two new profile likelihood-based approaches which overcome these problems. Another contribution is the extension of the benchmark dose determination to non full likelihood models, such as quasi-likelihood, generalized estimating equations, which are widely used in settings such as developmental toxicity where clustered data are encountered. This widening of the scope of application is possible by the use of (robust) score statistics. Benchmark dose methods are applied to a data set from a developmental toxicity study.     

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.     

Characterizing source-sink dynamics with genetic parentage assignments

Source-sink dynamics have been suggested to characterize the population structure of many species, but the prevalence of source-sink systems in nature is uncertain because of inherent challenges in estimating migration rates among populations. Migration rates are often difficult to estimate directly with demographic methods, and indirect genetic methods are subject to a variety of assumptions that are difficult to meet or to apply to evolutionary timescales. Furthermore, such methods cannot be rigorously applied to high-gene-flow species. Here, we employ genetic parentage assignments in conjunction with demographic simulations to infer the level of immigration into a putative sink population. We use individual-based demographic models to estimate expected distributions of parent-offspring dyads under competing sink and closed-population models. By comparing the actual number of parent-offspring dyads (identified from multilocus genetic profiles) in a random sample of individuals taken from a population to expectations under these two contrasting demographic models, it is possible to estimate the rate of immigration and test hypotheses related to the role of immigration on population processes on an ecological timescale. The difference in the expected number of parent-offspring dyads between the two population models was greatest when immigration into the sink population was high, indicating that unlike traditional population genetic inference models, the highest degree of statistical power is achieved for the approach presented here when migration rates are high. We used the proposed genetic parentage approach to demonstrate that a threatened population of Marbled Murrelets (Braclhyrarmphus marmotus) appears to be supplemented by a low level of immigration (approximately 2-6% annually) from other populations.     

A comparison of models for predicting population persistence

We consider a range of models that may be used to predict the future persistence of populations, particularly those based on discrete-state Markov processes. While the mathematical theory of such processes is very well-developed, they may be difficult to work with when attempting to estimate parameters or expected times to extinction. Hence, we focus on diffusion and other approximations to these models, presenting new and recent developments in parameter estimation for density dependent processes, and the calculation of extinction times for processes subject to catastrophes. We illustrate these and other methods using data from simulated and real time series. We give particular attention to a procedure, due to Ross et al. [Ross, J.V., Taimre, T., Pollett, P.K. On parameter estimation in population models, Theor. Popul. Biol., in press], for estimating the parameters of the stochastic SIS logistic model, and demonstrate ways in which these parameters may be used to estimate expected extinction times. Although the stochastic SIS logistic model is strictly density dependent and allows only for birth and death events, it nonetheless may be used to predict extinction times with some accuracy even for populations that are only weakly density dependent, or that are subject to catastrophes.     

A general random walk model for the leptokurtic distribution of organism movement: Theory and application

The movement of organisms is usually leptokurtic in which some individuals move long distances while the majority remains at or near the area they are released. There has been extensive research into the origin of such leptokurtic movement, but one important aspect that has been overlooked is that the foraging behaviour of most organisms is not Brownian as assumed in most existing models. In this paper we show that such non-Brownian foraging indeed gives rise to leptokurtic distribution. We first present a general random walk model to describe the organism movement by breaking the foraging of each individual into events of active movement and inactive stationary period; its foraging behaviour is therefore fully characterized by a joint probability of how far the individual can move in each active movement and the duration it remains stationary between two consecutive movements. The spatio-temporal distribution of the organism can be described by a generalized partial differential equation, and the leptokurtic distribution is a special case when the stationary period is not exponentially distributed. Empirical observations of some organisms living in different habitats indicated that their rest time shows a power-law distribution, and we speculate that this is general for other organisms. This leads to a fractional diffusion equation with three parameters to characterize the distributions of stationary period and movement distance. A method to estimate the parameters from empirical data is given, and we apply the model to simulate the movement of two organisms living in different habitats: a stream fish (Cyprinidae: Nocomis leptocephalus) in water, and a root-feeding weevil, Sitona lepidus in the soil. Comparison of the simulations with the measured data shows close agreement. This has an important implication in ecology that the leptokurtic distribution observed at population level does not necessarily mean population heterogeneity as most existing models suggested, in which the population consists of different phenotypes; instead, a homogeneous population moving in homogeneous habitat can also lead to leptokurtic distribution.     

Effective Population Size in Red-Cockaded Woodpeckers: Population and Model Differences

Loss of genetic variability in isolated populations is an important issue for conservation biology. Most studies involve only a single population of a given species and a single method of estimating rate of loss. Here we present analyses for three different Red-cockaded Woodpecker ( Picoides borealis ) populations from different geographic regions. We compare two different models for estimating the expected rate of loss of genetic variability, and test their sensitivity to model parameters. We found that the simpler model (Reed et al. 1988) consistently estimated a greater rate of loss of genetic variability from a population than did the Emigh and Pollak (1979) model. The ratio of effective population size (which describes the expected rate of loss of genetic variability) to breeder population size varied widely among Red-cockaded Woodpecker populations due to geographic variation in demography. For this species, estimates of effective size were extremely sensitive to survival parameters, but not to the probability of breeding or reproductive success. Sensitivity was sufficient that error in estimating survival rates in the field could easily mask true population differences in effective size. Our results indicate that accurate and precise demographic data are prerequisites to determining effective population size for this species using genetic models, and that a single estimate of rate of loss of genetic variability is not valid across populations.     

Applying the Patuxent Landscape Unit Model to human dominated ecosystems: the case of agriculture

Non-spatial dynamics are core to landscape simulations. Unit models simulate system interactions aggregated within one space unit of resolution used within a spatial model. For unit models to be applicable to spatial simulations they have to be formulated in a general enough way to simulate all habitat elements within the landscape. Within the Patuxent River watershed, human dominated land uses, such as agriculture and urban land, are already 50% of the current land use, while urban land is replacing forests, agriculture and wetlands at a rapid rate. The Patuxent Landscape Model (PLM) with the Patuxent General Unit Model as core (Pat-GEM) was developed as a predictive policy tool to estimate environmental impacts of such land use changes. The Pat-GEM is based on the General Ecosystem Model (GEM) developed by [Ecol. Modelling 88 1996 263]. Previous calibrations of the Pat-GEM for anthropogenic land uses have not been satisfactory due to the scarcity of appropriate data. This paper shows Pat-GEM simulations of biomass growth and nutrient uptake for crops typical within the Patuxent watershed. The Pat-GEM was expanded to include processes and fluxes that characterize agricultural land use. The most important extension was to include crop rotation into the model. Additionally, we refined the processes for planting, harvesting and fertilization by introducing specific growth parameters. Our revised Pat-GEM was calibrated against the results from Erosion Productivity Impact Calculator (EPIC) a widely used and calibrated agricultural model. We achieved high correlation between results generated with Pat-GEM and EPIC. The correlation coefficients (r²) varied between 0.87 and 0.98, with the simulation results for winter wheat showing the lowest correlation coefficients. Intercalibration using EPIC is a powerful method for calibrating the Pat-GEM model for agricultural land use. EPIC was able (a) to provide about 30% of the input data required for running the Pat-GEM model; and (b) to provide time series output data (with a daily time step) to calibrate the output variables biomass production and nutrient uptake.     

Quantifying parameter uncertainty in a coral reef model using Metropolis-Coupled Markov Chain Monte Carlo

Coral reefs are threatened ecosystems, so it is important to have predictive models of their dynamics. Most current models of coral reefs fall into two categories. The first is simple heuristic models which provide an abstract understanding of the possible behaviour of reefs in general, but do not describe real reefs. The second is complex simulations whose parameters are obtained from a range of sources such as literature estimates. We cannot estimate the parameters of these models from a single data set, and we have little idea of the uncertainty in their predictions.We have developed a compromise between these two extremes, which is complex enough to describe real reef data, but simple enough that we can estimate parameters for a specific reef from a time series. In previous work, we fitted this model to a long-term data set from Heron Island, Australia, using maximum likelihood methods. To evaluate predictions from this model, we need estimates of the uncertainty in our parameters. Here, we obtain such estimates using Bayesian Metropolis-Coupled Markov Chain Monte Carlo. We do this for versions of the model in which corals are aggregated into a single state variable (the three-state model), and in which corals are separated into four state variables (the six-state model), in order to determine the appropriate level of aggregation. We also estimate the posterior distribution of predicted trajectories in each case.In both cases, the fitted trajectories were close to the observed data, but we had doubts about the biological plausibility of some parameter estimates. We suggest that informative prior distributions incorporating expert knowledge may resolve this problem. In the six-state model, the posterior distribution of state frequencies after 40 years contained two divergent community types, one dominated by free space and soft corals, and one dominated by acroporid, pocilloporid, and massive corals. The three-state model predicts only a single community type. We conclude that the three-state model hides too much biological heterogeneity, but we need more data if we are to obtain reliable predictions from the six-state model. It is likely that there will be similarly large, but currently unevaluated, uncertainty in the predictions of other coral reef models, many of which are much more complex and harder to fit to real data.     

Regression analysis of forest damage by marginal models for correlated ordinal responses

Studies on forest damage generally cannot be carried out by common regression models, for two main reasons: Firstly, the response variable, damage state of trees, is usually observed in ordered categories. Secondly, responses are often correlated, either serially, as in a longitudinal study, or spatially, as in the application of this paper, where neighbourhood interactions exist between damage states of spruces determined from aerial pictures. Thus so-called marginal regression models for ordinal responses, taking into account dependence among observations, are appropriate for correct inference. To this end we extend the binary models of Liang and Zeger (1986) and develop an ordinal GEEI model, based on parametrizing association by global cross-ratios. The methods are applied to data from a survey conducted in Southern Germany. Due to the survey design, responses must be assumed to be spatially correlated. The results show that the proposed ordinal marginal regression models provide appropriate tools for analysing the influence of covariates, that characterize the stand, on the damage state of spruce.     

The value of information for integrated assessment models of climate change

We estimate the value of information (VOI) for three key parameters of climate integrated assessment models (IAMs): marginal damages at low temperature anomalies, marginal damages at high temperature anomalies, and equilibrium climate sensitivity. Most empirical studies of climate damages have examined temperature anomalies up to 3 °C, while some recent theoretical studies emphasize the risks of “climate catastrophes,” which depend on climate sensitivity and on marginal damages at higher temperature anomalies. We use a new IAM to estimate the VOI for each parameter over a range of assumed levels of study precision based on prior probability distributions calibrated using results from previous studies. We measure the VOI as the maximum fixed fraction of consumption that a social planner would be willing to pay to conduct a new study before setting a carbon tax. Our central results suggest that the VOI is greatest for marginal damages at high temperature anomalies.     

Plume rise and spread in buoyant releases from elevated sources in the lower atmosphere

This study focuses on the influence of emission conditions—velocity and temperature—on the dynamics of a buoyant gas release in the atmosphere. The investigations are performed by means of wind tunnel experiments and numerical simulations. The aim is to evaluate the reliability of a Lagrangian code to simulate the dispersion of a plume produced by pollutant emissions influenced by thermal and inertial phenomena. This numerical code implements the coupling between a Lagrangian stochastic model and an integral plume rise model being able to estimate the centroid trajectory. We verified the accuracy of the plume rise model and we investigated the ability of two Lagrangian models to evaluate the plume spread by means of comparisons between experiments and numerical solutions. A quantitative study of the performances of the models through some suitable statistical indices is presented and critically discussed. This analysis shows that an additional spread has to be introduced in the Lagrangian trajectory equation in order to account the dynamical and thermal effects induced by the source conditions.     

Using hidden Markov chains and empirical Bayes change-point estimation for transect data

Consider a lattice of locations in one dimension at which data are observed. We model the data as a random hierarchical process. The hidden process is assumed to have a (prior) distribution that is derived from a two-state Markov chain. The states correspond to the mean values (high and low) of the observed data. Conditional on the states, the observations are modelled, for example, as independent Gaussian random variables with identical variances. In this model, there are four free parameters: the Gaussian variance, the high and low mean values, and the transition probability in the Markov chain. A parametric empirical Bayes approach requires estimation of these four parameters from the marginal (unconditional) distribution of the data and we use the EM-algorithm to do this. From the posterior of the hidden process, we use simulated annealing to find the maximum a posteriori (MAP) estimate. Using a Gibbs sampler, we also obtain the maximum marginal posterior probability (MMPP) estimate of the hidden process. We use these methods to determine where change-points occur in spatial transects through grassland vegetation, a problem of considerable interest to plant ecologists.     

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.     

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.     

Comparison of some homogeneity tests in analysis of over-dispersed binomial data

This paper compares the procedures based on the extended quasi-likelihood, pseudo-likelihood and quasi-likelihood approaches for testing homogeneity of several proportions for over-dispersed binomial data. The type I error of the Wald tests using the model-based and robust variance estimates, the score test, and the extended quasi-likelihood ratio test (deviance reduction test) were examined by simulation. The extended quasi-likelihood method performs less well when mean responses are close to 1 or 0. The model-based Wald test based on the quasi-likelihood performs the best in maintaining the nominal level. The score test performs less well when the intracluster correlations are large or heterogeneous. In summary: (i) both the quasilikelihood and pseudo-likelihood methods appear to be acceptable but care must be taken when overfitting a variance function with small sample sizes; (ii) the extended quasi-likelihood approach is the least favourable method because its nominal level is much too high; and (iii) the robust variance estimator performs poorly, particularly when the sample size is small.