Probabilistic sensitivity analysis of complex models: a Bayesian approach

Summary.  In many areas of science and technology, mathematical models are built to simulate complex real world phenomena. Such models are typically implemented in large computer programs and are also very complex, such that the way that the model responds to changes in its inputs is not transparent. Sensitivity analysis is concerned with understanding how changes in the model inputs influence the outputs. This may be motivated simply by a wish to understand the implications of a complex model but often arises because there is uncertainty about the true values of the inputs that should be used for a particular application. A broad range of measures have been advocated in the literature to quantify and describe the sensitivity of a model's output to variation in its inputs. In practice the most commonly used measures are those that are based on formulating uncertainty in the model inputs by a joint probability distribution and then analysing the induced uncertainty in outputs, an approach which is known as probabilistic sensitivity analysis. We present a Bayesian framework which unifies the various tools of prob- abilistic sensitivity analysis. The Bayesian approach is computationally highly efficient. It allows effective sensitivity analysis to be achieved by using far smaller numbers of model runs than standard Monte Carlo methods. Furthermore, all measures of interest may be computed from a single set of runs.     

Bayesian emulation of complex multi-output and dynamic computer models

Computer models are widely used in scientific research to study and predict the behaviour of complex systems. The run times of computer-intensive simulators are often such that it is impractical to make the thousands of model runs that are conventionally required for sensitivity analysis, uncertainty analysis or calibration. In response to this problem, highly efficient techniques have recently been developed based on a statistical meta-model (the emulator) that is built to approximate the computer model. The approach, however, is less straightforward for dynamic simulators, designed to represent time-evolving systems. Generalisations of the established methodology to allow for dynamic emulation are here proposed and contrasted. Advantages and difficulties are discussed and illustrated with an application to the Sheffield Dynamic Global Vegetation Model, developed within the UK Centre for Terrestrial Carbon Dynamics.     

Estimating percentiles of uncertain computer code outputs

Summary.  A deterministic computer model is to be used in a situation where there is uncertainty about the values of some or all of the input parameters. This uncertainty induces uncertainty in the output of the model. We consider the problem of estimating a specific percentile of the distribution of this uncertain output. We also suppose that the computer code is computationally expensive, so we can run the model only at a small number of distinct inputs. This means that we must consider our uncertainty about the computer code itself at all untested inputs. We model the output, as a function of its inputs, as a Gaussian process, and after a few initial runs of the code use a simulation approach to choose further suitable design points and to make inferences about the percentile of interest itself. An example is given involving a model that is used in sewer design.     

Quantifying uncertainty in the biospheric carbon flux for England and Wales

Summary.  A crucial issue in the current global warming debate is the effect of vegetation and soils on carbon dioxide (CO₂) concentrations in the atmosphere. Vegetation can extract CO₂ through photosynthesis, but respiration, decay of soil organic matter and disturbance effects such as fire return it to the atmosphere. The balance of these processes is the net carbon flux. To estimate the biospheric carbon flux for England and Wales, we address the statistical problem of inference for the sum of multiple outputs from a complex deterministic computer code whose input parameters are uncertain. The code is a process model which simulates the carbon dynamics of vegetation and soils, including the amount of carbon that is stored as a result of photosynthesis and the amount that is returned to the atmosphere through respiration. The aggregation of outputs corresponding to multiple sites and types of vegetation in a region gives an estimate of the total carbon flux for that region over a period of time. Expert prior opinions are elicited for marginal uncertainty about the relevant input parameters and for correlations of inputs between sites. A Gaussian process model is used to build emulators of the multiple code outputs and Bayesian uncertainty analysis is then used to propagate uncertainty in the input parameters through to uncertainty on the aggregated output. Numerical results are presented for England and Wales in the year 2000. It is estimated that vegetation and soils in England and Wales constituted a net sink of 7.55 Mt C (1 Mt C = 10¹² g of carbon) in 2000, with standard deviation 0.56 Mt C resulting from the sources of uncertainty that are considered.     

Sequential design for response surface model fit in computer experiments using derivative information

This article considers the problem of response surface model fit in computer experiment. We propose a new sequential adaptive design through the “maximum expected improvement” approach. The new method defines the improvement by the first order approximation from the known design points using derivative information and sequentially seeks point in area with large curvature and variance. A version with distance penalty is also considered. We demonstrate their superiority over some existing methods by simulation.     

Reified Bayesian modelling and inference for physical systems

We describe an approach, termed reified analysis, for linking the behaviour of mathematical models with inferences about the physical systems which the models represent. We describe the logical basis for the approach, based on coherent assessment of the implications of deficiencies in the mathematical model. We show how the statistical analysis may be carried out by specifying stochastic relationships between the model that we have, improved versions of the model that we might construct, and the system itself. We illustrate our approach with an example concerning the potential shutdown of the Thermohaline circulation in the Atlantic Ocean.     

Bayesian Tobit quantile regression with single-index models

Based on the Bayesian framework of utilizing a Gaussian prior for the univariate nonparametric link function and an asymmetric Laplace distribution (ALD) for the residuals, we develop a Bayesian treatment for the Tobit quantile single-index regression model (TQSIM). With the location-scale mixture representation of the ALD, the posterior inferences of the latent variables and other parameters are achieved via the Markov Chain Monte Carlo computation method. TQSIM broadens the scope of applicability of the Tobit models by accommodating nonlinearity in the data. The proposed method is illustrated by two simulation examples and a labour supply dataset.     

Sequential tuning of complex computer models

We propose a method that uses a sequential design instead of a space filling design for estimating tuning parameters of a complex computer model. The goal is to bring the computer model output closer to the real system output. The method fits separate Gaussian process (GP) models to the available data from the physical experiment and the computer experiment and minimizes the discrepancy between the predictions from the GP models to obtain estimates of the tuning parameters. A criterion based on the discrepancy between the predictions from the two GP models and the standard error of prediction for the computer experiment output is then used to obtain a design point for the next run of the computer experiment. The tuning parameters are re-estimated using the augmented data set. The steps are repeated until the budget for the computer experiment data is exhausted. Simulation studies show that the proposed method performs better in bringing a computer model closer to the real system than methods that use a space filling design.     

Bayesian integrative analysis for multi-fidelity computer experiments

This paper proposes a Bayesian integrative analysis method for linking multi-fidelity computer experiments. Instead of assuming covariance structures of multivariate Gaussian process models, we handle the outputs from different levels of accuracy as independent processes and link them via a penalization method that controls the distance between their overall trends. Based on the priors induced by the penalty, we build Bayesian prediction models for the output at the highest accuracy. Simulated and real examples show that the proposed method is better than existing methods in terms of prediction accuracy for many cases.     

On asymptotic properties of Bayesian partially linear models

In this paper, we present large sample properties of a partially linear model from the Bayesian perspective, in which responses are explained by the semiparametric regression model with the additive form of the linear component and the nonparametric component. For this purpose, we investigate asymptotic behaviors of posterior distributions in terms of consistency. Specifically, we deal with a specific Bayesian partially linear regression model with additive noises in which the nonparametric component is modeled using Gaussian process priors. Under the Bayesian partially linear model using Gaussian process priors, we focus on consistency of posterior distribution and consistency of the Bayes factor, and extend these results to generalized additive regression models and study their asymptotic properties. In addition we illustrate the asymptotic properties based on empirical analysis through simulation studies.     

Estimation across multiple models with application to Bayesian computing and software development

Statistical models are sometimes incorporated into computer software for making predictions about future observations. When the computer model consists of a single statistical model this corresponds to estimation of a function of the model parameters. This paper is concerned with the case that the computer model implements multiple, individually-estimated statistical sub-models. This case frequently arises, for example, in models for medical decision making that derive parameter information from multiple clinical studies. We develop a method for calculating the posterior mean of a function of the parameter vectors of multiple statistical models that is easy to implement in computer software, has high asymptotic accuracy, and has a computational cost linear in the total number of model parameters. The formula is then used to derive a general result about posterior estimation across multiple models. The utility of the results is illustrated by application to clinical software that estimates the risk of fatal coronary disease in people with diabetes.     

Sensitivity Analysis in Gaussian Bayesian Networks Using a Divergence Measure

This article develops a method for computing the sensitivity analysis in a Gaussian Bayesian network. The measure presented is based on the Kullback–Leibler divergence and is useful to evaluate the impact of prior changes over the posterior marginal density of the target variable in the network. We find that some changes do not disturb the posterior marginal density of interest. Finally, we describe a method to compare different sensitivity measures obtained depending on where the inaccuracy was. An example is used to illustrate the concepts and methods presented.     

Bayesian Analysis of DSGE Models

This paper reviews Bayesian methods that have been developed in recent years to estimate and evaluate dynamic stochastic general equilibrium (DSGE) models. We consider the estimation of linearized DSGE models, the evaluation of models based on Bayesian model checking, posterior odds comparisons, and comparisons to vector autoregressions, as well as the non-linear estimation based on a second-order accurate model solution. These methods are applied to data generated from correctly specified and misspecified linearized DSGE models and a DSGE model that was solved with a second-order perturbation method.     

Screening Experiments with Many Factors

A new procedure is introduced for conducting screening experiments to find a small number of influential factors from among a large number of factors with negligible effects. It is intended for experiments in which the factors are easily controlled, as in simulation models. It adds observations sequentially after conducting a small initial experiment. The performance of the procedure is investigated using simulation, and evidence is presented that this and other procedures scale as the logarithm of the total number of factors if the number of influential factors is fixed. An investigation of the new procedure for 1–3 active factors shows that it compares favorably with competing methods, particularly when the size of the nonzero effects is 1–2 times the standard deviation. A limited look at the procedure for up to 6 active factors is also presented.     

Bayesian spatio‐temporal models based on discrete convolutions

Abstract: The authors consider a class of models for spatio‐temporal processes based on convolving independent processes with a discrete kernel that is represented by a lower triangular matrix. They study two families of models. In the first one, spatial Gaussian processes with isotropic correlations are convoluted with a kernel that provides temporal dependencies. In the second family, AR(p) processes are convoluted with a kernel providing spatial interactions. The covariance structures associated with these two families are quite rich. Their covariance functions that are stationary and separable in space and time as well as time dependent nonseparable and nonisotropic ones.     

Bayesian hierarchical models for food frequency assessment

The aim of this study is to assess the biases of a Food Frequency Questionnaire (FFQ) by comparing total energy intake (TEI) with total energy expenditure (TEE) obtained from doubly labelled water(DLW) biomarker after adjusting measurement errors in DLW. We develop several Bayesian hierarchical measurement error models of DLW with different distributional assumptions on TEI to obtain precise bias estimates of TEI. Inference is carried out by using MCMC simulation techniques in a fully Bayesian framework, and model comparisons are done via the mean square predictive error. Our results showed that the joint model with random effects under the Gamma distribution is the best fit model in terms of the MSPE and residual diagnostics, in which bias in TEI is not significant based on the 95% credible interval. The Canadian Journal of Statistics 38: 506–516; 2010 © 2010 Statistical Society of Canada     

Batch sequential designs for computer experiments

Computer models simulating a physical process are used in many areas of science. Due to the complex nature of these codes it is often necessary to approximate the code, which is typically done using a Gaussian process. In many situations the number of code runs available to build the Gaussian process approximation is limited. When the initial design is small or the underlying response surface is complicated this can lead to poor approximations of the code output. In order to improve the fit of the model, sequential design strategies must be employed. In this paper we introduce two simple distance based metrics that can be used to augment an initial design in a batch sequential manner. In addition we propose a sequential updating strategy to an orthogonal array based Latin hypercube sample. We show via various real and simulated examples that the distance metrics and the extension of the orthogonal array based Latin hypercubes work well in practice.     

Evaluation of Bayesian multiple stage estimation under spatial CAR model variants

In this study, an evaluation of Bayesian hierarchical models is made based on simulation scenarios to compare single-stage and multi-stage Bayesian estimations. Simulated datasets of lung cancer disease counts for men aged 65 and older across 44 wards in the London Health Authority were analysed using a range of spatially structured random effect components. The goals of this study are to determine which of these single-stage models perform best given a certain simulating model, how estimation methods (single- vs. multi-stage) compare in yielding posterior estimates of fixed effects in the presence of spatially structured random effects, and finally which of two spatial prior models – the Leroux or ICAR model, perform best in a multi-stage context under different assumptions concerning spatial correlation. Among the fitted single-stage models without covariates, we found that when there is low amount of variability in the distribution of disease counts, the BYM model is relatively robust to misspecification in terms of DIC, while the Leroux model is the least robust to misspecification. When these models were fit to data generated from models with covariates, we found that when there was one set of covariates – either spatially correlated or non-spatially correlated, changing the values of the fixed coefficients affected the ability of either the Leroux or ICAR model to fit the data well in terms of DIC. When there were multiple sets of spatially correlated covariates in the simulating model, however, we could not distinguish the goodness of fit to the data between these single-stage models. We found that the multi-stage modelling process via the Leroux and ICAR models generally reduced the variance of the posterior estimated fixed effects for data generated from models with covariates and a UH term compared to analogous single-stage models. Finally, we found the multi-stage Leroux model compares favourably to the multi-stage ICAR model in terms of DIC. We conclude that the mutli-stage Leroux model should be seriously considered in applications of Bayesian disease mapping when an investigator desires to fit a model with both fixed effects and spatially structured random effects to Poisson count data.     

Inverse Gaussian process models for bivariate degradation analysis: A Bayesian perspective

This article conducts a Bayesian analysis for bivariate degradation models based on the inverse Gaussian (IG) process. Assume that a product has two quality characteristics (QCs) and each of the QCs is governed by an IG process. The dependence of the QCs is described by a copula function. A bivariate simple IG process model and three bivariate IG process models with random effects are investigated by using Bayesian method. In addition, a simulation example is given to illustrate the effectiveness of the proposed methods. Finally, an example about heavy machine tools is presented to validate the proposed models.     

Convergence properties of the expected improvement algorithm with fixed mean and covariance functions

This paper deals with the convergence of the expected improvement algorithm, a popular global optimization algorithm based on a Gaussian process model of the function to be optimized. The first result is that under some mild hypotheses on the covariance function k of the Gaussian process, the expected improvement algorithm produces a dense sequence of evaluation points in the search domain, when the function to be optimized is in the reproducing kernel Hilbert space generated by k  . The second result states that the density property also holds for P-almost$\mathsf{P}\text{-almost}$ all continuous functions, where P is the (prior) probability distribution induced by the Gaussian process.