Areal Spatial Modelling (R)

spdep for autocorrelation, CARBayes for BYM, sf maps throughout

The R counterpart to areal_python.ipynb. Where the Python notebook builds the CAR sampler from scratch, this notebook uses the field-standard spatial stack: spdep for the spatial weights and Moran's I, CARBayes for the Bayesian BYM disease-mapping model, and sf + ggplot2 to draw the choropleth maps inline. Same data: the North Carolina SIDS counties.

options(repr.plot.width=8.2, repr.plot.height=3.4)
.libPaths(c("C:/Users/user/R/win-library/4.6", .libPaths()))
suppressMessages({library(sf); library(spdep); library(CARBayes); library(ggplot2)})
nc <- st_read(system.file("shape/nc.shp", package="sf"), quiet=TRUE)
nc$expected <- nc$BIR74 * sum(nc$SID74)/sum(nc$BIR74)          # expected SIDS given births
nc$SMR <- nc$SID74 / nc$expected
nc$nonwhite <- scale(nc$NWBIR74 / nc$BIR74)[,1]
cat(nrow(nc), "NC counties;", sum(nc$SID74), "SIDS deaths;", "raw SMR range", round(min(nc$SMR),1), "-", round(max(nc$SMR),1), "\n")

100 NC counties; 667 SIDS deaths; raw SMR range 0 - 4.7 

1. Raw SMR map and spatial weights

The standardised mortality ratio (observed / expected SIDS) mapped across counties — noisy in small counties. poly2nb builds the neighbour graph from the polygons and nb2listw the weights.

nb <- poly2nb(nc); lw <- nb2listw(nb, style="W")
cat("avg neighbours:", round(mean(card(nb)),1), "\n")
ggplot(nc) + geom_sf(aes(fill=SMR), color="white", linewidth=0.15) +
  scale_fill_distiller(palette="OrRd", direction=1) +
  labs(title="Raw SMR (observed / expected SIDS)", fill="SMR") + theme_void()

avg neighbours: 4.9 

2. Moran's I — spatial autocorrelation

moran.test gives the analytic test and moran.mc a permutation version; moran.plot draws the scatterplot of each county against its neighbours' mean.

mt <- moran.mc(nc$SMR, lw, nsim=999)
cat(sprintf("Moran's I = %.3f, permutation p = %.4f  -> SIDS rates are spatially autocorrelated\n", mt$statistic, mt$p.value))
# permutation null distribution: reshuffle the values over the map and recompute I
hist(mt$res, breaks=30, col="grey80", border="white", xlim=range(c(mt$res, mt$statistic)),
     xlab="Moran's I", main=sprintf("Permutation null (999 reshuffles),  p = %.3f", mt$p.value))
abline(v=mt$statistic, col="red", lwd=3); legend("topright", "observed I", col="red", lwd=3, bty="n")
moran.plot(nc$SMR, lw, xlab="county SMR", ylab="mean of neighbours' SMR", main="Moran scatterplot")

Moran's I = 0.231, permutation p = 0.0010  -> SIDS rates are spatially autocorrelated

3. BYM disease mapping with CARBayes

S.CARbym fits the Besag–York–Mollié Poisson model — a spatially structured plus an unstructured random effect — with our nonwhite-births covariate and the births-based offset. The fitted relative risk is the smoothed map.

W <- nb2mat(nb, style="B")
dat <- data.frame(SID74=nc$SID74, expected=nc$expected, nonwhite=nc$nonwhite)   # plain data.frame for the model
set.seed(1)
m <- S.CARbym(SID74 ~ offset(log(expected)) + nonwhite, family="poisson", data=dat,
              W=W, burnin=5000, n.sample=20000, thin=10, verbose=FALSE)
bet <- m$summary.results["nonwhite", c("Mean","2.5%","97.5%")]
cat(sprintf("covariate (nonwhite births): beta = %.2f  95%% CrI [%.2f, %.2f]\n", bet[1], bet[2], bet[3]))
nc$RR <- m$fitted.values / nc$expected
cat(sprintf("shrinkage: raw SMR sd %.2f -> BYM relative-risk sd %.2f\n", sd(nc$SMR), sd(nc$RR)))
p1<-ggplot(nc)+geom_sf(aes(fill=SMR),color="white",linewidth=0.15)+scale_fill_distiller(palette="OrRd",direction=1,limits=c(0,max(nc$SMR)))+labs(title="Raw SMR (noisy)",fill="")+theme_void()
p2<-ggplot(nc)+geom_sf(aes(fill=RR),color="white",linewidth=0.15)+scale_fill_distiller(palette="OrRd",direction=1,limits=c(0,max(nc$SMR)))+labs(title="BYM smoothed relative risk",fill="")+theme_void()
suppressMessages(library(patchwork)); print(p1 + p2)

covariate (nonwhite births): beta = 0.40  95% CrI [0.28, 0.54]

shrinkage: raw SMR sd 0.78 -> BYM relative-risk sd 0.44

The exceedance probability $P(\text{RR}>1.5)$ from the posterior samples flags the counties that are elevated with confidence — the decision map.

risk.samples <- m$samples$fitted / matrix(nc$expected, nrow=nrow(m$samples$fitted), ncol=nrow(nc), byrow=TRUE)
nc$exceed <- colMeans(risk.samples > 1.5)
ggplot(nc) + geom_sf(aes(fill=exceed), color="white", linewidth=0.15) +
  scale_fill_distiller(palette="Reds", direction=1, limits=c(0,1)) +
  labs(title="P(relative risk > 1.5) — exceedance", fill="P") + theme_void()
cat(sprintf("%d counties have P(RR>1.5) > 0.9 (elevated with high confidence)\n", sum(nc$exceed>0.9)))

7 counties have P(RR>1.5) > 0.9 (elevated with high confidence)

4. Summary

The standard R spatial stack reproduces the from-scratch results. spdep confirms significant spatial autocorrelation in the SIDS rates (Moran's I permutation test), CARBayes's S.CARbym fits the Besag–York–Mollié model — recovering the positive nonwhite-births covariate effect and shrinking the noisy raw SMR into a smooth relative-risk map — and sf/ggplot2 render every step as a choropleth, with an exceedance-probability map flagging the confidently-elevated counties.

spdep, CARBayes and sf are the field-standard areal-data tools and the package mirror of the Bayesian sampler in areal_python.ipynb; they are the same tools used for the Scotland lip-cancer BYM in Bayesian Hierarchical Spatial Poisson Model. This areal foundation underpins the spatial arc and maps directly onto regional disability/mortality surveillance. Next, spatial econometrics reuses this neighbour graph in a simultaneous-autoregressive model of spillovers.