Maximum likelihood from spatial random effects models via the stochastic approximation expectation maximization algorithm

We introduce a class of spatial random effects models that have Markov random fields (MRF) as latent processes. Calculating the maximum likelihood estimates of unknown parameters in SREs is extremely difficult, because the normalizing factors of MRFs and additional integrations from unobserved random effects are computationally prohibitive. We propose a stochastic approximation expectation-maximization (SAEM) algorithm to maximize the likelihood functions of spatial random effects models. The SAEM algorithm integrates recent improvements in stochastic approximation algorithms; it also includes components of the Newton-Raphson algorithm and the expectation-maximization (EM) gradient algorithm. The convergence of the SAEM algorithm is guaranteed under some mild conditions. We apply the SAEM algorithm to three examples that are representative of real-world applications: a state space model, a noisy Ising model, and segmenting magnetic resonance images (MRI) of the human brain. The SAEM algorithm gives satisfactory results in finding the maximum likelihood estimate of spatial random effects models in each of these instances.     

Quantitative convergence assessment for Markov chain Monte Carlo via cusums

Yu (1995) provides a novel convergence diagnostic for Markov chain Monte Carlo (MCMC) which provides a qualitative measure of mixing for Markov chains via a cusum path plot for univariate parameters of interest. The method is based upon the output of a single replication of an MCMC sampler and is therefore widely applicable and simple to use. One criticism of the method is that it is subjective in its interpretation, since it is based upon a graphical comparison of two cusum path plots. In this paper, we develop a quantitative measure of smoothness which we can associate with any given cusum path, and show how we can use this measure to obtain a quantitative measure of mixing. In particular, we derive the large sample distribution of this smoothness measure, so that objective inference is possible. In addition, we show how this quantitative measure may also be used to provide an estimate of the burn-in length for any given sampler. We discuss the utility of this quantitative approach, and highlight a problem which may occur if the chain is able to remain in any one state for some period of time. We provide a more general implementation of the method to overcome the problem in such cases.     

Bayesian inference in hidden Markov models through the reversible jump Markov chain Monte Carlo method

Hidden Markov models form an extension of mixture models which provides a flexible class of models exhibiting dependence and a possibly large degree of variability. We show how reversible jump Markov chain Monte Carlo techniques can be used to estimate the parameters as well as the number of components of a hidden Markov model in a Bayesian framework. We employ a mixture of zero-mean normal distributions as our main example and apply this model to three sets of data from finance, meteorology and geomagnetism.     

Parallel algorithms for Markov chain Monte Carlo methods in latent spatial Gaussian models

Markov chain Monte Carlo (MCMC) implementations of Bayesian inference for latent spatial Gaussian models are very computationally intensive, and restrictions on storage and computation time are limiting their application to large problems. Here we propose various parallel MCMC algorithms for such models. The algorithms' performance is discussed with respect to a simulation study, which demonstrates the increase in speed with which the algorithms explore the posterior distribution as a function of the number of processors. We also discuss how feasible problem size is increased by use of these algorithms.     

Exact inference in contingency tables via stochastic approximation Monte Carlo

Monte Carlo methods for the exact inference have received much attention recently in complete or incomplete contingency table analysis. However, conventional Markov chain Monte Carlo, such as the Metropolis–Hastings algorithm, and importance sampling methods sometimes generate the poor performance by failing to produce valid tables. In this paper, we apply an adaptive Monte Carlo algorithm, the stochastic approximation Monte Carlo algorithm (SAMC; Liang, Liu, & Carroll, 2007), to the exact test of the goodness-of-fit of the model in complete or incomplete contingency tables containing some structural zero cells. The numerical results are in favor of our method in terms of quality of estimates.     

Analyses of infectious disease data from household outbreaks by Markov chain Monte Carlo methods

The analysis of infectious disease data presents challenges arising from the dependence in the data and the fact that only part of the transmission process is observable. These difficulties are usually overcome by making simplifying assumptions. The paper explores the use of Markov chain Monte Carlo (MCMC) methods for the analysis of infectious disease data, with the hope that they will permit analyses to be made under more realistic assumptions. Two important kinds of data sets are considered, containing temporal and non-temporal information, from outbreaks of measles and influenza. Stochastic epidemic models are used to describe the processes that generate the data. MCMC methods are then employed to perform inference in a Bayesian context for the model parameters. The MCMC methods used include standard algorithms, such as the Metropolis–Hastings algorithm and the Gibbs sampler, as well as a new method that involves likelihood approximation. It is found that standard algorithms perform well in some situations but can exhibit serious convergence difficulties in others. The inferences that we obtain are in broad agreement with estimates obtained by other methods where they are available. However, we can also provide inferences for parameters which have not been reported in previous analyses.     

Discrete time modelling of disease incidence time series by using Markov chain Monte Carlo methods

Summary.  A stochastic discrete time version of the susceptible–infected–recovered model for infectious diseases is developed. Disease is transmitted within and between communities when infected and susceptible individuals interact. Markov chain Monte Carlo methods are used to make inference about these unobserved populations and the unknown parameters of interest. The algorithm is designed specifically for modelling time series of reported measles cases although it can be adapted for other infectious diseases with permanent immunity. The application to observed measles incidence series motivates extensions to incorporate age structure as well as spatial epidemic coupling between communities.     

Stochastic approximation Monte Carlo EM for change-point analysis

In the expectation–maximization (EM) algorithm for maximum likelihood estimation from incomplete data, Markov chain Monte Carlo (MCMC) methods have been used in change-point inference for a long time when the expectation step is intractable. However, the conventional MCMC algorithms tend to get trapped in local mode in simulating from the posterior distribution of change points. To overcome this problem, in this paper we propose a stochastic approximation Monte Carlo version of EM (SAMCEM), which is a combination of adaptive Markov chain Monte Carlo and EM utilizing a maximum likelihood method. SAMCEM is compared with the stochastic approximation version of EM and reversible jump Markov chain Monte Carlo version of EM on simulated and real datasets. The numerical results indicate that SAMCEM can outperform among the three methods by producing much more accurate parameter estimates and the ability to achieve change-point positions and estimates simultaneously.     

Bayesian texture segmentation of weed and crop images using reversible jump Markov chain Monte Carlo methods

Summary. A Bayesian method for segmenting weed and crop textures is described and implemented. The work forms part of a project to identify weeds and crops in images so that selective crop spraying can be carried out. An image is subdivided into blocks and each block is modelled as a single texture. The number of different textures in the image is assumed unknown. A hierarchical Bayesian procedure is used where the texture labels have a Potts model (colour Ising Markov random field) prior and the pixels within a block are distributed according to a Gaussian Markov random field, with the parameters dependent on the type of texture. We simulate from the posterior distribution by using a reversible jump Metropolis–Hastings algorithm, where the number of different texture components is allowed to vary. The methodology is applied to a simulated image and then we carry out texture segmentation on the weed and crop images that motivated the work.     

The use of simple reparameterizations to improve the efficiency of Markov chain Monte Carlo estimation for multilevel models with applications to discrete time survival models

Summary.  We consider the application of Markov chain Monte Carlo (MCMC) estimation methods to random-effects models and in particular the family of discrete time survival models. Survival models can be used in many situations in the medical and social sciences and we illustrate their use through two examples that differ in terms of both substantive area and data structure. A multilevel discrete time survival analysis involves expanding the data set so that the model can be cast as a standard multilevel binary response model. For such models it has been shown that MCMC methods have advantages in terms of reducing estimate bias. However, the data expansion results in very large data sets for which MCMC estimation is often slow and can produce chains that exhibit poor mixing. Any way of improving the mixing will result in both speeding up the methods and more confidence in the estimates that are produced. The MCMC methodological literature is full of alternative algorithms designed to improve mixing of chains and we describe three reparameterization techniques that are easy to implement in available software. We consider two examples of multilevel survival analysis: incidence of mastitis in dairy cattle and contraceptive use dynamics in Indonesia. For each application we show where the reparameterization techniques can be used and assess their performance.     

Markov chain Monte Carlo tests for designed experiments

We consider conditional exact tests of factor effects in designed experiments for discrete response variables. Similarly to the analysis of contingency tables, a Markov chain Monte Carlo method can be used for performing exact tests, when large-sample approximations are poor and the enumeration of the conditional sample space is infeasible. For designed experiments with a single observation for each run, we formulate log-linear or logistic models and consider a connected Markov chain over an appropriate sample space. In particular, we investigate fractional factorial designs with 2^p-q$2^{p-q}$ runs, noting correspondences to the models for 2^p-q$2^{p-q}$ contingency tables.     

Reversible jump Markov chain Monte Carlo algorithms for Bayesian variable selection in logistic mixed models

In this article, to reduce computational load in performing Bayesian variable selection, we used a variant of reversible jump Markov chain Monte Carlo methods, and the Holmes and Held (HH) algorithm, to sample model index variables in logistic mixed models involving a large number of explanatory variables. Furthermore, we proposed a simple proposal distribution for model index variables, and used a simulation study and real example to compare the performance of the HH algorithm with our proposed and existing proposal distributions. The results show that the HH algorithm with our proposed proposal distribution is a computationally efficient and reliable selection method.     

Bayesian image analysis with Markov chain Monte Carlo and coloured continuum triangulation models

It is now possible to carry out Bayesian image segmentation from a continuum parametric model with an unknown number of regions. However, few suitable parametric models exist. We set out to model processes which have realizations that are naturally described by coloured planar triangulations. Triangulations are already used, to represent image structure in machine vision, and in finite element analysis, for domain decomposition. However, no normalizable parametric model, with realizations that are coloured triangulations, has been specified to date. We show how this must be done, and in particular we prove that a normalizable measure on the space of triangulations in the interior of a fixed simple polygon derives from a Poisson point process of vertices. We show how such models may be analysed by using Markov chain Monte Carlo methods and we present two case-studies, including convergence analysis.     

Motor unit number estimation using reversible jump Markov chain Monte Carlo methods

Summary.  We present an application of reversible jump Markov chain Monte Carlo sampling from the field of neurophysiology where we seek to estimate the number of motor units within a single muscle. Such an estimate is needed for monitoring the progression of neuromuscular diseases such as amyotrophic lateral sclerosis. Our data consist of action potentials that were recorded from the surface of a muscle in response to stimuli of different intensities applied to the nerve supplying the muscle. During the gradual increase in intensity of the stimulus from the threshold to supramaximal, all motor units are progressively excited. However, at any given submaximal intensity of stimulus, the number of units that are excited is variable, because of random fluctuations in axonal excitability. Furthermore, the individual motor unit action potentials exhibit variability. To account for these biological properties, Ridall and co-workers developed a model of motor unit activation that is capable of describing the response where the number of motor units, N , is fixed. The purpose of this paper is to extend that model so that the possible number of motor units, N , is a stochastic variable. We illustrate the elements of our model, show that the results are reproducible and show that our model can measure the decline in motor unit numbers during the course of amyotrophic lateral sclerosis. Our method holds promise of being useful in the study of neurogenic diseases.     

Markov chain Monte Carlo methods for high dimensional inversion in remote sensing

Summary.  We discuss the inversion of the gas profiles (ozone, NO₃, NO₂, aerosols and neutral density) in the upper atmosphere from the spectral occultation measurements. The data are produced by the 'Global ozone monitoring of occultation of stars' instrument on board the Envisat satellite that was launched in March 2002. The instrument measures the attenuation of light spectra at various horizontal paths from about 100 km down to 10–20 km. The new feature is that these data allow the inversion of the gas concentration height profiles. A short introduction is given to the present operational data management procedure with examples of the first real data inversion. Several solution options for a more comprehensive statistical inversion are presented. A direct inversion leads to a non-linear model with hundreds of parameters to be estimated. The problem is solved with an adaptive single-step Markov chain Monte Carlo algorithm. Another approach is to divide the problem into several non-linear smaller dimensional problems, to run parallel adaptive Markov chain Monte Carlo chains for them and to solve the gas profiles in repetitive linear steps. The effect of grid size is discussed, and we present how the prior regularization takes the grid size into account in a way that effectively leads to a grid-independent inversion.     

Monte Carlo integration with Markov chain

There are two conceptually distinct tasks in Markov chain Monte Carlo (MCMC): a sampler is designed for simulating a Markov chain and then an estimator is constructed on the Markov chain for computing integrals and expectations. In this article, we aim to address the second task by extending the likelihood approach of Kong et al. for Monte Carlo integration. We consider a general Markov chain scheme and use partial likelihood for estimation. Basically, the Markov chain scheme is treated as a random design and a stratified estimator is defined for the baseline measure. Further, we propose useful techniques including subsampling, regulation, and amplification for achieving overall computational efficiency. Finally, we introduce approximate variance estimators for the point estimators. The method can yield substantially improved accuracy compared with Chib's estimator and the crude Monte Carlo estimator, as illustrated with three examples.     

Regenerative Markov Chain Monte Carlo for Any Distribution

While Markov chain Monte Carlo (MCMC) methods are frequently used for difficult calculations in a wide range of scientific disciplines, they suffer from a serious limitation: their samples are not independent and identically distributed. Consequently, estimates of expectations are biased if the initial value of the chain is not drawn from the target distribution. Regenerative simulation provides an elegant solution to this problem. In this article, we propose a simple regenerative MCMC algorithm to generate variates for any distribution.     

Markov chain Monte Carlo methods for parameter estimation of the modified Weibull distribution

In this paper, the Markov chain Monte Carlo (MCMC) method is used to estimate the parameters of a modified Weibull distribution based on a complete sample. While maximum-likelihood estimation (MLE) is the most used method for parameter estimation, MCMC has recently emerged as a good alternative. When applied to parameter estimation, MCMC methods have been shown to be easy to implement computationally, the estimates always exist and are statistically consistent, and their probability intervals are convenient to construct. Details of applying MCMC to parameter estimation for the modified Weibull model are elaborated and a numerical example is presented to illustrate the methods of inference discussed in this paper. To compare MCMC with MLE, a simulation study is provided, and the differences between the estimates obtained by the two algorithms are examined.     

Parallel Markov chain Monte Carlo for Bayesian hierarchical models with big data,in two stages

Due to the escalating growth of big data sets in recent years, new Bayesian Markov chain Monte Carlo (MCMC) parallel computing methods have been developed. These methods partition large data sets by observations into subsets. However, for Bayesian nested hierarchical models, typically only a few parameters are common for the full data set, with most parameters being group specific. Thus, parallel Bayesian MCMC methods that take into account the structure of the model and split the full data set by groups rather than by observations are a more natural approach for analysis. Here, we adapt and extend a recently introduced two-stage Bayesian hierarchical modeling approach, and we partition complete data sets by groups. In stage 1, the group-specific parameters are estimated independently in parallel. The stage 1 posteriors are used as proposal distributions in stage 2, where the target distribution is the full model. Using three-level and four-level models, we show in both simulation and real data studies that results of our method agree closely with the full data analysis, with greatly increased MCMC efficiency and greatly reduced computation times. The advantages of our method versus existing parallel MCMC computing methods are also described.     

Markov Chain Monte Carlo Methods for Fitting Spatiotemporal Stochastic Models in Plant Epidemiology

Strategies for controlling plant epidemics are investigated by fitting continuous time spatiotemporal stochastic models to data consisting of maps of disease incidence observed at discrete times. Markov chain Monte Carlo methods are used for fitting two such models to data describing the spread of citrus tristeza virus (CTV) in an orchard. The approach overcomes some of the difficulties encountered when fitting stochastic models to infrequent observations of a continuous process. The results of the analysis cast doubt on the effectiveness of a strategy identified from a previous spatial analysis of the CTV data. Extensions of the approaches to more general models and other problems are also considered.