Scotland lip cancer — the BYM spatial Poisson model¶

From-scratch Metropolis-within-Gibbs with an intrinsic CAR prior¶

Module: scotland_bym.py. The classic disease-mapping problem (Clayton & Kaldor 1987; Besag, York & Mollié 1991; WinBUGS/GeoBUGS "Lip"). 56 Scottish districts: we have observed lip-cancer cases $O_i$, age-standardised expected cases $E_i$ (the offset), and a covariate $x_i$ = % employed in agriculture/fishing/forestry (AFF — a proxy for outdoor sun exposure). The new ingredient over the hierarchical-Poisson projects is a spatially structured random effect: neighbouring districts borrow strength from each other, not just from a global mean.

Model — Besag-York-Mollié (BYM)¶

$$O_i\sim\text{Poisson}(E_i\,e^{\eta_i}),\qquad \eta_i=\alpha+\beta x_i+\theta_i+\phi_i.$$ Two area-level random effects, the heart of BYM:

  • $\theta_i\sim N(0,\sigma_\theta^2)$ — unstructured heterogeneity (exchangeable, as in the Epil GLMM): area-specific noise with no geography.
  • $\phi\sim\text{ICAR}(\sigma_\phi^2)$ — structured spatial effect via the intrinsic conditional autoregressive prior: $$\phi_i\mid\phi_{-i}\sim N\!\Big(\tfrac1{n_i}\sum_{j\sim i}\phi_j,\ \tfrac{\sigma_\phi^2}{n_i}\Big),$$ i.e. each district is pulled toward the average of its neighbours. Equivalently $\log p(\phi)=-\tfrac1{2\sigma_\phi^2}\sum_{i\sim j}(\phi_i-\phi_j)^2$ — a penalty on differences across the map's adjacency graph. It is improper (invariant to a global shift), so $\phi$ is re-centred to sum-to-zero with the level absorbed into $\alpha$.

The fitted relative risk of district $i$ is $e^{\eta_i}$ — a smoothed version of the raw SMR $O_i/E_i$.

Algorithm (Metropolis-within-Gibbs)¶

  1. $(\alpha,\beta)$ — block RW-Metropolis.
  2. $\theta_i$ — vectorised RW-Metropolis (the $\theta_i$ are conditionally independent), curvature-scaled step.
  3. $\phi_i$ — sequential RW-Metropolis (neighbours are coupled, so each update must use the current neighbour values), curvature-scaled.
  4. $\sigma_\theta^2$ Inverse-Gamma ($n/2$ d.f.); $\sigma_\phi^2$ Inverse-Gamma ($(n{-}1)/2$ d.f., since ICAR has rank $n-1$). Re-centre $\theta,\phi$ into $\alpha$ each sweep (location moves).

Where this sits¶

  • vs Epil/Pumps/Rugby (hierarchical Poisson): those pool areas toward a global distribution; BYM adds a local spatial term so nearby areas are similar — pooling with geography.
  • vs fish ZIP (mixture): different latent structure again — continuous spatial field, not a discrete class.
In [1]:
import numpy as np, matplotlib.pyplot as plt
from scotland_bym import scotland_data, bym_gibbs

O, E, x, adj, names = scotland_data('scotland_lip.csv'); n = len(O)
out = bym_gibbs(O, E, x, adj, R=20000, burn=5000, seed=1)
al, be = out['alpha'], out['beta']
print(f"acceptance  (alpha,beta) {out['accept_ab']:.2f}   theta {out['accept_theta']:.2f}   phi {out['accept_phi']:.2f}")
print(f"alpha = {al.mean():.3f}")
print(f"beta  = {be.mean():.3f}  95% CI [{np.percentile(be,2.5):.2f}, {np.percentile(be,97.5):.2f}]"
      f"   -> RR per +10% AFF = {np.exp(be.mean()):.2f}   P(beta>0) = {(be>0).mean():.3f}")
print(f"sd_theta (unstructured) = {out['sd_theta'].mean():.3f}")
print(f"sd_phi   (spatial)      = {out['sd_phi'].mean():.3f}")
frac = (out['sd_phi']**2 / (out['sd_phi']**2 + out['sd_theta']**2)).mean()
print(f"spatial fraction of area variance = {frac:.2%}  -> the heterogeneity is almost entirely SPATIAL")

acceptance  (alpha,beta) 0.39   theta 0.61   phi 0.61
alpha = -0.151
beta  = 0.302  95% CI [0.04, 0.55]   -> RR per +10% AFF = 1.35   P(beta>0) = 0.988
sd_theta (unstructured) = 0.032
sd_phi   (spatial)      = 0.728
spatial fraction of area variance = 99.58%  -> the heterogeneity is almost entirely SPATIAL
In [2]:
smr = O / E
rr  = np.exp(al.mean() + be.mean()*x + out['theta'].mean(0) + out['phi'].mean(0))
print('raw SMR range [%.2f, %.2f]   smoothed RR range [%.2f, %.2f]'%(smr.min(),smr.max(),rr.min(),rr.max()))
top = np.argsort(-rr)[:6]
print('\nhighest-risk districts:')
for i in top: print('  %-18s  O=%2d  E=%5.1f  SMR=%.2f  ->  RR=%.2f'%(names[i],int(O[i]),E[i],smr[i],rr[i]))

fig, ax = plt.subplots(1, 2, figsize=(12, 4.6))
# shrinkage: raw SMR -> smoothed RR, point size ∝ expected count (reliability)
ax[0].scatter(smr, rr, s=8+6*E, color='steelblue', alpha=.8, edgecolor='k', lw=.3, zorder=3)
lim=[0, max(smr.max(),rr.max())*1.05]; ax[0].plot(lim,lim,'k:',lw=1,label='SMR = RR')
ax[0].axhline(1,color='grey',lw=.7); ax[0].axvline(1,color='grey',lw=.7)
ax[0].set_xlabel('raw SMR  $O_i/E_i$'); ax[0].set_ylabel('smoothed RR  $e^{\\eta_i}$')
ax[0].set_title('Spatial smoothing (point size ∝ expected $E_i$)'); ax[0].legend(fontsize=8)
# beta posterior (AFF effect)
ax[1].hist(be, bins=50, density=True, color='seagreen', alpha=.8)
ax[1].axvline(be.mean(),color='k',ls='--',lw=1,label=f'mean {be.mean():.2f}')
ax[1].axvline(0,color='firebrick',ls=':',lw=1,label='no effect')
ax[1].set_xlabel('β  (AFF effect, per +10%)'); ax[1].set_yticks([])
ax[1].set_title(f'AFF covariate  (P(β>0) = {(be>0).mean():.3f})'); ax[1].legend(fontsize=8)
plt.tight_layout(); plt.savefig('scotland_bym.png', dpi=120, bbox_inches='tight'); plt.show()

raw SMR range [0.00, 6.52]   smoothed RR range [0.36, 4.51]

highest-risk districts:
  Skye-Lochalsh       O= 9  E=  1.4  SMR=6.52  ->  RR=4.51
  Banff-Buchan        O=39  E=  8.7  SMR=4.50  ->  RR=4.34
  Okney               O= 8  E=  2.4  SMR=3.33  ->  RR=3.67
  Caithness           O=11  E=  3.0  SMR=3.62  ->  RR=3.24
  Sutherland          O= 5  E=  1.8  SMR=2.79  ->  RR=3.24
  Ross-Cromarty       O=15  E=  4.3  SMR=3.52  ->  RR=3.15
No description has been provided for this image

Mapping the risk surface¶

The scatter above shows the smoothing abstractly; on the actual map of the 56 (pre-1996) Scottish districts it becomes geographic. The boundaries come from the SpatialEpi R package, exported to scotland_geo.csv by scotland_R.ipynb; the district order matches our data exactly, so each polygon takes our model's per-area value directly. Left = the raw SMR (noisy, with extreme values in tiny districts); right = the BYM-smoothed RR (a coherent risk surface — the speckle is replaced by a clear north-vs-central-belt gradient).

In [3]:
import pandas as pd, matplotlib as mpl
from matplotlib.collections import PolyCollection

# district boundaries exported from SpatialEpi by scotland_R.ipynb (same row order as our data)
geo = pd.read_csv('scotland_geo.csv')
polys, pareas = [], []
for (a, p), g in geo.groupby(['area', 'part']):
    polys.append(g[['x', 'y']].values); pareas.append(int(a))
pareas = np.array(pareas)

smr = O / E
rr  = np.exp(al.mean() + be.mean()*x + out['theta'].mean(0) + out['phi'].mean(0))
brks = [0, 0.5, 1, 1.5, 2, 3, 8]
cmap = mpl.colors.ListedColormap(['#ffffcc','#ffeda0','#fd8d3c','#f03b20','#bd0026','#800026'])
norm = mpl.colors.BoundaryNorm(brks, cmap.N)

fig, ax = plt.subplots(1, 2, figsize=(11, 6.5))
for k, (v, t) in enumerate([(smr, 'Raw SMR  (O/E)'), (rr, 'Smoothed RR  (BYM)')]):
    pc = PolyCollection(polys, array=v[pareas], cmap=cmap, norm=norm, edgecolor='grey', lw=.3)
    ax[k].add_collection(pc); ax[k].autoscale(); ax[k].set_aspect('equal'); ax[k].axis('off'); ax[k].set_title(t)
cb = fig.colorbar(pc, ax=ax, shrink=.6, ticks=brks); cb.set_label('rate ratio')
plt.savefig('scotland_maps.png', dpi=120, bbox_inches='tight'); plt.show()
No description has been provided for this image

Results & interpretation¶

  • Outdoor work raises lip-cancer risk: β ≈ 0.30 per +10% in AFF (RR ≈ 1.35), with P(β>0) ≈ 0.99 — consistent with the sun-exposure hypothesis behind this classic dataset.
  • The variation is almost entirely spatial: σ_φ ≈ 0.73 vs σ_θ ≈ 0.03, so ~99% of the area-level heterogeneity is structured — lip-cancer risk clusters geographically (high in the rural north — Skye-Lochalsh, Banff-Buchan, Caithness — low in the urban central belt), not as unstructured noise.
  • Smoothing shrinks the unreliable rates: the raw SMR ranges 0–6.5, but the smoothed RR is pulled to 0.36–4.5. Small-$E$ districts (tiny expected counts ⇒ noisy SMR) are pulled hardest toward their neighbours — exactly what the left panel shows (small points move farthest off the diagonal).

Takeaways¶

  • The CAR prior is "pooling with a map". Where the hierarchical-Poisson models borrowed strength from a global mean, the ICAR term borrows from adjacent areas — turning noisy small-area rates into a coherent risk surface.
  • Two effects, two roles. θ (unstructured) and φ (spatial) compete to explain the area heterogeneity; here the data overwhelmingly assign it to φ. (The θ/φ split is only weakly identified individually — the BYM identifiability issue — but their sum, the area effect, and the variance partition are well estimated.)
  • Cross-checks: PyMC (pm.ICAR) and R (CARBayes::S.CARbym) follow.