Bayesian spatially varying coefficient models in the spBayes R package

TL;DR

This paper introduces new spatially varying coefficient modeling capabilities in the spBayes R package, enabling more flexible spatial regression analysis with improved computational efficiency and parallelization.

Contribution

The paper extends spBayes to include spatially varying coefficients with Gaussian process effects, incorporating efficient MCMC algorithms and parallelization features.

Findings

Enhanced modeling flexibility for spatial data analysis.

Improved computational efficiency through OpenMP parallelization.

Successful application to PM10 pollution data in Europe.

Abstract

This paper describes and illustrates new functionality for fitting spatially varying coefficients models in the spBayes (version 0.4-2) R package. The new spSVC function uses a computationally efficient Markov chain Monte Carlo algorithm and extends current spBayes functions, that fit only space-varying intercept regression models, to fit independent or multivariate Gaussian process random effects for any set of columns in the regression design matrix. Newly added OpenMP parallelization options for spSVC are discussed and illustrated, as well as helper functions for joint and point-wise prediction and model fit diagnostics. The utility of the proposed models is illustrated using a PM10 analysis over central Europe.

Figures

Equations

$$y(\text{s})=({\beta }_1+{\delta }_1{w}_1(\text{s}))+\sum _{j=2}^p{x}_j(\text{s})\left\{{\beta }_j+{\delta }_j{w}_j(\text{s})\right\}+ϵ(\text{s}),$$

$$p(\mathbf{\tau },\mathbf{\theta })\times N(\mathbf{\beta }|{\mathbf{\mu }}_{\beta },{\mathbf{\Sigma }}_{\beta })\times N(\text{w}|\text{0},{\text{K}}_{\theta })\times N(\text{y}|\text{X}\mathbf{\beta }+\text{Z}\text{w},{\text{D}}_{\tau }),$$

$$p({\tau }^2,\mathbf{\theta })=IG({\tau }^2|a_{\tau },b_{\tau })\times IW(\text{A}{\text{A}}^{\top }|r_a,{\text{S}}_a)\times \prod _{j=1}^{r}p({\mathbf{\varphi }}_j),$$

$$\log p(\mathbf{\theta }|\text{y})=\text{mbox}const+\log p(\mathbf{\theta })-\frac{1}{2}\log |{\mathbf{\Sigma }}_{y|\theta }|-\frac{1}{2}Q(\mathbf{\theta }),$$

$$\log p(\mathbf{\theta }|\text{y})=\text{mbox}constant+\log p(\mathbf{\theta })-\frac{1}{2}\log |{\text{X}}^{\top }{\mathbf{\Sigma }}_{y|\beta ,\theta }^{-1}\text{X}|-\frac{1}{2}\log |{\mathbf{\Sigma }}_{y|\beta ,\theta }|-\frac{1}{2}Q(\mathbf{\theta }),$$

$$\log p(\mathbf{\theta })-\sum _{i=1}^p\log w_{i,i}-\sum _{i=1}^n\log l_{i,i}-\frac{1}{2}({\text{v}}^{\top }\text{v}-{\tilde{\text{b}}}^{\top }\tilde{\text{b}}),$$

$$\text{b}={\mathbf{\Sigma }}_{\beta }^{-1}{\mathbf{\mu }}_{\beta }+{\text{X}}^{\top }{\mathbf{\Sigma }}_{y|\beta ,\theta }^{-1}\text{y}\text{mbox}and\text{B}={\left({\mathbf{\Sigma }}_{\beta }^{-1}+{\text{X}}^{\top }{\mathbf{\Sigma }}_{y|\beta ,\theta }^{-1}\text{X}\right)}^{-1}.$$

$${\mathbf{\beta }}^{(k)}=\text{mbox}trsolve\left({\text{L}}_B^{\top },\text{mbox}trsolve({\text{L}}_B,\text{b})\right)+\text{mbox}trsolve({\text{L}}_B^{\top },\text{z}),$$

$$\text{b}={\text{Z}}^{\top }{\mathbf{\Sigma }}_{y|w,\theta }^{-1}(\text{y}-\text{X}{\mathbf{\mu }}_{\beta })\text{mbox}and\text{B}={\left({\text{K}}_{\theta }^{-1}+{\text{Z}}^{\top }{\mathbf{\Sigma }}_{y|w,\theta }^{-1}\text{Z}\right)}^{-1}.$$

$${\left({\text{K}}_{\theta }^{-1}+{\text{G}}_{\theta }^{-1}\right)}^{-1}={\text{G}}_{\theta }-{\text{G}}_{\theta }{\left({\text{K}}_{\theta }+{\text{G}}_{\theta }\right)}^{-1}{\text{G}}_{\theta }$$

$$\displaystyle\begin{bmatrix}\textbf{y}\\ \textbf{y}_{0}\end{bmatrix}\bigg{|}\boldsymbol{\beta},{\boldsymbol{\theta}}\sim N_{n_{0}+n}\left(\begin{bmatrix}\textbf{X}\\ \textbf{X}_{0}\end{bmatrix}\boldsymbol{\beta},\;\begin{bmatrix}\textbf{C}_{11}({\boldsymbol{\theta}})&\textbf{C}_{12}({\boldsymbol{\theta}})\\ \textbf{C}_{12}({\boldsymbol{\theta}})^{\top}&\textbf{C}_{22}({\boldsymbol{\theta}})\end{bmatrix}\right)\;,$$

$${\mathbf{\mu }}_p$$

$$y(\text{s})={\beta }_0+w_0(\text{s})+a(\text{s})\left\{{\beta }_a+w_a(\text{s})\right\}+b(\text{s})\left\{{\beta }_b+w_b(\text{s})\right\}+ϵ(\text{s}).$$

$$\left(\begin{array}{ccc}1& 0& 0\\ -1& 1& 0\\ 0& 1& 0.1\end{array}\right)\left(\begin{array}{ccc}\text{exp}(-4d_{i,j})& 0& 0\\ 0& \text{exp}(-6d_{i,j})& 0\\ 0& 0& \text{exp}(-6d_{i,j})\end{array}\right)\left(\begin{array}{ccc}1& -1& 0\\ 0& 1& 1\\ 0& 0& 0.1\end{array}\right),$$

$${\text{K}}_{\theta }=\text{A}{\text{A}}^{\top }\left(\begin{array}{ccc}1& 0& 0\\ -1& 1& 0\\ 0& 1& 0.1\end{array}\right)\left(\begin{array}{ccc}1& -1& 0\\ 0& 1& 1\\ 0& 0& 0.1\end{array}\right)=\left(\begin{array}{ccc}1& -1& 0\\ -1& 2& 1\\ 0& 1& 1.01\end{array}\right).$$

$${\text{PM}}_{10}(\text{s})={\beta }_0+w_0(\text{s})+\text{CTM}(\text{s})\left\{{\beta }_{CTM}+w_{CTM}(\text{s})\right\}+ϵ(\text{s}).$$

Full text

Bayesian spatially varying coefficient models in the \pkgspBayes \proglangR package

Andrew O. Finley

Michigan State University &Sudipto Banerjee

University of California

  

Los Angeles

finleya@msu.edu

sudipto@ucla.edu

\Plainauthor

Andrew O. Finley, Sudipto Banerjee

\PlaintitleBayesian spatially varying coefficient models in the spBayes R package \ShorttitleSpatially varying coefficient in \proglangR \Abstract This paper describes and illustrates new functionality for fitting spatially varying coefficients models in the \pkgspBayes (version 0.4-2) \proglangR package. The new \codespSVC function uses a computationally efficient Markov chain Monte Carlo algorithm and extends current \pkgspBayes functions, that fit only space-varying intercept regression models, to fit independent or multivariate Gaussian process random effects for any set of columns in the regression design matrix. Newly added \proglangOpenMP parallelization options for \codespSVC are discussed and illustrated, as well as helper functions for joint and point-wise prediction and model fit diagnostics. The utility of the proposed models is illustrated using a PM10 analysis over central Europe.

\KeywordsMCMC, multivariate Gaussian Process, kriging, \proglangR \PlainkeywordsMCMC, multivariate Gaussian Process, kriging, R \Address Andrew Finley

Department of Forestry

Michigan State University

Natural Resources Building

480 Wilson Road, Room 126

East Lansing, MI 48824-6402

E-mail:

URL: https://www.finley-lab.com

Sudipto Banerjee

Fielding School of Public Health

University of California, Los Angeles

650 Charles E. Young Dr. South

Los Angeles, CA 90095-1772

E-mail:

URL: https://ph.ucla.edu/faculty/banerjee

1 Introduction

In this paper we describe and illustrate extended functionality of a recent reformulation and rewrite of core functions in the \pkgspBayes (Finley et al., 2015) \proglangR (R Core Team, 2018) package. The \pkgspBayes package provides a suite of univariate and multivariate regression models for both Gaussian and non-Gaussian outcomes that are spatially indexed. There are, by now, many \proglangR packages that provide similar functionality. A recent read of the “Analysis of Spatial Data” CRAN Task View (Bivand, 2019) yielded $\sim$46 packages listed for geostatistical analysis—and this is not an exhaustive accounting of packages available for such analyses. Finley et al. (2015) focused on laying out computationally efficient and flexible MCMC algorithms for estimating an array of spatio-temporal Gaussian process (GP) models. However, while the proposed sampling algorithms were quite general, only a narrow set of models were implemented in \pkgspBayes. Specifically, users could only specify univariate or multivariate GPs on model intercepts. Now, the addition of the \codespSVC function to \pkgspBayes (version 0.4-2 available on CRAN 3/7/2019) aims to provide additional user options for placing univariate or multivariate GPs on any set of model regression coefficients.

Such functionality is not unique, there are several \proglangR packages capable of fitting spatially varying coefficient (SVC) models. Some are specifically designed to work with spatial or spatio-tempral data and others provide flexibility to allow coefficients to vary by some generic set of variables, which could be indexes in a coordinate system. Most of these packages employ some flavor of spline or kernel based regression method to allow varying impact of predictors. Hastie and Tibshirani (1993) and Fan and Zhang (2008) offer a general development of varying coefficient models and Gelfand et al. (2003a) provide treatment particular to spatial settings. Regarding implementation in \proglangR, the \pkgspgwr package (Bivand and Yu, 2017) implements geographically weighted regression as originally detailed in Fotheringham et al. (2002). Key spline-based package options include \pkgmgcv (Wood, 2017), \pkgsvcm (Heim, 2007), \pkgnp (Hayfield and Racine, 2008), and \pkgmboost (Hothorn et al., 2018). Bürgin and Ritschard (2017) recently developed a tree-based varying coefficient model (TVCM) algorithm and associated \pkgvcrpart package. The packages \pkgwalker (Helske, 2019; Vihola et al., 2017), \pkgspTimer Bakar and Sahu (2018), and \pkgspTDyn (Bakar et al., 2017, 2015a, 2015b) offer Bayesian time and space-time SVC models. Other Bayesian options include model development using more general software such as INLA (Rue et al., 2009; Lindgren and Rue, 2015; Bakka et al., 2018) and Stan (Carpenter et al., 2017; Stan Development Team, 2018), which can be called from their respective \proglangR packages.

The \codespSVC function offers Markov chain Monte Carlo (MCMC) based SVC inference using an efficient sampling algorithm. The algorithm’s efficiency derives from updates to only covariance parameters (i.e., regression coefficients and random effects are integrated out), computing parallelization, and use of tuned and/or multi-threaded matrix algebra libraries. Subsequent sections define the model and algorithm specifics, software features, and illustrative analyses of simulated and real data.

2 Models and software

Let $y(\textbf{s})$ be the dependent variable (response or outcome) at location s and consider the spatially varying regression model,

[TABLE]

where $x_{j}(\textbf{s})$, for each $j=2,\ldots,p$ with $p\geq 1$, is the known value of a predictor at location s, $\beta_{j}$ is the regression coefficient corresponding to $x_{j}(\textbf{s})$, $\beta_{1}$ is an intercept, and $\epsilon(\textbf{s})$ is a Gaussian measurement error process independently distributed for each s. The quantities $w_{1}(\textbf{s})$ and $w_{j}(\textbf{s})$ are spatial processes corresponding to the intercept and predictors, thereby yielding a spatially varying regression model. We further accommodate the possibility that not all the predictors will have spatially varying impact on the outcome. Thus, $\delta$’s in (1) are binary indicators assuming the value $1$ if the associated predictor has a spatially varying regression coefficient and [math] otherwise. For later convenience, when the respective $\delta=1$ we define $\tilde{\beta}_{1}(\textbf{s})=\beta_{1}+\delta_{1}w_{1}(\textbf{s})$ and $\tilde{\beta}_{j}(\textbf{s})=\beta_{j}+\delta_{j}w_{j}(\textbf{s})$ as the space-varying regression coefficients.

Let ${\cal S}=\{\textbf{s}_{1},\textbf{s}_{2},\ldots,\textbf{s}_{n}\}$ be the set of spatial locations from which $y(\textbf{s}_{i})$ and the predictors have been observed. Let w be the $nr\times 1$ vector obtained by stacking up $\textbf{w}(\textbf{s}_{i})$’s, where each $\textbf{w}(\textbf{s}_{i})$ is an $r\times 1$ vector with $j$-th entry $w_{j}(\textbf{s}_{i})$, $j=1,2,\ldots,r$ and $i=1,2,\ldots,n$. We treat $\textbf{w}(\textbf{s})$ as a multivariate Gaussian process (see, e.g., Banerjee et al., 2014) so the matrix $\textbf{K}_{\theta}$ is an $nr\times nr$ spatial covariance matrix constructed as a block matrix with $(i,j)$-th block obtained from the $r\times r$ cross-covariance matrix $\textbf{K}_{\theta}(\textbf{s}_{i},\textbf{s}_{j})$ specifying the multivariate spatial process $\textbf{w}(\textbf{s})$. In addition, $\boldsymbol{\beta}$ is the $p\times 1$ regression coefficient corresponding to X, and ${\boldsymbol{\theta}}$ and ${\boldsymbol{\tau}}$ are the parameters in $\textbf{K}_{\theta}$ and $\textbf{D}_{\tau}$, respectively.

Consider the Bayesian hierarchical model built from (1),

[TABLE]

where y is $n\times 1$ with $i$-th element $y(\textbf{s}_{i})$, X is $n\times p$ with the first column 1 and the remaining $p-1$ columns corresponding to the predictors $x_{j}(\textbf{s}_{i})$ in (1). The matrix Z is $n\times nr$, where $r\leq p$, with precisely those columns of X which have $\delta_{j}=1$.

We will offer users the option to scale and center the matrix X. Note that Z is constructed from X, and thus, for a scaled and centered X, the predictors used in Z will also be scaled and centered. Scaling and centering often improves numerical stability and provides more robust estimation of spatially varying regression models (see, e.g., Gelfand et al., 2003b).

Some further specifications are in order. In \codespSVC we will fix $\textbf{D}_{\tau}=\tau^{2}\textbf{I}_{n}$ so ${\boldsymbol{\tau}}=\{\tau^{2}\}$ is the scalar quantity representing the measurement error variance or “nugget” in geostatistics. The cross-covariance matrix $\textbf{K}_{\theta}(\textbf{s},\textbf{t})$, where t is a generic location, will most generally be modeled using the Linear Model of Coregionalization (LMC). Here, we will model $\textbf{K}_{\theta}(\textbf{s},\textbf{t})=\textbf{A}\Gamma(\textbf{s},\textbf{t})\textbf{A}^{\top}$, where A is an $r\times r$ lower triangular matrix and $\Gamma(\textbf{s},\textbf{t})$ is a diagonal matrix with $\rho_{j}(\textbf{s},\textbf{t})$ being the $j$-th diagonal element, where $\rho_{j}(\textbf{s},\textbf{t})$ is a spatial correlation function with parameters specific to $w_{j}(\textbf{s})$. Here, ${\boldsymbol{\theta}}$ in (2) corresponds to $\{\textbf{A},\{\boldsymbol{\phi}_{j}\}_{j=1}^{r}\}$ where each $\boldsymbol{\phi}_{j}$ is a collection of parameters in the spatial correlation function. For example, with the Matérn covariance function each $\boldsymbol{\phi}_{j}$ comprises a spatial decay parameter and a smoothness parameter.

The covariance structure for $\textbf{w}(\textbf{s})$ within any location s is captured by $\textbf{A}\textbf{A}^{\top}$, which identifies with the Cholesky decomposition for $\mbox{var}\{\textbf{w}(\textbf{s})\}$. In general, we will specify priors as

[TABLE]

where $IG$ is inverse-Gamma, $IW$ is inverse-Wishart, and each $p(\boldsymbol{\phi}_{j})$ can be one of the several distributions provided by \pkgspBayes. Another particular choice offered by \codespSVC specifies $\textbf{A}=\mbox{diag}(\sigma_{1},\ldots,\sigma_{r})$, so that $\textbf{K}_{\theta}(\textbf{s},\textbf{t})$ is diagonal with entries $\sigma_{j}^{2}\rho_{j}(\textbf{s},\textbf{t})$, in which case we assume $IG(\sigma_{j}^{2}\,|\,a_{\sigma},b_{\sigma})$ for $j=1,2,\ldots,r$. Choices for $\rho_{j}$ include any of the standard correlation functions offered by \pkgspBayes.

2.1 Parameter estimation and computational considerations

Bayesian inference for (1) involves sampling the parameters ${\boldsymbol{\theta}}$, $\boldsymbol{\beta}$ and w from their marginal posterior distributions. Such sampling algorithms require expensive operations on dense matrices such as decomposition and multiplication. Therefore, as we have outlined below, care is needed to use efficient numerical algorithms such as Cholesky factorizations, working with triangular systems, and avoiding redundant operations.

2.1.1 Sampling the process parameters

Sampling from (2) employs MCMC methods, in particular Gibbs sampling and random walk Metropolis steps (e.g., Robert and Casella 2004). For faster convergence, we integrate out $\boldsymbol{\beta}$ and w from the model and first sample from $\displaystyle p({\boldsymbol{\theta}}\,|\,\textbf{y})\propto p({\boldsymbol{\theta}})\times N(\textbf{y}\,|\,\textbf{X}\boldsymbol{\mu}_{\beta},\boldsymbol{\Sigma}_{y\,|\,\theta})$, where $\boldsymbol{\Sigma}_{y\,|\,\theta}=\textbf{X}\boldsymbol{\Sigma}_{\beta}\textbf{X}^{\top}+\textbf{Z}\textbf{K}_{\theta}\textbf{Z}^{\top}+\textbf{D}_{\tau}$. This matrix needs to be constructed for every update of ${\boldsymbol{\theta}}$. $\textbf{D}_{\tau}$ is diagonal and $\textbf{X}\boldsymbol{\Sigma}_{\beta}\textbf{X}^{\top}$ is fixed, so the computation involves the matrix $\textbf{Z}\textbf{K}_{\theta}\textbf{Z}^{\top}$ which requires $rn^{2}$ flops (floating point operations).

We adopt a random-walk Metropolis step with a multivariate normal proposal (same dimension as there are parameters in ${\boldsymbol{\theta}}$) after transforming parameters to have support over the entire real line. This involves evaluating

[TABLE]

where $Q({\boldsymbol{\theta}})=(\textbf{y}-\textbf{X}\boldsymbol{\mu}_{\beta})^{\top}\boldsymbol{\Sigma}^{-1}_{y\,|\,\theta}(\textbf{y}-\textbf{X}\boldsymbol{\mu}_{\beta})$. Generally, we compute $\textbf{L}=\mbox{{chol}}(\boldsymbol{\Sigma}_{y\,|\,\theta})$, where $\mbox{{chol}}(\boldsymbol{\Sigma}_{y\,|\,\theta})$ returns the lower-triangular Cholesky factor L of $\boldsymbol{\Sigma}_{y\,|\,{\boldsymbol{\theta}}}$. This involves $O(n^{3}/3)$ flops. Next, we obtain $\textbf{u}=\mbox{{trsolve}}(\textbf{L},\textbf{y}-\textbf{X}\boldsymbol{\mu}_{\beta})$, which solves the triangular system $\textbf{L}\textbf{u}=\textbf{y}-\textbf{X}\boldsymbol{\mu}_{\beta}$. This involves $O(n^{2})$ flops and $Q({\boldsymbol{\theta}})=\textbf{u}^{\top}\textbf{u}$ requires another $2n$ flops. The log-determinant in (4) is evaluated as $2\sum_{i=1}^{n}\log l_{i,i}$, where $l_{i,i}$ are the diagonal entries in L. Since L has already been obtained, the log-determinant requires another $n$ steps. Therefore, the Cholesky factorization dominates the work and computing (4) is achieved in $O(n^{3})$ flops.

If $\boldsymbol{\beta}$ is flat, i.e., $\boldsymbol{\Sigma}_{\beta}^{-1}=\textbf{O}$, the analogue of distribution (4) is

[TABLE]

where $\boldsymbol{\Sigma}_{y\,|\,\beta,\theta}=\textbf{Z}\textbf{K}_{\theta}\textbf{Z}^{\top}+\textbf{D}_{\tau}$ and $Q({\boldsymbol{\theta}})=\textbf{y}^{\top}\boldsymbol{\Sigma}_{y\,|\,\beta,\theta}^{-1}\textbf{y}-\textbf{b}^{\top}(\textbf{X}^{\top}\boldsymbol{\Sigma}_{y\,|\,\beta,\theta}^{-1}\textbf{X})^{-1}\textbf{b}$ and $\textbf{b}=\textbf{X}^{\top}\boldsymbol{\Sigma}_{y\,|\,\beta,\theta}^{-1}\textbf{y}$. Computations proceed similar to the above. We first evaluate $\textbf{L}=\mbox{{chol}}(\boldsymbol{\Sigma}_{y\,|\,\beta,\theta})$ and then obtain $[\textbf{v}:\textbf{U}]=\mbox{{trsolve}}(\textbf{L},[\textbf{y}:\textbf{X}])$, so $\textbf{L}\textbf{v}=\textbf{y}$ and $\textbf{L}\textbf{U}=\textbf{X}$. Next, we evaluate $\textbf{W}=\mbox{{chol}}(\textbf{U}^{\top}\textbf{U})$, $\textbf{b}=\textbf{U}^{\top}\textbf{v}$ and solve $\tilde{\textbf{b}}=\mbox{{trsolve}}(\textbf{W},\textbf{b})$. Finally, (5) is evaluated as

[TABLE]

where $w_{i,i}$’s and $l_{i,i}$’s are the diagonal elements in W and L respectively. The number of flops is again of cubic order in $n$.

Importantly, our strategy above avoids computing inverses. We use Cholesky factorizations and solve only triangular systems. If $n$ is not large, say \sim$$10^{2}, this strategy is feasible. As described in Section 2.2 and illustrated in Section 3, use of efficient and parallelized numerical linear algebra routines yields substantial gains in computing time.

2.1.2 Sampling the slope and the random effects

Once we have obtained marginal posterior samples ${\boldsymbol{\theta}}$ from $p({\boldsymbol{\theta}}\,|\,\textbf{y})$, we can draw posterior samples of $\boldsymbol{\beta}$ and w using composition sampling. Suppose $\{{\boldsymbol{\theta}}^{(1)},{\boldsymbol{\theta}}^{(2)},\ldots,{\boldsymbol{\theta}}^{(M)}\}$ are $M$ samples from $p({\boldsymbol{\theta}}\,|\,\textbf{y})$. Drawing $\boldsymbol{\beta}^{(k)}\sim p(\boldsymbol{\beta}\,|\,{\boldsymbol{\theta}}^{(k)},\textbf{y})$ and $\textbf{w}^{(k)}\sim p(\textbf{w}\,|\,{\boldsymbol{\theta}}^{(k)},\textbf{y})$ for $k=1,2,\ldots M$ results in $M$ samples from $p(\boldsymbol{\beta}\,|\,\textbf{y})$ and $p(\textbf{w}\,|\,\textbf{y})$ respectively. Only the samples of ${\boldsymbol{\theta}}$ obtained after convergence (i.e., post burn-in) of the MCMC algorithm need to be stored.

To elucidate further, note that $\boldsymbol{\beta}\,|\,{\boldsymbol{\theta}},\textbf{y}\sim N_{p}(\textbf{B}\textbf{b},\textbf{B})$ with mean B****b and variance-covariance matrix B, where

[TABLE]

For each $k=1,2,\ldots,M$, we compute B and b at the current value ${\boldsymbol{\theta}}^{(k)}$ and draw $\boldsymbol{\beta}^{(k)}\sim N_{p}(\textbf{B}\textbf{b},\textbf{B})$. This is achieved by computing $\textbf{b}=\boldsymbol{\Sigma}_{\beta}^{-1}\boldsymbol{\mu}_{\beta}+\textbf{U}^{\top}\textbf{v}$, where $\textbf{L}=\mbox{{chol}}(\boldsymbol{\Sigma}_{y\,|\,\beta,\theta^{(k)}})$ and $[\textbf{v}:\textbf{U}]=\mbox{{trsolve}}(\textbf{L},[\textbf{y}:\textbf{X}])$. Next, we generate $p$ independent standard normal variables, collect them into z and set

[TABLE]

where $\displaystyle\textbf{L}_{B}=\mbox{{chol}}\left(\boldsymbol{\Sigma}_{\beta}^{-1}+\textbf{U}^{\top}\textbf{U}\right)$. This completes the $k$-th iteration. After $M$ iterations, we obtain $\{\boldsymbol{\beta}^{(1)},\boldsymbol{\beta}^{(2)},\ldots,\boldsymbol{\beta}^{(M)}\}$, which are samples from $p(\boldsymbol{\beta}\,|\,\textbf{y})$.

Mapping point or interval estimates of spatial random effects is often helpful in identifying missing regressors and/or building a better understanding of model adequacy. $\boldsymbol{\Sigma}_{y\,|\,w,\theta}=\textbf{X}\boldsymbol{\Sigma}_{\beta}\textbf{X}^{\top}+\textbf{D}_{\tau}$ and note that $\textbf{w}\,|\,{\boldsymbol{\theta}},\textbf{y}\sim N(\textbf{B}\textbf{b},\textbf{B})$, where

[TABLE]

The vector b here is computed analogously as for $\boldsymbol{\beta}$. For each $k=1,2,\ldots,M$ we now evaluate $\textbf{L}=\mbox{{chol}}(\boldsymbol{\Sigma}_{y\,|\,w,\theta^{(k)}})$, $[\textbf{v}:\textbf{U}]=\mbox{{trsolve}}(\textbf{L},[\textbf{y}-\textbf{X}\boldsymbol{\mu}_{\beta}:\textbf{Z}({\boldsymbol{\theta}}^{(k)})])$ and set $\textbf{b}=\textbf{U}({\boldsymbol{\theta}}^{(k)})^{\top}\textbf{v}$. For computing B, one could proceed as for $\boldsymbol{\beta}$ but that would involve $\mbox{{chol}}(\textbf{K}({\boldsymbol{\theta}}))$, which may become numerically unstable for certain covariance functions (e.g., the Gaussian or the Matérn with large $\nu$). For robust software performance we define $\textbf{G}_{\theta}^{-1}=\textbf{Z}^{\prime}\boldsymbol{\Sigma}_{y\,|\,w,\theta}^{-1}\textbf{Z}$ and utilize the identity (Henderson and Searle, 1981)

[TABLE]

to devise a numerically stable algorithm. For each $k=1,2,\ldots,M$, we evaluate $\textbf{L}=\mbox{{chol}}(\textbf{K}_{\theta}^{(k)}+\textbf{G}_{\theta}^{(k)})$, $\textbf{W}=\mbox{{trsolve}}(\textbf{L},\textbf{G}_{\theta}^{(k)})$ and $\textbf{L}_{B}=\mbox{{chol}}(\textbf{G}_{\theta}^{(k)}-\textbf{W}^{\top}\textbf{W})$. If z is a $r\times 1$ vector of independent standard normal variables, then we set $\displaystyle\textbf{w}^{(k)}=\textbf{L}_{B}\textbf{L}_{B}^{\top}\textbf{b}+\textbf{L}_{B}\textbf{z}$. The resulting $\{\textbf{w}^{(1)},\textbf{w}^{(2)},\ldots,\textbf{w}^{(M)}\}$ are samples from $p(\textbf{w}\,|\,\textbf{y})$.

2.1.3 Spatial predictions

To predict a random $n_{0}\times 1$ vector $\textbf{y}_{0}$ associated with a $n_{0}\times p$ matrix of predictors, $\textbf{X}_{0}$, we assume that

[TABLE]

where $\textbf{C}_{11}({\boldsymbol{\theta}})=\boldsymbol{\Sigma}_{y\,|\,\beta,\theta}$, $\textbf{C}_{12}({\boldsymbol{\theta}})$ is the $n\times n_{0}$ cross-covariance matrix between y and $\textbf{y}_{0}$, and $\textbf{C}_{22}({\boldsymbol{\theta}})$ is the variance-covariance matrix for $\textbf{y}_{0}$. A valid joint distribution will supply a conditional distribution $p(\textbf{y}_{0}\,|\,\textbf{y},\boldsymbol{\beta},{\boldsymbol{\theta}})$, which is normal with mean and variance

[TABLE]

Bayesian prediction proceeds by sampling from the posterior predictive distribution $\displaystyle p(\textbf{y}_{0}\,|\,\textbf{y})=\int p(\textbf{y}_{0}\,|\,\textbf{y},\boldsymbol{\beta},{\boldsymbol{\theta}})p(\boldsymbol{\beta},{\boldsymbol{\theta}}\,|\,\textbf{y})d\boldsymbol{\beta}d{\boldsymbol{\theta}}$. For each posterior sample of $\{\boldsymbol{\beta},{\boldsymbol{\theta}}\}$, we draw a corresponding $\textbf{y}_{0}\sim N(\boldsymbol{\mu}_{p},\boldsymbol{\Sigma}_{p})$. This produces samples from the posterior predictive distribution.

The posterior predictive computations involve only the retained MCMC samples after convergence. For any posterior sample $\{\boldsymbol{\beta}^{(k)},{\boldsymbol{\theta}}^{(k)}\}$, we solve $[\textbf{u}:\textbf{V}]=\mbox{{trsolve}}(\textbf{L},[\textbf{y}-\textbf{X}\boldsymbol{\beta}^{(k)}:\textbf{C}_{12}({\boldsymbol{\theta}}^{(k)})])$, where $\textbf{L}=\mbox{{chol}}(\textbf{C}_{11}({\boldsymbol{\theta}}^{(k)}))$. Next, we set $\boldsymbol{\mu}_{p}^{(k)}=\textbf{X}_{0}\boldsymbol{\beta}^{(k)}+\textbf{V}^{\top}\textbf{u}$ and $\boldsymbol{\Sigma}_{p}^{(k)}=\textbf{C}_{22}({\boldsymbol{\theta}}^{(k)})-\textbf{V}^{\top}\textbf{V}$ and draw $\textbf{y}_{0}^{(k)}\sim N(\boldsymbol{\mu}_{p}^{(k)},\boldsymbol{\Sigma}_{p}^{(k)})$.

Updating $\textbf{y}_{0}^{(k)}$’s requires Cholesky factorization of $\boldsymbol{\Sigma}_{p}$, which is $n_{0}\times n_{0}$ and can be expensive if $n_{0}$ is large. In most practical settings, it is sufficient to take $n_{0}=1$ and perform point-wise predictions.

2.2 Software features

The \codespSVC function accommodates the \codespLM function in \pkgspBayes and offers additional user options to simplify analysis and inference. The list below highlights some of these new options.

Any set of predictors can receive either independent univariate GPs or a multivariate GP. 2. 2.

Prediction can be done by sampling from either the joint or point-wise (marginal) posterior predictive distribution. 3. 3.

\proglang

openMP (Dagum and Menon, 1998) support is available via the \coden.omp.threads argument for parameter estimation, composition sampling, model fit diagnostics, and prediction functions. 4. 4.

Matrix operation parallelization is available via multi-threaded implementations of Basic Linear Algebra Subprograms (BLAS; www.netlib.org/blas) and Linear Algebra Package (LAPACK; www.netlib.org/lapack). 5. 5.

Coordinate system used to index observed and prediction locations can be of arbitrary dimension—users were previously restricted to using 2-dimensional systems. 6. 6.

Univariate and multivariate random effect samples and space-varying coefficients are returned as lists with element names corresponding to the given predictor.

3 Illustrations

We consider two analyses to illustrate key features of \codespSVC along with supporting functions. The first analysis is of a simulated dataset and second is of an air pollution dataset that was previously analyzed in Hamm et al. (2015) and Datta et al. (2016).

3.1 Analysis of simulated data

The simulated data \codemvSVCData is available in \pkgspBayes and comprises $n$=500 observations distributed within a 2-dimensional unit square spatial domain. At generic location s the outcome was generated following

[TABLE]

The predictors $a(\textbf{s})$ and $b(\textbf{s})$ were drawn from independent normal distributions with mean zero and variance one. The regression coefficients $\beta_{0}$, $\beta_{a}$, and $\beta_{b}$ equaled 1, 10, and -10, respectively. The $r$=3 spatial random effects associated with the intercept and predictors were generated from a non-separable multivariate GP. The cross-covariance function used to construct the $(i,j)$-th $r\times r$ block in the multivariate GP’s $nr\times nr$ covariance matrix, i.e., $\textbf{K}_{\theta}(\textbf{s}_{i},\textbf{s}_{j})=\textbf{A}\Gamma(\textbf{s}_{i},\textbf{s}_{j})\textbf{A}^{\top}$, was

[TABLE]

where $d_{i,j}$ is the euclidean distance between location $\textbf{s}_{i}$ and $\textbf{s}_{j}$, and diagonal elements of $\Gamma(\textbf{s}_{i},\textbf{s}_{j})$ are the exponential correlation function $\text{exp}(-\phi_{k}d_{i,j})$ for $k=1,2,\ldots,r$. Figure 2(a)-2(c) display the realizations of $\textbf{w}_{0}$, $\textbf{w}_{a}$, and $\textbf{w}_{b}$. The residual term $\epsilon(\textbf{s})$ was simulated from a Normal distribution with mean zero and variance $\tau^{2}=0.1$.

The code below specifies the model covariance parameters’ prior distributions, and MCMC sampler starting and Metropolis proposal variance values. Here, we use a Uniform prior for the spatial decay parameters each with support from 1 to 10. The prior for the cross-covaraince matrix is an IW with degrees of freedom $r$ and identity scale matrix. The prior for the measurement error (or nugget variance) follows an IG with shape 2 and scale 1.

The parameter priors, starting values, and Metropolis sampler proposal variances are passed to \codespSVC via the \codepriors, \codestarting, and \codetuning arguments, respectively. The proposed model is specified via the \codeformula argument using syntax like that used in base \proglangR’s \codelm, with the addition of the \codesvc.cols argument that accepts a vector of either integer indexes or character names to indicate the space-varying design matrix columns (i.e., the columns of X with $\delta_{j}=1$). For example, in the call to \codespSVC below, the vector passed to \codesvc.cols indicates we want the intercept and columns labeled \codea and \codeb to follow a multivariate GP (or, equivalently, one could use the argument value \codec(1,2,3)).

\MakeFramed data(SVCMvData.dat)

r <- 3

n.ltr <- r*(r+1)/2

priors <- list("phi.Unif"=list(rep(1,r), rep(10,r)), "K.IW"=list(r, diag(rep(1,r))), "tau.sq.IG"=c(2, 1))

starting <- list("phi"=rep(3/0.5,r), "A"=c(1,0,0,1,0,1), "tau.sq"=1)

tuning <- list("phi"=rep(0.1,r), "A"=rep(0.01, n.ltr), "tau.sq"=0.01)

sim.m <- spSVC(y~a+b, coords=c("x.coords","y.coords"), data=SVCMvData.dat, starting=starting, svc.cols=c("(Intercept)","a","b"), tuning=tuning, priors=priors, cov.model="exponential", n.samples=10000, n.report=5000, n.omp.threads=4)

---

ΨGeneral model description

Model fit with 200 observations.

Number of covariates 3.

Number of space varying covariates 3.

Using the exponential spatial correlation model.

Number of MCMC samples 10000.

Priors and hyperpriors: Ψbeta flat. ΨK IW hyperpriors: Ψdf: 3.00000 ΨS: Ψ1.000Ψ0.000Ψ0.000Ψ Ψ0.000Ψ1.000Ψ0.000Ψ Ψ0.000Ψ0.000Ψ1.000Ψ

Ψphi Unif lower bound hyperpriors:Ψ1.000Ψ1.000Ψ1.000Ψ Ψphi Unif upper bound hyperpriors:Ψ10.000Ψ10.000Ψ10.000Ψ

Ψtau.sq IG hyperpriors shape=2.00000 and scale=1.00000

Source compiled with OpenMP, posterior sampling is using 4 thread(s).

ΨΨSampling

Sampled: 5000 of 10000, 50.00% Report interval Metrop. Acceptance rate: 34.64% Overall Metrop. Acceptance rate: 34.64%

Sampled: 10000 of 10000, 100.00% Report interval Metrop. Acceptance rate: 34.24% Overall Metrop. Acceptance rate: 34.44%

\endMakeFramed

As described in Section 2, \codespSVC computes and returns MCMC samples for only model covariance parameters. If \codeverbose=TRUE, basic model specifications are written to the terminal followed by updates on the sampler’s progress and Metropolis algorithm acceptance rate. The sampler progress report interval is controlled using the \coden.report argument. One should adjust the Metropolis sampler proposal variances to achieve an acceptance rate between $\sim$30-50% (see, e.g., Gelman et al., 2013, for model fitting best practices). If it proves difficult to maintain an acceptable acceptance rate, the \codeamcmc argument can be added to invoke an adaptive MCMC algorithm (Roberts and Rosenthal, 2009) that automatically adjusts the tuning to achieve a target acceptance rate (see the manual page for more details).

The \coden.omp.threads argument in \codespSVC call above requests that key \codefor loops within a given MCMC iteration use 4 threads via \proglangopenMP (Dagum and Menon, 1998). If the user’s \proglangR is set up to use a parallelized version of BLAS then \coden.omp.threads will also control the number of threads in some LAPACK matrix operations. Such parallelization can greatly reduce the sampler’s runtime.

The computer used to conduct this analysis has an Intel(R) Core(TM) i7-8550U CPU @ 1.80GHz chip with 4 cores and \proglangR compiled with \proglangopenMP, as confirmed in the “General model description” printed after calling \codespSVC, which notes Source compiled with OpenMP, posterior sampling is using 4 thread(s). \codespSVC will throw a warning if \proglangR was not compiled with \proglangopenMP support and \coden.omp.threads is set to a value greater than 1. In addition to \proglangopenMP support, the current implementation of \proglangR uses openBLAS (Zhang, 2016) which is a version of BLAS capable of exploiting multiple processors. Figure 1 shows the runtime needed to complete 10000 MCMC iterations across the number of available CPUs.

Following execution of \codespSVC, the \codesim.m object holds MCMC samples for covariance parameters (\codep.theta.samples) along with data and model fitting details. Using, possibly post burn-in and thinned, \codep.theta.samples, the \codespRecover function conducts composition sampling to generate samples from the regression coefficients $\boldsymbol{\beta}$ (\codep.beta.recover.samples), spatial random effects w (\codep.w.recover.samples), and model fitted values (\codep.y.samples). \codespRecover also returns the subset of \codep.theta.samples (\codep.theta.recover.samples) used in the composition sampling. Further, for convenience, \codespRecover returns samples of the space-varying regression coefficients $\tilde{\beta}(\textbf{s})$’s. \codespRecover appends these various composition sampling outputs to the \codespSVC input object, i.e., the \codesim.m object returned by \codespRecover below is identical to the \codesim.m object returned by \codespSVC except for the addition of the composition sampling results. In addition to providing posterior samples for all model parameters, a call to \codespRecover is necessary for subsequent prediction and model fit diagnostics, via \codespPredict and \codespDiag respectively. Like \codespSVC, \codespRecover takes advantage of multiple CPUs via \proglangopenMP when available.

\MakeFramed sim.m <- spRecover(sim.m, start=5000, thin=2, n.omp.threads=4, verbose=FALSE)

\endMakeFramed

\code

spSVC and \codespRecover return samples as \codecoda objects to simplify posterior summaries. Output below provides the post burn-in and thinned median with lower and upper 95% credible bounds for $\beta$’s, cross-covariance matrix used to construct $\textbf{K}_{\theta}$, spatial decays $\phi$’s, and $\tau^{2}$. These summaries show the model captures well the parameter values used to simulate the data.

\MakeFramed round(summary(sim.m$p.beta.recover.samples)$quantiles[,c(3,1,5)],2)

           50%   2.5% 97.5%

(Intercept) 0.20 -1.05 1.16 a 10.53 9.37 11.92 b -10.20 -11.08 -9.31

\endMakeFramed

Note, following the notation in Section 2 the cross-covaraince matrix used to simulated the data is

[TABLE]

The posterior summary of the covariance parameters is below. Observed versus estimated random effects are given Figure 2(d)-2(f).

\MakeFramed round(summary(sim.m$p.theta.recover.samples)$quantiles[,c(3,1,5)],2)

              50%  2.5% 97.5%

K[1,1] 1.11 0.58 2.78 K[2,1] -1.00 -2.49 -0.46 K[3,1] -0.06 -0.52 0.35 K[2,2] 1.89 1.17 3.60 K[3,2] 0.90 0.35 1.76 K[3,3] 1.14 0.63 2.03 tau.sq 0.18 0.10 0.33 phi.(Intercept) 4.16 1.57 8.65 phi.a 4.98 2.09 9.81 phi.b 6.13 1.55 9.89

\endMakeFramed

3.2 Analysis of air pollution data

Increases in human morbidity and mortality is a known outcome to airborne particulate matter (PM) exposure (Brunekreef and Holgate, 2002; Loomis et al., 2013; Hoek et al., 2013). In response, regulatory agencies have instigated monitor programs and regulate PM concentrations. One such regulation by the European Commission’s air quality standards limits PM10 (PM$<$10 $\mu$m in diameter) concentrations to 50 $\mu$g m*-3* average over 24 hours and 40 $\mu$g m*-3* over a year (European Commission, 2015).

Measurements made with instruments at monitoring stations are considered authoritative; however, these observations are often too sparse to deliver regional maps at sufficient resolution to assess progress with mitigation strategies and for monitoring compliance. One solution is to couple spatially sparse monitoring station observations with spatially complete chemistry transport model (CTM) output, (see, e.g., van de Kassteele and Stein, 2006; Denby et al., 2008; Candiani et al., 2013). In such settings, monitoring station observations serve as a regression model outcome with CTM output set as a predictor.

[TABLE]

This illustration draws on data and analyses presented in (Hamm et al., 2015; Datta et al., 2016). We consider April 6, 2010, PM10 measurements across central Europe with corresponding output from the LOTOS-EUROS (Schaap et al., 2008) CTM. Following (Hamm et al., 2015) we hypothesis a space-varying relationship between the PM10 measurements observed at monitoring stations and CTM output. In what follows, we compare fit metrics for three candidate models derived from (14): 1) a non-spatial regression; 2) space-varying intercept; 3) space-varying intercept and CTM output. Resulting model objects are called \codepm.1, \codepm.2, and \codepm.3, respectively. For brevity, code only for fitting \codepm.3 is shown. We then consider parameter estimates and associated plots of the spatial random effects from \codepm.3, followed by development of predictive maps of both the space-varying coefficients and PM10 prediction for a grid over the study area.

We begin by loading the data and separating it into a “model” set \codePM10.mod comprising locations where both PM10 measurements and CTM values are available, and a “prediction” set \codePM10.pred where only CTM values are available. Here too, we calculated the maximum distance between any two monitoring stations which will help with setting prior distributions for spatial decay parameters.

\MakeFramed data(PM10.dat)

PM10.mod <- PM10.dat[!is.na(PM10.dat$pm10.obs),] PM10.pred <- PM10.dat[**is.na**(PM10.dat$pm10.obs),]

d.max <- max(iDist(PM10.mod[,c("x.coord","y.coord")])) d.max #km

[1] 2929.193

\endMakeFramed

The code below specifies the model covariance parameters’ prior distributions, and MCMC sampler starting and Metropolis proposal variance values. Unlike the simulated data analysis, here we demonstrate placing independent GPs on the intercept and CTM predictor. This requires priors for a process specific spatial decay parameter $\phi$ and variance $\sigma^{2}$. We again use a Uniform prior for the process’ decay parameters that provides support for an effective spatial range between $\sim$ 3 and 2197 km, given an exponential covariance function. The two spatial variances and single observational variance $\tau^{2}$ each are assumed to follow an IG with shape 2 and scale 1. We center the IG’s on 1, because it is approximately equal to the residual variance from the first candidate model, i.e., the non-spatial regression. One should generally do careful exploratory data analysis to arrive at a robust set of prior distributions and hyperparameters

\MakeFramed r <- 2

priors <- list("phi.Unif"=list(rep(3/(0.75d.max), r), rep(3/(0.001d.max), r)), "sigma.sq.IG"=list(rep(2, r), rep(1, r)), "tau.sq.IG"=c(2, 1))

starting <- list("phi"=rep(3/(0.1*d.max), r), "sigma.sq"=rep(1, r), "tau.sq"=1)

tuning <- list("phi"=rep(0.1, r), "sigma.sq"=rep(0.05, r), "tau.sq"=0.1)

n.samples <- 10000

m.3 <- spSVC(pm10.obs ~ pm10.ctm, coords=c("x.coord","y.coord"), data=PM10.mod, starting=starting, svc.cols=c(1,2), tuning=tuning, priors=priors, cov.model="exponential", n.samples=n.samples, n.report=5000, n.omp.threads=4)

---

ΨGeneral model description

Model fit with 256 observations.

Number of covariates 2.

Number of space varying covariates 2.

Using the exponential spatial correlation model.

Number of MCMC samples 10000.

Priors and hyperpriors: Ψbeta flat. ΨDiag(K) IG hyperpriors ΨΨparameterΨshapeΨΨscale ΨΨK[1,1]ΨΨ2.000000Ψ1.000000 ΨΨK[2,2]ΨΨ2.000000Ψ1.000000

Ψphi Unif lower bound hyperpriors:Ψ0.001Ψ0.001Ψ Ψphi Unif upper bound hyperpriors:Ψ1.024Ψ1.024Ψ

Ψtau.sq IG hyperpriors shape=2.00000 and scale=1.00000

Source compiled with OpenMP, posterior sampling is using 4 thread(s).

ΨΨSampling

Sampled: 5000 of 10000, 50.00% Report interval Metrop. Acceptance rate: 36.84% Overall Metrop. Acceptance rate: 36.84%

Sampled: 10000 of 10000, 100.00% Report interval Metrop. Acceptance rate: 36.08% Overall Metrop. Acceptance rate: 36.46%

\endMakeFramed

We again pass the \codespSVC object to \codespRecover for composition sampling of the remaining model parameters needed for posterior summaries, model assessment, and subsequent prediction.

\MakeFramed m.3 <- spRecover(m.3, start=floor(0.75*n.samples), thin=2, n.omp.threads=4, verbose=FALSE)

\endMakeFramed

Passing the \codespRecover object to \codespDiag yields several popular model fit diagnostics, two of which are summarized in Tables 1 and 2. Table 1 shows the deviance information criterion (DIC) and associated effective number of parameters pD (Spiegelhalter et al., 2001), while Table 2 presents a posterior predictive loss metric D = G+P proposed by (Gelfand and Ghosh, 1998), where G measures goodness of fit and P penalizes complexity. Models with lower values of DIC or D are preferred over those with higher values. Both metrics favor Model 3 which allows both the intercept and CTM predictor to vary spatially over the study area.

[TABLE]

[TABLE]

Again, passing \codespRecover’s \codecoda objects to \codesummary provides posterior summaries of regression coefficients and covariance parameters.

\MakeFramed

round(summary(m.3$p.beta.recover.samples)$quantiles[,c(3,1,5)],3)

          50%   2.5% 97.5%

(Intercept) 3.189 2.105 4.286 pm10.ctm 0.324 -0.087 0.726

round(summary(m.3$p.theta.recover.samples)$quantiles[,c(3,1,5)],3)

                   50%  2.5% 97.5%

sigma.sq.(Intercept) 0.278 0.146 0.480 sigma.sq.pm10.ctm 0.103 0.066 0.153 tau.sq 0.286 0.137 0.463 phi.(Intercept) 0.426 0.072 0.909 phi.pm10.ctm 0.001 0.001 0.002

\endMakeFramed

Given the spatial decay parameter estimates, the corresponding effective spatial range (defined as the distance at which the correlation drops to 0.05) posterior median and 95% CI for the intercept and CTM processes are approximately 7.03 (3.3, 41.82) km and 2005.4 (1330, 2183.12) km, respectively. While the CTM predictor does have a long spatial range relative to the size of the study area, model fit metrics and the magnitude of its process variance \codesigma.sq.pm10.ctm estimates relative to the intercept process and nugget variance, offer evidence for a space-varying relationship with the outcome variable. This conclusion is further reinforced by Figure 3(b), which shows the posterior median for the CTM predictor regression coefficient, $\tilde{\beta}_{CTM}(\textbf{s})$’s, over observed monitoring locations. These posterior samples, along with those of the space-varying intercept, $\tilde{\beta}_{0}(\textbf{s})$’s, are extracted from \codem.3 and summarized in the code below (\codetilde.beta.0 and \codetilde.beta.ctm are displayed in Figure 3(a)-3(b)).

\MakeFramed tilde.beta.0 <- apply( m.3$p.tilde.beta.recover.samples[["tilde.beta.(Intercept)"]], 1, median)

tilde.beta.ctm <- apply( m.3$p.tilde.beta.recover.samples[["tilde.beta.pm10.ctm"]], 1, median)

\endMakeFramed

We next turn to prediction over the grid of 2336 CTM output locations via a call to \codespPredict. As illustrated below, this call uses samples from a \codespRecover object along with the prediction locations (\codepred.coords) and associated design matrix (\codepred.covars). The argument \codejoint specifies if posterior predictive samples should be drawn from the joint or point-wise distribution.

\MakeFramed m.3.pred <- spPredict(m.3, pred.covars=cbind(1, PM10.pred$pm10.ctm), pred.coords=PM10.pred[,1:2], thin=25, joint=TRUE, n.omp.threads=4, verbose=FALSE)

\endMakeFramed

If the number of prediction locations is large, joint prediction can be prohibitively expensive. Even here with 2336 locations, 51 samples, and using 4 CPUs, joint posterior sampling takes 3.53 minutes verses 0.72 minutes for point-wise sampling.

Joint prediction results are given in the bottom row of Figure 3. The posterior predictive distribution median map (Figure 3(c)) shows three distinct zones of high PM10 values over Central Europe. A compelling quality of MCMC-based inference is access to the posterior predictive distribution. This access facilitates summaries like that given in Figure 3(d) which identifies the probability that a given location will exceed a PM10 value 50 $\mu$g m*-3* (as further explored in Hamm et al. (2015) and Datta et al. (2016)).

4 Summary

The new \codespSVC function more fully implements the computationally efficient MCMC algorithm detailed in Finley et al. (2015) and provides a flexible software tool for fitting spatially varying coefficient models. While other software, some of which are noted in Section 1, offer similar spatially adaptive regression, few provide both univariate and multivariate GP specifications and the computational efficiency delivered by the proposed sampling algorithm and use of \proglangOpenMP parallelization in combination with optional calls to multi-lower-level BLAS and LAPACK multi-threaded matrix algebra libraries. Future work will focus on extending this function to accommodate non-Guassian and multivariate outcomes, as well as for settings where the number of locations precluded the use of full-rank spatial GPs.

Acknowledgments

Finley was supported by National Science Foundation (NSF) EF-1253225 and DMS-1916395, and National Aeronautics and Space Administration’s Carbon Monitoring System project. Banerjee was supported by NSF DMS-1513654, IIS-1562303, and DMS-1916349.