Spatiotemporal Modelling in R — CARBayesST¶

The space-time CAR with the package that defines it¶

The Python notebook built the space-time ANOVA from scratch and cross-checked it in PyMC. Here we fit the same model with CARBayesST — the R package purpose-built for spatiotemporal areal data — on the identical simulated Georgia panel (real 159-county geography, real population, a simulated 12-year process with known truth).

ST.CARanova fits exactly the Knorr-Held ANOVA decomposition of the Python notebook: $$y_{it}\sim\text{Poisson}\!\big(E_{it}\,e^{\alpha+x_i'\beta+\phi_i+\delta_t\,(+\,\psi_{it})}\big),$$ a spatial main effect $\phi_i$ (a Leroux CAR field), a temporal main effect $\delta_t$ (also CAR-smoothed over time), and an optional space-time interaction $\psi_{it}$. We fit the additive form (interaction=FALSE) to match the from-scratch model, recover the truth, and map it. The sibling ST.CARar (spatial CAR × AR(1) in time) is the autoregressive alternative.

In [1]:
suppressWarnings(suppressMessages({
  .libPaths(c("C:/Users/user/R/win-library/4.6", .libPaths()))
  library(CARBayesST); library(sf); library(ggplot2)
}))
options(repr.plot.width=9.8, repr.plot.height=3.9)
w   <- read.csv("ga_st_counts.csv"); N <- nrow(w); Tt <- 12
Wm  <- as.matrix(read.csv("ga_adj.csv")); dimnames(Wm) <- NULL
phi_true   <- read.csv("ga_st_truth_spatial.csv")$phi_true
gamma_true <- read.csv("ga_st_truth_temporal.csv")$gamma_true
# long format: spatial units nested within time periods (CARBayesST ordering) = column-major stack
long <- data.frame(
  Y = as.vector(as.matrix(w[, paste0("y", 0:11)])),
  x = rep(w$x, times = Tt),
  E = rep(w$pop / 1000, times = Tt),
  fips = rep(w$fips, times = Tt),
  year = rep(0:11, each = N))
cat(sprintf("Georgia space-time panel: %d counties x %d years = %d rows; %d events\n",
            N, Tt, nrow(long), sum(long$Y)))
cat("W:", nrow(Wm), "x", ncol(Wm), " spatial adjacency (symmetric:", isSymmetric(Wm), ")\n")

Georgia space-time panel: 159 counties x 12 years = 1908 rows; 59114 events

W: 159 x 159  spatial adjacency (symmetric: TRUE )

1. Fit the space-time ANOVA¶

ST.CARanova with a Poisson likelihood, the covariate x, a log(E) offset, and the 159×159 spatial weight matrix. interaction=FALSE gives the additive spatial + temporal model matching the from-scratch engine.

In [2]:
set.seed(1)
m <- ST.CARanova(Y ~ x + offset(log(E)), family="poisson", data=long, W=Wm,
                 interaction=FALSE, burnin=2000, n.sample=7000, thin=5, verbose=FALSE)
cat("samples available:", paste(names(m$samples), collapse=", "), "\n\n")
beta_x <- mean(m$samples$beta[, 2])
bq <- quantile(m$samples$beta[, 2], c(.025, .975))
phi_e <- colMeans(m$samples$phi)          # spatial main effect (length N)
del_e <- colMeans(m$samples$delta)        # temporal main effect (length Tt)
cat(sprintf("covariate beta:  CARBayesST %.2f  [%.2f, %.2f]   (true 0.30)\n", beta_x, bq[1], bq[2]))
cat(sprintf("spatial field recovery:  corr(estimated, true) = %.2f\n", cor(phi_e, phi_true)))
cat(sprintf("temporal wave recovery:  corr(estimated, true) = %.2f\n", cor(del_e, gamma_true)))

samples available: beta, phi, delta, gamma, tau2, rho, fitted, Y 


covariate beta:  CARBayesST 0.29  [0.25, 0.34]   (true 0.30)

spatial field recovery:  corr(estimated, true) = 0.99

temporal wave recovery:  corr(estimated, true) = 1.00

2. Map the space-time structure¶

The estimated spatial field (which counties run high, averaged over time) beside the truth, and the temporal trend — CARBayesST's smoothed epidemic wave against the simulated one. Maps rendered inline from the Georgia county geometry.

In [3]:
geo <- st_read("ga_counties.geojson", quiet=TRUE)
geo$phi_e <- phi_e[match(geo$id, w$fips)]
geo$phi_t <- phi_true[match(geo$id, w$fips)]
lim <- max(abs(c(geo$phi_e, geo$phi_t)))
suppressWarnings(suppressMessages(library(gridExtra)))
g1 <- ggplot(geo) + geom_sf(aes(fill=phi_t), color="white", size=0.05) +
  scale_fill_gradient2(limits=c(-lim,lim), name=expression(phi)) +
  ggtitle("TRUE spatial field") + theme_void() + theme(plot.title=element_text(size=11))
g2 <- ggplot(geo) + geom_sf(aes(fill=phi_e), color="white", size=0.05) +
  scale_fill_gradient2(limits=c(-lim,lim), name=expression(phi)) +
  ggtitle("ESTIMATED spatial field (CARBayesST)") + theme_void() + theme(plot.title=element_text(size=11))
grid.arrange(g1, g2, ncol=2)
No description has been provided for this image
In [4]:
options(repr.plot.width=6, repr.plot.height=3.4)
td <- data.frame(year=0:11, est=del_e, true=gamma_true)
ggplot(td, aes(year)) +
  geom_hline(yintercept=0, linewidth=0.3, color="grey60") +
  geom_line(aes(y=true, color="true"), linewidth=1) + geom_point(aes(y=true, color="true"), size=2) +
  geom_line(aes(y=est, color="estimated"), linewidth=1, linetype=2) + geom_point(aes(y=est, color="estimated"), size=2) +
  scale_color_manual(values=c(true="#2f855a", estimated="#c53030"), name="") +
  labs(x="year", y=expression(delta[t]), title="Temporal main effect: CARBayesST vs truth") +
  theme_minimal() + theme(legend.position="top")
No description has been provided for this image

3. Summary¶

CARBayesST fits the space-time ANOVA directly and recovers the simulated structure — the covariate near its true 0.30, the spatial field and temporal wave both closely correlated with the truth — reproducing the from-scratch and PyMC results of the Python notebook on the identical Georgia panel. The three engines agree.

The package's interaction=TRUE adds the space-time interaction $\psi_{it}$ (counties departing from the shared wave), and ST.CARar swaps the temporal CAR for an AR(1) process — the natural next steps once a real multi-year panel is in hand. That real application is the disability belt over time: SSA's yearly county receipt panel run through exactly this model would show whether the belt is expanding, stable, or shifting — the spatiotemporal continuation of The Disability Belt — a Spatial Application, and the capstone of the spatial arc (areal → econometric → geostatistical → point-process → spatiotemporal).