Bayesian Hierarchical Spatial Poisson Model

Model

The spatial extension of the hierarchical Poisson models: the classic Besag–York–Mollié (BYM) disease-mapping model (Clayton & Kaldor 1987; Besag, York & Mollié 1991; the WinBUGS/GeoBUGS "Lip" example). The new ingredient is a spatially structured random effect — neighbouring areas borrow strength from each other, not just from a global mean. Applied to 56 Scottish districts with observed lip-cancer cases $O_i$, age-standardised expected cases $E_i$ (the offset), and an outdoor-work covariate (agriculture/fishing/forestry — a sun-exposure proxy).

• $O_i,\, E_i$ — observed and age-standardised expected cases (offset) in district $i$
• $\theta_i \sim N(0, \sigma_\theta^2)$ — unstructured (exchangeable) heterogeneity
• $\phi \sim \text{ICAR}(\sigma_\phi^2)$ — structured spatial effect over the adjacency graph

Two random effects — unstructured + spatial

BYM splits the area heterogeneity into two competing random effects: $\theta_i \sim N(0, \sigma_\theta^2)$ is unstructured (exchangeable noise with no geography, as in the Epil GLMM), while $\phi$ is the structured spatial effect via the intrinsic conditional autoregressive (ICAR) prior — each district pulled toward the average of its neighbours, $\phi_i \mid \phi_{-i} \sim N(\frac{1}{n_i}\sum_{j\sim i}\phi_j,\ \sigma_\phi^2/n_i)$. Equivalently the ICAR penalises differences across the map's adjacency graph. It is improper (invariant to a global shift), so $\phi$ is re-centred to sum-to-zero each sweep with the level absorbed into $\alpha$. The fitted relative risk $e^{\eta_i}$ is a smoothed version of the raw SMR $O_i/E_i$ — the standardised morbidity ratio, observed cases divided by the number expected once the district's age profile is accounted for, so 1.0 is the national average and 3.0 is three times it. Raw SMRs are notoriously unstable in small districts, where one extra case can move the ratio enormously, which is the problem the model exists to solve.

Sampler — Metropolis-within-Gibbs

The ICAR coupling forces a structural difference from the earlier GLMMs: the unstructured $\theta_i$ are conditionally independent and updated vectorised, but the spatial $\phi_i$ are coupled to their neighbours and must be updated sequentially (each step using the current neighbour values). Variances are conjugate Inverse-Gamma — with $\sigma_\phi^2$ carrying $(n-1)/2$ degrees of freedom since the ICAR has rank $n-1$.

Re-centre $\theta, \phi$ into $\alpha$ each sweep (location moves; the ICAR level is unidentified).

Notebooks

The fit recovers the classic findings: outdoor work raises risk ($\beta \approx 0.3\text{–}0.37$ per +10% in the AFF covariate, RR $\approx 1.35\text{–}1.45$, $P(\beta>0)\approx 0.99$) — the sun-exposure hypothesis behind this dataset. And the variation is almost entirely spatial: $\sigma_\phi \gg \sigma_\theta$, so ~99% of the area heterogeneity is structured — lip-cancer risk clusters geographically (high in the rural north — Skye-Lochalsh, Banff-Buchan, Caithness — low in the urban central belt). The CAR prior is "pooling with a map": where the hierarchical-Poisson models borrowed strength from a global mean, the ICAR borrows from adjacent areas, turning noisy small-area rates (raw SMR 0–6.5) into a coherent risk surface (smoothed RR $\approx 0.36\text{–}4.5$). The $\theta$/$\phi$ split is only weakly identified individually (the BYM identifiability issue), but the covariate effect, spatial dominance, and district risk map agree across all three engines — from-scratch Metropolis-within-Gibbs, PyMC (pm.ICAR, non-centered NUTS), and R CARBayes::S.CARbym.

Downloads

References

• Besag, J., York, J. & Mollié, A. (1991). Bayesian image restoration, with two applications in spatial statistics. Annals of the Institute of Statistical Mathematics 43(1), 1–20. — the BYM model (unstructured + ICAR-structured random effects)
• Clayton, D. & Kaldor, J. (1987). Empirical Bayes estimates of age-standardized relative risks for use in disease mapping. Biometrics 43(3), 671–681. — the Scotland lip-cancer data and the disease-mapping problem
• Besag, J. (1974). Spatial interaction and the statistical analysis of lattice systems. Journal of the Royal Statistical Society: Series B 36(2), 192–236. — the conditional autoregressive (CAR) construction underlying the ICAR prior
• Riebler, A., Sørbye, S. H., Simpson, D. & Rue, H. (2016). An intuitive Bayesian spatial model for disease mapping that accounts for scaling. Statistical Methods in Medical Research 25(4), 1145–1165. — the modern (BYM2) reparameterisation and the θ/φ identifiability issue
• Lee, D. (2013). CARBayes: an R package for Bayesian spatial modeling with conditional autoregressive priors. Journal of Statistical Software 55(13), 1–24. — the CARBayes::S.CARbym implementation used in the R cross-check

BYM Sampler — Source Code

"""
Scotland lip cancer -- Besag-York-Mollie (BYM) spatial Poisson, from-scratch Metropolis-within-Gibbs.
The classic disease-mapping model (Clayton & Kaldor 1987; Besag, York & Mollie 1991), WinBUGS/GeoBUGS
"Lip" example. 56 Scottish districts; the new ingredient vs the hierarchical-Poisson projects is a
STRUCTURED spatial random effect with an intrinsic CAR prior over the map's adjacency graph.

Model
-----
  O_i ~ Poisson(E_i * exp(eta_i)),   eta_i = alpha + beta * x_i + theta_i + phi_i
    O_i = observed lip-cancer cases,  E_i = age-standardised EXPECTED cases (offset),
    x_i = AFF/10  (% in agriculture/fishing/forestry -- outdoor sun exposure), divided by 10.
  theta_i ~ N(0, sigma_theta^2)                     -- UNSTRUCTURED heterogeneity (exchangeable)
  phi    ~ ICAR(sigma_phi^2)                          -- STRUCTURED spatial effect:
           phi_i | phi_{-i} ~ N( mean_{j~i} phi_j ,  sigma_phi^2 / n_i )
  alpha, beta ~ N(0, 100^2);  precisions 1/sigma^2 ~ Gamma(0.5, 0.0005)  (WinBUGS vague priors).

The ICAR prior penalises differences between neighbours: log p(phi) = -1/(2 sigma_phi^2) * sum_{i~j}
(phi_i - phi_j)^2. It is improper (invariant to adding a constant) -> phi is re-centred to sum-to-zero
each sweep, with the level absorbed into alpha (likewise theta). The fitted relative risk of area i is
exp(eta_i); the smoothed SMR.

Sampler (Metropolis-within-Gibbs)
---------------------------------
  1. (alpha, beta)        -- block RW-Metropolis.
  2. theta_i              -- VECTORISED RW-Metropolis (theta_i conditionally independent given the rest).
  3. phi_i                -- SEQUENTIAL RW-Metropolis (neighbours are coupled: must use current values).
  +  centre theta, phi -> alpha (location moves; phi's ICAR level is unidentified).
  4. sigma_theta^2        -- Inverse-Gamma (n/2 df).
     sigma_phi^2          -- Inverse-Gamma ((n-1)/2 df; ICAR has rank n-1).
"""
import numpy as np
import pandas as pd
import ast


def scotland_data(csv='scotland_lip.csv'):
    d = pd.read_csv(csv)
    O = d['CANCER'].to_numpy(float); E = d['CEXP'].to_numpy(float)
    x = d['AFF'].to_numpy(float) / 10.0                       # per 10% in AFF
    adj = [np.array(ast.literal_eval(a), int) - 1 for a in d['ADJ']]   # 0-indexed neighbour lists
    return O, E, x, adj, d['NAME'].tolist()


def bym_gibbs(O, E, x, adj, R=20000, burn=5000, seed=0, prior_fixed_sd=100.0,
              g_a=0.5, g_b=0.0005, s_ab=0.07, s_th=1.4, s_ph=1.4):
    rng = np.random.default_rng(seed); n = len(O)
    ni = np.array([len(a) for a in adj], float)
    alpha = np.log(O.sum() / E.sum()); beta = 0.0
    theta = np.zeros(n); phi = np.zeros(n); s2_th = 0.1; s2_ph = 0.1
    keep = R - burn
    A = np.zeros(keep); B = np.zeros(keep); STH = np.zeros(keep); SPH = np.zeros(keep)
    TH = np.zeros((keep, n)); PH = np.zeros((keep, n))
    acc_ab = 0; acc_th = 0.0; acc_ph = 0.0
    for it in range(R):
        # 1. (alpha, beta) block RW-Metropolis
        e0 = alpha + beta * x + theta + phi; mu = E * np.exp(e0)
        ll = np.sum(O * e0 - mu) - 0.5 * (alpha**2 + beta**2) / prior_fixed_sd**2
        ap = alpha + s_ab * rng.standard_normal(); bp = beta + s_ab * rng.standard_normal()
        ep = ap + bp * x + theta + phi; mup = E * np.exp(ep)
        llp = np.sum(O * ep - mup) - 0.5 * (ap**2 + bp**2) / prior_fixed_sd**2
        if np.log(rng.random()) < llp - ll:
            alpha, beta = ap, bp; acc_ab += 1
        # 2. theta_i  (vectorised; conditionally independent; curvature-scaled step)
        e0 = alpha + beta * x + theta + phi; mu = E * np.exp(e0)
        curv = mu + 1.0 / s2_th                              # Poisson curvature + prior precision
        prop = theta + s_th * rng.standard_normal(n) / np.sqrt(curv)
        mup = E * np.exp(e0 + (prop - theta))
        d = O * (prop - theta) - (mup - mu) - 0.5 * (prop**2 - theta**2) / s2_th
        ok = np.log(rng.random(n)) < d; theta = np.where(ok, prop, theta); acc_th += ok.mean()
        # 3. phi_i  (sequential; ICAR conditional N(neighbour mean, s2_ph/n_i); curvature-scaled)
        a_ph = 0
        for i in range(n):
            nbar = phi[adj[i]].mean(); prec = ni[i] / s2_ph
            cur = alpha + beta * x[i] + theta[i] + phi[i]; mui = E[i] * np.exp(cur)
            step = s_ph / np.sqrt(mui + prec)
            pp = phi[i] + step * rng.standard_normal(); muip = E[i] * np.exp(cur + (pp - phi[i]))
            d = O[i] * (pp - phi[i]) - (muip - mui) - 0.5 * prec * ((pp - nbar)**2 - (phi[i] - nbar)**2)
            if np.log(rng.random()) < d:
                phi[i] = pp; a_ph += 1
        acc_ph += a_ph / n
        # centre theta, phi -> alpha (location moves)
        m = theta.mean(); theta -= m; alpha += m
        m = phi.mean();   phi -= m;   alpha += m
        # 4. variances (Inverse-Gamma; precision ~ Gamma(g_a, g_b))
        s2_th = 1.0 / rng.gamma(g_a + n / 2.0, 1.0 / (g_b + 0.5 * np.sum(theta**2)))
        S = 0.5 * sum(np.sum((phi[i] - phi[adj[i]])**2) for i in range(n))   # sum over edges
        s2_ph = 1.0 / rng.gamma(g_a + (n - 1) / 2.0, 1.0 / (g_b + 0.5 * S))
        if it >= burn:
            k = it - burn
            A[k] = alpha; B[k] = beta; STH[k] = np.sqrt(s2_th); SPH[k] = np.sqrt(s2_ph)
            TH[k] = theta; PH[k] = phi
    return dict(alpha=A, beta=B, sd_theta=STH, sd_phi=SPH, theta=TH, phi=PH,
                accept_ab=acc_ab / R, accept_theta=acc_th / R, accept_phi=acc_ph / R)


def summary(draws, names=None):
    d = draws.reshape(draws.shape[0], -1); q = d.shape[1]
    names = names or [f'p{j}' for j in range(q)]
    return [(nm, d[:, j].mean(), d[:, j].std(), *np.percentile(d[:, j], [2.5, 97.5]))
            for j, nm in enumerate(names)]