Medical, Pharma, Engineering, Science, Technology and Business

^{1}Department of Statistics and CIS, Zicklin School of Business, Baruch College, The City University of New York, New York, NY, USA

^{2}Department of Quantitative Health Sciences / Biostatistics Section, Cleveland Clinic Lerner Research Institute, Cleveland, OH, USA

- *Corresponding Author:
- Xiao-Feng Wang

Department of Quantitative Health Sciences

Cleveland Clinic Lerner Research Institute

9500 Euclid Avenue/JJN3, Cleveland, OH 44195, USA

**Tel:**216-445-7737

**Fax:**216-444-8023

**E-mail:**[email protected]

**Received date:** April 22, 2014; **Accepted date:** April 24, 2014; **Published date:** April 28, 2014

**Citation:** Yue YR, Wang XF (2014) Spatial Gaussian Markov Random Fields: Modelling, Applications and Efficient Computations. J Biomet Biostat 5:e128. doi:10.4172/2155-6180.1000e128

**Copyright:** © 2014 Yue YR, et al. This is an open-access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are are credited.

**Visit for more related articles at** Journal of Biometrics & Biostatistics

A powerful modelling tool for spatial data is the framework of Gaussian Markov random fields (GMRFs), which are discrete domain Gaussian random fields equipped with a Markov property. GMRFs allow us to combine the analytical results for the Gaussian distribution as well as Markov properties, thus allow for the development of computationally efficient algorithms. Here we briefly review popular spatial GMRFs, show how to construct them, and outline their recent developments and possible future work.

Gaussian Markov random fields; Markov chain Monte Carlo; spatial statistics; Cholesky factorization; Integrated nested laplace approximation

Gaussian Markov random fields (GMRFs) are powerful and important tools for modeling spatial data. They have been widely used in different areas of spatial statistics including disease mapping, spatialtemporal modeling and image analysis. Constructing a GMRF is straightforward: it is just a finite-dimensional random vector following a multivariate Gaussian distribution with additional conditional independence properties, hence termed as Markov. It is convenient and invaluable to combine the analytical results for the Gaussian distribution and the Markov properties, which enables us to solve a large class of statistical models. Historically, the most common method to make inference for the parameters in GMRFs has been maximum likelihood [1,2]. The behavior of maximum likelihood estimator is asymptotic in nature and their small sample behavior is often unknown. On the other hand, the Markov property has become a requirement for constructing efficient Markov chain Monte Carlo (MCMC) algorithms for GMRFs. Rue [3] showed that the Markov property makes it possible to apply numerical methods on sparse matrices. He proposed fast algorithms for sampling and evaluating the log-density of a GMRF, and conducted efficient MCMC-based inferences. Rue and Held [4] provides a comprehensive account of the main properties of GMRFs, emphasizes the strong connection between GMRFs and numerical methods for sparse matrices, and outlines various applications of GMRFs for statistical inference (e.g., spatial statistics, time-series analysis, graphical models).

More specifically, a GMRF is a Gaussian random vector x = (x_{1},… ,x_{n})’with Markov properties: for some i≠j, x_{i} and x_{j} are independent conditional on other variables x_{-ij}. It is defined over a set of discrete indexed locations connected by a graph labelled by G, which shows the conditional independence properties of x. We say x is a GMRF with respect to G with mean μ and precision (inverse covariance) matrix Q if its density has the form

(1)

The Markov properties in x are conveniently encoded in matrix Q: Q_{ij}= 0 if and only if x_{i} and x_{j} conditionally independent. The pattern of zero and non-zero elements in such a matrix is called its sparsity structure. The total number of non-zero elements divided by the total number of elements is called the density of the matrix. In most cases only O(n) of the n^{2} entries of Q are not zeroes, so Q has a very low density and is a sparse matrix. This allows for a fast Cholesky decomposition of Q as LL0, where L is the lower triangular matrix that inherits the sparseness of Q. Therefore, only non-zero terms in L are computed and the nodes can be reordered to decrease the number of the non-zero terms. The typical cost of this decomposition depends on the dimension of the GMRF: it is O(n) for one dimension, O(n^{3/2}) for two dimensions and O(n^{2}) for three dimensions. Using the Cholesky triangle, it is easy to produce random samples from a GMRF, and compute log-density of (1) and marginal variances; see Rue and Held [4] for more technical details.

GMRFs have been used in a wide range of common statistical models such as generalized linear (mixed) models, generalized additive models, dynamic linear models, and spatial and spatiotemporal models, among others. Those models can be written as hierarchical models which assume an-dimensional latent GMRF with a sparse precision matrix to be point-wise observed through n_{d} conditional independent data y. The estimation of the models takes advantage of modern techniques for sparse matrices. We here focus on the (two dimensional) spatial GMRFs that have been extensively used in the hierarchical analysis for spatial data [5]. We briefly review popular spatial GMRFs, show how to construct them, and outline their recent developments and possible future work.

From definition (1), we can see that a GMRF is completely determined by its precision matrix **Q:** different **Q**’s bring different Markov properties to the field. To build matrix **Q**, one has to specify a particular neighbourhood structure, corresponding to the type of the spatial data. Below we show how to construct a few GMRFs that have been widely used in spatial statistics.

**GMRFs on irregular lattices**

An important type of spatial data is so called areal data, where the locations of observations are geographic regions (e.g., the states of the US) with adjacency information. For areal spatial data, the firstorder GMRFs on irregular lattices are often preferred. To construct neighbourhood structure, we may define two regions are neighbors if they share a common border. Other ways to define neighbors are also possible. Between neighboring regions i and j, a Gaussian increment is defined as x_{i}-x_{j}~N(0,1/w_{ij}τ), where w_{ij} are positive and symmetric weights. We can let w_{ij}= 1 if we believe region i equally depends on its neighbors, or let w_{ij} be, for example, the inverse Euclidean distance between region centroids if we think the neighbors somehow contribute differently. Assuming the increments are independent, the density of this GMRF is given by:

(2)

Where i ~ j denotes the set of all unordered pairs of neighbors. The corresponding precision matrix **Q** has entries

(3)

Where the summation over neighbors of node i on a lattice. Since the sum of each row is zero, **Q** is singular with rank n-1. It is easy to see that the conditional distribution of x_{i} is

(4)

where the conditional mean of xi depends on its neighboring nodes x j through weights wi j and its conditional variance depends on weight sum w_{i+}.

**GMRFs on regular lattices**

Another important type of spatial data are point-referenced data, where the spatial locations are points with known coordinates. When the coordinates constitute a regular lattice, one can construct first and second-order GMRFs as described below.

For a regular lattice with n=n_{1}n_{2} nodes, let (i,j) denote the node in the ith row and jth column. In the interior, we can define the nearest four nodes of (i,j) as neighbors, i.e., (i+1,j), (i-1,j), (i,j+1), (i,j-1). Along the boundaries, we define the neighbors of (i,j) to be the two or three adjacent nodes. For example, the neighbors of (i,j) are (i,j+1) and (i+1, j) if (i,j) is the upper left corner. Without further weights, the corresponding precision matrix and the full conditionals of x_{i} are given in (3) and (4) with w_{ij}=1, respectively. By weighting the horizontal and vertical neighbors differently, this GMRF can be extended to an anisotropic model [4]. In practice the first-order GMRF models may not provide sufficient spatial smoothing as we need. To increase the smoothness of the field, we can use higher-order neighborhood structures. One way to build a second-order neighborhood structure for xijin the interior is based on the increment

(x_{i+1,j}+x_{i-1,j}+x_{i,j+1}+x_{i,j-1})-4x_{i,j}

which is the sum of second-order differences in vertical and horizontal directions. For the nodes on or near the boundaries we need different increments; see Yue and Speckman [6] for details. By assuming all the increments follow independent Gaussian distributions, we build a second-order GMRF on a regular lattice. The co-efficients of the corresponding **Q** can be found by expanding the quadratic terms of those increments. Using graphical notation, the conditional distribution of xij in the interior has mean and precision

(5)

where the locations denoted by a ‘ • ’ represent the neighbors that xi j depends on and the number in front of each grid denotes the weight given to the corresponding ‘ ○ ’ locations. This second-order GMRF is closely related to thin-plate spline [7], which has been widely used as a spatial smoother [1,8]. Yue and Speckman [6] showed that this GMRF can also be derived by discretizing the penalty function of the thinplate spline.

**Continuous indexed GMRFs**

The GMRFs mentioned above are all discrete indexed in nature and thus only work appropriately for the spatial data on lattices. If the observations are measured at irregularly-spaced locations, one often bins the locations to a regular lattice first and then apply the GMRF to the summary statistics calculated for the grids. As a result, we lose information from binning process and the the spatial resolution of the lattice significantly affects the inference. Another problem is that a regular GMRF can only capture the smoothness of a spatial field, but not its structure of spatial correlations, because it defines the precision matrix (not covariance matrix). Thus, we cannot make inference regarding the correlation between two locations.

To address the issues mentioned above, Lindgren et al. [9] derived a new class of continuous indexed GMRFs, which are explicit mappings of Mat´ern Gaussian fields that have been extensively used in statistical modeling of spatial data. Letting x(u) be a realization of spatial field x at location , such GMRFs are derived by solving a stochastic partial differential equation (SPDE)

(6)

where is the two-dimensional Laplacian operator, κ>0 is the spatial scale, α> 0 controls the smoothness of the realizations, τ>0 controls the variance, and wis the spatial Gaussian white noise. The link to the Mat´ern smoothness n and variance σ^{2} is ν= α-d/2 and σ^{2}=Γ(ν)(Γ(α)(4π)^{d/2}κ^{2ν}τ^{2})^{-1}, where d is the spatial dimension. A measure of the spatial range can be empirically derived as Lindgren et al. [9] solve the SPDE (6) using finite element method.They first approximate x(u) with piecewise linear basis functions defined on a triangular domain, and then turn (6) into a system of linear equations to derive the GMRF solution.

As shown in Lindgren et al. [9], the derived GMRF is the best piecewise linear approximation to the continuous solution to the SPDE given a triangulation. Since it is a GMRF representation of Mat´ern fields, it allows us to capture both spatial correlation and spatial smoothness in a spatial process. Another exciting aspect of such models is their flexibility. There is no conceptual or computational barrier to extending them to being, for example, non-stationary, multivariate and spatial-temporal GMRFs. It is even possible to construct them on the sphere and other manifolds.

**Adaptive GMRFs**

The smoothness of a GMRF is determined by the scale of increments, which is often invariant across the space. It means that the GMRF provides the same amount of smoothing at every location. This, however, limits the GMRF to stationary spatial processes only. We here present a few adaptive GMRFs that are recently developed to deal with non-stationary spatial processes.

Brezger et al. [10] extended the first-order GMRF in (2) to being adaptive by letting weights wi j vary with locations. More specifically, they take independent gamma priors on each w_{ij}, that is w_{ij}~ Gamma (1/2,1/2), making the marginal distribution of difference (x_{i}-x_{j}) a Cauchy distribution. This approach allows varying strength of interactions between neighboring sites i ~ j, but without spatial prior structure.

Similarly, Yue and Speckman [6] extended the second-order GMRF in (5) to being non-stationary by spatially adaptive modeling its conditional variance, i.e., Var (x_{ij}|x_{-ij}) = (20τ_{ij})^{-1}. Thus, a small value of τ_{ij} (large variance) leads to less smoothing, appropriate when x_{ij} shows increased local variation. Then, let τ_{ij} = τγ_{ij}, so that τ is the global scale parameter and γ_{ij} controls the local smoothing of the field. Finally, a first-order spatial GMRF is taken as a prior on log(γ_{ij}) to make the smoothing vary over the space. The resulting GMRF is able to capture both local and global features of a spatial process while retains the nice Markov properties for computation. Yue et al. [11] proposed a similar adaptive GMRF but with independent gamma priors on γ_{ij}.

The SPDE models can be generalized to be non-stationary in several ways Lindgren et al. [9]. Let spatial scale κ and precision τ in (6) to depend on the coordinate u, and allow them vary slowly through log linear models. Therefore, one can analyze how the spatial nonstationarity depends on certain covariates. Bolin and Lindgren [9] used a nested SPDE approach to develop a large class of non-stationary covariance functions, such as oscillating covariance functions. Fuglstad et al. [12] introduced a new class of non-stationary spatial GMRFs with varying local anisotropy by adding a diffusion matrix that varies with position to the Δ operator in (6).

The adaptive GMRF models mentioned above have been successfully applied to the function magnetic resonance imaging (fMRI) data [10,13], the precipitation data [6,14], the global ozone data [15] and the neuroimaging meta-analysis [11].

**Other extensions**

Spatial-temporal GMRF models are extensions of spatial GMRF models to account for additional temporal variation. Think about a sequence of T graphs in time and let x_{it} denote the ith node in tth graph. A common extension to a spatial-temporal GMRF is to take into account temporal neighbors in addition to spatial neighbors. The temporal neighbors can be the same nodes in the previous and next graphs, that is, x_{it-1} and x_{it+1}. The precision matrix of the corresponding GMRF can be written as where Q_{T} and Q_{S} are precision matrices in time and space, respectively. We can use conditional autoregressive process (Besag, 1974) [16] to model temporal correlations and the spatial GMRFs presented above to model spatial correlations. Then, both Q_{T} and Q_{S} are sparse matrices, making Q sparse as well. Cameletti et al. [17] employed a model of this type to perform spatial-temporal analysis on particulate matter concentration data.

It is possible to build a multivariate GMRF by extending the SPDE framework considered above [18]. The idea is to replace a single SPDE with a system of SPDEs. The co-variances matrices constructed with this approach are automatically symmetric positive definite. The sparse precision matrix of this multivariate GMRF facilitates its implementation to large data sets. Hu et al. [19] have shown that these models are related (but not equivalent) to the multivariate Mat´ern fields constructed by.

The GMRFs are often implemented in a class of structured additive regression models, which are quite flexible and have been extensively used [20]. In these models, the response variable y_{i} is assumed to belong to an exponential family, where the mean μ_{i} is linked to a structured additive predictor η_{i} through a link function g(•), so that g(μ_{i}) = η_{i}. The predictor hi accounts for various covariate effects in an additive way:

(7)

Where fj’s are unknown functions of the covariates μ_{i}, β_{k}’s represent the linear effect of z_{k}, and the ε_{i}’s are error terms. The structured additive regression models cover a large class of spatial and spatiotemporal models, geostatistical and geoadditive models. The spatial GMRFs can be used to model f_{j} when the covariate u_{j} are spatial locations. The model (7) can handle various types of responses (e.g., continuous, binary and count data) and very different forms that the unknown f_{j} can take (e.g., random, nonlinear and spatial effects). As a matter of fact, most Bayesian models used in spatial statistics are of this form [1].

To fit the spatial models of form (7), one can perform standard MCMC methods to simulate the posterior distributions. Using GMRFs as priors in the Bayesian hierarchical models, it is feasible to construct accurate GMRF approximations to the full conditionals (if they are not Gaussian), based on which efficient blockwise Metropolis-Hastings algorithms can be used to explore the posteriors [4]. One is able to take advantage of the sparseness of GMRFs and implement the fast Cholesky factorization algorithms to speed up the MCMC algorithms. This can now be easily done using an R package named spam developed by Furrer and Sain [21]. Recently, Rue et al. [22] has introduced a novel Bayesian inference tool based on integrated nested Laplace approximations (INLA). The INLA method can directly compute very accurate approximations to the posterior marginal distributions for latent Gaussian models. It is much faster than MCMC and can be easily implemented using R-INLA package (http://www.r-inla.org). The package contains a variety of popular GMRFs that are readily to be used.

- CressieN (1993) Statistics for Spatial Data. New York: Wiley-Interscience.
- Richardson S,Guihenneuc C, Lasserre V (1992)Spatial linear models with autocorrelated errorstructure. The Statistician 41: 539-557.
- RueH (2001) Fast sampling of Gaussian Markov random fields. Journal of the Royal StatisticalSociety Series B 65: 325–338.
- Rue H, Held L (2005) Gaussian Markov Random Fields: Theory and Applications, vol. 104 of Monographson Statistics and Applied Probability. London: Chapman & Hall.
- BanerjeeS, Gelfand AE, Carlin BP (2004) Hierarchical modeling and analysis for spatialdata. CRC Press.
- YueY, Speckman PL(2010). Nonstationary spatial Gaussian Markov random fields.Journal of Computational and Graphical Statistics 19: 96–116.
- WahbaG (1990) Spline Models for Observational Data. Philadelphia: SIAM [Society forIndustrial and Applied Mathematics].
- NychkaDW(2000) Spatial-process estimates as smoothers. In Smoothing and regression:approaches, computation, and application. MGaSchimek, ed. John Wiley &Sons.
- LindgrenF, Rue H, Lindstr¨Om J (2011)An explicit link between gaussian fields and Gaussianmarkov random fields: the stochastic partial differential equation approach(with discussion). Journal of the Royal Statistical Society: Series B(Statistical Methodology) 73: 423–498.
- BrezgerA, Fahrmeir L, Hennerfeind A (2007) Adaptive Gaussian Markov random fields withapplications in human brain mapping. Journal of the Royal Statistical Society: SeriesC (Applied Statistics) 56: 327–345.
- Yue YR,Lindquist MA, Loh JM (2012) Meta-analysis of functional neuroimagingdatausingbayesian nonparametric binary regression. Annals of Applied Statistics6: 697–718.
- Fuglstad GA, Lindgren F, Simpso D,Rue H (2014). Exploring a new class of nonstationary spatial Gaussian randomfields with varying local anisotropy. StatisticaSinica.
- YueY, Loh JM, Lindquist MA (2010). Adaptive spatial smoothing of fMRI images.Statistics and Its Interface 3: 3–13.
- Fuglstad GA, Lindgren F, Simpson D,Rue H (2013) Non-stationary spatial modelling with applications to spatialprediction of precipitation. Tech. rep., Norwegian University of Science andTechnology.
- Bolin D, Lindgren F(2011).Spatial models generated by nested stochastic partial differential equations,with an application to global ozone mapping. The Annals of Applied Statistics5: 523–550.
- BesagJ(1974). Spatial interaction and the statistical analysis of lattice systems(with discussion). Journal of the Royal Statistical Society, Series B:Methodological 36: 192–236.
- CamelettiM, Lindgren F, Simpson D, Rue H(2013)Spatio-temporal modeling of particulatematter concentration through the SPDE approach. AStA Advances in StatisticalAnalysis 97: 109–131.
- Hu X, Simpson D, Lindgren F, Rue H(2013)Multivariate Gaussian Random Fields Using Systems of Stochastic PartialDifferential Equations. Tech. rep., Norwegian University of Science andTechnology.
- GneitingT, Kleiber W, Schlather M(2010)Mat´ern cross-covariance functions formultivariate random fields. Journal of the American Statistical Association105: 1167 – 1177.
- FahrmeirL, Tutz G(2001) Multivariate Statistical Modeling based on Generalized LinearModels. Berlin: Springer.
- Furrer R, Sain SR(2010) spam: A sparse matrix r package with emphasis on mcmc methods forgaussianmarkov random fields. Journal of Statistical Software 36: 1–25.
- RueH, Martino S, Chopin N (2009) Approximate Bayesian inference for latentGaussian models using integrated nested Laplace approximations (withdiscussion). Journal of the Royal Statistical Society, Series B: StatisticalMethodology 71: 319–392.

Select your language of interest to view the total content in your interested language

- Adomian Decomposition Method
- Algebra
- Algebraic Geometry
- Algorithm
- Analytical Geometry
- Applied Mathematics
- Artificial Intelligence Studies
- Axioms
- Balance Law
- Behaviometrics
- Big Data Analytics
- Big data
- Binary and Non-normal Continuous Data
- Binomial Regression
- Bioinformatics Modeling
- Biometrics
- Biostatistics methods
- Biostatistics: Current Trends
- Clinical Trail
- Cloud Computation
- Combinatorics
- Complex Analysis
- Computational Model
- Computational Sciences
- Computer Science
- Computer-aided design (CAD)
- Convection Diffusion Equations
- Cross-Covariance and Cross-Correlation
- Data Mining Current Research
- Deformations Theory
- Differential Equations
- Differential Transform Method
- Findings on Machine Learning
- Fourier Analysis
- Fuzzy Boundary Value
- Fuzzy Environments
- Fuzzy Quasi-Metric Space
- Genetic Linkage
- Geometry
- Hamilton Mechanics
- Harmonic Analysis
- Homological Algebra
- Homotopical Algebra
- Hypothesis Testing
- Integrated Analysis
- Integration
- Large-scale Survey Data
- Latin Squares
- Lie Algebra
- Lie Superalgebra
- Lie Theory
- Lie Triple Systems
- Loop Algebra
- Mathematical Modeling
- Matrix
- Microarray Studies
- Mixed Initial-boundary Value
- Molecular Modelling
- Multivariate-Normal Model
- Neural Network
- Noether's theorem
- Non rigid Image Registration
- Nonlinear Differential Equations
- Number Theory
- Numerical Solutions
- Operad Theory
- Physical Mathematics
- Quantum Group
- Quantum Mechanics
- Quantum electrodynamics
- Quasi-Group
- Quasilinear Hyperbolic Systems
- Regressions
- Relativity
- Representation theory
- Riemannian Geometry
- Robotics Research
- Robust Method
- Semi Analytical-Solution
- Sensitivity Analysis
- Smooth Complexities
- Soft Computing
- Soft biometrics
- Spatial Gaussian Markov Random Fields
- Statistical Methods
- Studies on Computational Biology
- Super Algebras
- Symmetric Spaces
- Systems Biology
- Theoretical Physics
- Theory of Mathematical Modeling
- Three Dimensional Steady State
- Topologies
- Topology
- mirror symmetry
- vector bundle

- Total views:
**12276** - [From(publication date):

April-2014 - Jan 22, 2018] - Breakdown by view type
- HTML page views :
**8415** - PDF downloads :
**3861**

Peer Reviewed Journals

International Conferences
2018-19