Markov Chain Monte Carlo Software Research Articles

Monte Carlo Simulation on the Stiefel Manifold via Polar Expansion

Motivated by applications to Bayesian inference for statistical models with orthogonal matrix parameters, we present a general approach to Monte Carlo simulation from probability distributions on the Stiefel manifold. To bypass many of the well-established challenges of simulating from the distribution of a random orthogonal matrix we construct a distribution for an unconstrained random matrix X such that the orthogonal component of the polar decomposition of is equal in distribution to The distribution of X is amenable to Markov chain Monte Carlo (MCMC) simulation using standard methods, and an approximation to the distribution of Q can be recovered from a Markov chain on the unconstrained space. When combined with modern MCMC software, polar expansion allows for routine and flexible posterior inference in models with orthogonal matrix parameters. We find that polar expansion with adaptive Hamiltonian Monte Carlo is an order of magnitude more efficient than competing MCMC approaches in a benchmark protein interaction network application. We also propose a new approach to Bayesian functional principal component analysis which we illustrate in a meteorological time series application. Supplementary materials for this article are available online.

Post-processing of Markov chain Monte Carlo output in Bayesian latent variable models with application to multidimensional scaling

In Bayesian analysis of multidimensional scaling model with MCMC algorithm, we encounter the indeterminacy of rotation, reflection and translation of the parameter matrix of interest. This type of indeterminacy may be seen in other multivariate latent variable models as well. In this paper, we propose to address this indeterminacy problem with a novel, offline post-processing method that is easily implemented using easy-to-use Markov chain Monte Carlo (MCMC) software. Specifically, we propose a post-processing method based on the generalized extended Procrustes analysis to address this problem. The proposed method is compared with four existing methods to deal with indeterminacy thorough analyses of artificial as well as real datasets. The proposed method achieved at least as good a performance as the best existing method. The benefit of the offline processing approach in the era of easy-to-use MCMC software is discussed.

Independent sampling for Bayesian normal conditional autoregressive models with OpenCL acceleration

Abstract A new computational strategy produces independent samples from the joint posterior distribution for a broad class of Bayesian spatial and spatiotemporal conditional autoregressive models. The method is based on reparameterization and marginalization of the posterior distribution and massive parallelization of rejection sampling using graphical processing units (GPUs) or other accelerators. It enables very fast sampling for small to moderate-sized datasets (up to approximately 10,000 observations) and feasible sampling for much larger datasets. Even using a mid-range GPU and a high-end CPU, the GPU-based implementation is up to 30 times faster than the same algorithm run serially on a single CPU, and the numbers of effective samples per second are orders of magnitude higher than those obtained with popular Markov chain Monte Carlo software. The method has been implemented in the R package CARrampsOcl. This work provides both a practical computing strategy for fitting a popular class of Bayesian models and a proof of concept that GPU acceleration can make independent sampling from Bayesian joint posterior densities feasible.

A Linear Regression Framework for the Verification of Bayesian Model Calibration Algorithms

We describe a framework for the verification of Bayesian model calibration routines. The framework is based on linear regression and can be configured to verify calibration to data with a range of observation error characteristics. The framework is designed for efficient implementation and is suitable for verifying code intended for large-scale problems. We propose an approach for using the framework to verify Markov chain Monte Carlo (MCMC) software by combining it with a nonparametric test for distribution equality based on the energy statistic. Our matlab-based reference implementation of the framework is shown to correctly distinguish between output obtained from correctly and incorrectly implemented MCMC routines. Since correctness of output from an MCMC software depends on choosing settings appropriate for the problem-of-interest, our framework can potentially be used for verifying such settings.

Finding binaries from phase modulation of pulsating stars with Kepler – IV. Detection limits and radial velocity verification

We explore the detection limits of the phase modulation (PM) method of finding binary systems among multi-periodic pulsating stars. The method is an attractive way of finding non-transiting planets in the habitable zones of intermediate mass stars, whose rapid rotation inhibits detections via the radial velocity (RV) method. While oscillation amplitudes of a few mmag are required to find planets, many $\delta$ Scuti stars have these amplitudes. In sub-optimal cases where the signal-to-noise of the oscillations is lower, low-mass brown dwarfs ($\sim$13 M$_{\rm Jup}$) are detectable at orbital periods longer than about 1 yr, and the lowest mass main-sequence stars (0.1-0.2 M$_{\odot}$) are detectable at all orbital periods where the PM method can be applied. We use purpose-written Markov chain Monte Carlo (MCMC) software for the calculation of the PM orbits, which offers robust uncertainties for comparison with RV solutions. Using Kepler data and ground-based RVs, we verify that these two methods are in agreement, even at short orbital periods where the PM method undersamples the orbit. We develop new theory to account for the undersampling of the time delays, which is also necessary for the inclusion of RVs as observational data in the MCMC software. We show that combining RVs with time delays substantially refines the orbits because of the complementarity of working in both the spatial (PM) and velocity (RV) domains simultaneously. Software outputs were tested through an extensive hare and hounds exercise, covering a wide range of orbital configurations including binaries containing two pulsators.

Phylogenetics, likelihood, evolution and complexity

Phylogenetics, likelihood, evolution and complexity (PLEX) is a flexible and fast Bayesian Markov chain Monte Carlo software program for large-scale analysis of nucleotide and amino acid data using complex evolutionary models in a phylogenetic framework. The program gains large speed improvements over standard approaches by implementing 'partial sampling of substitution histories', a data augmentation approach that can reduce data analysis times from months to minutes on large comparative datasets. A variety of nucleotide and amino acid substitution models are currently implemented, including non-reversible and site-heterogeneous mixture models. Due to efficient algorithms that scale well with data size and model complexity, PLEX can be used to make inferences from hundreds to thousands of taxa in only minutes on a desktop computer. It also performs probabilistic ancestral sequence reconstruction. Future versions will support detection of co-evolutionary interactions between sites, probabilistic tests of convergent evolution and rigorous testing of evolutionary hypotheses in a Bayesian framework. PLEX v1.0 is licensed under GPL. Source code and documentation will be available for download at www.evolutionarygenomics.com/ProgramsData/PLEX. PLEX is implemented in C++ and supported on Linux, Mac OS X and other platforms supporting standard C++ compilers. Example data, control files, documentation and accessory Perl scripts are available from the website. David.Pollock@UCDenver.edu. Supplementary data are available at Bioinformatics online.

Automating network meta-analysis.

Mixed treatment comparison (MTC) (also called network meta-analysis) is an extension of traditional meta-analysis to allow the simultaneous pooling of data from clinical trials comparing more than two treatment options. Typically, MTCs are performed using general-purpose Markov chain Monte Carlo software such as WinBUGS, requiring a model and data to be specified using a specific syntax. It would be preferable if, for the most common cases, both could be derived from a well-structured data file that can be easily checked for errors. Automation is particularly valuable for simulation studies in which the large number of MTCs that have to be estimated may preclude manual model specification and analysis. Moreover, automated model generation raises issues that provide additional insight into the nature of MTC. We present a method for the automated generation of Bayesian homogeneous variance random effects consistency models, including the choice of basic parameters and trial baselines, priors, and starting values for the Markov chain(s). We validate our method against the results of five published MTCs. The method is implemented in freely available open source software. This means that performing an MTC no longer requires manually writing a statistical model. This reduces time and effort, and facilitates error checking of the dataset. Copyright © 2012 John Wiley & Sons, Ltd.

A statistical model for under- or overdispersed clustered and longitudinal count data

We propose a likelihood-based model for correlated count data that display under- or overdispersion within units (e.g. subjects). The model is capable of handling correlation due to clustering and/or serial correlation, in the presence of unbalanced, missing or unequally spaced data. A family of distributions based on birth-event processes is used to model within-subject underdispersion. A computational approach is given to overcome a parameterization difficulty with this family, and this allows use of common Markov Chain Monte Carlo software (e.g. WinBUGS) for estimation. Application of the model to daily counts of asthma inhaler use by children shows substantial within-subject underdispersion, between-subject heterogeneity and correlation due to both clustering of measurements within subjects and serial correlation of longitudinal measurements. The model provides a major improvement over Poisson longitudinal models, and diagnostics show that the model fits well.

GSEVM v.2: MCMC software to analyze genetically structured environmental variance models

This note provides a description of software that allows to fit Bayesian genetically structured variance models using Markov chain Monte Carlo (MCMC). The gsevm v.2 program was written in Fortran 90. The DOS and Unix executable programs, the user's guide, and some example files are freely available for research purposes at http://www.bdporc.irta.es/estudis.jsp. The main feature of the program is to compute Monte Carlo estimates of marginal posterior distributions of parameters of interest. The program is quite flexible, allowing the user to fit a variety of linear models at the level of the mean and the logvariance.

Towards quantitative uncertainty assessment for cancer risks: Central estimates and probability distributions of risk in dose–response modeling

Regulatory agencies and the scientific community have been engaged in a long-term effort to strengthen health risk assessment procedures. Recently the momentum of this effort has accelerated to increasing biological information for a variety of toxic compounds and emphasis on the policy goal of broader characterization of scientific uncertainty (in contrast to providing only a single risk estimate). For example, the OMB Regulatory Analysis Guidelines [OMB, 2003. Office of Management and Budget. Circular A-4. Available from: <http://www.whitehouse.gov/omb/circulars/a004/a-4.html/>] suggest that a formal quantitative uncertainty analysis be performed for economic assessments in support of major regulatory analyses, a process that can utilize both expected values and probability distributions for risk estimates. Some efforts have been made in the past to provide probability distributions of risk estimates. In this article, we examine a procedure for constructing probability distributions and expected values of risk estimates using a Bayesian framework. This approach has the advantage of mathematical soundness and computational feasibility, given the Markov chain Monte Carlo software tools that are available today. Importantly, the Bayesian framework can serve as a unifying platform for uncertainty analysis in cancer risk assessment. This paper provides some initial applications of Bayesian methods in quantitative analysis of uncertainty in cancer risk assessment, including implementation with cancer dose–response data sets for two chemicals. The Bayesian expected risk calculations provide an approach to generating a central estimate of risk that does not have the instability problems that have often limited utility of MLE risk estimates.

Data cloning: easy maximum likelihood estimation for complex ecological models using Bayesian Markov chain Monte Carlo methods

We introduce a new statistical computing method, called data cloning, to calculate maximum likelihood estimates and their standard errors for complex ecological models. Although the method uses the Bayesian framework and exploits the computational simplicity of the Markov chain Monte Carlo (MCMC) algorithms, it provides valid frequentist inferences such as the maximum likelihood estimates and their standard errors. The inferences are completely invariant to the choice of the prior distributions and therefore avoid the inherent subjectivity of the Bayesian approach. The data cloning method is easily implemented using standard MCMC software. Data cloning is particularly useful for analysing ecological situations in which hierarchical statistical models, such as state-space models and mixed effects models, are appropriate. We illustrate the method by fitting two nonlinear population dynamics models to data in the presence of process and observation noise.

Assessing the convergence of Markov Chain Monte Carlo methods: An example from evaluation of diagnostic tests in absence of a gold standard

The accessibility of Markov Chain Monte Carlo (MCMC) methods for statistical inference have improved with the advent of general purpose software. This enables researchers with limited statistical skills to perform Bayesian analysis. Using MCMC sampling to do statistical inference requires convergence of the MCMC chain to its stationary distribution. There is no certain way to prove convergence; it is only possible to ascertain when convergence definitely has not been achieved. These methods are rather subjective and not implemented as automatic safeguards in general MCMC software. This paper considers a pragmatic approach towards assessing the convergence of MCMC methods illustrated by a Bayesian analysis of the Hui–Walter model for evaluating diagnostic tests in the absence of a gold standard. The Hui–Walter model has two optimal solutions, a property which causes problems with convergence when the solutions are sufficiently close in the parameter space. Using simulated data we demonstrate tools to assess the convergence and mixing of MCMC chains using examples with and without convergence. Suggestions to remedy the situation when the MCMC sampler fails to converge are given. The epidemiological implications of the two solutions of the Hui–Walter model are discussed.

Assessing Evidence Inconsistency in Mixed Treatment Comparisons

Randomized comparisons among several treatments give rise to an incomplete-blocks structure known as mixed treatment comparisons (MTCs). To analyze such data structures, it is crucial to assess whether the disparate evidence sources provide consistent information about the treatment contrasts. In this article we propose a general method for assessing evidence inconsistency in the framework of Bayesian hierarchical models. We begin with the distinction between basic parameters, which have prior distributions, and functional parameters, which are defined in terms of basic parameters. Based on a graphical analysis of MTC structures, evidence inconsistency is defined as a relation between a functional parameter and at least two basic parameters, supported by at least three evidence sources. The inconsistency degrees of freedom (ICDF) is the number of such inconsistencies. We represent evidence consistency as a set of linear relations between effect parameters on the log odds ratio scale, then relax these relations to allow for inconsistency by adding to the model random inconsistency factors (ICFs). The number of ICFs is determined by the ICDF. The overall consistency between evidence sources can be assessed by comparing models with and without ICFs, whereas their posterior distribution reflects the extent of inconsistency in particular evidence cycles. The methods are elucidated using two published datasets, implemented with standard Markov chain Monte Carlo software.