Calculating spatial correlations in a MRF

coupling.R

library(mgcv)
library(spdep)
library(sf)
library(ggplot2)

# Columbus crime data from mgcv/spdep
example(columbus)

# Make a neighborhood object compatible with mgcv's MRF
nb <- poly2nb(columbus)
names(nb) <- attr(nb, "region.id")
columbus$ID <- as.factor(names(nb))

# Fit a model with an MRF to represent spatial structure but nothing else for now
b <- gam(CRIME ~ s(ID, bs="mrf", xt=list(nb=nb)), data=columbus, method="REML")

# Calculate predictions and their covariance
predmat <-  predict(b, type="lpmatrix")
V <- vcov(b, unconditional = TRUE)
C <- coef(b)
preds <- predmat %*% C
pred_vcov <- predmat %*% V %*% t(predmat)  # <-- YO IS THIS RIGHT DAVE? (DAVE HAS SINCE SAID THAT YES LOOKS RIGHT)
pred_corr <- cov2cor(pred_vcov)

columbus$pred <- preds

# Convert correlations to an sf network of lines with weights
nb_corr <- pred_corr * nb2mat(nb, style = "B")
centroids <- st_coordinates(st_centroid(columbus))
nb_net <-  nb2lines(nb = nb, coords = centroids, as_sf = TRUE)
for (z in seq_len(nrow(nb_net))) {
  nb_net$wt[z] <- pred_corr[nb_net$i[z], nb_net$j[z]]
}

myplot <- ggplot() +
  geom_sf(data = columbus, aes(fill = preds)) +
  geom_sf(data = nb_net, aes(col = wt, size = abs(wt))) +
  geom_point(data = as.data.frame(centroids), aes(x=X, y=Y)) +
  scale_fill_continuous(name = "Predicted value") +
  scale_color_distiller(type = "div", limits = c(-1,1)*max(abs(nb_net$wt)), name = "correlation") +
  scale_size_continuous(range = c(0.25, 2), name = "abs(correlation)") +
  theme_minimal() +
  ggtitle(paste0("\"Spatial Coupling\" / Correlation of Predictions in a Markov Random Field (k = ", b$smooth[[1]]$bs.dim, ")"))
myplot

ggsave("coupling.png", plot = myplot)