Latent Gaussian Markov Random Fields Research Articles

Lung Cancer Prevalence in Virginia: A Spatial Zipcode-Level Analysis via INLA.

Examining lung cancer (LC) cases in Virginia (VA) is essential due to its significant public health implications. By studying demographic, environmental, and socioeconomic variables, this paper aims to provide insights into the underlying drivers of LC prevalence in the state adjusted for spatial associations at the zipcode level. We model the available VA zipcode-level LC counts via (spatial) Poisson and negative binomial regression models, taking into account missing covariate data, zipcode-level spatial association and allow for overdispersion. Under latent Gaussian Markov Random Field (GMRF) assumptions, our Bayesian hierarchical model powered by Integrated Nested Laplace Approximation (INLA) considers simultaneous (spatial) imputation of all missing covariates through elegant prediction. The spatial random effect across zip codes follows a Conditional Autoregressive (CAR) prior. Zip codes with elevated smoking indices demonstrated a corresponding increase in LC counts, underscoring the well-established connection between smoking and LC. Additionally, we observed a notable correlation between higher Social Deprivation Index (SDI) scores and increased LC counts, aligning with the prevalent pattern of heightened LC prevalence in regions characterized by lower income and education levels. On the demographic level, our findings indicated higher LC counts in zip codes with larger White and Black populations (with Whites having higher prevalence than Blacks), lower counts in zip codes with higher Hispanic populations (compared to non-Hispanics), and higher prevalence among women compared to men. Furthermore, zip codes with a larger population of elderly people (age ≥ 65 years) exhibited higher LC prevalence, consistent with established national patterns. This comprehensive analysis contributes to our understanding of the complex interplay of demographic and socioeconomic factors influencing LC disparities in VA at the zip code level, providing valuable information for targeted public health interventions and resource allocation. Implementation code is available at GitHub.

Fitting double hierarchical models with the integrated nested Laplace approximation

Double hierarchical generalized linear models (DHGLM) are a family of models that are flexible enough as to model hierarchically the mean and scale parameters. In a Bayesian framework, fitting highly parameterized hierarchical models is challenging when this problem is addressed using typical Markov chain Monte Carlo (MCMC) methods due to the potential high correlation between different parameters and effects in the model. The integrated nested Laplace approximation (INLA) could be considered instead to avoid dealing with these problems. However, DHGLM do not fit within the latent Gaussian Markov random field (GMRF) models that INLA can fit. In this paper, we show how to fit DHGLM with INLA by combining INLA and importance sampling (IS) algorithms. In particular, we will illustrate how to split DHGLM into submodels that can be fitted with INLA so that the remainder of the parameters are fit using adaptive multiple IS (AMIS) with the aid of the graphical representation of the hierarchical model. This is illustrated using a simulation study on three different types of models and two real data examples.

Reconstruction of past human land use from pollen data and anthropogenic land cover changes

AbstractAccurate maps of past land cover and human land use are necessary for studying the impact of anthropogenic land‐cover changes, such as deforestation, on the climate. The maps of past land cover should ideally be separated into naturally occurring vegetation and human‐induced changes, thereby enabling the quantification of the effect of human land use on the past climate. We developed a Bayesian hierarchical model that combines fossil pollen‐based reconstructions of actual land cover with estimates of past human land use. The model interpolates the fractions of unforested land as well as coniferous and broadleaved forest from the pollen data, and uses the human land‐use estimates to decompose the unforested land into natural vegetation and human deforestation. This results in maps of both natural and human‐induced vegetation, which can be used by climate modelers to quantify the influence of deforestation on the past climate. The model was applied to five time periods from 1900 CE to 4000 BCE over Europe. The model uses a latent Gaussian Markov random field (GMRF) for the interpolation and Markov chain Monte Carlo for the estimation. The sparse precision matrix of the GMRF, together with an adaptive Metropolis‐adjusted Langevin step, allows for rapid inference.

Reduced-Dimensional Monte Carlo Maximum Likelihood for Latent Gaussian Random Field Models

Monte Carlo maximum likelihood (MCML) provides an elegant approach to find maximum likelihood estimators (MLEs) for latent variable models. However, MCML algorithms are computationally expensive when the latent variables are high-dimensional and correlated, as is the case for latent Gaussian random field models. Latent Gaussian random field models are widely used, for example, in building flexible regression models and in the interpolation of spatially dependent data in many research areas such as analyzing count data in disease modeling and presence-absence satellite images of ice sheets. We propose a computationally efficient MCML algorithm by using a projection-based approach to reduce the dimensions of the random effects. We develop an iterative method for finding an effective importance function; this is generally a challenging problem and is crucial for the MCML algorithm to be computationally feasible. We find that our method is applicable to both continuous (latent Gaussian process) and discrete domain (latent Gaussian Markov random field) models. We illustrate the application of our methods to challenging simulated and real data examples for which maximum likelihood estimation would otherwise be very challenging. Furthermore, we study an often overlooked challenge in MCML approaches to latent variable models: practical issues in calculating standard errors of the resulting estimates, and assessing whether resulting confidence intervals provide nominal coverage. Our study therefore provides useful insights into the details of implementing MCML algorithms for high-dimensional latent variable models. Supplementary materials for this article are available online.

Bayesian Model Averaging with the Integrated Nested Laplace Approximation

The integrated nested Laplace approximation (INLA) for Bayesian inference is an efficient approach to estimate the posterior marginal distributions of the parameters and latent effects of Bayesian hierarchical models that can be expressed as latent Gaussian Markov random fields (GMRF). The representation as a GMRF allows the associated software R-INLA to estimate the posterior marginals in a fraction of the time as typical Markov chain Monte Carlo algorithms. INLA can be extended by means of Bayesian model averaging (BMA) to increase the number of models that it can fit to conditional latent GMRF. In this paper, we review the use of BMA with INLA and propose a new example on spatial econometrics models.

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.

Spatial Bayesian modeling of GLCM with application to malignant lesion characterization

ABSTRACTThe emerging field of cancer radiomics endeavors to characterize intrinsic patterns of tumor phenotypes and surrogate markers of response by transforming medical images into objects that yield quantifiable summary statistics to which regression and machine learning algorithms may be applied for statistical interrogation. Recent literature has identified clinicopathological association based on textural features deriving from gray-level co-occurrence matrices (GLCM) which facilitate evaluations of gray-level spatial dependence within a delineated region of interest. GLCM-derived features, however, tend to contribute highly redundant information. Moreover, when reporting selected feature sets, investigators often fail to adjust for multiplicities and commonly fail to convey the predictive power of their findings. This article presents a Bayesian probabilistic modeling framework for the GLCM as a multivariate object as well as describes its application within a cancer detection context based on computed tomography. The methodology, which circumvents processing steps and avoids evaluations of reductive and highly correlated feature sets, uses latent Gaussian Markov random field structure to characterize spatial dependencies among GLCM cells and facilitates classification via predictive probability. Correctly predicting the underlying pathology of 81% of the adrenal lesions in our case study, the proposed method outperformed current practices which achieved a maximum accuracy of only 59%. Simulations and theory are presented to further elucidate this comparison as well as ascertain the utility of applying multivariate Gaussian spatial processes to GLCM objects.

Stochastic Modelling of the Spatial Spread of Influenza in Germany

In geographical epidemiology, disease counts are typically available in discrete spatial units and at discrete time-points. For example, surveillance data on infectious diseases usually consists of weekly counts of new infections in pre-defined geographical areas. Similarly, but on a different timescale, cancer registries typically report yearly incidence or mortality counts in administrative regions.A major methodological challenge lies in building realistic models for spacetime interactions on discrete irregular spatial graphs. In this paper we will discuss an observation-driven approach, where past observed counts in neighboring areas enter directly as explanatory variables, in contrast to the parameterdriven approach through latent Gaussian Markov random fields (Rue and Held, 2005) with spatio-temporal structure. The main focus will lie on the demonstration of the spread of influenza in Germany, obtained through the design and simulation of a spatial extension of the classical SIR model (Hufnagel et al., 2004).

Bayesian motion estimation for dust aerosols

Dust storms in the earth’s major desert regions significantly influence microphysical weather processes, the CO$_{2}$-cycle and the global climate in general. Recent increases in the spatio-temporal resolution of remote sensing instruments have created new opportunities to understand these phenomena. However, the scale of the data collected and the inherent stochasticity of the underlying process pose significant challenges, requiring a careful combination of image processing and statistical techniques. Using satellite imagery data, we develop a statistical model of atmospheric transport that relies on a latent Gaussian Markov random field (GMRF) for inference. In doing so, we make a link between the optical flow method of Horn and Schunck and the formulation of the transport process as a latent field in a generalized linear model. We critically extend this framework to satisfy the integrated continuity equation, thereby incorporating a flow field with nonzero divergence, and show that such an approach dramatically improves performance while remaining computationally feasible. Effects such as air compressibility and satellite column projection hence become intrinsic parts of this model. We conclude with a study of the dynamics of dust storms formed over Saharan Africa and show that our methodology is able to accurately and coherently track storm movement, a critical problem in this field.

Classification of brain activation via spatial Bayesian variable selection in fMRI regression

Functional magnetic resonance imaging (fMRI) is the most popular technique in human brain mapping, with statistical parametric mapping (SPM) as a classical benchmark tool for detecting brain activity. Smith and Fahrmeir (J Am Stat Assoc 102(478):417---431, 2007) proposed a competing method based on a spatial Bayesian variable selection in voxelwise linear regressions, with an Ising prior for latent activation indicators. In this article, we alternatively link activation probabilities to two types of latent Gaussian Markov random fields (GMRFs) via a probit model. Statistical inference in resulting high-dimensional hierarchical models is based on Markov chain Monte Carlo approaches, providing posterior estimates of activation probabilities and enhancing formation of activation clusters. Three algorithms are proposed depending on GMRF type and update scheme. An application to an active acoustic oddball experiment and a simulation study show a substantial increase in sensitivity compared to existing fMRI activation detection methods like classical SPM and the Ising model.

Skew‐elliptical spatial random effect modeling for areal data with application to mapping health utilization rates

Mixed models incorporating spatially correlated random effects are often used for the analysis of areal data. In this setting, spatial smoothing is introduced at the second stage of a hierarchical framework, and this smoothing is often based on a latent Gaussian Markov random field. The Markov random field provides a computationally convenient framework for modeling spatial dependence; however, the Gaussian assumption underlying commonly used models can be overly restrictive in some applications. This can be a problem in the presence of outliers or discontinuities in the underlying spatial surface, and in such settings, models based on non-Gaussian spatial random effects are useful. Motivated by a study examining geographic variation in the treatment of acute coronary syndrome, we develop a robust model for smoothing small-area health service utilization rates. The model incorporates non-Gaussian spatial random effects, and we develop a formulation for skew-elliptical areal spatial models. We generalize the Gaussian conditional autoregressive model to the non-Gaussian case, allowing for asymmetric skew-elliptical marginal distributions having flexible tail behavior. The resulting new models are flexible, computationally manageable, and can be implemented in the standard Bayesian software WinBUGS. We demonstrate performance of the proposed methods and comparisons with other commonly used Gaussian and non-Gaussian spatial prior formulations through simulation and analysis in our motivating application, mapping rates of revascularization for patients diagnosed with acute coronary syndrome in Quebec, Canada.

Space‐Time Data fusion Under Error in Computer Model Output: An Application to Modeling Air Quality

We provide methods that can be used to obtain more accurate environmental exposure assessment. In particular, we propose two modeling approaches to combine monitoring data at point level with numerical model output at grid cell level, yielding improved prediction of ambient exposure at point level. Extending our earlier downscaler model (Berrocal, V. J., Gelfand, A. E., and Holland, D. M. (2010b). A spatio-temporal downscaler for outputs from numerical models. Journal of Agricultural, Biological and Environmental Statistics 15, 176-197), these new models are intended to address two potential concerns with the model output. One recognizes that there may be useful information in the outputs for grid cells that are neighbors of the one in which the location lies. The second acknowledges potential spatial misalignment between a station and its putatively associated grid cell. The first model is a Gaussian Markov random field smoothed downscaler that relates monitoring station data and computer model output via the introduction of a latent Gaussian Markov random field linked to both sources of data. The second model is a smoothed downscaler with spatially varying random weights defined through a latent Gaussian process and an exponential kernel function, that yields, at each site, a new variable on which the monitoring station data is regressed with a spatial linear model. We applied both methods to daily ozone concentration data for the Eastern US during the summer months of June, July and August 2001, obtaining, respectively, a 5% and a 15% predictive gain in overall predictive mean square error over our earlier downscaler model (Berrocal et al., 2010b). Perhaps more importantly, the predictive gain is greater at hold-out sites that are far from monitoring sites.

Bayesian multiscale feature detection of log-spectral densities

A fully-automatic Bayesian visualization tool to identify periodic components of evenly sampled stationary time series, is presented. The given method applies the multiscale ideas of the SiZer-methodology to the log-spectral density of a given series. The idea is to detect significant peaks in the true underlying curve viewed at different resolutions or scales. The results are displayed in significance maps, illustrating for which scales and for which frequencies, peaks in the log-spectral density are detected as significant. The inference involved in producing the significance maps is performed using the recently developed simplified Laplace approximation. This is a Bayesian deterministic approach used to get accurate estimates of posterior marginals for latent Gaussian Markov random fields at a low computational cost, avoiding the use of Markov chain Monte Carlo techniques. Application of the given exploratory tool is illustrated analyzing both synthetic and real time series.

Estimation of Binary Markov Random Fields Using Markov chain Monte Carlo

This article compares three binary Markov random fields (MRFs) which are popular Bayesian priors for spatial smoothing. These are the Ising prior and two priors based on latent Gaussian MRFs. Concern is given to the selection of a suitable Markov chain Monte Carlo (MCMC) sampling scheme for each prior. The properties of the three priors and sampling schemes are investigated in the context of three empirical examples. The first is a simulated dataset, the second involves a confocal fluorescence microscopy dataset, while the third is based on the analysis of functional magnetic resonance imaging (fMRI) data. In the case of the Ising prior, single site and multi-site Swendsen-Wang sampling schemes are both considered. The single site scheme is shown to work consistently well, while it is shown that the Swendsen-Wang algorithm can have convergence problems. The sampling schemes for the priors are extended to generate the smoothing parameters, so that estimation becomes fully automatic. Although this works well, it is found that for highly contiguous images fixing smoothing parameters to very high values can improve results by injecting additional prior information concerning the level of contiguity in the image. The relative properties of the three binary MRFs are investigated, and it is shown how the Ising prior in particular defines sharp edges and encourages clustering. In addition, one of the latent Gaussian MRF priors is shown to be unable to distinguish between higher levels of smoothing. In the context of the fMRI example we also undertake a simulation study.

Parallel exact sampling and evaluation of Gaussian Markov random fields

Markov chain Monte Carlo algorithms are computationally expensive for large models. Especially, the so-called one-block Metropolis–Hastings (M–H) algorithm demands large computational resources, and parallel computing seems appealing. A parallel one-block M–H algorithm for latent Gaussian Markov random field (GMRF) models is introduced. Important parts of this algorithm are parallel exact sampling and evaluation of GMRFs. Parallelisation is achieved with parallel algorithms from linear algebra for sparse symmetric positive definite matrices. The parallel GMRF sampler is tested for GMRFs on lattices and irregular graphs, and gives both good speed-up and good scalability. The parallel one-block M–H algorithm is used to make inference for a geostatistical GMRF model with a latent spatial field of 31,500 variables.