Hierarchical Poisson — PyMC cross-check

pm.Poisson with a non-centered patient random intercept

Companion to hpois_python.ipynb (from-scratch Gibbs) and hpois_R.ipynb (lme4::glmer). Same Poisson GLMM, fit with PyMC/NUTS on the same data (hpois_synth.csv, epil.csv).

Non-centered parameterisation. A direct b_j ~ N(0, σ_b²) creates the funnel geometry that cripples NUTS. We instead sample standard-normal $z_j$ and set $b_j = \sigma_b z_j$ — the standard fix for hierarchical models, giving a smooth posterior NUTS can explore efficiently.

import os
_lib = r"C:\Users\user\anaconda3\envs\pymc-env\Library\bin"
if _lib not in os.environ.get("PATH", ""): os.environ["PATH"] = _lib + os.pathsep + os.environ.get("PATH", "")
os.environ["KMP_DUPLICATE_LIB_OK"] = "TRUE"
os.environ.setdefault("PYTENSOR_FLAGS", "cxx=C:/Users/user/anaconda3/envs/pymc-env/Library/bin/g++.exe")
import numpy as np, pandas as pd, matplotlib.pyplot as plt, pymc as pm, arviz as az
import pymc.sampling.mcmc as _mcmc
class _Noop:
    def __init__(self,*a,**k): pass
    def __enter__(self): return self
    def __exit__(self,*a): pass
_mcmc.threadpool_limits = _Noop

def hpois_model(X, y, grp, J):
    with pm.Model() as m:
        beta = pm.Normal('beta', 0., 10., shape=X.shape[1])
        sigma_b = pm.HalfNormal('sigma_b', 1.0)
        z = pm.Normal('z', 0., 1., shape=J)                 # non-centered random intercepts
        b = pm.Deterministic('b', z * sigma_b)
        pm.Poisson('y', mu=pm.math.exp(pm.math.dot(X, beta) + b[grp]), observed=y)
    return m
print('ready')

ready

---

1. Synthetic validation

from hpois_gibbs import simulate_hpois
y, X, g, b_true = simulate_hpois(J=120, ni=6, beta=(0.5, 0.8, -0.4), sigma_b=0.7, seed=1)
with hpois_model(X, y, g, 120):
    idata_s = pm.sample(1000, tune=1000, chains=2, cores=1, target_accept=0.9, random_seed=1, progressbar=False)
print(az.summary(idata_s, var_names=['beta', 'sigma_b'])[['mean','sd','r_hat']].to_string())
print("\n(truth beta 0.5/0.8/-0.4, sigma_b 0.7;  from-scratch 0.525/0.811/-0.417, 0.679;  glmer 0.512/0.807/-0.415, 0.675)")

Initializing NUTS using jitter+adapt_diag...

Sequential sampling (2 chains in 1 job)

NUTS: [beta, sigma_b, z]

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

We recommend running at least 4 chains for robust computation of convergence diagnostics

           mean      sd r_hat
beta[0]   0.503   0.068  1.01
beta[1]   0.807  0.0274  1.00
beta[2]  -0.415  0.0242  1.00
sigma_b   0.685   0.057  1.00

(truth beta 0.5/0.8/-0.4, sigma_b 0.7;  from-scratch 0.525/0.811/-0.417, 0.679;  glmer 0.512/0.807/-0.415, 0.675)

1. three engines agree

PyMC: β = 0.503 / 0.807 / −0.415, σ_b = 0.685 (r̂ ≈ 1.00) — identical to the from-scratch Gibbs and lme4::glmer. The non-centered NUTS samples cleanly.

---

2. Epil: the treatment effect on the significance boundary

ep = pd.read_csv('epil.csv'); ye = ep['seizures'].values.astype(float)
Xe = np.column_stack([np.ones(len(ep)), np.log(ep['Base']/4), ep['Trt'], np.log(ep['Age']), ep['V4']])
ge = ep['patient'].values - 1
with hpois_model(Xe, ye, ge, 59):
    idata_e = pm.sample(1500, tune=1500, chains=2, cores=1, target_accept=0.95, random_seed=2, progressbar=False)
print(az.summary(idata_e, var_names=['beta', 'sigma_b'])[['mean','sd','r_hat']].to_string())
B = idata_e.posterior['beta'].values.reshape(-1, 5); tr = B[:, 2]
print(f"\nTrt: RR {np.exp(tr.mean()):.2f}  95% CI [{np.exp(np.percentile(tr,2.5)):.2f}, {np.exp(np.percentile(tr,97.5)):.2f}]"
      f"   P(RR<1) = {(tr<0).mean():.3f}")

Initializing NUTS using jitter+adapt_diag...

Sequential sampling (2 chains in 1 job)

NUTS: [beta, sigma_b, z]

C:\Users\user\anaconda3\envs\pymc-env\Lib\site-packages\pymc\step_methods\hmc\quadpotential.py:321: RuntimeWarning: overflow encountered in dot
  return 0.5 * np.dot(x, v_out)

Sampling 2 chains for 1_500 tune and 1_500 draw iterations (3_000 + 3_000 draws total) took 29 seconds.

We recommend running at least 4 chains for robust computation of convergence diagnostics

           mean     sd r_hat
beta[0]    -0.9   1.27  1.00
beta[1]   1.025  0.105  1.00
beta[2]   -0.32   0.16  1.00
beta[3]    0.29   0.36  1.00
beta[4]  -0.163  0.052  1.00
sigma_b   0.553  0.067  1.00

Trt: RR 0.73  95% CI [0.52, 0.99]   P(RR<1) = 0.978

# The hierarchical treatment effect across the three engines -- it straddles RR = 1
engines = [('from-scratch Gibbs', -0.297, 0.180),
           ('PyMC (NUTS)',        tr.mean(), tr.std()),
           ('lme4::glmer',        -0.315, 0.150)]
fig, ax = plt.subplots(1, 2, figsize=(12, 4.2))
for i, (lab, est, se) in enumerate(engines):
    incl = (est-1.96*se) < 0 < (est+1.96*se)
    ax[0].errorbar(np.exp(est), i, xerr=[[np.exp(est)-np.exp(est-1.96*se)], [np.exp(est+1.96*se)-np.exp(est)]],
                   fmt='o', capsize=4, color='firebrick' if incl else 'steelblue')
ax[0].axvline(1.0, color='gray', ls='--'); ax[0].set_yticks(range(3)); ax[0].set_yticklabels([e[0] for e in engines])
ax[0].invert_yaxis(); ax[0].set_xlabel('treatment rate ratio (95% interval)')
ax[0].set_title('Hierarchical Trt effect straddles RR=1 (borderline)')
# data vs fitted, by arm (PyMC hierarchical)
bpost = B.mean(0); bj = idata_e.posterior['b'].values.reshape(-1, 59).mean(0)
ep['fit'] = np.exp(Xe @ bpost + bj[ge])
obs = ep.groupby(['visit','Trt'])['seizures'].mean().unstack(1); fitm = ep.groupby(['visit','Trt'])['fit'].mean().unstack(1)
sem = ep.groupby(['visit','Trt'])['seizures'].sem().unstack(1)
v = obs.index.values; cols = {0:'steelblue', 1:'firebrick'}; lab = {0:'placebo', 1:'progabide'}
for t in (0, 1):
    ax[1].errorbar(v, obs[t], yerr=sem[t], fmt='o', color=cols[t], capsize=3, label=f'{lab[t]} obs')
    ax[1].plot(v, fitm[t], '--', color=cols[t], lw=2, label=f'{lab[t]} fitted')
ax[1].set_xlabel('visit'); ax[1].set_ylabel('mean seizures / 2 weeks'); ax[1].set_xticks(v)
ax[1].set_title('Data vs fitted (PyMC hierarchical), by arm'); ax[1].legend(fontsize=8)
plt.tight_layout(); plt.savefig('hpois_pymc.png', dpi=120, bbox_inches='tight'); plt.show()

2. Results: three engines, one borderline verdict

PyMC hierarchical fixed effects (matching from-scratch / glmer): log(Base/4) ≈ 1.03, Trt ≈ −0.32, log(Age) ≈ 0.29, V4 ≈ −0.16, σ_b ≈ 0.55.

The hierarchical treatment effect, across all three engines, lands right on the RR = 1 boundary:

engine	RR	95% interval	verdict
from-scratch Gibbs (Bayes)	0.74	[0.53, 1.07]	just includes 1
PyMC NUTS (Bayes)	0.73	[0.52, 0.99]	just excludes 1
lme4::glmer (freq., p=0.036)	0.73	[0.54, 0.98]	just excludes 1

All three agree on the point estimate (RR ≈ 0.73) and the inflated SE (~0.15–0.18 vs the single-level NegBin's 0.10) — but they fall on both sides of the 5% line. That is the honest conclusion: once within-patient correlation is modelled, the progabide effect is borderline / fragile, not the clear effect the single-level analysis reported. PyMC's P(RR<1) ≈ 0.98 says the same thing probabilistically — it just clears the conventional one-sided 0.975 threshold, consistent with the 95% interval that just excludes 1.

Conclusion

PyMC (non-centered NUTS), the from-scratch Gibbs, and lme4::glmer give the same hierarchical Poisson fit (β, σ_b, the data-vs-fit by arm). The hierarchical model's headline — correlation widens the treatment uncertainty to the significance boundary — is reproduced by all three; the fact that they split on the verdict is the lesson that the effect is genuinely marginal.

Probability of benefit (PyMC posterior)

The clinically meaningful summary: the overlap of the placebo vs progabide rate posteriors, and P(rate ratio < 1) — the probability progabide lowers seizures.

from scipy.stats import gaussian_kde
mlb = np.log(ep['Base']/4).mean(); mla = np.log(ep['Age']).mean()
p0 = np.exp(B[:,0] + B[:,1]*mlb + B[:,3]*mla); p1 = p0 * np.exp(B[:,2]); rr = np.exp(B[:,2])
P_red = (B[:,2] < 0).mean(); P_20 = (rr < 0.8).mean()
fig, ax = plt.subplots(1, 2, figsize=(12, 4.2))
xs = np.linspace(min(p0.min(),p1.min()), max(p0.max(),p1.max()), 300)
ax[0].fill_between(xs, gaussian_kde(p0)(xs), alpha=.45, color='steelblue', label='placebo')
ax[0].fill_between(xs, gaussian_kde(p1)(xs), alpha=.45, color='firebrick', label='progabide')
ax[0].set_xlabel('adjusted mean seizures / 2 weeks'); ax[0].set_yticks([])
ax[0].set_title('Posterior rate: placebo vs progabide (overlap)'); ax[0].legend()
xr = np.linspace(rr.min(), rr.max(), 300); dr = gaussian_kde(rr)(xr)
ax[1].plot(xr, dr, 'k', lw=1.5); ax[1].fill_between(xr, dr, where=xr<1, color='seagreen', alpha=.5, label=f'RR<1: P = {P_red:.2f}')
ax[1].axvline(1.0, color='gray', ls='--'); ax[1].set_xlabel('treatment rate ratio  exp(Trt)'); ax[1].set_yticks([])
ax[1].set_title('P(progabide lowers seizures)'); ax[1].legend()
plt.tight_layout(); plt.savefig('hpois_pymc_prob.png', dpi=120, bbox_inches='tight'); plt.show()
print(f'P(RR<1) = {P_red:.3f}   P(RR<0.8) = {P_20:.3f}   (from-scratch: 0.943 / 0.647)')

P(RR<1) = 0.978   P(RR<0.8) = 0.715   (from-scratch: 0.943 / 0.647)

PyMC gives P(RR<1) ≈ 0.98 (vs the from-scratch 0.94) — both say progabide is very likely beneficial even though the 95% interval grazes 1. The small method gap (0.94 vs 0.98) is the same boundary sensitivity that flipped the significance verdict; the probability-of-benefit framing makes the agreement obvious where the binary verdict made them look contradictory.