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

Defining the MRA model

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.

Simulating data from the model

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

Plotting the simulated data

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))
}

plot of chunk plot-precision-matrices

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()