Sandwich algorithms for Bayesian variable selection

Bayesian model updating provides a rigorous framework to account for uncertainty induced by lack of knowledge about engineering systems in their respective mathematical models through updates of the joint probability density function (PDF), the so-called posterior PDF, of the unknown model parameters. The Markov chain Monte Carlo (MCMC) methods are currently the most popular approaches for generating samples from the posterior PDF. However, these methods often found wanting when sampling from difficult distributions (e.g., high-dimensional PDFs, PDFs with flat manifolds, multimodal PDFs, and very peaked PDFs). This paper introduces a new multi-level sampling approach for Bayesian model updating, called Sequential Gauss-Newton algorithm, which is inspired by the Transitional Markov chain Monte Carlo (TMCMC) algorithm. The Sequential Gauss-Newton algorithm improves two aspects of TMCMC to make an efficient and effective MCMC algorithm for drawing samples from difficult posterior PDFs. First, the statistical efficiency of the algorithm is enhanced by use of the systematic resampling scheme. Second, a new MCMC algorithm, called Gauss-Newton MCMC algorithm, is proposed which is essentially an M-H algorithm with a Gaussian proposal PDF tailored to the posterior PDF using the gradient and Hessian information of the negative log posterior. The effectiveness of the proposed algorithm for solving the Bayesian model updating problem is illustrated using three examples with irregularly shaped posterior PDFs.

In an effort to build seismic models that are the most consistent with multiple data sets we have applied a new probabilistic inverse technique. This method uses a Markov chain Monte Carlo (MCMC) algorithm to sample models from a prior distribution and test them against multiple data types to generate a posterior distribution. While computationally expensive, this approach has several advantages over deterministic models, notably the seamless reconciliation of different data types that constrain the model, the proper handling of both data and model uncertainties, and the ability to easily incorporate a variety of prior information, all in a straightforward, natural fashion. A real advantage of the technique is that it provides a more complete picture of the solution space. By mapping out the posterior probability density function, we can avoid simplistic assumptions about the model space and allow alternative solutions to be identified, compared, and ranked. Here we use this method to determine the crust and upper mantle structure of the Yellow Sea and Korean Peninsula region. The model is parameterized as a series of seven layers in a regular latitude‐longitude grid, each of which is characterized by thickness and seismic parameters (Vp, Vs, and density). We use surface wave dispersion and body wave traveltime data to drive the model. We find that when properly tuned (i.e., the Markov chains have had adequate time to fully sample the model space and the inversion has converged), the technique behaves as expected. The posterior model reflects the prior information at the edge of the model where there is little or no data to constrain adjustments, but the range of acceptable models is significantly reduced in data‐rich regions, producing values of sediment thickness, crustal thickness, and upper mantle velocities consistent with expectations based on knowledge of the regional tectonic setting.

Sampling across tree space is one of the major challenges in Bayesian phylogenetic inference using Markov chain Monte Carlo (MCMC) algorithms. Standard MCMC tree moves consider small random perturbations of the topology, and select from candidate trees at random or based on the distance between the old and new topologies. MCMC algorithms using such moves tend to get trapped in tree space, making them slow in finding the globally most probable trees (known as “convergence”) and in estimating the correct proportions of the different types of them (known as “mixing”). Here, we introduce a new class of moves, which propose trees based on their parsimony scores. The proposal distribution derived from the parsimony scores is a quickly computable albeit rough approximation of the conditional posterior distribution over candidate trees. We demonstrate with simulations that parsimony-guided moves correctly sample the uniform distribution of topologies from the prior. We then evaluate their performance against standard moves using six challenging empirical data sets, for which we were able to obtain accurate reference estimates of the posterior using long MCMC runs, a mix of topology proposals, and Metropolis coupling. On these data sets, ranging in size from 357 to 934 taxa and from 1740 to 5681 sites, we find that single chains using parsimony-guided moves usually converge an order of magnitude faster than chains using standard moves. They also exhibit better mixing, that is, they cover the most probable trees more quickly. Our results show that tree moves based on quick and dirty estimates of the posterior probability can significantly outperform standard moves. Future research will have to show to what extent the performance of such moves can be improved further by finding better ways of approximating the posterior probability, taking the trade-off between accuracy and speed into account. [Bayesian phylogenetic inference; MCMC; parsimony; tree proposal.]

Variable selection in regression models has been well studied in the literature, with many non‐Bayesian and Bayesian methods available in this regard. An important class of regression models is generalized linear models, which involve situations where the response variable is discrete. To add more flexibility, generalized additive partial linear models can be considered, where some predictors can have a non‐linear effect while some predictors have a strictly linear effect. We consider Bayesian variable selection in these models. The functions in the non‐parametric additive part of the model are expanded in a B‐spline basis and multivariate Laplace prior put on the coefficients with point mass at zero. The coefficients corresponding to the strictly linear components are assigned a univariate Laplace prior with point mass at zero. The prior times the likelihood is mathematically intractable, but we find an approximation by expansion around the posterior mode, which is the group lasso solution in generalized linear model setting for the choice of prior. We thus completely avoid Markov chain Monte Carlo methods, which are extremely slow and unreliable in high‐dimensional models. We evaluate the performance of the Bayesian method by conducting simulation studies and real data analyses. Copyright © 2014 John Wiley & Sons, Ltd.

Recent progress using geometry in the design of efficient Markov chain Monte Carlo (MCMC) algorithms have shown the effectiveness of the Fisher Riemannian structure. Furthermore, the theory of the underlying geometry of spaces of statistical models has made an important breakthrough by extending the classical theory on exponential families to their closures, the so-called extended exponential families. This paper looks at the underlying geometry of the Fisher information, in particular its limiting behaviour near boundaries, which illuminates the excellent behaviour of the corresponding geometric MCMC algorithms. Further, the paper shows how Fisher geodesics in extended exponential families smoothly attach the boundaries of extended exponential families to their relative interior. We conjecture that this behaviour could be exploited for trans-dimensional MCMC algorithms.

Analyses of multivariate ordinal probit models typically use data augmentation to link the observed (discrete) data to latent (continuous) data via a censoring mechanism defined by a collection of “cutpoints.” Most standard models, for which effective Markov chain Monte Carlo (MCMC) sampling algorithms have been developed, use a separate (and independent) set of cutpoints for each element of the multivariate response. Motivated by the analysis of ratings data, we describe a particular class of multivariate ordinal probit models where it is desirable to use a common set of cutpoints. While this approach is attractive from a data-analytic perspective, we show that the existing efficient MCMC algorithms can no longer be accurately applied. Moreover, we show that attempts to implement these algorithms by numerically approximating required multivariate normal integrals over high-dimensional rectangular regions can result in severely degraded estimates of the posterior distribution. We propose a new data augmentation that is based on a covariance decomposition and that admits a simple and accurate MCMC algorithm. Our data augmentation requires only that univariate normal integrals be evaluated, which can be done quickly and with high accuracy. We provide theoretical results that suggest optimal decompositions within this class of data augmentations, and, based on the theory, recommend default decompositions that we demonstrate work well in practice. This article has supplementary material online.

Markov Chain Monte Carlo (MCMC) is widely applied to fit complex, practical models in various disciplines, such as engineering, physics, agriculture, etc. Indeed, when the target density $f$ is not available in closed form, a Markov chain is constructed with stationary density $f$ for forming Monte Carlo estimators of means with respect to $f$. Although a Harris ergodic Markov chain will converge to the target distribution as the Monte Carlo sample size tends to infinity, and the population means can be consistently estimated by the corresponding MCMC estimators, in practice, it is important to estimate the standard errors associated with these Monte Carlo estimates due to the availability of only a finite sample. Therefore, the geometric ergodicity of a Markov Chain is highly desired as it is the most standard method for guaranteeing a central limit theorem (CLT) for MCMC estimators and finding its standard errors. Among the different MCMC algorithms, Metropolis-Hastings (MH) algorithms are predominant. This thesis establishes conditions under which MH algorithms with a position-dependent proposal covariance matrix will or will not have the geometric convergence rate. Indeed, we provide these results for MH algorithms where the normal proposal density has a general mean function $c(x)$ and a covariance matrix $G(x)$ depending on the current position $x$ of the Markov chain. For modern variants of the popular Metropolis adjusted Langevin algorithms (MALA), like the manifold MALA (MMALA) that adapt to the geometry of the target distributions, the proposal covariance matrix changes in every iteration. Thus, we provide conditions for the geometric ergodicity of different variations of the Langevin algorithms, e.g., MMALA, position-dependent MALA (PMALA), and the pre-conditioned MALA (PCMALA). We also provide conditions guaranteeing geometric ergodicity of the pre-conditioned unadjusted Langevin algorithms (PCULA). Generalized linear mixed models (GLMMs), where the response variable is allowed to follow a non-Gaussian distribution, are popular statistical models, although their likelihood functions are not available in closed form. GLMMs can handle the over-dispersion often present in real-world data, and MCMC methods are generally used to fit frequentist and Bayesian GLMMs. Since the likelihood function for GLMMs is not available in closed form, Monte Carlo EM or Monte Carlo maximum likelihood methods are often used for making inferences from GLMMs. Simulation from the conditional (on the data) distribution of the random effects is, in turn, needed for implementing the Monte Carlo EM or Monte Carlo maximum likelihood methods. In this thesis, our general conditions for geometric ergodicity of variants of Langevin algorithms are also verified in the context of conditional simulation from the two most popular GLMMs, namely the binomial GLMM with the logit link and the Poisson GLMM with the log link. These results have important practical implications as users can now confidently fit these models by attaching asymptotically valid standard errors with their Markov chain-based Monte Carlo estimates. Spatial generalized linear mixed models (SGLMMs) are GLMMs where the random effects form a spatial process and are useful for modeling spatially correlated data. We compare the performance of different variants of MALA mentioned above as well as other popular MH algorithms like Random Walk Metropolis in the context of conditional simulation from SGLMMs. Next, we consider MCMC algorithms for full Bayesian analysis of GLMMs. Blocking and reparameterization are techniques used to improve the efficiency and convergence properties of MCMC algorithms. We combine these strategies with pre-conditioning to develop efficient Langevin algorithms for Bayesian GLMMs and SGLMMs. We illustrate our algorithms using extensive real and simulated data analyses.

One of the most common types of models that helps us to understand neuron behavior is based on the Hodgkin–Huxley ion channel formulation (HH model). A major challenge with inferring parameters in HH models is non-uniqueness: many different sets of ion channel parameter values produce similar outputs for the same input stimulus. Such phenomena result in an objective function that exhibits multiple modes (i.e., multiple local minima). This non-uniqueness of local optimality poses challenges for parameter estimation with many algorithmic optimization techniques. HH models additionally have severe non-linearities resulting in further challenges for inferring parameters in an algorithmic fashion. To address these challenges with a tractable method in high-dimensional parameter spaces, we propose using a particular Markov chain Monte Carlo (MCMC) algorithm, which has the advantage of inferring parameters in a Bayesian framework. The Bayesian approach is designed to be suitable for multimodal solutions to inverse problems. We introduce and demonstrate the method using a three-channel HH model. We then focus on the inference of nine parameters in an eight-channel HH model, which we analyze in detail. We explore how the MCMC algorithm can uncover complex relationships between inferred parameters using five injected current levels. The MCMC method provides as a result a nine-dimensional posterior distribution, which we analyze visually with solution maps or landscapes of the possible parameter sets. The visualized solution maps show new complex structures of the multimodal posteriors, and they allow for selection of locally and globally optimal value sets, and they visually expose parameter sensitivities and regions of higher model robustness. We envision these solution maps as enabling experimentalists to improve the design of future experiments, increase scientific productivity and improve on model structure and ideation when the MCMC algorithm is applied to experimental data.

Traditional full-waveform inversion (FWI) methods only render a “best-fit” model that cannot account for uncertainties of the ill-posed inverse problem. Additionally, local optimization-based FWI methods cannot always converge to a geologically meaningful solution unless the inversion starts with an accurate background model. We seek the solution for FWI in the Bayesian inference framework to address those two issues. In Bayesian inference, the model space is directly probed by sampling methods such that we obtain a reliable uncertainty appraisal, determine optimal models, and avoid entrapment in a small local region of the model space. The solution of such a statistical inverse method is completely described by the posterior distribution, which quantifies the distributions for parameters and inversion uncertainties. To efficiently sample the posterior distribution, we introduce a sampling algorithm in which the proposal distribution is constructed by the local gradient and the diagonal approximate Hessian of the local log posterior. Our algorithm is called the gradient-based Markov chain Monte Carlo (GMCMC) method. The GMCMC FWI method can quantify inversion uncertainties with estimated posterior distribution given sufficiently long Markov chains. By directly sampling the posterior distribution, we obtain a global view of the model space. Theoretically speaking, statistical assessments do not depend on starting models. Our method is applied to the 2D Marmousi model with the frequency-domain FWI setting. Numerical results suggest that our method can be readily applied to 2D cases with affordable computational efforts.

In high-throughput genetics studies, an important aim is to identify gene–environment interactions associated with the clinical outcomes. Recently, multiple marginal penalization methods have been developed and shown to be effective in G×E studies. However, within the Bayesian framework, marginal variable selection has not received much attention. In this study, we propose a novel marginal Bayesian variable selection method for G×E studies. In particular, our marginal Bayesian method is robust to data contamination and outliers in the outcome variables. With the incorporation of spike-and-slab priors, we have implemented the Gibbs sampler based on Markov Chain Monte Carlo (MCMC). The proposed method outperforms a number of alternatives in extensive simulation studies. The utility of the marginal robust Bayesian variable selection method has been further demonstrated in the case studies using data from the Nurse Health Study (NHS). Some of the identified main and interaction effects from the real data analysis have important biological implications.

Selecting a small number of relevant genes for classification has received a great deal of attention in microarray data analysis. While the development of methods for microarray data with only two classes is relevant, developing more efficient algorithms for classification with any number of classes is important. In this paper, we propose a Bayesian stochastic search variable selection approach for multi-class classification, which can identify relevant genes by assessing sets of genes jointly. We consider a multinomial probit model with a generalized $$g$$g-prior for the regression coefficients. An efficient algorithm using simulation-based MCMC methods are developed for simulating parameters from the posterior distribution. This algorithm is robust to the choice of initial value, and produces posterior probabilities of relevant genes for biological interpretation. We demonstrate the performance of the approach with two well-known gene expression profiling data: leukemia data, lymphoma data, SRBCTs data and NCI60 data. Compared with other classification approaches, our approach selects smaller numbers of relevant genes and obtains competitive classification accuracy based on obtained results.

The paper presents a novel method for rapid quantification of uncertainty in history matching reservoir models using a two-stage Markov Chain Monte Carlo (MCMC) method. Our approach is based on a combination of fast linearized approximation to the dynamic data and the MCMC algorithm. In the first stage, we use streamline-derived sensitivities to obtain an analytical approximation in a small neighborhood of the previously computed dynamic data. The sensitivities can be conveniently obtained using either a finite-difference or streamline simulator. The approximation of the dynamic data is then used to modify the instrumental proposal distribution during MCMC. In the second stage, those proposals that pass the first stage are assessed by running full flow simulations to assure rigorousness in sampling. The uncertainty analysis is carried out by analyzing multiple models sampled from the posterior distribution in the Bayesian formulation for history matching. We demonstrate that the two-stage approach increases the acceptance rate, and significantly reduces the computational cost compared to conventional MCMC sampling without sacrificing accuracy. Finally, both two-dimensional synthetic and three-dimensional field examples demonstrate the power and utility of the two-stage MCMC method for history matching and uncertainty analysis.

Bayesian model updating has received considerable attention and has been extensively used in structural damage detection. It provides a rigorous statistical framework for realizing structural system identification and characterizing uncertainties associated with modeling and measurements. The Markov Chain Monte Carlo (MCMC) is a promising tool for inferring the posterior distribution of model parameters to avoid the intractable evaluation of multi-dimensional integration. However, the efficacy of most MCMC techniques suffers from the curse of parameter dimension, which restricts the application of Bayesian model updating to the damage detection of large-scale systems. In addition, there are several MCMC techniques that require users to properly choose application-specific models, based on the understanding of algorithm mechanisms and limitations. As seen in the literature, there is a lack of comprehensive work that investigates the performances of various MCMC algorithms in their application of structural damage detection. In this study, the Differential Evolutionary Adaptive Metropolis (DREAM), a multi-chain MCMC, is explored and adapted to Bayesian model updating. This paper illustrates how DREAM is used for model updating with many uncertainty parameters (i.e., 40 parameters). Furthermore, the study provides a tutorial to users who may be less experienced with Bayesian model updating and MCMC. Two advanced single-chain MCMC algorithms, namely, the Delayed Rejection Adaptive Metropolis (DRAM) and Transitional Markov Chain Monte Carlo (TMCMC), and DREAM are elaborately introduced to allow practitioners to understand better the concepts and practical implementations. Their performances in model updating and damage detection are compared through three different engineering applications with increased complexity, e.g., a forty-story shear building, a two-span continuous steel beam, and a large-scale steel pedestrian bridge.

The Bayesian approach to inverse problems, in which the posterior probability distribution on an unknown field is sampled for the purposes of computing posterior expectations of quantities of interest, is starting to become computationally feasible for partial differential equation (PDE) inverse problems. Balancing the sources of error arising from finite-dimensional approximation of the unknown field, the PDE forward solution map and the sampling of the probability space under the posterior distribution are essential for the design of efficient computational Bayesian methods for PDE inverse problems. We study Bayesian inversion for a model elliptic PDE with an unknown diffusion coefficient. We provide complexity analyses of several Markov chain Monte Carlo (MCMC) methods for the efficient numerical evaluation of expectations under the Bayesian posterior distribution, given data δ. Particular attention is given to bounds on the overall work required to achieve a prescribed error level ε. Specifically, we first bound the computational complexity of ‘plain’ MCMC, based on combining MCMC sampling with linear complexity multi-level solvers for elliptic PDE. Our (new) work versus accuracy bounds show that the complexity of this approach can be quite prohibitive. Two strategies for reducing the computational complexity are then proposed and analyzed: first, a sparse, parametric and deterministic generalized polynomial chaos (gpc) ‘surrogate’ representation of the forward response map of the PDE over the entire parameter space, and, second, a novel multi-level Markov chain Monte Carlo strategy which utilizes sampling from a multi-level discretization of the posterior and the forward PDE. For both of these strategies, we derive asymptotic bounds on work versus accuracy, and hence asymptotic bounds on the computational complexity of the algorithms. In particular, we provide sufficient conditions on the regularity of the unknown coefficients of the PDE and on the approximation methods used, in order for the accelerations of MCMC resulting from these strategies to lead to complexity reductions over ‘plain’ MCMC algorithms for the Bayesian inversion of PDEs.

We study the computational complexity of Markov chain Monte Carlo (MCMC) methods for high-dimensional Bayesian linear regression under sparsity constraints. We first show that a Bayesian approach can achieve variable-selection consistency under relatively mild conditions on the design matrix. We then demonstrate that the statistical criterion of posterior concentration need not imply the computational desideratum of rapid mixing of the MCMC algorithm. By introducing a truncated sparsity prior for variable selection, we provide a set of conditions that guarantee both variable-selection consistency and rapid mixing of a particular Metropolis-Hastings algorithm. The mixing time is linear in the number of covariates up to a logarithmic factor. Our proof controls the spectral gap of the Markov chain by constructing a canonical path ensemble that is inspired by the steps taken by greedy algorithms for variable selection.