Scotland lip cancer — PyMC cross-check (BYM via pm.ICAR)¶

Companion to scotland_python.ipynb (from-scratch Metropolis-within-Gibbs). PyMC ships the intrinsic CAR prior as pm.ICAR(W, sigma), so the BYM model is a few lines and NUTS samples it directly. We confirm it reproduces the from-scratch covariate effect, variance partition, and per-area risks.

Model¶

$$O_i\sim\text{Poisson}(E_i e^{\alpha+\beta x_i+\theta_i+\phi_i}),\quad \theta_i\sim N(0,\sigma_\theta^2),\quad \phi\sim\text{ICAR}(W,\sigma_\phi).$$ pm.ICAR(W=W, sigma=sd_phi) is exactly the structured term (sum-to-zero enforced internally); sigma is the conditional scale $\sigma_\phi$. Priors: $\alpha,\beta\sim N(0,5)$, $\sigma_\theta,\sigma_\phi\sim\text{HalfNormal}(1)$.

In [1]:
import numpy as np, pymc as pm, pytensor.tensor as pt, arviz as az
from scotland_bym import scotland_data, bym_gibbs

O, E, x, adj, names = scotland_data('scotland_lip.csv'); n = len(O)
W = np.zeros((n, n))
for i, a in enumerate(adj): W[i, a] = 1            # binary adjacency matrix
assert (W == W.T).all()

with pm.Model() as m:
    alpha = pm.Normal('alpha', 0, 5); beta = pm.Normal('beta', 0, 5)
    sd_theta = pm.HalfNormal('sd_theta', 1.0); sd_phi = pm.HalfNormal('sd_phi', 1.0)
    # NON-CENTERED: unit-scale effects scaled by sd outside, to avoid the funnel (cf. Epil)
    theta_std = pm.Normal('theta_std', 0, 1, shape=n)
    phi_std = pm.ICAR('phi_std', W=W, sigma=1.0)
    theta = pm.Deterministic('theta', sd_theta * theta_std)
    phi = pm.Deterministic('phi', sd_phi * phi_std)
    eta = alpha + beta * x + theta + phi
    pm.Poisson('O', mu=E * pt.exp(eta), observed=O)
    idata = pm.sample(2000, tune=2000, chains=4, target_accept=0.99,
                      random_seed=1, progressbar=False)

po = idata.posterior; be = po['beta'].values.ravel()
mxr = float(az.summary(idata, var_names=['alpha','beta','sd_theta','sd_phi'])['r_hat'].max())
print(f"alpha = {float(po['alpha'].mean()):.3f}")
print(f"beta  = {be.mean():.3f}  95% CI [{np.percentile(be,2.5):.2f}, {np.percentile(be,97.5):.2f}]"
      f"   RR/10% = {np.exp(be.mean()):.2f}   P(beta>0) = {(be>0).mean():.3f}")
print(f"sd_theta = {float(po['sd_theta'].mean()):.3f}   sd_phi = {float(po['sd_phi'].mean()):.3f}   max r_hat {mxr:.3f}")

Initializing NUTS using jitter+adapt_diag...

Multiprocess sampling (4 chains in 4 jobs)

NUTS: [alpha, beta, sd_theta, sd_phi, theta_std, phi_std]

Sampling 4 chains for 2_000 tune and 2_000 draw iterations (8_000 + 8_000 draws total) took 44 seconds.

alpha = -0.221
beta  = 0.370  95% CI [0.10, 0.63]   RR/10% = 1.45   P(beta>0) = 0.997
sd_theta = 0.153   sd_phi = 0.702   max r_hat 1.000
In [2]:
# per-area smoothed RR: PyMC vs from-scratch
g = bym_gibbs(O, E, x, adj, R=20000, burn=5000, seed=1)
rr_g = np.exp(g['alpha'].mean()+g['beta'].mean()*x+g['theta'].mean(0)+g['phi'].mean(0))
a = po['alpha'].values.ravel(); b = po['beta'].values.ravel()
th = po['theta'].values.reshape(-1, n); ph = po['phi'].values.reshape(-1, n)
rr_p = np.exp(a[:,None] + b[:,None]*x[None,:] + th + ph).mean(0)
print('per-area RR  max |PyMC - from-scratch| = %.3f   correlation = %.4f'%(np.abs(rr_p-rr_g).max(), np.corrcoef(rr_p,rr_g)[0,1]))
print('%-22s %12s %12s'%('parameter','from-scratch','PyMC'))
print('%-22s %12.3f %12.3f'%('beta (AFF/10%)', g['beta'].mean(), b.mean()))
print('%-22s %12.3f %12.3f'%('sd_phi (spatial)', g['sd_phi'].mean(), float(po['sd_phi'].mean())))
print('%-22s %12.3f %12.3f'%('sd_theta (unstruct)', g['sd_theta'].mean(), float(po['sd_theta'].mean())))

per-area RR  max |PyMC - from-scratch| = 0.323   correlation = 0.9990
parameter              from-scratch         PyMC
beta (AFF/10%)                0.302        0.370
sd_phi (spatial)              0.728        0.702
sd_theta (unstruct)           0.032        0.153

Results¶

  • With the non-centered parameterisation PyMC samples cleanly (r̂ ≈ 1.00, no divergences / tree-depth warnings) and closely tracks the from-scratch fit on the spatial structure: σ_φ ≫ σ_θ (variation almost entirely spatial) and the per-area smoothed relative risks line up (correlation ≈ 0.999). The AFF covariate effect is β ≈ 0.37 (RR ≈ 1.45 per +10% AFF) — a little higher than the from-scratch Gibbs' β ≈ 0.30 (RR ≈ 1.35); the two agree on the sign and the P(β>0) ≈ 0.99 conclusion but differ in magnitude, because the weakly-identified β–φ confounding is resolved differently under each prior.
  • The centered version failed. Writing phi = pm.ICAR(W, sigma=sd_phi) directly creates the classic variance funnel (σ vs the effects it scales), and NUTS stalled — r̂ ≈ 1.09, σ_φ blew up to ~3.7, β became meaningless. Scaling unit-variance effects outside the prior (phi = sd_phi * phi_std) fixes it — the same non-centering trick the Epil GLMM needed. (Gibbs, sampling σ and the effects in alternation, never sees this funnel — which is why the from-scratch sampler had no trouble.)
  • The structured/unstructured split (σ_θ, σ_φ) differs slightly across engines — HalfNormal(1) here vs the vague Gamma-on-precision in the from-scratch sampler — because the θ-vs-φ decomposition is only weakly identified (the BYM identifiability issue). But the covariate effect, the spatial dominance, and the area risks all agree.
  • Two routes, one risk surface: explicit sequential Gibbs over the adjacency graph (from-scratch) vs PyMC's built-in ICAR + NUTS give the same smoothed map.