Fast sampling of Gaussian Markov random fields

This paper demonstrates how Gaussian Markov random fields (conditional autoregressions) can be sampled quickly by using numerical techniques for sparse matrices. The algorithm is general and efficient, and expands easily to various forms for conditional simulation and evaluation of normalization constants. We demonstrate its use by constructing efficient block updates in Markov chain Monte Carlo algorithms for disease mapping.     

Approximating hidden Gaussian Markov random fields

Summary.  Gaussian Markov random-field (GMRF) models are frequently used in a wide variety of applications. In most cases parts of the GMRF are observed through mutually independent data; hence the full conditional of the GMRF, a hidden GMRF (HGMRF), is of interest. We are concerned with the case where the likelihood is non-Gaussian, leading to non-Gaussian HGMRF models. Several researchers have constructed block sampling Markov chain Monte Carlo schemes based on approximations of the HGMRF by a GMRF, using a second-order expansion of the log-density at or near the mode. This is possible as the GMRF approximation can be sampled exactly with a known normalizing constant. The Markov property of the GMRF approximation yields computational efficiency.The main contribution in the paper is to go beyond the GMRF approximation and to construct a class of non-Gaussian approximations which adapt automatically to the particular HGMRF that is under study. The accuracy can be tuned by intuitive parameters to nearly any precision. These non-Gaussian approximations share the same computational complexity as those which are based on GMRFs and can be sampled exactly with computable normalizing constants. We apply our approximations in spatial disease mapping and model-based geostatistical models with different likelihoods, obtain procedures for block updating and construct Metropolized independence samplers.     

Gaussian Markov random field spatial models in GAMLSS

This paper describes the modelling and fitting of Gaussian Markov random field spatial components within a Generalized AdditiveModel for Location, Scale and Shape (GAMLSS) model. This allows modelling of any or all the parameters of the distribution for the response variable using explanatory variables and spatial effects. The response variable distribution is allowed to be a non-exponential family distribution. A new package developed in R to achieve this is presented. We use Gaussian Markov random fields to model the spatial effect in Munich rent data and explore some features and characteristics of the data. The potential of using spatial analysis within GAMLSS is discussed. We argue that the flexibility of parametric distributions, ability to model all the parameters of the distribution and diagnostic tools of GAMLSS provide an ideal environment for modelling spatial features of data.     

Bayesian inference for generalized additive mixed models based on Markov random field priors

Most regression problems in practice require flexible semiparametric forms of the predictor for modelling the dependence of responses on covariates. Moreover, it is often necessary to add random effects accounting for overdispersion caused by unobserved heterogeneity or for correlation in longitudinal or spatial data. We present a unified approach for Bayesian inference via Markov chain Monte Carlo simulation in generalized additive and semiparametric mixed models. Different types of covariates, such as the usual covariates with fixed effects, metrical covariates with non-linear effects, unstructured random effects, trend and seasonal components in longitudinal data and spatial covariates, are all treated within the same general framework by assigning appropriate Markov random field priors with different forms and degrees of smoothness. We applied the approach in several case-studies and consulting cases, showing that the methods are also computationally feasible in problems with many covariates and large data sets. In this paper, we choose two typical applications.     

Modeling material stress using integrated Gaussian Markov random fields

The equations of a physical constitutive model for material stress within tantalum grains were solved numerically using a tetrahedrally meshed volume. The resulting output included a scalar vonMises stress for each of the more than 94,000 tetrahedra within the finite element discretization. In this paper, we define an intricate statistical model for the spatial field of vonMises stress which uses the given grain geometry in a fundamental way. Our model relates the three-dimensional field to integrals of latent stochastic processes defined on the vertices of the one- and two-dimensional grain boundaries. An intuitive neighborhood structure of the said boundary nodes suggested the use of a latent Gaussian Markov random field (GMRF). However, despite the potential for computational gains afforded by GMRFs, the integral nature of our model and the sheer number of data points pose substantial challenges for a full Bayesian analysis. To overcome these problems and encourage efficient exploration of the posterior distribution, a number of techniques are now combined: parallel computing, sparse matrix methods, and a modification of a block update strategy within the sampling routine. In addition, we use an auxiliary variables approach to accommodate the presence of outliers in the data.     

Selection of the neighborhood structure for space-time Markov random field models

A space-time, univariate dataset is assumed to have been sampled from a 3-dimensional Markov Random Field where the data dependence structure is modeled through pairwise interaction parameters. The likelihood function depends upon (1) an undirected, 3-dimensional graph, where edges connect observation points, and (2) the parameter dimension that captures possible space-time anisotropy of data interaction. Automatic model selection to discriminate both the graph and the model dimension is suggested on the basis of a penalized Pseudo-likelihood function. In most cases, the procedure can be implemented using standard statistical packages capable of GLM estimation. Weak consistency of the criterion is shown to hold under mild and easily verifiable sufficient conditions. Its performance in small samples is studied providing simulation results.     

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.     

On Block Updating in Markov Random Field Models for Disease Mapping

Gaussian Markov random field (GMRF) models are commonly used to model spatial correlation in disease mapping applications. For Bayesian inference by MCMC, so far mainly single-site updating algorithms have been considered. However, convergence and mixing properties of such algorithms can be extremely poor due to strong dependencies of parameters in the posterior distribution. In this paper, we propose various block sampling algorithms in order to improve the MCMC performance. The methodology is rather general, allows for non-standard full conditionals, and can be applied in a modular fashion in a large number of different scenarios. For illustration we consider three different applications: two formulations for spatial modelling of a single disease (with and without additional unstructured parameters respectively), and one formulation for the joint analysis of two diseases. The results indicate that the largest benefits are obtained if parameters and the corresponding hyperparameter are updated jointly in one large block. Implementation of such block algorithms is relatively easy using methods for fast sampling of Gaussian Markov random fields ( Rue, 2001 ). By comparison, Monte Carlo estimates based on single-site updating can be rather misleading, even for very long runs. Our results may have wider relevance for efficient MCMC simulation in hierarchical models with Markov random field components.     

Singular Gaussian graphical models: Structure learning

ABSTRACT

The goal of this article is to introduce singular Gaussian graphical models and their conditional independence properties. In fact, we extend the concept of Gaussian Markov Random Field to the case of a multivariate normally distributed vector with a singular covariance matrix. We construct, then, the associated graph’s structure from the covariance matrix’s pseudo-inverse on the basis of a characterization of the pairwise conditional independence. The proposed approach can also be used when the covariance matrix is ill-conditioned, through projecting data on a smaller subspace. In this case, our method ensures numerical stability and consistency of the constructed graph and significantly reduces the inference problem’s complexity. These aspects are illustrated using numerical experiments.     

A comparison of Bayesian model selection based on MCMC with an application to GARCH-type models

This paper presents a comprehensive review and comparison of five computational methods for Bayesian model selection, based on MCMC simulations from posterior model parameter distributions. We apply these methods to a well-known and important class of models in financial time series analysis, namely GARCH and GARCH-t models for conditional return distributions (assuming normal and t-distributions). We compare their performance with the more common maximum likelihood-based model selection for simulated and real market data. All five MCMC methods proved reliable in the simulation study, although differing in their computational demands. Results on simulated data also show that for large degrees of freedom (where the t-distribution becomes more similar to a normal one), Bayesian model selection results in better decisions in favor of the true model than maximum likelihood. Results on market data show the instability of the harmonic mean estimator and reliability of the advanced model selection methods.     

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.     

Longitudinal network models and permutation-uniform Markov chains

Consider longitudinal networks whose edges turn on and off according to a discrete-time Markov chain with exponential-family transition probabilities. We characterize when their joint distributions are also exponential families with the same parameter, improving data reduction. Further we show that the permutation-uniform subclass of these chains permit interpretation as an independent, identically distributed sequence on the same state space. We then apply these ideas to temporal exponential random graph models, for which permutation uniformity is well suited, and discuss mean-parameter convergence, dyadic independence, and exchangeability. Our framework facilitates our introducing a new network model; simplifies analysis of some network and autoregressive models from the literature, including by permitting closed-form expressions for maximum likelihood estimates for some models; and facilitates applying standard tools to longitudinal-network Markov chains from either asymptotics or single-observation exponential random graph models.     

Ecological regression analysis of environmental benzene exposure and childhood leukaemia: sensitivity to data inaccuracies, geographical scale and ecological bias

Benzene is classified as a group 1 human carcinogen by the International Agency for Research on Cancer, and it is now accepted that occupational exposure is associated with an increased risk of various leukaemias. However, occupational exposure accounts for less than 1% of all benzene exposures, the major sources being cigarette smoking and vehicle exhaust emissions. Whether such low level exposures to environmental benzene are also associated with the risk of leukaemia is currently not known. In this study, we investigate the relationship between benzene emissions arising from outdoor sources (predominantly road traffic and petrol stations) and the incidence of childhood leukaemia in Greater London. An ecological design was used because of the rarity of the disease, the difficulty of obtaining individual level measurements of benzene exposure and the availability of data. However, some methodological difficulties were encountered, including problems of case registration errors, the choice of geographical areas for analysis, exposure measurement errors and ecological bias. We use a Bayesian hierarchical modelling framework to address these issues, and we investigate the sensitivity of our inference to various modelling assumptions.     

Comparison of three-dimensional profiles over time

In this paper, we describe an analysis for data collected on a three-dimensional spatial lattice with treatments applied at the horizontal lattice points. Spatial correlation is accounted for using a conditional autoregressive model. Observations are defined as neighbours only if they are at the same depth. This allows the corresponding variance components to vary by depth. We use the Markov chain Monte Carlo method with block updating, together with Krylov subspace methods, for efficient estimation of the model. The method is applicable to both regular and irregular horizontal lattices and hence to data collected at any set of horizontal sites for a set of depths or heights, for example, water column or soil profile data. The model for the three-dimensional data is applied to agricultural trial data for five separate days taken roughly six months apart in order to determine possible relationships over time. The purpose of the trial is to determine a form of cropping that leads to less moist soils in the root zone and beyond. We estimate moisture for each date, depth and treatment accounting for spatial correlation and determine relationships of these and other parameters over time.     

Propriety of posteriors in structured additive regression models: Theory and empirical evidence

Structured additive regression comprises many semiparametric regression models such as generalized additive (mixed) models, geoadditive models, and hazard regression models within a unified framework. In a Bayesian formulation, non-parametric functions, spatial effects and further model components are specified in terms of multivariate Gaussian priors for high-dimensional vectors of regression coefficients. For several model terms, such as penalized splines or Markov random fields, these Gaussian prior distributions involve rank-deficient precision matrices, yielding partially improper priors. Moreover, hyperpriors for the variances (corresponding to inverse smoothing parameters) may also be specified as improper, e.g. corresponding to Jeffreys prior or a flat prior for the standard deviation. Hence, propriety of the joint posterior is a crucial issue for full Bayesian inference in particular if based on Markov chain Monte Carlo simulations. We establish theoretical results providing sufficient (and sometimes necessary) conditions for propriety and provide empirical evidence through several accompanying simulation studies.     

Bayesian analysis of longitudinal ordered data with flexible random effects using McMC: application to diabetic macular Edema data

In the analysis of correlated ordered data, mixed-effect models are frequently used to control the subject heterogeneity effects. A common assumption in fitting these models is the normality of random effects. In many cases, this is unrealistic, making the estimation results unreliable. This paper considers several flexible models for random effects and investigates their properties in the model fitting. We adopt a proportional odds logistic regression model and incorporate the skewed version of the normal, Student's t and slash distributions for the effects. Stochastic representations for various flexible distributions are proposed afterwards based on the mixing strategy approach. This reduces the computational burden being performed by the McMC technique. Furthermore, this paper addresses the identifiability restrictions and suggests a procedure to handle this issue. We analyze a real data set taken from an ophthalmic clinical trial. Model selection is performed by suitable Bayesian model selection criteria.     

Bayesian comparison of private and common values in structural second-price auctions

Private and common values (CVs) are the two main competing valuation models in auction theory and empirical work. In the framework of second-price auctions, we compare the empirical performance of the independent private value (IPV) model to the CV model on a number of different dimensions, both on real data from eBay coin auctions and on simulated data. Both models fit the eBay data well with a slight edge for the CV model. However, the differences between the fit of the models seem to depend to some extent on the complexity of the models. According to log predictive score the IPV model predicts auction prices slightly better in most auctions, while the more robust CV model is much better at predicting auction prices in more unusual auctions. In terms of posterior odds, the CV model is clearly more supported by the eBay data.     

Complete convergence for weighted sums of END random variables and its application to nonparametric regression models

In this article, the complete convergence for weighted sums of extended negatively dependent (END, for short) random variables is investigated. Some sufficient conditions for the complete convergence are provided. In addition, the Marcinkiewicz–Zygmund type strong law of large numbers for weighted sums of END random variables is obtained. The results obtained in the article generalise and improve the corresponding one of Wang et al. [(2014b), ‘On Complete Convergence for an Extended Negatively Dependent Sequence’, Communications in Statistics-Theory and Methods, 43, 2923–2937]. As an application, the complete consistency for the estimator of nonparametric regression model is established.     

Bayesian semi-parametric analysis of Poisson change-point regression models: application to policy-making in Cali,Colombia

A Poisson regression model with an offset assumes a constant baseline rate after accounting for measured covariates, which may lead to biased estimates of coefficients in an inhomogeneous Poisson process. To correctly estimate the effect of time-dependent covariates, we propose a Poisson change-point regression model with an offset that allows a time-varying baseline rate. When the non-constant pattern of a log baseline rate is modeled with a non-parametric step function, the resulting semi-parametric model involves a model component of varying dimensions and thus requires a sophisticated varying-dimensional inference to obtain the correct estimates of model parameters of a fixed dimension. To fit the proposed varying-dimensional model, we devise a state-of-the-art Markov chain Monte Carlo-type algorithm based on partial collapse. The proposed model and methods are used to investigate the association between the daily homicide rates in Cali, Colombia, and the policies that restrict the hours during which the legal sale of alcoholic beverages is permitted. While simultaneously identifying the latent changes in the baseline homicide rate which correspond to the incidence of sociopolitical events, we explore the effect of policies governing the sale of alcohol on homicide rates and seek a policy that balances the economic and cultural dependencies on alcohol sales to the health of the public.     

Bivariate beta regression models: joint modeling of the mean,dispersion and association parameters

In this paper a bivariate beta regression model with joint modeling of the mean and dispersion parameters is proposed, defining the bivariate beta distribution from Farlie–Gumbel–Morgenstern (FGM) copulas. This model, that can be generalized using other copulas, is a good alternative to analyze non-independent pairs of proportions and can be fitted applying standard Markov chain Monte Carlo methods. Results of two applications of the proposed model in the analysis of structural and real data set are included.