Hierarchical Poisson regression from scratch — Poisson GLMM¶
Group random intercepts; the Breslow–Clayton / BUGS "Epil" model¶
From-scratch module: hpois_gibbs.py. The hierarchical extension of count_reg_mcmc.py: a Poisson model with a patient random intercept that absorbs within-patient correlation and extra-Poisson heterogeneity. Companion engines (PyMC, R glmer) follow in separate notebooks.
| Section | Content |
|---|---|
| Model & algorithm | Poisson GLMM; Metropolis-within-Gibbs |
| Section 1 | Synthetic validation — recover β, σ_b, and the random effects |
| Section 2 | Epil — how the patient random effect changes the treatment inference |
Model¶
$$y_i \sim \text{Poisson}(\mu_i),\qquad \log\mu_i = x_i'\beta + b_{g_i},\qquad b_j \sim N(0,\sigma_b^2)$$
for groups (patients) $j=1,\dots,J$ with $g_i$ the group of observation $i$. The random intercept $b_j$ lets each unit have its own baseline rate; it soaks up the within-group correlation of repeated counts and the residual over-dispersion. Statistically it is a log-normal mixing of the Poisson rate — the hierarchical cousin of the Negative Binomial's gamma mixing (both explain over-dispersion; NB marginally, the GLMM through explicit per-unit effects).
Priors: $\beta \sim N(0,\sigma_p^2 I)$, $\sigma_b^2 \sim \text{IG}(a_0,b_0)$.
MCMC algorithm — Metropolis-within-Gibbs¶
The Poisson log-likelihood is non-conjugate and (unlike the logit) has no Pólya–Gamma augmentation, so we cycle through blocks, sampling each full conditional exactly where conjugate and by Metropolis otherwise. One sweep:
- $\beta \mid b$ — block random-walk Metropolis. Conditional on $b$, β is a Poisson regression with offsets $b_{g_i}$. Propose $\beta^*=\beta+s\,Lz$ with $L=\text{chol}[(\mathcal I_{\text{pool}}+\Sigma_p^{-1})^{-1}]$ (the pooled Poisson Fisher information, computed once) and $s=2.38/\sqrt{k}$; accept on the Poisson log-likelihood + prior. (acc ≈ 0.3)
- $b_j \mid \beta,\, \sigma_b$ — per-group RW-Metropolis, vectorised. Each $b_j$'s conditional ∝ (group-$j$ Poisson likelihood)×$N(0,\sigma_b^2)$. Propose $b_j^*=b_j+(s_b/\sqrt{c_j})z_j$ with curvature $c_j=\sum_{i\in j}\mu_i+1/\sigma_b^2$; accept/reject all $J$ at once. (acc ≈ 0.57)
- location move — exact Gibbs. The likelihood depends only on $\beta_0+b_j$, so the level mixes slowly between them. Draw $\delta\sim N(\bar b,\sigma_b^2/J)$ and set $\beta_0{+}{=}\delta$, $b{-}{=}\delta$ (hierarchical centering); keeps $b$ mean-zero and de-correlates $\beta_0$ from the random effects.
- $\sigma_b^2 \mid b$ — exact Gibbs (conjugate): $\sigma_b^2\sim\text{IG}(a_0+J/2,\ b_0+\tfrac12\sum_j b_j^2)$.
import numpy as np, pandas as pd, matplotlib.pyplot as plt
from hpois_gibbs import simulate_hpois, hpois_gibbs
from count_reg_mcmc import poisson_mle, negbin_mle # single-level fits for comparison
print('hpois_gibbs loaded')
hpois_gibbs loaded
1. Synthetic validation¶
$J=120$ groups × 6 observations, $\beta=(0.5,0.8,-0.4)$, random-intercept SD $\sigma_b=0.7$. We check recovery of β, σ_b, and the per-group effects $b_j$.
y, X, g, b_true = simulate_hpois(J=120, ni=6, beta=(0.5, 0.8, -0.4), sigma_b=0.7, seed=1)
fit = hpois_gibbs(y, X, g, R=8000, burn=3000, seed=2)
bm, bsd = fit['beta'].mean(0), fit['beta'].std(0)
print(f"{'param':>8}{'truth':>8}{'post.mean':>11}{'post.sd':>9}")
for j in range(3):
print(f"{'b'+str(j):>8}{(0.5,0.8,-0.4)[j]:>8.2f}{bm[j]:>11.3f}{bsd[j]:>9.3f}")
print(f"{'sigma_b':>8}{0.7:>8.2f}{fit['sigma_b'].mean():>11.3f}{fit['sigma_b'].std():>9.3f}")
print(f"\nrandom-effect recovery: corr(true b_j, estimated b_j) = {np.corrcoef(b_true, fit['b'])[0,1]:.2f}")
print(f"acceptance: beta {fit['accept_beta']:.2f}, b {fit['accept_b']:.2f}")
param truth post.mean post.sd
b0 0.50 0.525 0.071
b1 0.80 0.811 0.028
b2 -0.40 -0.417 0.026
sigma_b 0.70 0.679 0.052
random-effect recovery: corr(true b_j, estimated b_j) = 0.92
acceptance: beta 0.33, b 0.57
fig, ax = plt.subplots(1, 2, figsize=(12, 4.2))
truth = [0.5, 0.8, -0.4]
for j in range(3):
ax[0].hist(fit['beta'][:, j], bins=50, density=True, alpha=0.6, label=f'b{j}')
ax[0].axvline(truth[j], color='k', ls='--', lw=1)
ax[0].set_title('β posteriors (dashed = truth)'); ax[0].legend(fontsize=8); ax[0].set_yticks([])
ax[1].scatter(b_true, fit['b'], s=14, color='steelblue')
lo, hi = b_true.min(), b_true.max(); ax[1].plot([lo, hi], [lo, hi], 'k:', lw=1)
ax[1].set_xlabel('true random effect $b_j$'); ax[1].set_ylabel('posterior mean $\\hat b_j$')
ax[1].set_title(f'Random-effect recovery (r={np.corrcoef(b_true, fit["b"])[0,1]:.2f})')
plt.tight_layout(); plt.savefig('hpois_synth.png', dpi=120, bbox_inches='tight'); plt.show()
1. Results¶
| param | truth | post. mean (sd) |
|---|---|---|
| b0 | 0.50 | 0.525 (0.071) |
| b1 | 0.80 | 0.811 (0.028) |
| b2 | −0.40 | −0.417 (0.026) |
| σ_b | 0.70 | 0.679 (0.05) |
The sampler recovers the fixed effects, the random-intercept SD (σ_b ≈ 0.68 vs 0.70), and the 120 per-group effects (correlation 0.92 with truth). Acceptance ≈ 0.33 (β) and 0.57 (b) — healthy. The location move keeps the intercept unbiased (without it β₀ and the $b_j$ trade off the level and mix badly).
2. Epil: does the patient random effect change the treatment conclusion?¶
Same model as count_reg_python.ipynb Section 3 — seizures ~ log(Base/4) + Trt + log(Age) + V4 — but now with a patient random intercept $b_j$. The single-level NegBin treated the 236 visits as independent; the random effect properly accounts for the 4 correlated visits per patient.
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
enames = ['intercept', 'log(Base/4)', 'Trt', 'log(Age)', 'V4']
he = hpois_gibbs(ye, Xe, ge, R=10000, burn=4000, seed=3)
print(f"{'coef':>12}{'post.mean':>11}{'post.sd':>9}")
for j, nm in enumerate(enames):
print(f"{nm:>12}{he['beta'].mean(0)[j]:>11.3f}{he['beta'].std(0)[j]:>9.3f}")
sb = he['sigma_b']
print(f"{'sigma_b':>12}{sb.mean():>11.3f}{sb.std():>9.3f} 95% CI [{np.percentile(sb,2.5):.2f}, {np.percentile(sb,97.5):.2f}]")
tr = he['beta'][:, 2]
print(f"\nTREATMENT (progabide): {tr.mean():+.3f} SD {tr.std():.3f} "
f"RR {np.exp(tr.mean()):.2f} 95% CI [{np.exp(np.percentile(tr,2.5)):.2f}, {np.exp(np.percentile(tr,97.5)):.2f}]")
print(f"between-patient heterogeneity: exp(sigma_b) = {np.exp(sb.mean()):.2f}x multiplicative SD")
coef post.mean post.sd
intercept -0.836 1.509
log(Base/4) 1.133 0.100
Trt -0.297 0.180
log(Age) 0.195 0.445
V4 -0.159 0.055
sigma_b 0.576 0.065 95% CI [0.46, 0.72]
TREATMENT (progabide): -0.297 SD 0.180 RR 0.74 95% CI [0.53, 1.07]
between-patient heterogeneity: exp(sigma_b) = 1.78x multiplicative SD
# Treatment effect across the three models (single-level Poisson, single-level NegBin, hierarchical)
bp, sep, _ = poisson_mle(Xe, ye)
bn, lrn, sen = negbin_mle(Xe, ye)
models = [('Poisson (single-level)', bp[2], sep[2]),
('NegBin (single-level)', bn[2], sen[2]),
('Poisson GLMM (hierarchical)', tr.mean(), tr.std())]
fig, ax = plt.subplots(1, 2, figsize=(12, 4.2))
yy = np.arange(3)
for i, (lab, est, se) in enumerate(models):
ax[0].errorbar(est, i, xerr=1.96*se, fmt='o', capsize=4,
color='firebrick' if 'hier' in lab else 'steelblue')
ax[0].axvline(0, color='gray', ls=':'); ax[0].set_yticks(yy); ax[0].set_yticklabels([m[0] for m in models])
ax[0].invert_yaxis(); ax[0].set_xlabel('Trt coefficient (95% interval)')
ax[0].set_title('Treatment effect: SE widens, CI crosses 0')
# data vs fitted, by arm (hierarchical fitted mean uses beta + b_j)
ep['fit'] = np.exp(Xe @ he['beta'].mean(0) + he['b'][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 (hierarchical), by arm'); ax[1].legend(fontsize=8)
plt.tight_layout(); plt.savefig('hpois_epil.png', dpi=120, bbox_inches='tight'); plt.show()
print('single-level Trt: Poisson %.3f(%.3f) NegBin %.3f(%.3f)' % (bp[2],sep[2],bn[2],sen[2]))
print('hierarchical Trt: %.3f(%.3f)' % (tr.mean(), tr.std()))
single-level Trt: Poisson -0.017(0.048) NegBin -0.233(0.101) hierarchical Trt: -0.297(0.180)
2. Results — the hierarchical model overturns the single-level conclusion¶
| coef | post. mean (sd) |
|---|---|
| intercept | −0.84 (1.5) |
| log(Base/4) | 1.13 (0.10) |
| Trt | −0.30 (0.18) |
| log(Age) | 0.20 (0.45) |
| V4 | −0.16 (0.06) |
| σ_b | 0.58 [0.46, 0.72] |
The treatment effect across the three models (Trt coefficient, 95% interval):
| model | Trt | rate ratio | 95% CI (RR) | significant? |
|---|---|---|---|---|
| Poisson (single-level) | −0.02 | 0.98 | — | no (but mis-specified) |
| NegBin (single-level) | −0.23 | 0.79 | [0.65, 0.96] | yes |
| Poisson GLMM (hierarchical) | −0.30 | 0.74 | [0.53, 1.07] | borderline — CI just includes 1 |
This is the key result. Accounting for the within-patient correlation (the 4 visits per patient) via the random intercept widens the treatment SE from 0.10 to ~0.18, and the 95% credible interval for the rate ratio now spans 1.0 — so the progabide effect is no longer clearly significant. It is borderline: a frequentist glmer puts it right at p ≈ 0.04 (hpois_R.ipynb), while the fuller Bayesian uncertainty tips the CI just past 1. Either way the clear single-level result is gone. The point estimate is still a ~26% reduction, but the data are equally consistent with a 47% reduction or a 7% increase. The single-level NegBin's "significant" effect was an artefact of pseudo-replication (treating 236 correlated visits as 236 independent observations). This reproduces the classic Breslow & Clayton (1993) finding.
The random effect is large and real: σ_b ≈ 0.58 ⇒ patients differ by a factor of $e^{0.58}\approx1.8$ in baseline rate beyond what Base/Age explain — this is the over-dispersion the NegBin captured as a single dispersion parameter, here resolved into per-patient effects. The data-vs-fitted plot shows the hierarchical fit tracks each arm closely once these effects are included.
Takeaways¶
- The Poisson GLMM is the right model for clustered counts: it absorbs over-dispersion and correlation, where single-level Poisson/NegBin handle at most one.
- Modelling correlation can change conclusions, not just standard errors — here it pushes a "clearly significant" single-level effect onto the significance boundary (borderline either way).
- Cross-checks (PyMC
glmer-style model, Rlme4::glmer) and the connection to the Negative Binomial follow in the companion notebooks.
2. Probability of benefit (a better summary than "significant")¶
Rather than the binary 5% verdict, the posterior answers the clinical question directly: how probable is it that progabide lowers seizures, and by how much? Below: (left) the posterior densities of the adjusted mean rate for a typical patient under placebo vs progabide — their overlap is the uncertainty; (right) the posterior of the rate ratio $e^{\text{Trt}}$, with the area below 1 = P(reduction).
from scipy.stats import gaussian_kde
B = he['beta'] # from-scratch Epil posterior draws
mlb = np.log(ep['Base']/4).mean(); mla = np.log(ep['Age']).mean()
p0 = np.exp(B[:,0] + B[:,1]*mlb + B[:,3]*mla) # placebo mean rate (typical patient)
p1 = p0 * np.exp(B[:,2]) # progabide
rr = np.exp(B[:,2]) # rate ratio
P_red = (B[:,2] < 0).mean() # P(RR < 1)
P_20 = (rr < 0.8).mean() # P(>=20% reduction)
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_epil_prob.png', dpi=120, bbox_inches='tight'); plt.show()
print(f'P(progabide reduces seizures) P(RR<1) = {P_red:.3f}')
print(f'P(at least a 20% reduction) P(RR<0.8) = {P_20:.3f}')
P(progabide reduces seizures) P(RR<1) = 0.943 P(at least a 20% reduction) P(RR<0.8) = 0.647
Reading it. Even though the 95% interval grazes RR = 1, the posterior puts ~94% probability on some reduction (P(RR<1) ≈ 0.94) and a solid probability on a clinically meaningful (≥20%) reduction. This is the constructive way to report a borderline effect: not "fails significance," but "likely beneficial, with non-trivial uncertainty" — and it makes the placebo/treatment overlap explicit. (The left panel's two clouds clearly overlap, which is exactly why the verdict is borderline.)