Bayesian Analysis of Spatial Generalized Linear Mixed Models with Laplace Random Fields

TL;DR

This paper explores replacing Gaussian random fields with Laplace moving averages in spatial generalized linear mixed models, showing improved prediction for data with localized spikes while maintaining similar inference and computational efficiency.

Contribution

It introduces a novel discrete space LMA model for irregular lattices and develops conjugate samplers, expanding Bayesian spatial modeling capabilities.

Findings

LMAs improve predictive accuracy for spike-like spatial data

Parameter inference remains consistent with Gaussian models

Computational methods are comparable to existing Gaussian SGLMMs

Abstract

Gaussian random field (GRF) models are widely used in spatial statistics to capture spatially correlated error. We investigate the results of replacing Gaussian processes with Laplace moving averages (LMAs) in spatial generalized linear mixed models (SGLMMs). We demonstrate that LMAs offer improved predictive power when the data exhibits localized spikes in the response. SGLMMs with LMAs are shown to maintain analogous parameter inference and similar computing to Gaussian SGLMMs. We propose a novel discrete space LMA model for irregular lattices and construct conjugate samplers for LMAs with georeferenced and areal support. We provide a Bayesian analysis of SGLMMs with LMAs and GRFs over multiple data support and response types.

Figures

Tables

Table 1. Table 1: BCVS and ESS for ten-fold cross validation on the Slovenia stomach cancer outbreak dataset.

Model	BCVS	ESS
GRF	560.61	4.12
LMA	$563.04$	3.81

Table 2. Table 2: BCVS from ten-fold cross validation and ESS for Columbus Ohio Crime dataset.

Model	BCVS	ESS
GRF	$341.58$	18.13
LRF	$329.85$	14.07

Table 3. Table 3: Parameter estimates for discrete space data analysis examples of Sections ( 5.1 – 5.2 )

Slovenia Stomach Outbreak
Predictor	Parameter	GRF Estimate	95% CI	LMA Estimate	95% CI
Intercept	${\beta }_0$	$0.097$	$(-0.036,0.202)$	$0.115$	$(-0.010,0.339)$
SEc	${\beta }_1$	$-0.078$	$(-0.170,0.031)$	$-0.068$	$(-0.162,0.040)$
	$\xi$	$0.898$	$(0.492,1.561)$	NA	NA
	$\sigma$	$0.290$	$(0.238,0.349)$	$0.282$	$(0.230,0.339)$
	$\kappa$	1.427	(0.616, 2.523)	1.608	(0.814, 2.543)
	$\lambda$	NA	NA	1.127	(0.718, 1.701)
Columbus Crime Data
Predictor	Parameter	GRF Estimate	95% CI	LMA Estimate	95% CI
Intercept	${\beta }_0$	35.107	(-27.661, 76.053)	29.379	(-39.597, 83.936)
Avg. Inc	${\beta }_1$	-0.321	(-0.391, -0.252)	-0.235	(-0.327, -0.149)
Avg. Value	${\beta }_2$	-0.981	(-1.247, -0.716)	-1.057	(-1.351, -0.755)
	$\xi$	4.125	(3.536, 4.676)	NA	NA
	$\sigma$	3.370	(3.191, 3.550)	2.803	(2.541, 3.050)
	$\kappa$	0.162	(0.051, 0.298)	0.219	(0.081, 0.373)
	$\lambda$	NA	NA	4.582	(4.086, 5.040)

Table 4. Table 4: BCVS from ten-fold cross validation and ESS for malaria in the Gambia, Africa.

Model	BCVS	ESS
GRF	1292.20	2.66
LMA	1284.08	2.90

Table 5. Table 5: BCVS for ten fold cross validation and ESS for the LAGOS dataset.

Model	BCVS	ESS
GRF	5733.309	0.35
LMA	5727.993	0.39

Table 6. Table 6: Parameter estimates for continuous space data analysis examples of Sections ( 5.3 – 5.4 )

Malaria of the Gambia
Predictor	Parameter	GRF Estimate	95% CI	LMA Estimate	95% CI
Intercept	${\beta }_0$	0.123	(-2.408, 2.677)	-0.975	(-4.563, 2.554)
Age	${\beta }_1$	0.163	(0.111, 0.215)	0.178	(0.124, 0.232)
Bed Net Use	${\beta }_2$	-0.262	(-0.413, -0.109)	-0.218	(-0.398, -0.036)
Insecticide Use	${\beta }_3$	-0.104	(-0.270, 0.062)	-0.211	(-0.417, -0.007)
$\log \left(\text{NDVI}\right)$	${\beta }_4$	-0.007	(-0.659, 0.637)	0.283	(-0.633, 1.220)
Health	${\beta }_5$	-0.253	(-0.397, -0.109)	-0.089	(-0.297, 0.119)
	$\xi$	2.759	(1.612, 4.203)	NA	NA
	$\kappa$	1.413	(0.399, 2.433)	1.192	(0.518, 1.759)
	$\tau$	NA	NA	13.272	(12.074, 14.540)
	$\lambda$	NA	NA	4.730	(3.506, 5.824)
LAGOS
Predictor	Parameter	GRF Estimate	95% CI	LMA Estimate	95% CI
Intercept	${\beta }_0$	3.93	(3.773,4.087)	3.925	(3.771,4.081)
Area	${\beta }_1$	6.33e-05	(5.00e-05,7.71e-05)	6.33e-05	(4.99e-05,7.69e-05)
Max depth	${\beta }_2$	-0.036	(-0.038,-0.034)	-0.036	(-0.038,-0.034)
IWS Urban	${\beta }_3$	1.296	(1.103,1.489)	1.302	(1.105,1.498)
IWS Ag	${\beta }_4$	1.594	(1.485,1.701)	1.598	(1.491,1.707)
IWS Wetland	${\beta }_5$	0.635	(0.448,0.822)	0.639	(0.454,0.826)
Road Density	${\beta }_6$	-0.002	(-0.003,-0.001)	-0.002	(-0.003,-0.001)
Baseflow	${\beta }_7$	-0.012	(-0.014,-0.01)	-0.012	(-0.014,-0.01)
Runoff	${\beta }_8$	-0.026	(-0.03,-0.021)	-0.026	(-0.03,-0.021)
	$\xi$	0.383	( 0.220 , 0.580 )	NA	NA
	${\sigma }^2$	0.464	( 0.444 , 0.485 )	0.464	( 0.444 , 0.485 )
	$\kappa$	0.772	( 0.397 , 1.178 )	0.724	( 0.477 , 0.988 )
	$\tau$	NA	NA	1.751	( 1.322 , 2.263 )
	$\lambda$	NA	NA	0.266	( 0.184 , 0.38 )

Equations

$$\displaystyle g(\mu_{i})=\boldsymbol{x}_{i}^{\prime}\boldsymbol{\beta}+\eta(\boldsymbol{u}_{i})+\mbox{$\epsilon$}(\boldsymbol{u}_{i}).$$

$$({\kappa }^2-△{)}^{\alpha /2}\eta (\mathbf{u})=\xi \mathcal{W}(\mathbf{u}),\mathbf{u}\in {\mathbb{R}}^d,\alpha =\nu +d/2$$

$$C(\mathbf{u},\mathbf{v})=\frac{{\varphi }^2}{2^{\nu -1}\Gamma (\nu )}(\kappa ||\mathbf{v}-\mathbf{u}||{)}^{\nu }{K}_{\nu }(\kappa ||\mathbf{v}-\mathbf{u}||),\mathbf{u},\mathbf{v}\in {\mathbb{R}}^d$$

$$\dot{\Lambda }=\sum _{k=1}^{\infty }\left({\Gamma }_k+\sqrt{{\Gamma }_k}{G}_k\right){\delta }_{{\mathbf{u}}_k}$$

$$({\kappa }^2-△{)}^{\alpha /2}\eta (\mathbf{u})=\dot{\Lambda }(\mathbf{u}),\alpha =\nu +d/2.$$

$$({\kappa }^2-△{)}^{\alpha /2}\eta (\psi )=\xi \mathcal{W}(\psi ),\alpha =\nu +d/2$$

$$\eta (\mathbf{u})=\sum _{i=1}^n{\varphi }_i(\mathbf{u})w_i,\mathbf{u}\in \Omega$$

$$\mathbf{w}|\xi ,\kappa \sim \mathcal{N}\left(0,{\xi }^2{\mathbf{Q}}_{\alpha }^{-2}\right)$$

$${\mathbf{Q}}_{\alpha }=\{\begin{array}{ll}{\mathbf{Q}}_1=\mathbf{L},& \alpha =1\\ {\mathbf{Q}}_2=\mathbf{L}{\mathbf{C}}^{-1}\mathbf{L},& \alpha =2\\ {\mathbf{Q}}_{\alpha }=\mathbf{L}{\mathbf{C}}^{-1}{\mathbf{Q}}_{(\alpha -2)}{\mathbf{C}}^{-1}\mathbf{L},& \alpha \ge 3\end{array}$$

$$C_{ij}$$

$$G_{ij}$$

$$C_{ii}={\int }_{\Omega }{\varphi }_i(\mathbf{u})d\mathbf{u}$$

$$({\kappa }^2-△{)}^{\alpha /2}\eta (\psi )=\lambda \dot{\Lambda }(\psi )$$

$${\mathbf{K}}_{\alpha }\mathbf{w}|\mathbf{\Gamma }$$

$${\mathbf{K}}_{\alpha }=\{\begin{array}{ll}{\mathbf{K}}_2=\mathbf{L},& \alpha =2\\ {\mathbf{K}}_{\alpha }=\mathbf{L}{\mathbf{C}}^{-1}{\mathbf{K}}_{\alpha -2},& \alpha =4,6,8,...\end{array}$$

$$\left[\mathbf{w}|\mathbf{y},{\sigma }^2,\mathbf{\beta },\mathbf{\Gamma }\right]\sim \mathcal{N}\left({\left[{\sigma }^2\mathbf{L}{\mathbf{\Gamma }}^{-1}\mathbf{L}+{\mathbf{A}}^{\prime }\mathbf{A}\right]}^{-1}\left[\mathbf{y}-\mathbf{A}\mathbf{w}\right],{\left[\mathbf{L}{\mathbf{\Gamma }}^{-1}\mathbf{L}+\left(\frac{1}{{\sigma }^2}\right){\mathbf{A}}^{\prime }\mathbf{A}\right]}^{-1}\right).$$

$$f(x;p,a,b)\propto x^{p-1}\exp \left(-\frac{ax+b/x}{2}\right).$$

$$\hat{\mathbf{\beta }}=\underset{\mathbf{\beta }\in {\mathbb{R}}^n}{\mathrm{argmin}}\frac{1}{2}||\mathbf{y}-\mathbf{\beta }|{|}_2^2+\lambda ||{△}^{(k+1)}\mathbf{\beta }|{|}_1$$

$${△}^{(k+1)}=\{\begin{array}{ll}{\left({△}^{(1)}\right)}^{\prime }{△}^{(k)},& \text{for odd }k\\ {△}^{(1)}{△}^{(k)},& \text{for even }k\end{array}.$$

$$A_{ij}=\{\begin{array}{ll}-1,& \text{i is neighboring j}\\ {\sum }_k|A_{ik}|,& i=j\\ 0,& \text{else}\end{array}.$$

$${\mathbf{Q}}_k={\left({△}^{(k)}\right)}^{\prime }{△}^{(k)}$$

$${△}^{(k)}=\{\begin{array}{ll}{\mathbf{L}}^{\frac{(k+1)}{2}},& \text{for odd k}\\ \mathbf{D}{\mathbf{L}}^{\frac{k}{2}},& \text{for even k }\end{array}$$

$${△}^{(k)}\mathbf{\eta }|{\xi }^2,{\kappa }^2\sim \mathcal{N}\left(0,{\xi }^2\mathbf{I}\right)$$

$$Cov\left(\mathbf{\eta }|{\xi }^2,{\kappa }^2\right)={\xi }^2{\mathbf{L}}^{-(k+1)}={\xi }^2{\mathbf{Q}}_k^{-1}$$

$${{△}_l^{\prime }}^{(k)}\mathbf{\eta }|\lambda ,{\kappa }^2\stackrel{iid}{\sim }\mathcal{L}(\lambda ),l=1,...,n$$

$$Cov\left(\mathbf{\eta }|\lambda ,{\kappa }^2\right)=\frac{2}{\lambda }{\mathbf{L}}^{-(k+1)}=\frac{2}{\lambda }{\mathbf{Q}}_k^{-1}$$

$$\log [\mathbf{\eta }|\mathbf{\beta },\mathbf{y},{\mathbf{\theta }}_G]\approx \log [\mathbf{y}|\mathbf{\beta },\mathbf{\eta },{\mathbf{\theta }}_G]-\frac{1}{2{\xi }^2}\sum _{i=1}^n({\kappa }^2{\eta }_i+\sum _{j\in N(i)}({\eta }_i-{\eta }_j){)}^2+Const.$$

$$\log [\mathbf{\eta }|\mathbf{\beta },\mathbf{y},{\mathbf{\theta }}_L]\approx \log [\mathbf{y}|\mathbf{\beta },\mathbf{\eta },{\mathbf{\theta }}_L]-\frac{\lambda }{2}\sum _{i=1}^n|{\kappa }^2{\eta }_i+\sum _{j\in N(i)}({\eta }_i-{\eta }_j)|+Const.$$

$${△}^{(k)}\mathbf{\eta }|\mathbf{S},{\kappa }^2\sim N(0,\mathbf{S}),\mathbf{S}=diag(S_1,..,S_n)$$

$$f(x;a,b)={\left(\frac{b}{2\pi x^3}\right)}^{1/2}\exp \left(-\frac{b(x-a{)}^2}{2b^{2}x}\right)$$

$$BCVS=-\sum _{k=1}^{10}\log \left(\frac{{\sum }_{t=1}^T[{\mathbf{y}}_k|{\mathbf{y}}_{-k},{\mathbf{\theta }}^{(t)}]}{T}\right)$$

Full text

††journal: Journal of Computational Statistics and Data Analysis

1 Introduction

Gaussian processes are the most common method for describing spatially and/or temporally correlated errors. The Gaussian random field (GRF) possesses an intuitive dependence structure, as well as a globally flexible fit. These desirable features have popularized the use of Gaussian processes in spatial and temporal statistics, as well as design of experiments and other fields. Despite the GRF’s flexible nature, Gaussian processes can over-smooth in the presence of local spikes [17]. In this work we consider the use of Laplace moving average models (LMAs) in place of traditional Gaussian processes in spatial generalized linear mixed models (SGLMMs)

LMAs have received sporadic attention over the past decade as alternatives to GRFs [1, 27, 16]. However, there has been no systematic comparison of LMAs and GRFs for spatial generalized linear mixed models (SGLMMs) in the Bayesian framework. Our contributions in this work include

The development of a novel discrete space (areal) LMA model for irregular lattices. 2. 2.

The construction of conjugate samplers for both continuous (point-referenced) and discrete (areal) SGLMMs with LMAs. 3. 3.

A Bayesian analysis comparing the predictive power and computational efficiency of LMAs and GRFs over a range of scenarios, including continuous, binary, and count data collected both in discrete (areal) and point-referenced (geostatistical) spatial support.

Whittle [29] demonstrated that continuous space GRFs with Matérn covariance arise as solutions to a stochastic partial differential equation (SPDE). Lindgren et al. [15] constructed a sparse finite element representation of this Gaussian Matérn SPDE. As a result, a sparse form of the multivariate normal distribution can be used to fit Matérn GRF models in a computationally efficient manner [21]. Bolin [6] extended the finite element approximation of Lindgren et al. [15] to the case of Type-G Matérn random fields of which the LMA is a special case. In Section 3, we provide a summary of this extension for the symmetric LMA. Following Bolin [6], the LMA can similarly be expressed as a conditionally sparse Gaussian random field through the introduction of auxiliary data. We also provide insights for handling the computational issues associated with MCMC implementation.

Wallin and Bolin [27] explored the LMA for geostatistical data with Gaussian responses. Though the discrete space model was claimed to be analogous, no further exploration was considered. Faulkner and Minin [9] provided a Bayesian implementation of the graph trend filtering (GTF) estimates of Wang et al. [28] for temporal data. Both works, found that replacing traditional Gaussian priors for Laplace priors provided a model with better adaptivity in the presence of local “spikes” in the response. We develop a novel discrete space analog to the continuous space LMA model that is an extension of Wang et al. [28] for SGLMMs. Our model can easily be implemented in place of any Gaussian conditionally autoregressive (CAR) or simultaneously autoregressive (SAR) model. In Section 4, we propose a novel MCMC implementation of our Bayesian hierarchical model for discrete space SGLMMs.

The LMA models of this paper offer an intuitive alternative to traditional GRF SGLMMs. The discrete space and continuous space LMA models are constructed based on sparse matrix operations making for fast and efficient fitting. We provide Bayesian analyses based on our novel MCMC implementation of the LMA models over four separate data sets. Our MCMC implementation and model construction allows for Bayesian inference with LMAs that is just as interpretable as with GRF models. In some cases, the LMA is shown to provide a better fit than the Gaussian model. Given the ease of implementation, and familiarity of inference, LMA models can be useful alternatives to GRF models.

The paper is organized as follows: In Section 2 we provide background material needed to develop our models. In Section 3 we discuss finite element approximations for continuous space LMAs. We also provide details related to the numerical issues involved with fitting LMAs via MCMC. In Section 4 we detail our discrete space LMA model and its relation to the GTF estimates of Wang et al. [28]. In Section 5 we consider four datasets on which the LMA model is compared to its GRF counterpart. We conclude with a discussion in Section 6.

2 Spatial Models with Matérn Random Fields

In this section we provide background information to assist in developing the framework of the hierarchical spatial models considered in this work. We begin by describing the SGLMM, and follow with a discussion of the Matérn random field as a solution to a stochastic partial differential equation (SPDE).

2.1 SGLMMs

Generalized linear models (GLMs) model the mean $\boldsymbol{\mu}$ of a distribution $f$ with linear predictors through an invertible link function $g(\cdot)$. Generalized linear mixed models (GLMMs) are GLMs that allow for the inclusion of random effects in the linear predictor. The spatial GLMM (SGLMM) attempts to capture an unobserved spatially varying trend by imposing a dependence structure in the random effect.

Let $\{\boldsymbol{u}_{i}\}_{i=1}^{N}$ be a collection of locations observed in some spatial domain $\Omega\subset\mathbb{R}^{d}$. Consider $Y(\boldsymbol{u}_{i})\sim f(y_{i})$ such that the mean, \mu_{i}=E\left(Y(\boldsymbol{u}_{i})\mid\boldsymbol{\beta},\eta(\boldsymbol{u}_{i}),\mbox{\epsilon}(\boldsymbol{u}_{i})\right), is modeled through the link function

[TABLE]

In some cases an uncorrelated random effect, \mbox{\epsilon}(\boldsymbol{u})|\sigma^{2}\sim\mathcal{N}(0,\sigma^{2}), is included. In the Gaussian response case, $\sigma^{2}$, is thought of as homogeneous measurement error. The spatially varying random effect $\eta(\boldsymbol{u})|\boldsymbol{\theta}$ is included in (1) to capture spatial dependence. The hyper-parameters ($\boldsymbol{\theta}$) govern the mean and covariance structure of the spatial random effect. The most common distribution assumed for $\eta(\boldsymbol{u})|\boldsymbol{\theta}$ is Gaussian. This work considers replacing the traditional Gaussian prior for a less common Laplace moving average (LMA). We provide a thorough comparison of the two prior choices over four datasets in Section 5.

2.2 Random Fields with Matérn Covariance

Stationary Matérn random fields arise as stationary solutions to the SPDE

[TABLE]

where $\triangle=\sum_{i=1}^{d}\frac{\partial^{2}}{\partial u_{i}^{2}}$ is the Laplacian operator, $\kappa$ is a spatial scale parameter, $\alpha$ controls the smoothness of the realized fields, and $\xi$ is a variance parameter. The SPDE in (2) is driven by Gaussian white noise, $\mathcal{W}(\boldsymbol{u})$. Whittle [29, 30] showed that stationary solutions to (2) are GRFs with Matérn covariance

[TABLE]

where $||\cdot||$ denotes the Euclidean norm in $\mathbb{R}^{d}$. The marginal variance is $\phi^{2}=\left(\xi^{2}\Gamma(\nu)\right)/\left(\Gamma{(\alpha)}(4\pi)^{d/2}\kappa^{2\nu}\right)$ and the approximation of the effective range is given by $\rho=\sqrt{8\nu}/\kappa$ [15]. Lindgren et al. [15] provided a finite element representation of the SPDE in (2). In Section 3.2 we detail Lindgren’s approximation. The result is a sparse Gauss Markov random field (GMRF) approximation to a GRF with Matérn covariance given in (3).

3 LMA Models for Point Referenced Data

GRF modeling remains appealing to scientists due to familiarity of implementation and the intuitive dependence structure of Gaussian processes. The Laplace distribution offers a wider tailed, sharper peaked, alternative to the Gaussian distribution (see Figure 1). In this Section we provide an overview of the results of Bolin [6], which demonstrate that the LMA can be expressed as the solution to an SPDE similar to its Gaussian counterpart. We then lay out the finite element approximations of Lindgren et al. [15] and Bolin [6] which provide sparse representations of the analytic solutions of SPDEs driven by Gaussian and Laplace noise respectively. We conclude the Section by providing a novel exploration of model fitting via MCMC.

3.1 Laplace Moving Average Models as SPDEs

Gaussian priors often produce marginal distributions with light tails. Ȧberg and Podgórski [1] suggested the use of LMAs to obtain asymmetric and heavier tailed marginals. Ȧberg and Podgórski [1] showed that the LMA can be expressed as a convolution of a kernal with Laplace noise. Bolin [6] showed that LMA models with Matérn covariance can equivalently be expressed as the solution to an SPDE by replacing the Gaussian white noise, $\mathcal{W}(\boldsymbol{u})$, in (2) with Laplace noise. Laplace noise on a compact set $\Omega\subset\mathbb{R}^{d}$ can be expressed as

[TABLE]

where $G_{k}\stackrel{{\scriptstyle iid}}{{\sim}}\mathcal{N}(0,1)$, $\Gamma_{k}\sim e^{-\nu\gamma_{k}}W_{k}$, $W_{k}$ are iid standard exponentials, $\gamma_{k}$ are arrival times of a Poisson process with intensity one, and $\delta_{\boldsymbol{u}_{k}}$ is a Dirac distribution centered at $\boldsymbol{u}_{k}$ with each $\boldsymbol{u}_{k}$ uniformly distributed in $\Omega$ (see Bolin [6] for details).

The SPDE defining an LMA process, $\eta(\boldsymbol{u})$, with Matérn covariance is given by

[TABLE]

Both the SPDEs driven by Gaussian and Laplace noise produce random fields with Matérn covariance. However, the SPDE driven by Laplace noise provides a random field with the ability to capture “spikey” spatial behavior better than GRFs. Unlike its Gaussian counterpart, there does not exist a closed form solution for the Laplace driven SPDE in (5). Bolin [6] used a finite element approximation of the SPDE in (5) to fit LMA models using an EM algorithm. In Sections (3.2–3.3) we detail the finite element approximations to the Gaussian and LMA SPDEs proposed by Lindgren et al. [15] and Bolin [6] respectively.

3.2 Finite Element Approximations to Matérn GRFs

The finite element representation of the LMA model proposed by Bolin [6] is an extension of the results in the Gaussian case. For the sake of comparison, we first detail the finite element approximation of Lindgren et al. [15].

In Section 2.2 we expressed GRFs with stationary Matérn covariances as analytic solutions to SPDEs. Though the analytic solution provides useful insights, model fitting and parameter estimation are often facilitated by considering a numerical approximation. Lindgren et al. [15] proposed the use of a finite element approximation to the stochastic weak formulation of the SPDE

[TABLE]

where $\{\psi\}$ is a set of test functions. The finite element method (FEM) solution begins by expressing the solution, $\eta(\boldsymbol{u})$, as a basis expansion

[TABLE]

where $\{\phi_{i}(\boldsymbol{u})\}_{i=1}^{n}$ is a set of basis functions on $\Omega$. The solution in (6) is only required to hold for a finite collection of $\psi_{i}$. The Galerkin method approximate solution is obtained by setting $\{\psi_{i}\}_{i=1}^{n}=\{\phi_{i}\}_{i=1}^{n}$.

Lindgren et al. [15] formulated an FEM approximation by considering $\{\phi_{i}(\boldsymbol{u})\}_{i=1}^{n}$ to be piecewise triangular basis functions. The basis functions are constructed by partitioning the spatial region of interest, $\Omega\subset\mathbb{R}^{d}$, into non-overlapping triangular regions. The corners of the triangles, referred to as vertices, are assigned $n$ Gaussian weights, denoted $w_{i}$. Each $\phi_{i}$ is defined to be 1 at vertex i and 0 at all other vertices. Lindgren et al. [15] derived the distribution of the weights

[TABLE]

where $\boldsymbol{Q}_{\alpha}$ is determined by the choice of $\alpha$ in (6).

Following Lindgren et al. [15], the precision matrix in (8) is defined as

[TABLE]

where $\boldsymbol{L}=\kappa^{2}\boldsymbol{C}+\boldsymbol{G}$. The matrices used to define L, are given by

[TABLE]

Note that C as defined in (10) is sparse, but its inverse, $\boldsymbol{C}^{-1}$, which is required in $\boldsymbol{Q}_{\alpha}$ for $\alpha\geq 2$, may not be. In turn, the resulting precision matrix, $\boldsymbol{Q}_{\alpha}$, may not be sparse. To ensure sparsity in $\boldsymbol{Q}_{\alpha}$, Lindgren et al. [15] provided a GMRF approximation to the GRF representing the numerical solution to the SPDE in (6) by replacing C in (10) with the diagonal matrix,

[TABLE]

Under lattice refinement, the sparse representation of C as defined in (12) converges to the same solution given by C in (10) (see Lindgren et al. [15] Appendix C.5). C in (12) is now a diagonal matrix relating to the volume of the regions produced by the mesh. G is a sparse matrix with zero entries on the diagonal, describing the connectivity of the mesh nodes.

3.3 Finite Element Approximations to Matérn LMAs

Bolin [6] extended the results of Lindgren et al. [15] to the case of Type-G Matérn random fields. We consider the FEM approximation for the symmetric LMA model with Matérn covariance. Similar to the Gaussian case, the FEM approximation begins with a stochastic weak formulation of the SPDE in (5) given by

[TABLE]

Following Section 3.2 we construct piecewise linear basis functions $\{\phi_{i}\}_{i=1}^{n}$. Let $\Gamma_{i}\sim Gamma\left(\tau C_{ii},\lambda^{2}\right)$ where $C_{ii}$ is the $i^{th}$ element of C in (12). Define $\boldsymbol{\Gamma}=diag(\Gamma_{1},...,\Gamma_{n})$. Bolin [6] showed that the distribution of the basis expansion weights given by the Galerkin method is

[TABLE]

where

[TABLE]

with $\boldsymbol{L}=\kappa^{2}\boldsymbol{C}+\boldsymbol{G}$ defined as in the Gaussian case. We note that this Section contains all the information needed to construct and fit an FEM solution to LMA models. We refer the reader to Bolin [6] for a mathematically rigorous construction of the FEM solution.

There currently exists no extension of the LMA approximation for odd valued $\alpha$’s. We also point out that (14) results in a sparse precision matrix for small values of $\alpha$. For $\alpha=2$, the precision matrix, $\boldsymbol{K}_{2}\boldsymbol{\Gamma}^{-1}\boldsymbol{K}_{2}$, defines a sparse representation for the roughest covariance function offered by the finite element approximation to the LMA. This corresponds to Matérn covariance with $\nu=1$ for spatial models in $\mathbb{R}^{2}$ and $\nu=1.5$ for $\mathbb{R}$.

In summary, we see that the LMA model can be expressed as a sparse conditionally Gaussian distribution conditioned on Gamma-distributed auxiliary variables. We discuss the numerical issues and computational considerations implied by this approximation in Section 3.4.

3.4 Model Fitting: Continuous Space

Parameter estimation of LMAs is difficult since no closed form of the likelihood exists. Ȧberg and Podgórski [19] proposed a method of moments-based estimation for LMA models. In Section 3.3 we showed that expressing the LMA as an SPDE allows for inference in the likelihood framework. Bolin [6] constructed an EM algorithm for parameter estimation following the FEM approximation described in Section 3.3. Wallin and Bolin [27] considered the use of an MCEM algorithm in order to provide a computationally efficient extension to SGLMMs. Persistent numerical issues contributed to difficult parameter estimation, in both the works of Bolin [6] and Wallin and Bolin [27]. The issues stem from the fact that $\text{E}\left[\Gamma_{i}^{-1}|\boldsymbol{\cdot}\right]$ is unbounded for small $\underset{1\leq i\leq n}{\min}|\tau C_{ii}-1/2|$. Bolin [6] suggested truncating the expected value at 1000 to avoid numerical instabilities.

To our knowledge, there has been no systematic Bayesian comparison of SGLMMs with LMAs and GRFs fit via MCMC. Samplers for GRFs require samples from the full-conditionals of the $n$-dimensional basis expansion weights $\boldsymbol{w}$, variance parameter $\xi$, spatial scale parameter $\kappa^{2}$, fixed effects $\boldsymbol{\beta}$, and a homogeneous random effect variance $\sigma^{2}$ (if applicable). Samplers for the continuous space LMA require n additional auxiliary variables ($\Gamma_{i}^{\prime}s)$ and a shape parameter $\tau$.

For the sake of illustration, consider a continuous response point referenced model with N unique observed locations denoted $\{\boldsymbol{u}_{i}\}_{i=1}^{N}$. Define the N by n projection matrix $(\boldsymbol{A})$ with entries $A_{ij}=\phi_{j}(\boldsymbol{u}_{i})$, where $\{\phi_{l}(\boldsymbol{u})\}_{l=1}^{n}$ are triangular basis functions formed as described in Section 3.2. Using the basis expansion of $\eta(\boldsymbol{u})$ (see (7) Section 3.2) and assuming homogeneous error measurement \mbox{\epsilon}(\boldsymbol{u}_{i})\stackrel{{\scriptstyle iid}}{{\sim}}\mathcal{N}(0,\sigma^{2}), we can write the discretized likelihood as $[\boldsymbol{y}|\boldsymbol{\beta},\boldsymbol{w},\sigma^{2}]\sim\mathcal{N}\left(\boldsymbol{X}\boldsymbol{\beta}+\boldsymbol{A}\boldsymbol{w},\sigma^{2}\boldsymbol{I}\right)$.

Conditioned on n auxiliary gamma random variables, $\Gamma_{i}$, we express the LMA as a scale mixture of normals with gamma variance. For $\alpha=2$ it was shown in Section (3.3) equation (14) that $\boldsymbol{K}_{2}\boldsymbol{w}|\boldsymbol{\Gamma}\sim\mathcal{N}(\boldsymbol{0},\boldsymbol{\Gamma})$, where $\boldsymbol{\Gamma}=diag(\Gamma_{1},...,\Gamma_{n})$. The conjugate full-conditional for the weights of the finite element approximation of the LMA model are

[TABLE]

Let $\boldsymbol{t}=\boldsymbol{K}_{\alpha}\boldsymbol{w}|\boldsymbol{\Gamma}\sim N(\boldsymbol{0},\boldsymbol{\Gamma})$ as defined in equation (14) of Section 3.3. This leads to conjugate generalized inverse Gaussian (GIG) full conditionals for each gamma random variable, $\Gamma_{i}|t_{i},\tau,\lambda\sim GIG\left(\tau C_{ii}-1/2,2/\lambda^{2},t_{i}^{2}\right)$, where the GIG(p,a,b) density is given by

[TABLE]

Thus conjugate updates are available for the FEM weights $\boldsymbol{w}$ and the auxiliary variables $\Gamma_{i}$ in the LMA model.

In practice we found that Gibbs updates for the conjugate $\Gamma_{i}^{\prime}s$ resulted in poor-mixing (see Sections 5.3–5.4). This is not an uncommon issue, as conjugacy does not always produce an efficient sampler [22]. We found that the overall mixing of the Markov chains were improved by using independent one-at-a-time adaptively tuned Metropolis Hastings updates for each $\Gamma_{i}$.

It is possible that samples from $\Gamma_{i}|t_{i},\tau,\lambda$ are close to zero. This results in a numerically negative-definite precision matrix $\left[\boldsymbol{L}\boldsymbol{\Gamma}^{-1}\boldsymbol{L}+\left(\frac{1}{\sigma^{2}}\right)\boldsymbol{A}^{\prime}\boldsymbol{A}\right]^{-1}$ for the full-conditional distribution in (16). To overcome this issue, we re-sampled the $\Gamma_{i}^{\prime}s$ if $\left[\boldsymbol{L}\boldsymbol{\Gamma}^{-1}\boldsymbol{L}+\left(\frac{1}{\sigma^{2}}\right)\boldsymbol{A}^{\prime}\boldsymbol{A}\right]^{-1}$ was numerically rank-deficient. This amounts to a Metropolis Hastings update for $\boldsymbol{w}$ and $\boldsymbol{\Gamma}$ with the constraint that $\left[\boldsymbol{L}\boldsymbol{\Gamma}^{-1}\boldsymbol{L}+\left(\frac{1}{\sigma^{2}}\right)\boldsymbol{A}^{\prime}\boldsymbol{A}\right]^{-1}$ is positive-definite.

In binary error response SGLMMs, Gaussian full-conditionals can be constructed using data augmentation [2]. The LMA model can then be fit via MCMC exactly as described above. For Poisson error response distributions, the full-conditionals of $\boldsymbol{w}$ are not Gaussian. In this case, we suggest the use of conditionally independent one-at-a-time Metropolis Hastings updates for each $w_{i}|\boldsymbol{w}_{-i}$ (see Appendix A.3). We note that one-at-a-time conditionally independent block sampling is applicable for Gaussian and binary responses as well. However, in practice we found that block proposing $\boldsymbol{w}$ and $\boldsymbol{\Gamma}$ produced a faster and more efficient sampler for Gaussian and binary error responses. With the above numerical considerations we were able to fit the LMA models with a sampler that is a simple extension of traditional Gaussian samplers for each of point-referenced data analyses considered in Sections (5.3–5.4).

In summary, we have shown how LMA models with Matérn covariance can be expressed as an SPDE. We demonstrated that the FEM approximations for the GRF and LMA models result in sparse conditionally Gaussian representations. Following the numerical cautions detailed in this Section, we were able to fit the LMA using MCMC. This allows for Bayesian analyses familiar to traditional SGLMM models with GRFs.

4 LMA Models in Discrete Space

In the previous Section, we detailed the construction of the LMA model and provided a method for parameter estimation via MCMC. Wallin and Bolin [27] acknowledged the potential use for LMAs in discrete space, but no further investigation was pursued. In this Section, we present a novel Bayesian hierarchical formulation for discrete space SGLMMs with LMAs. Our model is shown to be an adaptation of the widely recognized graph trend filtering (GTF) estimates proposed by Wang et al. [28]. Similar to the continuous space model, we introduce auxiliary variables to express the discrete space LMA as a conditionally Gaussian distribution. The resulting model exhibits computing and inference similar to its Gaussian analogue.

4.1 Graph Trend Filtering

Heavier tailed alternatives to GRFs have been shown to provide improved predictive power [25, 9, 28] for discrete space models. Wang et al. [28] extended the univariate trend filter of Kim et al. [14] to irregular graphs. We provide a Bayesian extension of the graph trend filer (GTF) that is analogous to the LMA model for discrete space SGLMMs.

To motivate our discrete space model, we provide a brief overview of GTF estimates. Let G=(V,E) be a graph with vertices $i=\{1,...,n\}$ and undirected edges $\{e_{1},...,e_{m}\}$. Suppose that $\boldsymbol{y}=(y_{1},...,y_{n})$ are observed at the vertices. The GTF estimates, $\hat{\boldsymbol{\beta}}=(\hat{\beta}_{1},...,\hat{\beta}_{n})$, are the solution to

[TABLE]

where $\triangle^{(k+1)}$ is a recursive graph difference operator of order k, and $\lambda$ is a shrinkage parameter determined by cross-validation. For the case of $k=0$, $\triangle^{(1)}$ is defined such that $||\triangle^{(1)}\boldsymbol{\beta}||_{1}=\sum_{\begin{subarray}{c}(i,j)\in E\end{subarray}}|\beta_{i}-\beta_{j}|$ produces a first order difference penalty over G. The $k^{th}$ order GTF differencing matrix seen in (18) is recursively defined as

[TABLE]

Higher orders of $k$ in (19) correspond to higher order differencing penalties. For a chain graph, or a graph over a one-dimensional line, the GTF estimates reduce to the trend filter estimates of Kim et al. [14].

In this Section 4.2, we propose a discrete space version of the LMA for SGLMMs. The GTF estimates are shown to be a special case of our model. In turn, our model specification serves as a generalized extension of the GTF to SGLMMs in the Bayesian framework.

4.2 LMAs in Discrete Space

Areal datasets are composed of aggregated responses. Examples of common spatial aggregations include cumulative measurements for cities, states, or countries. The areal units at which responses were recorded determine the discretization of space. The spatial discretization can be summarized by a graph with nodes at each areal unit and undirected edges defined by the spatial neighborhood structure of the areal units.

The graph detailing the spatial neighborhood relationships can be expressed as a graph Laplacian matrix, A, of the form

[TABLE]

The graph Laplacian (A) serves as the precision matrix of the popular ICAR model for areal spatial random effect models [3]. A is positive semi definite with rank n-1, where n is the number of observations. It is common to add a positive value to the diagonal to ensure that A is diagonally dominant. Define $\boldsymbol{L}=\kappa^{2}\boldsymbol{I}+\boldsymbol{A}$. Let D be an upper triangular matrix such that $\boldsymbol{L}=\boldsymbol{D}^{\prime}\boldsymbol{D}$ (i.e., D could be a Cholesky decomposition). We consider the discrete space precision matrices given by,

[TABLE]

where

[TABLE]

Traditional simultaneously autoregressive (SAR) models (see Appendix A.2) assume a Gaussian prior for the n-dimensional random effect. We assume a Gaussian prior on the weighted sum of differences

[TABLE]

which implies that

[TABLE]

The discrete space LMA is obtained by placing an iid Laplace prior on each sum of weighted differences $\triangle^{(k)}\boldsymbol{\eta}|\lambda,\kappa^{2}$. That is,

[TABLE]

which implies that

[TABLE]

We observe that the discrete space LMA model is obtained by replacing the Gaussian prior in (24) for an iid Laplace prior in (26).

4.3 L2 vs. L1 Penalization

In the proceeding Section we demonstrated that the discrete space LMA is obtained by replacing the Gaussian prior on the weighted sum of differences in (23) with iid Laplace priors in (25). Here we detail the penalization implications that result from altering the prior. The contrast in penalizations is presented to provide intuition as to why the LMA model should be considered in place of the the GRF model.

Let $N(i):=\{j:j\text{ is a neighbor of }i\}$ and $\boldsymbol{\theta}_{G}=(\xi^{2},\kappa^{2},\sigma^{2})$. The order $k=1$ discrete space GRF has log full-conditional distribution

[TABLE]

Notice that (27) resembles penalized regression with a squared penalty term placed on the weighted sum of differenced nodes ($\eta_{i}$). If $\kappa^{2}>0$ in (27), $\boldsymbol{\eta}$ is penalized based on the sum of the localized differences relative to the magnitude of each node. We recover the intrinsic conditionally auto-regressive model (ICAR) by taking $\kappa^{2}=1$ and $k=0$.

Now define $\boldsymbol{\theta}_{L}=(\lambda,\kappa^{2},\sigma^{2})$. The order $k=1$ discrete space LMA log full-conditional is given by

[TABLE]

Our inclusion of $\kappa^{2}$ allows for a full rank precision matrix, as well as an L-1 penalty based on the magnitude of each node, $\eta_{i}$. The Bayesian equivalent to the $k^{th}$ order GTF estimates proposed by Wang et al. [28] are obtained by setting $\kappa^{2}=0$ in (28). We also point out that for $\kappa^{2}=1$ and $k=0$, we recover the Laplace analog to the ICAR model.

We note that the LMA model resembles a LASSO style (L-1) penalty while the GRF model resembles a ridge like (L-2) penalty on the sum of differenced nodes ($\eta_{i}$). There appears to be an analogous extension to be made for penalization in the continuous space. For a discrete space model on a regular grid, lattice refinement results in convergence to the continuous space LMA.

4.4 Model Fitting: Discrete Space

We aim to preserve the familiarity of the Gaussian fit for our discrete space model. To do so, we first recognize the Laplace distribution as a scale mixture of normals. This can be seen by taking a Gaussian variable, $Z\sim N(0,S_{i})$, where $S_{i}$ is an independent exponential random variable with rate $\lambda^{2}/2$. Marginalizing over $S_{i}$, it follows that $Z|\lambda\sim\mathcal{L}(\lambda)$.

Recall, the discrete space LMA places an $\mathcal{L}(\lambda)$ prior on the $l^{th}$ weighted difference, $\triangle^{(k)}_{l}\boldsymbol{\eta}|\lambda,\kappa^{2}\stackrel{{\scriptstyle iid}}{{\sim}}\mathcal{L}(\lambda)$. By introducing $n$ auxiliary variables $S_{i}\stackrel{{\scriptstyle iid}}{{\sim}}Exp(\lambda^{2}/2)$, we can express the prior as

[TABLE]

Park and Casella [18] showed that $S_{i}^{-1}$ are conjugate inverse Gaussian ($InvGauss$). To see this, first define $\boldsymbol{t}=\triangle^{(k)}\boldsymbol{\eta}|\boldsymbol{S},\kappa^{2}$. Then $S_{i}^{-1}|t_{i},\lambda\sim InvGauss\left(\lambda^{2},\sqrt{\frac{\lambda^{2}}{t_{i}^{2}}}\right)$, where the $InvGauss(a,b)$ has density

[TABLE]

We note that the InvGauss is a special case of the GIG with $p=-1/2$. We fit the discrete space LMA and GRF via MCMC as well. From (29), we observe that the discrete space LMA model can be expressed as a conditionally sparse Gaussian distribution. Additionally, the full-conditionals for the auxiliary random variables are conjugate. In turn, Gibbs updates can be used for the full-conditionals of the n auxiliary random variables. In Section 3.4 we discussed the numerical issues induced by using conjugate updates for the auxiliary mixing variables. This is not an issue in the discrete space sampler, as the InvGauss does not require estimation of a shape parameter.

In summary, we have provided a novel discrete space Bayesian hierarchical LMA. The LMA is conditionally Gaussian, as in continuous space, again allowing for inference and computing familiar to GRF models. In Section 5 we demonstrate that the discrete space LMA offers improved out-of-sample predictive power in the presence of localized trends.

5 Example Analyses

In this section we compare the LMA and GRF models over four datasets; including Gaussian, Poisson, and binary responses observed both on continuous (point referenced) and discrete (areal) support. We perform 10-fold cross validation by randomly splitting the data set into 10 roughly equal sized groups. We withhold a test set $\boldsymbol{y}_{k}$ and train the model on the remaining observations $\boldsymbol{y}_{-k}$. We use the Bayesian cross validation scoring criterion (BCVS) of Hooten and Hobbs [13] given by

[TABLE]

where $T$ is the total number of stored MCMC iterations, and $\boldsymbol{\theta}^{(t)}$ is the $t^{th}$ sample of the model’s parameters. Note that a smaller value implies a better model fit.

We used effective samples per second (ESS) to compare the computational performance of each model fit. All MCMC iterations had a burn-in phase in which the normal proposals were adaptively tuned according to Roberts and Rosenthal [20]. Every data analysis involved 50,000 post-burn-in MCMC states. The number of burn-in iterations for the GRF and LMA model were held constant for each individual data analysis. The ESS was computed based on the likelihood evaluated at the stored 50,000 states.

5.1 Discrete Space Examples: Slovenia Stomach Cancer

We compare models trained on areal count data with the Slovenian stomach cancer outbreak dataset. The dataset consists of 194 responses of aggregated stomach cancer counts in each municipality of Slovenia collected from 1994 to 2001. The model proposed by Hodges and Reich [12] was fit to investigate the relationship between standardized socioeconomic status, $SEc_{i}$, and the occurrence of stomach cancer, $y_{i}$.

We fit the Poisson spatial model considered by Hodges and Reich [12] assuming stomach cancer counts are Poisson, $y_{i}\sim Poisson(\mu_{i})$. The mean $\mu_{i}$ is modeled through the log-link function

[TABLE]

The offset term $(o_{i})$ in (32), is the expected number of stomach cancer counts for observation i. The covariate matrix $(\boldsymbol{X})$ contains an intercept and the standardized socioeconomic status ( SEc ).

We use an order $k=1$ differencing matrix $\triangle^{(1)}$ (see equation (22) in Section 4.2) to define the priors for the spatially correlated random effects in the GRF and LMA models defined in equations (23) and (25) of Section 4.2. The fixed effects $(\boldsymbol{\beta})$ are assigned iid normal priors with variance $10^{6}$. The variance of the spatially homogeneous random effect $(\sigma^{2})$ is assigned an inverse gamma prior with shape and scale of one. The scale parameter for the GRF model $(\xi)$, the scale parameter for the LMA model $(\lambda)$, and the spatial scale parameter $(\kappa)$ are all assigned independent half-normal priors with scale one. The full-conditionals for the Poisson response data are detailed in Appendix A.4.

Inference on the fixed effects is similar between the LMA and GRF models (see Table 3). Table 1 shows that the GRF provides a better fit than the LMA. This is perhaps due to the areal dataset being smooth. We do note that the LMA sampler produces roughly the same ESS.

5.2 Discrete Space Examples: Columbus Crime Data

The Columbus crime data are found in the “spdep” R package [5]. The dataset provides a spatial map with crime rates ($y_{i}$) for each county of Ohio. Ver Hoef et al. [26] suggested modeling the data as an intercept only model with Gaussian response. We include average household income and average household value, as well as an intercept as covariates. Figure 3 shows a plot of the crime rates in each county of Ohio. There appears to be a few localized spikes in the data, namely the counties that appear in white. This is sufficient reason to suspect that the LMA model should provide an improved model fit.

The model is of form

[TABLE]

We use the adjacency matrix (A) provided from the “spdep” package to construct the first order differencing matrix $\triangle^{(1)}$ (see equation (22) in Section 4.2). The precision matrices for the priors of the GRF and LMA models are constructed from $\triangle^{(1)}$ as in Section 4.2. We assumed iid normal variance $10^{3}$ priors for the fixed effects. For both models $\kappa^{2}$ was assigned a half-normal scale one prior. The variance parameter for the GRF model $(\xi^{2})$ and the variance parameter for the LMA model $(\lambda^{2})$ were assigned half normal scale 10 priors. The full-conditionals for all GRF and LMA model parameters can be found in Appendix A.5.

In Table 2, we see that the LMA model is the clear favorite in terms of BCVS. Inference on the fixed effects $\boldsymbol{\beta}$ is similar in both cases (see Table 3). The ESS for the LMA model seems substantially smaller than the GRF model. We note that the data set is only of size 49, so this difference may be a bit inflated. The ESS for both models is quite large for an MCMC sampler.

5.3 Continuous Space Examples: Malaria in the Gambia, Africa

A model comparison for the binary continuous case is illustrated with the use of presence absence data of malaria in the Gambia, Africa. The data was made publicly available by Diggle and Riberio [8] in the “geoR” package of R [4]. The dataset consists of 2035 children records recorded at 65 village locations denoted $\{\boldsymbol{u}_{i}\}_{i=1}^{65}$. Each village, i, has $n_{i}$, respondents. We consider a model similar to that of Hanks et al. [11]. Let $y^{(i)}_{j}$ be the indicator for the presence ($y^{(i)}_{j}=1$) of malaria in the $j^{th}$ child at village location $i$. The covariates considered for each child are composed of the child’s age, an indicator of whether or not a bed net was used, an indicator for whether or not an insecticide was applied to the bed net, the log normalized difference vegetation index (NDVI) at each village, and an indicator of presence or absence of a health center in the village.

We construct a triangular mesh with $n=288$ nodes containing all 65 village locations on a node (see Figure 4). We specify a binary probit model with $y^{(i)}_{j}\sim\text{Binom}\left(1,p_{j}^{(i)}\right)$, where the probability of malaria presence in the $j^{th}$ child at village location $i$ is linked through the probit function (standard normal CDF) and auxiliary data

[TABLE]

The superscripts for each variable denote the village location $\boldsymbol{u}_{i}$. We fix $\alpha=2$ leading to covariance matrices for the GRF and LMA priors for $\eta(\boldsymbol{u})$ in (34) given by (9) and (15) respectively. We assume $iid$ normal variance 10 priors for all fixed effects. We follow the data augmentation approach of Albert and Chib [2] for binary probit GLMs to obtain Gibbs updates for the fixed and random effects. The priors for hyper-parameters $\xi,\kappa,\tau,\lambda$ are all iid half normal scale one. Further details pertaining to prior specifications and model fitting are included in Appendix A.6.

For this dataset, we randomly split the dataset into 10 groups of village locations. All observations associated with a given test set were withheld for validation. Table 4 shows that the LMA model provides the better fit based on BCVS. In this case the LMA also provided a higher ESS. This may be due to the use of a data augmentation approach suggested by of Albert and Chib [2]. Estimates for all model parameters can be found in Table 6.

5.4 Continuous Space Example: LAGOS

Lastly, we compare the performance of the LMA for the continuous response LAGOS lake dataset. The Lake multi-scaled geospatial and temporal database (LAGOS) is a publicly accessible US lake water quality database [24]. The dataset used in this paper contains records for 5526 unique lakes observed over Iowa, Missouri, and Illinois at locations $\{\boldsymbol{u}_{i}\}_{i=1}^{5526}$. We are interested in modeling the log total phosphorus recorded in each lake. First we reduce the many covariates recorded for each lake in the LAGOS database by performing a step-wise regression assuming uncorrelated residuals. The covariates selected by AIC in this stepwise procedure are: an intercept, lake area (in hectares), max depth of the lake (meters), mean runoff (ground-water discharge into streams), the average annual runoff (in/yr), inter lake water shed (IWS) measurements (the area of land that drains directly into a lake) for urban, agricultural, road density, and total wetland regions. We then used these selected covariates to model log total phosphorus ($y(\boldsymbol{u})$) as

[TABLE]

We fit the LMA model and GRF model by fixing $\alpha=2$. Again this corresponds to Matérn smoothness parameter $\nu=1$. We construct a triangular mesh with $n=671$ nodes according to Section 3.2. The fixed effects $\boldsymbol{\beta}$ are assigned a normal prior with variance $10^{3}$. The nugget term $(\sigma^{2})$ and spatial scale parameter $(\kappa)$ are assumed to follow independent half-normal distribution with scale one. The variance term of the GRF $(\xi)$, as well as the variance $(\lambda)$ and scale $(\tau)$ parameters of the LMA are assigned independent half-normal scale one priors.

Table 5 shows that the LMA model provides a better fit than the GRF in this case. The ESS is better for the LMA as well, further suggesting the LMA model should be favored. As in all other cases, the fixed effect estimates remain roughly similar. All model parameter estimates can be found in Appendix A.7 Table 6. Details of the LAGOS data analysis and all full-conditionals are included in Appendix A.7.

6 Discussion

In this work we proposed a novel discrete space LMA model for irregular lattices and constructed Metropolis Hastings samplers for continuous and discrete space SGLMMs with LMAs. Using the Metropolis Hastings samplers, we provided a broad Bayesian analysis of SGLMMs with LMAs for continuous, binary and Poisson error responses. We found that the LMA offered a better fit than the GRF for datasets which exhibit “spikes” in the response. Provided our MCMC implementation, we saw that the LMA offers similar computational performance to the GRF.

In this work we restricted our attention to the symmetric LMA. Bolin [7] considered more general continuous space LMA which allowed for asymmetric posterior distributions. We note that the asymmetry parameter can be estimated in the MCMC framework, however we elected to compare symmetric models only. We also note that the extension of the asymmetric LMA to discrete space models on irregular graphs could be considered as well. We leave this for future work, but note that the extension should be straightforward.

The choice of half-normal priors for each of the complexity parameters in the GRF and LMA models, given by $\xi$, $\kappa$, $\tau$, and $\lambda$, were made for the sake of comparison. Fuglstad et al. [10] proposed the use of a joint prior on the variance $\xi$ and scaling parameter $\kappa^{2}$ in the Gaussian case following the work of Simpson et al. [23]. We note that their prior choice was motivated by the desire to provide a weakly informative prior and deal with a partial identifiability issue associated with the two parameters. Though identifiability issues with the complexity parameters persist, one does not usually concern themselves with the value of the parameter. We found that independent half-normal priors resulted in the best mixing.

In summary, we have developed a novel discrete space SGLMM with LMAs. We have proposd the use of Metropolis Hastings samplers to fit the LMA models as a simple extension of the GRF models. Through our extensive data analyses, we have provided evidence of cases in which the LMA model outperforms the traditional GRF SGLMM, while maintaining similar computation efficiency. Due to the comparable computation times and similarity of implementation, we recommend LMA models be considered when modeling correlated error structures.

Appendix A Appendix

A.1 CAR Models

Definition: (Conditionally Autoregressive Model) A Conditionally Autoregressive Model (CAR) takes on the form

[TABLE]

where C is the spatial dependence matrix with $C_{ii}=0$, and M is a diagonal matrix with entries $M_{ii}$. The conditional mean of each $\eta_{i}$ is determined by a weighted sum of neighboring $\eta_{j}$’s. Note that each marginal variance, $M_{ii}$, varies, so M is often non stationary.

The CAR model in (36) was shown to lead to the full distribution of $\boldsymbol{\eta}$ by Besag (1974). For positive definite $\boldsymbol{Q}^{-1}=(\boldsymbol{I}-\boldsymbol{C})^{-1}\boldsymbol{M}$, (36) leads to the full distribution $\boldsymbol{\eta}\sim N(0,\boldsymbol{Q}^{-1})$. Matrices M and C are defined from Q as follows: Write, $\boldsymbol{Q}=\boldsymbol{D}-\boldsymbol{R}$ with

[TABLE]

This gives $\boldsymbol{M}=\boldsymbol{D}^{-1}$ and $\boldsymbol{C}=\boldsymbol{DR}$ in (36).

A.2 SAR Models

Ver Hoef et al. [26] summarized the relationships between SAR and CAR models. CAR and SAR models are widely used in both temporal and spatial statistics due to their intuitive dependence structures. We provide a brief summary of the SAR model for the unfamiliar audience.

Consider a collection of random variables at n spatial locations or graph nodes, $\boldsymbol{Y}=(Y_{1}...,Y_{n})$. Let $\boldsymbol{\Lambda}$ be a positive diagonal matrix. A SAR model imposes an explicit spatial dependence structure,

[TABLE]

where B is a spatial dependency matrix that is not necessarily symmetric. Note that B relates $\boldsymbol{Y}$ to itself, and no site can depend on itself, so $B_{ii}$ must be zero for all i. Solving for $\boldsymbol{Y}$ in (37) we have, $(\boldsymbol{I}-\boldsymbol{B})\boldsymbol{Y}\sim\mathcal{N}(0,\boldsymbol{\Lambda})$. The covariance of the SAR model can then be written as $(\boldsymbol{I}-\boldsymbol{B})^{-1}\boldsymbol{\Lambda}[(\boldsymbol{I}-\boldsymbol{B})^{\prime}]^{-1}$, provided $(\boldsymbol{I}-\boldsymbol{B})$ is invertible. For a thorough comparison of SAR and CAR models see Ver Hoef et al. [26].

A.3 Conditionally Independent Block Proposals

Consider $\boldsymbol{\eta}\sim N(\boldsymbol{0},\boldsymbol{Q}^{-1})$, where $\boldsymbol{Q}$ is a GMRF. Define $N(i):=\{j:Q_{ij}\neq 0\}$. $N(i)$ is the collection of indices j, such that i and j are neighboring points in the spatial lattice. $\eta$ can be expressed as a CAR model with sparse Q

[TABLE]

where the procedure to obtain matrices M and C are described in Appendix A.1.

To produce one at a time Metropolis Hastings samples, we consider grouping subsets of conditionally independent $\eta_{i}$’s into blocks. Let $A_{k}$ be the collection of indices such that, for all $i,j\in A_{k}$

[TABLE]

We can now perform one at a time Metropolis Hastings updates for each individual $\eta_{i}$ within each block $A_{k}$.

A.4 Details of Slovenia Data Analysis

We consider the Poisson SGLMM of form

[TABLE]

where $o_{i}$ is an offset for individual $i$.

We use an order $k=1$ differencing matrix for the GRF and LMA model (see equation (22) of Section 4.2). We use the conditional distribution form of the CAR model to perform one at a time block Metropolis Hastings updates following the results of Appendix A.1–A.3. Matrices M and C of Appendix A.1 are defined for the GRF and LMA models from sparse $\boldsymbol{Q}_{1}$ of equation (21). Denote ${\eta}_{\mu_{i}}=\sum_{C_{ij}\neq 0}C_{ij}\eta_{j}$. The full-conditionals for $\boldsymbol{\eta}$ are given by

[TABLE]

In (39), $\boldsymbol{\theta}=(\kappa^{2},\xi)$ for the GRF model and $\boldsymbol{\theta}=(\kappa^{2},\lambda,\boldsymbol{S})$ for the LMA.

A normal prior with variance $\sigma_{\beta}^{2}=10^{6}$ is assumed for the fixed effects giving log full-conditionals

[TABLE]

The prior for the homogeneous spatial random effect is assumed to be iid $\mathcal{N}(0,\sigma^{2})$, giving log full-conditionals

[TABLE]

An inverse gamma prior with scale and shape one is assumed for $\sigma^{2}$ giving conjugate full-conditional

[TABLE]

The priors for the variance parameter $(\xi)$ and $(\kappa)$ for the GRF model are assumed to be independent scale one half-normals. The log full-conditionals are

[TABLE]

where $U_{ii}$ in (41) is the $i^{th}$ diagonal entry of the Cholesky decomposition of $\boldsymbol{L}$. The scale parameter, $\lambda$, and $\kappa$ of the LMA have log full-conditionals

[TABLE]

A.5 Details of the Columbus Crime Dataset Analysis

The model is of form

[TABLE]

where the covariates and response are detailed in Section 5.2. The fixed effects were assigned a normal prior with variance $10^{6}$ giving conjugate full-conditionals.

[TABLE]

The variance of the spatially homogeneous random effect $(\sigma^{2})$ is given a half-normal scale one prior. The log full-conditional is

[TABLE]

We use an order $k=1$ differencing matrix to define the covariance structure of the GRF and LMA model (see equation (22) of Section 4.2). The conjugate full-conditionals for the GRF random effects are

[TABLE]

The scale parameter for the GRF $(\xi)$ is given a half-normal scale 10 prior while $\kappa^{2}$ is given an independent half-normal scale one prior leading to log full-conditionals

[TABLE]

where $U_{ii}$ is the Cholesky decomposition of $\boldsymbol{L}=\triangle^{(1)}$. The conjugate full-conditionals for the LMA random effects are

[TABLE]

The scale parameter $(\lambda)$ for the LMA is given a half-normal scale 10 prior and $\kappa^{2}$ is given an independent scale one half-normal prior. The log full-conditionals for LMA model parameters are

[TABLE]

We observed spatial confounding among the random effects and the intercept. This is not uncommon, however to assess convergence we analyzed the trace plots of $\beta_{0}\boldsymbol{1}+\boldsymbol{\eta}$.

A.6 Details of Malaria Data Analyses

We follow the auxiliary data approach of Albert and Chib [2]. Let $\Phi(\cdot)$ denote the standard normal CDF. Consider the continuous space binary response model $y^{(i)}_{j}\sim Bernoulli(p^{(i)}_{j})$ where $p_{j}^{(i)}$ is the probability that the $j^{th}$ child in the $i^{th}$ village has malaria. We model $p_{j}^{(i)}$ through the probit link function by introducing auxiliary data $z^{(i)}_{j}$ as follows

[TABLE]

The covariates are as described in Section 5.3. Define $A_{ij}=\phi_{j}(\boldsymbol{u}_{i})$, where $\{\phi_{l}(\boldsymbol{u})\}_{l=1}^{n}$ are the basis functions corresponding to the triangular mesh with $n=288$ mesh nodes formed in Section 5.3. Define the 2035 by 65 matrix B such that $n_{i}$ entries of column $\boldsymbol{b}_{i}$ corresponding to responses $y^{(i)}_{j}$ are 1, and the remaining entries are 0. The auxiliary variables can be equivalently expressed in matrix form as

[TABLE]

The priors for $\boldsymbol{w}$ are constructed following Sections (3.2-3.3) with $\alpha=2$. We place a normal prior on the fixed effects with variance 100. In both models the spatial scale parameter $\kappa$ was assigned a half-normal prior with scale one. Additionally, in the Gaussian model, the scale parameter $(\xi)$ was assigned a half-normal scale one prior as well. We used joint MH proposals for $\kappa$ and $\xi$ in the Gaussian model. For the LMA, shape parameter $(\tau)$ and scale parameter $(\lambda)$ were jointly proposed with independent scale one half-normal priors.

**Full-Conditionals

**For both models, we have conjugate truncated normal ($TN_{(a,b)}$) updates for the auxiliary variables,

[TABLE]

and conjugate normal updates for the fixed effects,

[TABLE]

The weights of the basis expansion for the GRF are given by

[TABLE]

The full-conditional for the spatial scale $(\kappa)$ and variance parameter $(\xi)$ are given by

[TABLE]

where $U_{ii}$ denotes the $i^{th}$ diagonal entry of the Cholesky of $\boldsymbol{L}=\kappa^{2}\boldsymbol{C}+\boldsymbol{G}$. The weights of the basis expansion for the LMA are given by

[TABLE]

The spatial scale $(\kappa)$, shape $(\tau)$, and variance parameter $(\lambda)$, have log full-conditionals

[TABLE]

A.7 Details of LAGOS Analysis

We fit a continuous response point referenced model with 5526 unique lake locations denoted $\{\boldsymbol{u}_{i}\}_{i=1}^{5526}$. We form a mesh with $n=671$ nodes. The model considered is of the form

[TABLE]

Define the 5526 by 671 projection matrix $(\boldsymbol{A})$ with entries $A_{ij}=\phi_{j}(\boldsymbol{u}_{i})$, where $\{\phi_{l}(\boldsymbol{u})\}_{l=1}^{671}$ are the basis functions corresponding to the mesh formed in Section 5.3. Using the resulting basis expansion of $\eta(\boldsymbol{u})$ (see (7) Section 3.2) and assuming \mbox{\epsilon}(\boldsymbol{u}_{i})\stackrel{{\scriptstyle iid}}{{\sim}}\mathcal{N}(0,\sigma^{2}), we can re-write the discretized likelihood as follows

[TABLE]

The priors for $\boldsymbol{w}$ are constructed following Sections (3.2-3.3) with $\alpha=2$. We have assumed half-normal scale one priors for $\sigma^{2},\kappa,\xi,\tau$, and $\lambda$.

Full-Conditionals

The fixed effects for both models have conjugate full-conditionals

[TABLE]

The conjugate full-conditionals for the weights of the GRF are

[TABLE]

The log full-conditionals for $\kappa$ and $\xi$ of the GRF model are as seen in equation (44) of Appendix A.6. The conjugate full-conditionals for the weights of the LMA are

[TABLE]

The log full-conditionals for the parameters $\kappa,\tau$ and $\lambda$ are as seen in (44) of Appendix A.6. For model fitting we used normal proposals for all parameters. $\tau$ and $\lambda$ were jointly proposed for the LMA model, while $\kappa$ and $\xi$ were jointly proposed for the GRF.