MGPs for univariate spatial gridded data

library(magrittr)
library(dplyr)
library(ggplot2)
library(meshed)
set.seed(2021)

SS <- 50 # coord values for jth dimension 
dd <- 2 # spatial dimension
n <- SS^2 # number of locations
p <- 3 # number of covariates

xlocs <- seq(0, 1, length.out=SS)
coords <- expand.grid(list(xlocs, xlocs))

sigmasq <- 1.5
nu <- 0.5
phi <- 10
tausq <- .1

# covariance at coordinates
CC <- sigmasq * exp(-phi * as.matrix(dist(coords)))
# cholesky of covariance matrix
LC <- t(chol(CC))
# spatial latent effects are a GP
w <- LC %*% rnorm(n)
# measurement errors
eps <- tausq^.5 * rnorm(n)

# covariates and coefficients 
X <- matrix(rnorm(n*p), ncol=p)
Beta <- matrix(rnorm(p), ncol=1)

# univariate outcome, fully observed
y_full <- X %*% Beta + w + eps

# now introduce some NA values in the outcomes
y <- y_full %>% matrix(ncol=1)
y[sample(1:n, n/5, replace=FALSE), ] <- NA

simdata <- coords %>%
  cbind(data.frame(Outcome_full= y_full, 
                   Outcome_obs = y,
                   w = w)) 

simdata %>% ggplot(aes(Var1, Var2, fill=w)) +
    geom_raster() + 
    scale_fill_viridis_c() +
  theme_minimal() + ggtitle("w: Latent GP")

simdata %>% ggplot(aes(Var1, Var2, fill=y)) +
    geom_raster() + 
    scale_fill_viridis_c() + 
  theme_minimal() + ggtitle("y: Observed outcomes")

We now fit the following spatial regression model \[ y = X \beta + w + \epsilon, \] where $w \sim MGP$ are spatial random effects, and $\epsilon \sim N(0, \tau^2)$. For spatial data, an exponential covariance function is used by default: $C(h) = \sigma^2 \exp( -\phi h )$ where $h$ is the spatial distance.

The prior for the spatial decay $\phi$ is $U[l,u]$ and the values of $l$ and $u$ must be specified. We choose $l=1$, $u=30$ for this dataset.¹

Setting verbose tells spmeshed how many intermediate messages to output while running MCMC. We set up MCMC to run for 3000 iterations, of which we discard the first third. We use the argument block_size=25 to specify the coarseness of domain partitioning. In this case, the same result can be achieved by setting axis_partition=c(10, 10).

Finally, we note that NA values for the outcome are OK since the full dataset is on a grid. spmeshed will figure this out and use settings optimized for partly observed lattices, and automatically predict the outcomes at missing locations. On the other hand, X values are assumed to not be missing.

mcmc_keep <- 200 # too small! this is just a vignette.
mcmc_burn <- 400
mcmc_thin <- 2

mesh_total_time <- system.time({
  meshout <- spmeshed(y, X, coords,
                      #axis_partition=c(10,10), #same as block_size=25
                      block_size = 25,
                      n_samples = mcmc_keep, n_burn = mcmc_burn, n_thin = mcmc_thin, 
                      n_threads = 2,
                      verbose = 0,
                      prior=list(phi=c(1,30))
  )})

We can now do some postprocessing of the results. We extract posterior marginal summaries for $\sigma^2$, $\phi$, $\tau^2$, and $\beta_2$. The model that spmeshed targets is a slight reparametrization of the above:² \[ y = X \beta + \lambda w + \epsilon, \] where $w\sim MGP$ has unitary variance. This model is equivalent to the previous one and in fact we find $\sigma^2=\lambda^2$.

summary(meshout$lambda_mcmc[1,1,]^2)
#>    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
#>   1.310   1.833   1.976   2.140   2.319   4.758
summary(meshout$theta_mcmc[1,1,])
#>    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
#>   2.839   5.933   7.305   7.137   7.949  12.564
summary(meshout$tausq_mcmc[1,])
#>    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
#> 0.08837 0.11358 0.12587 0.12465 0.13383 0.17332
summary(meshout$beta_mcmc[1,1,])
#>    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
#> -0.8278 -0.8049 -0.7970 -0.7958 -0.7860 -0.7515

We proceed to plot predictions across the domain along with the recovered latent effects.

# process means
wmesh <- data.frame(w_mgp = meshout$w_mcmc %>% summary_list_mean())
# predictions
ymesh <- data.frame(y_mgp = meshout$yhat_mcmc %>% summary_list_mean())

outdf <- meshout$coordsdata %>% 
  cbind(wmesh, ymesh) %>%
  left_join(simdata)

outdf %>% 
    ggplot(aes(Var1, Var2, fill=w_mgp)) +
    geom_raster() +
    scale_fill_viridis_c() +
    theme_minimal() + ggtitle("w: recovered")


outdf %>% 
    ggplot(aes(Var1, Var2, fill=y_mgp)) +
    geom_raster() +
    scale_fill_viridis_c() +
    theme_minimal() + ggtitle("y: predictions")

spmeshed implements a model which can be interpreted as assigning $\sigma^2$ a folded-t via parameter expansion if settings$ps=TRUE, or an inverse Gamma with parameters $a=2$ and $b=1$ if settings$ps=FALSE, which cannot be changed at this time. $\tau^2$ is assigned an exponential prior.↩︎
At its core, spmeshed implements the spatial factor model $Y(s) = X(s)\beta + \Lambda v(s) + \epsilon(s)$ where $w(s) = \Lambda v(s)$ is modeled via linear coregionalization.↩︎

MGPs for univariate spatial gridded data

M Peruzzi

2022-09-19