Joint Disease Mapping — Multivariate CAR (R)¶

CARBayes::MVS.CARleroux — the two outcomes with a between-outcome covariance¶

The R counterpart to mcar_python.ipynb. The Python notebook builds a shared-component CAR (a common field plus period-specific fields) from scratch. Here we fit the alternative, package-standard multivariate CAR: CARBayes::MVS.CARleroux models the two outcomes jointly through a Leroux CAR with a $2\times2$ between-outcome covariance $\Sigma$, whose off-diagonal gives the same cross-period correlation from the other direction. Same data: North Carolina SIDS in 1974–78 and 1979–84.

In [1]:
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); library(patchwork)})
nc <- st_read(system.file("shape/nc.shp", package="sf"), quiet=TRUE)
nc$E1 <- nc$BIR74*sum(nc$SID74)/sum(nc$BIR74); nc$E2 <- nc$BIR79*sum(nc$SID79)/sum(nc$BIR79)
nc$SMR1 <- nc$SID74/nc$E1; nc$SMR2 <- nc$SID79/nc$E2
cat("NC SIDS: 1974-78 =", sum(nc$SID74), "deaths, 1979-84 =", sum(nc$SID79), "deaths;", nrow(nc), "counties\n")

NC SIDS: 1974-78 = 667 deaths, 1979-84 = 836 deaths; 100 counties

1. The two raw maps¶

Raw SMR for each period, side by side — broadly similar but noisy.

In [2]:
vmax <- max(nc$SMR1, nc$SMR2)
p1<-ggplot(nc)+geom_sf(aes(fill=SMR1),color="white",linewidth=0.15)+scale_fill_distiller(palette="OrRd",direction=1,limits=c(0,vmax))+labs(title="Raw SMR 1974-78",fill="")+theme_void()
p2<-ggplot(nc)+geom_sf(aes(fill=SMR2),color="white",linewidth=0.15)+scale_fill_distiller(palette="OrRd",direction=1,limits=c(0,vmax))+labs(title="Raw SMR 1979-84",fill="")+theme_void()
print(p1 + p2)
cat("raw-SMR correlation between periods:", round(cor(nc$SMR1,nc$SMR2),2), "(contaminated by small-area noise)\n")

raw-SMR correlation between periods: 0.21 (contaminated by small-area noise)
No description has been provided for this image

2. Multivariate CAR with MVS.CARleroux¶

The response is a two-column matrix (the two periods); the two-column offset carries each period's expected counts. MVS.CARleroux fits a joint Leroux CAR and estimates the between-outcome covariance $\Sigma$; the correlation $\Sigma_{12}/\sqrt{\Sigma_{11}\Sigma_{22}}$ is the cross-period spatial correlation.

In [3]:
nb <- poly2nb(nc); W <- nb2mat(nb, style="B")
Y <- cbind(nc$SID74, nc$SID79); logE <- cbind(log(nc$E1), log(nc$E2))
set.seed(1)
m <- MVS.CARleroux(Y ~ offset(logE), family="poisson", W=W, burnin=5000, n.sample=25000, thin=10, verbose=FALSE)
S <- m$samples$Sigma                                   # (draws, 2, 2)
cors <- S[,1,2]/sqrt(S[,1,1]*S[,2,2])
cat(sprintf("between-period spatial correlation: %.2f  95%% CrI [%.2f, %.2f]\n", mean(cors), quantile(cors,.025), quantile(cors,.975)))
cat(sprintf("Leroux spatial-dependence rho: %.2f\n", mean(m$samples$rho)))
hist(cors, breaks=40, col="#2f855a", border="white", xlab="between-period correlation", main="MVS.CARleroux: cross-period spatial correlation")
abline(v=mean(cors), col="red", lwd=2)
cat("The multivariate model reaches the same verdict as the from-scratch shared-component model: a moderate-to-\n")
cat("substantial positive correlation -- the two periods' SIDS geographies are related but not identical.\n")

between-period spatial correlation: 0.53  95% CrI [0.11, 0.86]

Leroux spatial-dependence rho: 0.55

The multivariate model reaches the same verdict as the from-scratch shared-component model: a moderate-to-

substantial positive correlation -- the two periods' SIDS geographies are related but not identical.
No description has been provided for this image

3. The two smoothed maps¶

The fitted relative risks for the two periods. Estimated jointly, they borrow strength across periods, so the shared high-risk core is rendered more cleanly than in either raw map.

In [4]:
fit <- m$fitted.values
if (is.matrix(fit)) { mu1<-fit[,1]; mu2<-fit[,2] } else { mu1<-fit[1:nrow(nc)]; mu2<-fit[(nrow(nc)+1):(2*nrow(nc))] }
nc$RR1 <- mu1/nc$E1; nc$RR2 <- mu2/nc$E2
vmax2 <- max(nc$RR1, nc$RR2)
q1<-ggplot(nc)+geom_sf(aes(fill=RR1),color="white",linewidth=0.15)+scale_fill_distiller(palette="OrRd",direction=1,limits=c(0,vmax2))+labs(title="Smoothed RR 1974-78",fill="")+theme_void()
q2<-ggplot(nc)+geom_sf(aes(fill=RR2),color="white",linewidth=0.15)+scale_fill_distiller(palette="OrRd",direction=1,limits=c(0,vmax2))+labs(title="Smoothed RR 1979-84",fill="")+theme_void()
print(q1 + q2)
cat(sprintf("smoothed RR correlation between periods: %.2f (vs raw %.2f)\n", cor(nc$RR1,nc$RR2), cor(nc$SMR1,nc$SMR2)))

smoothed RR correlation between periods: 0.78 (vs raw 0.21)
No description has been provided for this image

4. Summary¶

MVS.CARleroux fits the joint (multivariate) CAR directly: a two-column response, a two-column offset, and a between-outcome covariance $\Sigma$ whose correlation is the cross-period spatial correlation. On NC SIDS it agrees with the from-scratch shared-component model in mcar_python.ipynb — a moderate-to-substantial positive correlation, so the two periods' geographies are related but not identical — and it renders the two smoothed maps more cleanly than the raw SMRs by borrowing strength across periods.

The two approaches are complementary parametrisations of the same idea: the shared-component model (Python) makes the common field explicit and decomposes variance into shared vs specific; the multivariate-CAR model (MVS.CARleroux) works through a between-outcome covariance and reads the correlation from $\Sigma$. CARBayes (with sf/ggplot2 for the maps) is the standard R route, extending the areal foundation to multiple outcomes — the tool for jointly mapping related regional measures. Next the arc changes philosophy rather than outcome type: Spatial Econometrics — Lag, Error, and Spillovers keeps this same neighbour graph but replaces the conditional autoregression with a simultaneous one, where a region's outcome depends on its neighbours' outcomes directly.