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}")

g++ not available, if using conda: `conda install gxx`

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 43 seconds.

alpha = -0.219
beta  = 0.367  95% CI [0.10, 0.63]   RR/10% = 1.44   P(beta>0) = 0.996
sd_theta = 0.147   sd_phi = 0.713   max r_hat 1.000

# 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.317   correlation = 0.9991
parameter              from-scratch         PyMC
beta (AFF/10%)                0.302        0.367
sd_phi (spatial)              0.728        0.713
sd_theta (unstruct)           0.032        0.147

# ── Cross-engine comparison, drawn ────────────────────────────────────────────
import matplotlib.pyplot as plt
fig, ax = plt.subplots(1, 3, figsize=(16, 4.4))
GREY = '#cbd5e0'

# (1) per-area smoothed relative risk, engine against engine
lim = [min(rr_g.min(), rr_p.min()) * .9, max(rr_g.max(), rr_p.max()) * 1.05]
ax[0].plot(lim, lim, color=GREY, lw=2, zorder=1)
ax[0].scatter(rr_g, rr_p, s=26, color='steelblue', alpha=.75, lw=0, zorder=3)
worst = int(np.argmax(np.abs(rr_p - rr_g)))
ax[0].annotate(names[worst], (rr_g[worst], rr_p[worst]), fontsize=8, color='firebrick',
               xytext=(6, -10), textcoords='offset points')
ax[0].set_xscale('log'); ax[0].set_yscale('log'); ax[0].set_xlim(lim); ax[0].set_ylim(lim)
ax[0].set_xlabel('from-scratch Metropolis-within-Gibbs'); ax[0].set_ylabel('PyMC (pm.ICAR + NUTS)')
ax[0].set_title(f'Smoothed RR per district (r = {np.corrcoef(rr_p, rr_g)[0,1]:.4f})', fontsize=10)

# (2) the covariate effect: the thing both engines are asked for
ax[1].hist(g['beta'], bins=45, density=True, color='steelblue', alpha=.6, label='from-scratch Gibbs')
ax[1].hist(b, bins=45, density=True, color='firebrick', alpha=.5, label='PyMC NUTS')
ax[1].axvline(0, color='k', lw=1)
ax[1].set_xlabel(r'$\beta$  (AFF covariate, per +10%)'); ax[1].set_yticks([])
ax[1].legend(fontsize=8); ax[1].set_title('The covariate effect: agreement', fontsize=10)

# (3) the variance split: where they do NOT agree, and why that is expected
ax[2].hist(g['sd_phi'], bins=45, density=True, color='steelblue', alpha=.6, label=r'from-scratch $\sigma_\phi$')
ax[2].hist(po['sd_phi'].values.ravel(), bins=45, density=True, color='firebrick', alpha=.5,
           label=r'PyMC $\sigma_\phi$')
ax[2].hist(g['sd_theta'], bins=45, density=True, color='steelblue', alpha=.6, histtype='step',
           lw=1.8, label=r'from-scratch $\sigma_\theta$')
ax[2].hist(po['sd_theta'].values.ravel(), bins=45, density=True, color='firebrick', alpha=.9,
           histtype='step', lw=1.8, ls='--', label=r'PyMC $\sigma_\theta$')
ax[2].set_xlabel('posterior standard deviation'); ax[2].set_yticks([])
ax[2].legend(fontsize=7.5); ax[2].set_title(r'$\sigma_\phi$ agrees, $\sigma_\theta$ does not', fontsize=10)

fig.suptitle('BYM: two engines, one map — and one parameter they are entitled to disagree about', y=1.02)
plt.tight_layout(); plt.show()

print(f'{"":24s}{"from-scratch":>14s}{"PyMC":>10s}')
print('-' * 48)
for lab, gv, pv in [('beta (AFF per +10%)', g['beta'].mean(), b.mean()),
                    ('sd_phi (spatial)',    g['sd_phi'].mean(), float(po['sd_phi'].mean())),
                    ('sd_theta (unstruct.)',g['sd_theta'].mean(), float(po['sd_theta'].mean()))]:
    print(f'{lab:24s}{gv:>14.3f}{pv:>10.3f}')
print(f'\nper-area smoothed RR: correlation {np.corrcoef(rr_p, rr_g)[0,1]:.4f}, '
      f'largest gap {np.abs(rr_p-rr_g).max():.3f} (at {names[worst]})')
print(f'\nThe first two panels are the cross-check passing. The district risk surface is identical for')
print(f'practical purposes -- correlation {np.corrcoef(rr_p, rr_g)[0,1]:.4f} across all {n} districts. The covariate')
print(f'posteriors overlap heavily and both sit firmly above zero, but they are not on top of each')
print(f'other: PyMC is centred about {b.mean()-g["beta"].mean():.2f} higher. Against a posterior SD of roughly {b.std():.2f} that is')
print(f'well under half a standard deviation -- the ordinary gap between two independently written')
print(f'samplers with slightly different priors -- but it is a real offset, not the same number twice.')
print('The third panel is the interesting one, and it is not a failure. BYM splits the area variation')
print('into a spatial part and an unstructured part, and the data identify their SUM far better than')
print('either piece -- a district that is high can be high because its neighbours are, or high on its')
print('own, and the observed count cannot distinguish those. So sigma_phi, which carries almost all')
print('the signal, agrees across engines, while sigma_theta -- small, weakly identified, and pinned')
print('near zero from one direction by its prior -- comes out several times larger under NUTS than')
print('under the from-scratch sampler. Both agree it is negligible next to sigma_phi; they simply')
print('disagree about how negligible. That is the BYM identifiability issue, and it is why the')
print('conclusions reported from this model are about beta and the risk map rather than about the')
print('split itself.')

                          from-scratch      PyMC
------------------------------------------------
beta (AFF per +10%)              0.302     0.367
sd_phi (spatial)                 0.728     0.713
sd_theta (unstruct.)             0.032     0.147

per-area smoothed RR: correlation 0.9991, largest gap 0.317 (at Skye-Lochalsh)

The first two panels are the cross-check passing. The district risk surface is identical for
practical purposes -- correlation 0.9991 across all 56 districts. The covariate
posteriors overlap heavily and both sit firmly above zero, but they are not on top of each
other: PyMC is centred about 0.07 higher. Against a posterior SD of roughly 0.13 that is
well under half a standard deviation -- the ordinary gap between two independently written
samplers with slightly different priors -- but it is a real offset, not the same number twice.
The third panel is the interesting one, and it is not a failure. BYM splits the area variation
into a spatial part and an unstructured part, and the data identify their SUM far better than
either piece -- a district that is high can be high because its neighbours are, or high on its
own, and the observed count cannot distinguish those. So sigma_phi, which carries almost all
the signal, agrees across engines, while sigma_theta -- small, weakly identified, and pinned
near zero from one direction by its prior -- comes out several times larger under NUTS than
under the from-scratch sampler. Both agree it is negligible next to sigma_phi; they simply
disagree about how negligible. That is the BYM identifiability issue, and it is why the
conclusions reported from this model are about beta and the risk map rather than about the
split itself.

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

Model¶

Cross-engine comparison¶

Results¶

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

Model¶

Cross-engine comparison¶

Results¶

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