On Estimation and Selection of Autologistic Regression Models via Penalized Pseudolikelihood

Autologistic regression models are suitable for relating spatial binary responses in ecology to covariates such as environmental factors. For big ecological data, pseudolikelihood estimation is appealing due to its ease of computation, but at least two challenges remain. Although an important issue, it is unclear how model selection may be carried out under pseudolikelihood. In addition, for assessing the variation of pseudolikelihood estimates, parametric bootstrap using Monte Carlo simulation is often used but may be infeasible for very large data sizes. Here both these issues are addressed by developing a penalized pseudolikelihood estimation method and an approximation of the variance of the parameter estimates. A simulation study is conducted to evaluate the performance of the proposed method, followed by a data example in a study of land cover in relation to land ownership characteristics. Extension of these models and methods to spatial-temporal binary data is further discussed. This article has supplementary material online.     

Autologistic regression analysis of spatial-temporal binary data via Monte Carlo maximum likelihood

This article considers logistic regression analysis of binary data that are measured on a spatial lattice and repeatedly over discrete time points. We propose a spatial-temporal autologistic regression model and draw statistical inference via maximum likelihood. Due to an unknown normalizing constant in the likelihood function, we use Monte Carlo to obtain maximum likelihood estimates of the model parameters and predictive distributions at future time points. We also use path sampling to estimate the unknown normalizing constant and approximate an information criterion for model assessment. The methodology is illustrated by the analysis of a dataset of mountain pine beetle outbreaks in western Canada.     

A linear functional relationship model for circular data with an application to the assessment of ocean wave measurements

This article suggests a linear functional relationship model for comparing two sets of circular data subject to unobservable errors. Unlike the corresponding and relatively well-studied model for linear data, maximum likelihood estimation for this model is very complicated and no explicit solutions are possible. Using a numerical approximation, we are able to solve the likelihood equations approximately, and to obtain good approximations to the likelihood estimates of the parameters. The quality of our estimates and the feasibility of the estimation method are illustrated via simulation. By establishing a parallel with the model for linear data, we are able to explain the various problems occurring in the process of estimation and to substantiate our numerical results. The interest in the model arose in connection with the study of ocean wave data; an application to such data is also given.     

Using the truncated auto-Poisson model for spatially correlated counts of vegetation

With vegetation data there are often physical reasons for believing that the response of neighbors has a direct influence on the response at a particular location. In terms of modeling such scenarios the family of auto-models or Markov random fields is a useful choice. If the observed responses are counts, the auto-Poisson model can be used. There are different ways to formulate the auto-Poisson model, depending on the biological context. A drawback of this model is that for positive autocorrelation the likelihood of the auto-Poisson model is not available in closed form. We investigate how this restriction can be avoided by right truncating the distribution. We review different parameter estimation techniques which apply to auto-models in general and compare them in a simulation study. Results suggest that the method which is most easily implemented via standard statistics software, maximum pseudo-likelihood, gives unbiased point estimates, but its variance estimates are biased. An alternative method, Monte Carlo maximum likelihood, works well but is computer-intensive and not available in standard software. We illustrate the methodology and techniques for model checking with clover leaf counts and seed count data from an agricultural experiment.     

Capture-Recapture Smooth Estimation of Age-Specific Survival Probabilities in Animal Populations

Yielding sound estimates of survival according to age in wild populations where senescence or other age-related variations may occur is very important to management decision makers, and remains challenging. This paper proposes to use penalized maximum likelihood to obtain smooth estimates of annual survival probabilities across age in populations of wild animals followed by capture–recapture. We propose to use two different types of smoothing penalties, and we use ν-fold cross-validation to select the best value of the tuning parameter for the intensity of smoothing. We then assess the accuracy of the method by a simulation study with two different shapes of the relationship between age and survival, and we conclude that a careful use of this method provides reliable noise-free estimates of age-specific annual survival. We apply this procedure to the motivating data from a population of roe deer known to exhibit a marked decrease of survival with age, and we compare our results with those previously published on this population.     

Bayesian Methods for Hierarchical Distance Sampling Models

The few distance sampling studies that use Bayesian methods typically consider only line transect sampling with a half-normal detection function. We present a Bayesian approach to analyse distance sampling data applicable to line and point transects, exact and interval distance data and any detection function possibly including covariates affecting detection probabilities. We use an integrated likelihood which combines the detection and density models. For the latter, densities are related to covariates in a log-linear mixed effect Poisson model which accommodates correlated counts. We use a Metropolis-Hastings algorithm for updating parameters and a reversible jump algorithm to include model selection for both the detection function and density models. The approach is applied to a large-scale experimental design study of northern bobwhite coveys where the interest was to assess the effect of establishing herbaceous buffers around agricultural fields in several states in the US on bird densities. Results were compared with those from an existing maximum likelihood approach that analyses the detection and density models in two stages. Both methods revealed an increase of covey densities on buffered fields. Our approach gave estimates with higher precision even though it does not condition on a known detection function for the density model.     

Bayesian estimation of survival from mark-recapture data

An understanding of survival patterns is a fundamental component of animal population biology. Mark-recapture models are often used in the estimation of animal survival rates. Maximum likelihood estimation, via either analytic solution or numerical approximation, has typically been used for inference in these models throughout the literature. In this article, a Bayesian approach is outlined and an easily applicable implementation via Markov chain Monte Carlo is described. The method is illustrated using 13 years of mark-recapture data for fulmar petrels on an island in Orkney. Point estimates of survival are similar to the maximum likelihood estimates (MLEs), but the posterior variances are smaller than the corresponding asymptotic variances of the MLEs. The Bayesian approach yields point estimates of 0.9328 for the average annual survival probability and 14.37 years for the expected lifetime of the fulmar petrels. A simple modification that accounts for missing data is also described. The approach is easier to apply than augmentation methods in this case, and simulations indicate that the performance of the estimators is not significantly diminished by the missing data.     

Modeling spatial-temporal binary data using Markov random fields

An autologistic regression model consists of a logistic regression of a response variable on explanatory variables and an autoregression on responses at neighboring locations on a lattice. It is a Markov random field with pairwise spatial dependence and is a popular tool for modeling spatial binary responses. In this article, we add a temporal component to the autologistic regression model for spatial-temporal binary data. The spatial-temporal autologistic regression model captures the relationship between a binary response and potential explanatory variables, and adjusts for both spatial dependence and temporal dependence simultaneously by a space-time Markov random field. We estimate the model parameters by maximum pseudo-likelihood and obtain optimal prediction of future responses on the lattice by a Gibbs sampler. For illustration, the method is applied to study the outbreaks of southern pine bettle in North Carolina. We also discuss the generality of our approach for modeling other types of spatial-temporal lattice data.     

A resource selection model for analyzing pseudoreplicated data due to grouping behavior of animals

Quantifying resource selection is of primary interest in animal ecology. Many analyses of resource selection assume spatial and temporal independence of the sampling unit. Autocorrelation between observations, which is a general property of ecological variables, causes difficulties for most standard statistical procedures of resource selection because autocorrelated data violate the assumption of independence. To overcome this problem, we develop a mixed-effects model to estimate resource selection functions from data that are autocorrelated because of unobserved grouping behavior by animals. In the application of the expectation-maximization (EM) algorithm, the computation of the conditional expectation of the complete-data log-likelihood function does not have a closed-form solution requiring numerical integration. A Monte Carlo EM algorithm with Gibbs sampling can be used effectively in such situations to find exact maximum likelihood estimates. We propose a simple automated Monte Carlo EM algorithm with multiple sequences in which the Monte Carlo sample size is increased when the EM step is swamped by Monte Carlo errors.We demonstrate that the model can detect inherent autocorrelation and provide reasonable variance estimates when applied to nocturnal bird migration data. This approach could also be applied to ecological processes with various types of spatially and temporally autocorrelated data, circumventing serious problems caused by dangerous pseudoreplication.     

Estimating Parameters in Delay Differential Equation Models

Delay differential equations (DDEs) are widely used in ecology, physiology and many other areas of applied science. Although the form of the DDE model is usually proposed based on scientific understanding of the dynamic system, parameters in the DDE model are often unknown. Thus it is of great interest to estimate DDE parameters from noisy data. Since the DDE model does not usually have an analytic solution, and the numeric solution requires knowing the history of the dynamic process, the traditional likelihood method cannot be directly applied. We propose a semiparametric method to estimate DDE parameters. The key feature of the semiparametric method is the use of a flexible nonparametric function to represent the dynamic process. The nonparametric function is estimated by maximizing the DDE-defined penalized likelihood function. Simulation studies show that the semiparametric method gives satisfactory estimates of DDE parameters. The semiparametric method is demonstrated by estimating a DDE model from Nicholson’s blowfly population data.     

Likelihood-based modeling and analysis of possum trapping data

Possums are a major environmental threat in New Zealand. There is no simple way to estimate possum numbers directly, and most estimates are based on an index of the proportion of leg-hold traps that catch possums. In this article, possum trapping data are used in conjunction with a plausible stochastic model and maximum likelihood estimation to construct a direct estimate of the number of possums in the vicinity of the trap lines. The model assumes that possums are caught in the traps in a Poisson process at a rate that is proportional to the product of the declining density of possums and the number of free traps, with the constant of proportionality log-linearly dependent on the accumulating number of trapped possums and the number of possums caughton previousnights. Numerical solutions of the differential equations for the probabilities associated with this stochastic process were used to construct a full likelihood of the data and hence maximum likelihood estimation of all parameters specifying the model. Based on the likelihood ratio statistic, strong serial dependence of successive nightly numbers of trapped possums was found, together with a weaker attractive effect whereby trapped possums tended to attract other possums into neighboring traps. Additionally, a maximum likelihood estimate of the local number of possums present in the vicinity of the trap lines was determined, with confidence intervals constructed from the profile log likelihood.     

A hierarchical model for compositional data analysis

This article introduces a hierarchical model for compositional analysis. Our approach models both source and mixture data simultaneously, and accounts for several different types of variation: these include measurement error on both the mixture and source data; variability in the sample from the source distributions; and variability in the mixing proportions themselves, generally of main interest. The method is an improvement on some existing methods in that estimates of mixing proportions (including their interval estimates) are sure to lie in the range [0, 1]; in addition, it is shown that our model can help in situations where identification of appropriate source data is difficult, especially when we extend our model to include a covariate. We first study the likelihood surface of a base model for a simple example, and then include prior distributions to create a Bayesian model that allows analysis of more complex situations via Markov chain Monte Carlo sampling from the likelihood. Application of the model is illustrated with two examples using real data: one concerning chemical markers in plants, and another on water chemistry.     

Dose Response Models with Natural Mortality and Random Effects

When fitting dose–response models to entomological data it is often necessary to take account of natural mortality and/or overdispersion. The standard approach to handle natural mortality is to use Abbott’s formula, which allows for a constant underlying mortality rate. Commonly used overdispersion models include the beta-binomial model, logistic-normal, and discrete mixtures. Here we extend the standard natural mortality model by including a random effect to account for overdispersion. Parameter estimation is based on a combined EM Newton–Raphson algorithm, which provides a simple framework for maximum likelihood estimation of the natural mortality model. We consider the application of this model to data from an experiment on the use of a virus (PhopGV) for the biological control of worm larvae (Phthorimaea operculella) in potatoes. For this natural mortality model with a random effect we introduce the likelihood ratio test, effective dose, and the use of a simulated residual envelope for model checking. Comparisons are made with an equivalent beta-binomial model. The procedures are implemented in the R system.     

Polar bears in the Beaufort Sea: A 30-year mark-recapture case history

Knowledge of population size and trend is necessary to manage anthropogenic risks to polar bears (Ursus maritimus). Despite capturing over 1,025 females between 1967 and 1998, previously calculated estimates of the size of the southern Beaufort Sea (SBS) population have been unreliable. We improved estimates of numbers of polar bears by modeling heterogeneity in capture probability with covariates. Important covariates referred to the year of the study, age of the bear, capture effort, and geographic location. Our choice of best approximating model was based on the inverse relationship between variance in parameter estimates and likelihood of the fit and suggested a growth from ≈ 500 to over 1,000 females during this study. The mean coefficient of variation on estimates for the last decade of the study was 0.16—the smallest yet derived. A similar model selection approach is recommended for other projects where a best model is not identified by likelihood criteria alone.     

An Efficient Approach to Spatiotemporal Analysis and Modeling of Air Pollution Data

A statistically efficient approach is adopted for modeling spatial time-series of large data sets. The estimation of the main diagnostic tool such as the likelihood function in Gaussian spatiotemporal models is a cumbersome task when using extended spatial time-series such as air pollution. Here, using the Innovation Algorithm, we manage to compute it for many spatiotemporal specifications. These specifications refer to the spatial periodic-trend, the spatial autoregressive moving average, the spatial autoregressive integrated and fractionally integrated moving average Gaussian models. Our method is applied to daily pollutants over a large metropolitan area like Athens. In the applied part of our paper, we first diagnose temporal and spatial structures of data using non-likelihood based criteria, such as the empirical autocorrelation and covariance functions. Second, we use likelihood and non-likelihood based criteria to select a spatiotemporal model among various specifications. Finally, using kriging we regionalize the resulting parameter estimates of the best-fitted model in space at any unmonitored location in the Athens region. The results show that a specific autoregressive integrated moving average spatiotemporal model can optimally perform in within and out of spatial sample estimation. Supplemental materials for this article are available from the journal website.     

Estimating variograms of soil properties by the method-of-moments and maximum likelihood

Variograms of soil properties are usually obtained by estimating the variogram for distinct lag classes by the method‐of‐moments and fitting an appropriate model to the estimates. An alternative is to fit a model by maximum likelihood to data on the assumption that they are a realization of a multivariate Gaussian process. This paper compares the two using both simulation and real data. The method‐of‐moments and maximum likelihood were used to estimate the variograms of data simulated from stationary Gaussian processes. In one example, where the simulated field was sampled at different intensities, maximum likelihood estimation was consistently more efficient than the method‐of‐moments, but this result was not general and the relative performance of the methods depends on the form of the variogram. Where the nugget variance was relatively small and the correlation range of the data was large the method‐of‐moments was at an advantage and likewise in the presence of data from a contaminating distribution. When fields were simulated with positive skew this affected the results of both the method‐of‐moments and maximum likelihood. The two methods were used to estimate variograms from actual metal concentrations in topsoil in the Swiss Jura, and the variograms were used for kriging. Both estimators were susceptible to sampling problems which resulted in over‐ or underestimation of the variance of three of the metals by kriging. For four other metals the results for kriging using the variogram obtained by maximum likelihood were consistently closer to the theoretical expectation than the results for kriging with the variogram obtained by the method‐of‐moments, although the differences between the results using the two approaches were not significantly different from each other or from expectation. Soil scientists should use both procedures in their analysis and compare the results.     

Spatial designs and properties of spatial correlation: Effects on covariance estimation

In a spatial regression context, scientists are often interested in a physical interpretation of components of the parametric covariance function. For example, spatial covariance parameter estimates in ecological settings have been interpreted to describe spatial heterogeneity or “patchiness” in a landscape that cannot be explained by measured covariates. In this article, we investigate the influence of the strength of spatial dependence on maximum likelihood (ML) and restricted maximum likelihood (REML) estimates of covariance parameters in an exponential-with-nugget model, and we also examine these influences under different sampling designs—specifically, lattice designs and more realistic random and cluster designs—at differing intensities of sampling (n=144 and 361). We find that neither ML nor REML estimates perform well when the range parameter and/or the nugget-to-sill ratio is large—ML tends to underestimate the autocorrelation function and REML produces highly variable estimates of the autocorrelation function. The best estimates of both the covariance parameters and the autocorrelation function come under the cluster sampling design and large sample sizes. As a motivating example, we consider a spatial model for stream sulfate concentration.     

Efficient profile-likelihood confidence intervals for capture-recapture models

In a capture-recapture analysis, uncertainty in the parameter estimates is usually expressed by presenting classical Wald-type confidence intervals. This approach involves (1) the assumption that the maximum likelihood estimates are asymptotically normal and (2) numerical computation of the variance-covariance matrix of these estimates. When the sample size is small or when the estimates are on the boundary of their domain, a Wald confidence interval often performs badly. A natural alternative is to use profile-likelihood confidence intervals. In general, these intervals require a greater amount of computation. We propose a new implementation of this approach that is efficient, both in reducing the amount of computation and in coping with boundary estimates. We also show how profile-likelihood confidence intervals can be adjusted for overdispersion. Simulations were used to check whether nominal coverage levels were attained, and allowed us to compare this approach with the classical Wald procedure. We illustrate this work by considering a multi-state model for a sooty shearwater (Puffinus griseus) population.     

Estimating species number under an inconvenient abundance model

Estimating the number of species in a biological community based on a multinomial sample of individual organisms is a classical problem in statistical ecology. A central issue in parametric estimation is the specification of a model of the relative abundances of species given their number. A common approach to this problem is to assume that relative abundances follow a symmetric Dirichlet distribution. This is mathematically convenient but is unconnected to work by ecologists on abundance distributions in real communities. In this article we describe ML estimation based on the sequential broken stick model that has been proposed for abundances. This model is defined mechanistically, requiring that the likelihood be approximated numerically. For this to be feasible, the likelihood must be based on a small number of summary statistics. We present simulation results that show that the observed number of species and the observed number of species represented by a single individual is a reasonable set of summary statistics on which to base estimation. We apply the method to two published data sets, one involving insect species on Mount Kenya and the other involving spider species in an Appalachian forest.     

Improving Estimates of Abundance by Aggregating Sparse Capture-Recapture Data

Inferences about abundance often are based on unadjusted counts of individuals observed, in part, because of the large amount of data required to generate reliable estimates of abundance. Where capture-recapture data are sparse, aggregating data across multiple sample elements by pooling species, locations, and sampling periods increases the information available for modeling detection probability, a necessary step for estimating abundance reliably. The process of aggregating sample elements involves balancing trade-offs related to the number of aggregated elements; although larger aggregates increase the amount of information available for estimation, they often require more complex models. We describe a heuristic approach for aggregating data for studies with multiple sample elements, use simulated data to evaluate the efficacy of aggregation, and illustrate the approach using data from a field study. Aggregating data systematically improved reliability of model selection and increased accuracy of abundance estimates while still providing estimates of abundance for each original sample unit, an important benefit necessary to maintain the design and sampling structure of a study. Within the framework of capture-recapture sampling, aggregating data improves estimates of abundance and increases the reliability of subsequent inferences made from sparse data. Additional tables and datasets may be found in the online supplements.