`mcmc_mra()`

In this Vignette, we demonstrate how to fit Bayesian spatial models using the multi-resolution model using Markov Chain Monte Carlo (MCMC).

This package provides a Bayesian implementation of the multi-resolution Gaussian process model presented by Nychka et al. (2015). This package is similar to the R package `LatticeKrig::LatticeKrig`

(Nychka et al. (2016)) and provides an implementation that is suitable for use in MCMC samplers.

The model that the function `mcmc_mra()`

fits is defined below. Let \(\mathbf{s}\) be a spatial location in a domain of interest \(\mathcal{D}\) and let \(y(\mathbf{s})\) be the observation of the spatial process at location \(\mathbf{s}\). Then, we can write the observation as

\[\begin{align*} y(\mathbf{s}) & = \mathbf{x}(\mathbf{s})' \boldsymbol{\beta} + \eta(\mathbf{s}) + \varepsilon(\mathbf{s}), \end{align*}\]

where \(\mathbf{x}(s)\) is a vector of covariates at site \(\mathbf{s}\), \(\boldsymbol{\beta}\) are regression coefficients, \(\boldsymbol{\eta}(\mathbf{s})\) is the realization of a correlated Gaussian process at location \(\mathbf{s}\) and \(\varepsilon(\mathbf{s})\) is the realization of an uncorrelated Gaussian error process. The model can be written in matrix form where

\[\begin{align*} \mathbf{y} & = \mathbf{X} \boldsymbol{\beta} + \boldsymbol{\eta} + \boldsymbol{\varepsilon}. \end{align*}\]

In matrix form, the spatially correlated random process has the distribution \(\boldsymbol{\eta} \sim \operatorname{N}(\mathbf{0}, \mathbf{C})\) where the \(\mathbf{C}_{ij}\) (the \(i\)th row and \(j\)th column of \(\mathbf{C}\)) is defined as the output of the covariance function \(cov(\mathbf{s}_i, \mathbf{s}_j)\) at sites \(\mathbf{s}_i\) and \(\mathbf{s}_j\). The residual process \(\boldsymbol{\varepsilon} \sim \operatorname{N}(\mathbf{0}, \sigma^2 \mathbf{I})\) models measurement error, commonly referred to as the nugget in the spatial statistics literature.

Rather than specifying a parametric form for the covariance function \(cov(\mathbf{s}_i, \mathbf{s}_j)\), we approximate the covariance function effect with a multi-resolution approximation. This approximation gives \(\boldsymbol{\eta} \approx \sum_{m=1}^M \mathbf{W}_m \boldsymbol{\alpha}_m\) where \(\mathbf{W}_m\) is the Wendland basis expansion at the \(m\)th resolution and \(\boldsymbol{\alpha}_m\) are the \(m\)th resolution random effect.

First, we simulate some data from the process of interest. We simulate a spatially explicit set of fired effects X that represent spatial variables like elevation, latitude, etc. Next, the spatial process is initialized

```
set.seed(11)
N <- 40^2
## setup the spatial process
locs <- as.matrix(
expand.grid(
seq(0, 1, length.out = sqrt(N)),
seq(0, 1, length.out = sqrt(N))
)
)
D <- fields::rdist(locs)
## fixed effects include intercept, elevation, and latitude
X <- cbind(1, as.vector(mvnfast::rmvn(1, rep(0, N), 3 * exp(-D / 20))), locs[, 2])
p <- ncol(X)
beta <- rnorm(ncol(X))
## MRA spatio-temporal random effect
M <- 3
n_coarse <- 10
MRA <- mra_wendland_2d(locs, M = M, n_coarse = n_coarse, use_spam = TRUE)
# MRA <- mra_wendland_2d(locs, M = M, n_max_fine_grid = 2^8, use_spam = TRUE)
W_list <- MRA$W
n_dims <- rep(NA, length(W_list))
dims_idx <- NULL
for (i in 1:M) {
n_dims[i] <- ncol(W_list[[i]])
dims_idx <- c(dims_idx, rep(i, n_dims[i]))
}
W <- do.call(cbind, W_list)
Q_alpha <- make_Q_alpha_2d(sqrt(n_dims), rep(0.999, length(n_dims)), prec_model = "CAR")
tau2 <- 10 * 2^(2 * (1:M - 1))
Q_alpha_tau2 <- make_Q_alpha_tau2(Q_alpha, tau2)
## initialize the random effect
## set up a linear constraint so that each resolution sums to one
A_constraint <- sapply(1:M, function(i){
tmp = rep(0, sum(n_dims))
tmp[dims_idx == i] <- 1
return(tmp)
})
a_constraint <- rep(0, M)
alpha <- as.vector(rmvnorm.prec.const(n = 1, mu = rep(0, sum(n_dims)), Q = Q_alpha_tau2, A = t(A_constraint), a = a_constraint))
sigma2 <- runif(1, 0.25, 0.5)
y <- X %*% beta + W %*% alpha + rnorm(N, 0, sqrt(sigma2))
```

Next, the sum-to-one constraint is verified by checking that within each dimension the sum of the \(\boldsymbol{\alpha}_m\) parameters is 0.

`sum(alpha)`

`## [1] 1.98782e-14`

`sum(alpha[dims_idx == 1])`

`## [1] 1.450229e-15`

`sum(alpha[dims_idx == 2])`

`## [1] 4.22995e-14`

`sum(alpha[dims_idx == 3])`

`## [1] -2.387153e-14`

To explore the structure of the precision matrix of the spatial random effects, we plot the pattern of nonzero elements in the precision matrices \(\mathbf{Q}_{m}\) for each resolution \(m\). The sparse, banded structure in the precision matrices at each resolution are what enable efficient computation.

```
## plot Q_alpha_tau2 structure
## note: the display function gives a warning about the default cex setting
layout(matrix(1:4, 2, 2, byrow = TRUE))
for (m in 1:M) {
display(Q_alpha_tau2[which(dims_idx == m), which(dims_idx == m)])
title(main = paste("Nonzero elements of \n precision matrix; resolution", m))
}
```

The MRA approach requires a set of grid point at each resolution. The set of grid points at each resolution extends beyond the spatial domain that is highlighted in the black box.

```
data.frame(
x = unlist(sapply(1:M, function(i) MRA$locs_grid[[i]][, 1])),
y = unlist(sapply(1:M, function(i) MRA$locs_grid[[i]][, 2])),
resolution = factor(unlist(sapply(1:M, function(i) rep(i, each = nrow(MRA$locs_grid[[i]])))))
) %>%
ggplot(aes(x = x, y = y, color = resolution)) +
geom_point(alpha = 0.75, size = 0.75) +
ggtitle("MRA grid points") +
theme_bw() +
scale_colour_viridis_d(end = 0.8) +
geom_rect(
data = NULL,
aes(xmin = min(locs[, 1]),
xmax = max(locs[, 1]),
ymin = min(locs[, 2]),
ymax = max(locs[, 2])
),
fill = NA,
color = "black",
alpha = 0.5
) +
theme_custom()
```