Scotland lip cancer — R cross-check (CARBayes::S.CARbym)¶

Companion to scotland_python.ipynb / scotland_pymc.ipynb. CARBayes is the standard R package for areal disease mapping; S.CARbym fits exactly the Besag-York-Mollié convolution model (unstructured + spatial CAR random effects) by MCMC. We pass the same adjacency matrix and check it returns the same AFF effect and smoothed risks.

In [1]:
ok <- suppressWarnings(suppressMessages(require(CARBayes)))
if (!ok) { install.packages('CARBayes', repos='https://cloud.r-project.org'); library(CARBayes) }
d <- read.csv('scotland_lip.csv'); n <- nrow(d)
# build adjacency matrix W from the ADJ column (1-indexed neighbour lists)
W <- matrix(0, n, n)
for (i in 1:n) { a <- as.integer(strsplit(gsub('\\[|\\]','',d$ADJ[i]), ',')[[1]]); W[i, a] <- 1 }
stopifnot(isSymmetric(W))
df <- data.frame(O = d$CANCER, E = d$CEXP, x = d$AFF/10)
set.seed(1)
fit <- S.CARbym(O ~ offset(log(E)) + x, family='poisson', data=df, W=W,
                burnin=5000, n.sample=30000, thin=5, verbose=FALSE)
print(round(fit$summary.results, 3))

              Mean   2.5% 97.5% n.sample % accept n.effective Geweke.diag
(Intercept) -0.221 -0.457 0.028     5000     34.6      1118.4         1.9
x            0.374  0.106 0.621     5000     34.6       845.7        -2.0
tau2         0.471  0.191 0.904     5000    100.0       917.4         0.7
sigma2       0.013  0.002 0.067     5000    100.0        73.3         0.0
In [2]:
sr <- fit$summary.results
bx <- sr['x','Mean']
cat(sprintf('\nAFF effect beta = %.3f  95%% CI [%.2f, %.2f]   RR per +10%% = %.2f\n',
            bx, sr['x','2.5%'], sr['x','97.5%'], exp(bx)))
cat(sprintf('spatial variance tau2 = %.3f (sd %.3f)   unstructured sigma2 = %.3f (sd %.3f)\n',
            sr['tau2','Mean'], sqrt(sr['tau2','Mean']), sr['sigma2','Mean'], sqrt(sr['sigma2','Mean'])))

# smoothed relative risks RR_i = fitted / E ; show the highest-risk districts
rr <- fitted(fit) / df$E; smr <- df$O / df$E
ord <- order(-rr)[1:6]
cat('\nhighest-risk districts:\n')
for (i in ord) cat(sprintf('  %-18s O=%2d E=%5.1f  SMR=%.2f -> RR=%.2f\n',
                           d$NAME[i], df$O[i], df$E[i], smr[i], rr[i]))
cat(sprintf('\nsmoothed RR range [%.2f, %.2f]  (raw SMR [%.2f, %.2f])\n',
            min(rr), max(rr), min(smr), max(smr)))

AFF effect beta = 0.374  95% CI [0.11, 0.62]   RR per +10% = 1.45

spatial variance tau2 = 0.471 (sd 0.686)   unstructured sigma2 = 0.013 (sd 0.115)

highest-risk districts:

  Skye-Lochalsh      O= 9 E=  1.4  SMR=6.52 -> RR=4.72
  Banff-Buchan       O=39 E=  8.7  SMR=4.50 -> RR=4.35
  Okney              O= 8 E=  2.4  SMR=3.33 -> RR=3.86
  Sutherland         O= 5 E=  1.8  SMR=2.79 -> RR=3.35
  Caithness          O=11 E=  3.0  SMR=3.62 -> RR=3.24
  Ross-Cromarty      O=15 E=  4.3  SMR=3.52 -> RR=3.21

smoothed RR range [0.36, 4.72]  (raw SMR [0.00, 6.52])

Mapping the risk surface (SpatialEpi boundaries)¶

SpatialEpi::scotland bundles the 56 district polygons together with the same case/expected data, in the same row order as our analysis (verified by stopifnot below) — so we can colour each district directly by its raw SMR and its BYM-smoothed RR. The smoothing replaces the noisy raw map with a coherent risk surface. This cell also exports the geometry to scotland_geo.csv for the Python notebook to render the identical maps.

In [3]:
ok <- suppressWarnings(suppressMessages(require(SpatialEpi)))
if (!ok) { install.packages('SpatialEpi', repos='https://cloud.r-project.org'); library(SpatialEpi) }
data(scotland); sp <- scotland$spatial.polygon       # 56 district polygons, SAME row order as our data
stopifnot(all(scotland$data$cases == df$O))          # confirm alignment

smr <- df$O / df$E; rr <- fitted(fit) / df$E         # raw vs BYM-smoothed
brks <- c(0, 0.5, 1, 1.5, 2, 3, 8)
pal  <- c('#ffffcc','#ffeda0','#fd8d3c','#f03b20','#bd0026','#800026')
lab  <- levels(cut(smr, brks, include.lowest=TRUE))
options(repr.plot.width = 10, repr.plot.height = 6)
par(mfrow = c(1, 2), mar = c(1, 1, 2, 1))
plot(sp, col = pal[cut(smr, brks, include.lowest=TRUE)], border='grey50'); title('Raw SMR  (O/E)')
legend('topleft', legend = lab, fill = pal, bty='n', cex=0.7, title='rate ratio')
plot(sp, col = pal[cut(rr,  brks, include.lowest=TRUE)], border='grey50'); title('Smoothed RR  (BYM)')

# export polygon geometry for the Python notebook (area 0-indexed; one row per vertex)
rows <- list(); k <- 1
for (i in seq_along(sp@polygons)) for (j in seq_along(sp@polygons[[i]]@Polygons)) {
  co <- sp@polygons[[i]]@Polygons[[j]]@coords
  rows[[k]] <- data.frame(area = i-1, part = j-1, x = co[,1], y = co[,2]); k <- k + 1
}
write.csv(do.call(rbind, rows), 'scotland_geo.csv', row.names = FALSE); cat('wrote scotland_geo.csv\n')

Installing package into 'C:/Users/user/R/win-library/4.6'
(as 'lib' is unspecified)


package 'SpatialEpi' successfully unpacked and MD5 sums checked

The downloaded binary packages are in
	C:\Users\user\AppData\Local\Temp\Rtmp2B1tyB\downloaded_packages

Loading required package: sp


wrote scotland_geo.csv
No description has been provided for this image

Results¶

  • CARBayes::S.CARbym matches the PyMC fit closely: β ≈ 0.37 (RR ≈ 1.45 per +10% AFF), CI excluding 0, and a spatial variance (tau2, sd ≈ 0.69) far larger than the unstructured (sigma2, sd ≈ 0.12) — the same "risk is geographically clustered" conclusion (σ_phi ≈ 0.69–0.73 across all three engines).
  • The same northern rural districts (Skye-Lochalsh, Banff-Buchan, Caithness, …) top the smoothed-risk list, and the raw SMR range (0–6.5) is pulled in toward neighbours (≈ 0.36–4.7) just as in the Python/PyMC fits.
  • Three engines, one spatial model: from-scratch Metropolis-within-Gibbs ≈ PyMC (pm.ICAR + NUTS, non-centered) ≈ R CARBayes. The exact tau2/sigma2 split varies a little with each package's variance priors (and the BYM θ-vs-φ identifiability), but the spatial dominance and the district risk map agree closely across all three; on the covariate effect PyMC and CARBayes land at RR ≈ 1.45 (β ≈ 0.37) while the from-scratch Gibbs is a little lower at RR ≈ 1.35 (β ≈ 0.30) — same sign and significance.