Bayesian Hierarchical Poisson

Model

The hierarchical extension of the count regression example: a Poisson GLMM in which each group $j$ carries its own random effect on the log-rate, absorbing within-group correlation and extra-Poisson heterogeneity. The random intercept is a log-normal mixing of the Poisson rate — the hierarchical cousin of the Negative Binomial's gamma mixing (both explain overdispersion; NB marginally, the GLMM through explicit per-unit effects). It is the right model for clustered counts (repeated measures per unit), and the model behind the Breslow–Clayton (1993) / BUGS Epil analysis.

• $b_j$ — group $j$'s random intercept; absorbs within-group correlation and overdispersion
• $g_i$ — group (patient / pump / team) of observation $i$
• $\sigma_b^2$ — random-effect variance; Inverse-Gamma prior

Three applications across the conjugacy spectrum

Three applications span the conjugacy spectrum — the same hierarchical-Poisson idea, three different random-effect structures and samplers:

Pumps

Ten pumps at a nuclear power station, each with a failure count and an exposure time — the canonical BUGS example. Putting a Gamma prior on every pump's failure rate keeps the model fully conjugate: the per-pump rates and the population scale are drawn exactly, and only the Gamma shape needs a Metropolis step. It is the cleanest illustration of shrinkage — the noisiest, lowest-exposure pumps borrow the most strength from the population, while the high-exposure pumps are left near their raw rates.

Epil

A randomised trial of the anticonvulsant progabide (Thall & Vail, 1990): 59 patients, each contributing four repeated seizure counts. A per-patient random intercept absorbs the correlation between a patient's own visits — the feature that makes this the headline case, because once that correlation is modelled the treatment effect moves onto the significance boundary that the naïve single-level fit had comfortably cleared. Non-conjugate (a Normal on the log-rate, needed for covariates), so the coefficients and random effects use Metropolis-within-Gibbs.

Rugby

The 2014 Six Nations — 15 matches (Baio & Blangiardo, 2010). Every team carries two crossed random effects, an attack strength and a defence strength, which enter each match through different team indices, so the effects are non-nested. It shows the same machinery handling crossed rather than strictly hierarchical grouping, and recovers the final tournament standings — home advantage, England's attack, Ireland's defence — from just 15 games.

Sampler — Metropolis-within-Gibbs

Because the Poisson log-likelihood is non-conjugate and (unlike the logit) has no Pólya–Gamma augmentation, the general GLMM uses Metropolis-within-Gibbs: block RW-Metropolis for the fixed effects $\beta$, vectorised per-group RW-Metropolis for the random effects, an exact location/translation move to de-confound the level ($\beta_0$ vs. the $b_j$), and a conjugate Inverse-Gamma draw for $\sigma_b^2$. The Pumps model is the exception: putting the random effect as a gamma on the rate (not a normal on the log-rate) restores conjugacy, so 2 of its 3 blocks are exact Gamma draws and only the shape $\alpha$ needs Metropolis.

Notebooks

Epil (the headline). Refitting the epilepsy trial with a patient random intercept properly accounts for the 4 correlated visits per patient. This overturns the single-level conclusion: the treatment SE widens from 0.10 to ~0.15–0.18 and the progabide rate-ratio interval moves onto the significance boundary (RR $\approx 0.73$, 95% interval grazing 1) — the frequentist glmer lands just inside ($p=0.036$), the full-Bayes posterior just outside, with P(reduction) $\approx 0.94\text{–}0.98$. Modelling correlation changed the conclusion, not just the standard error — the classic Breslow–Clayton lesson; the single-level "significant" effect was pseudo-replication. Pumps (the conjugate counterpoint). The BUGS gamma–Poisson example (10 pumps): $\theta_i \sim \text{Gamma}(\alpha,\beta)$ on the rate makes the per-pump rates and the scale exact Gibbs draws — only the shape needs one Metropolis step. Low-exposure pumps shrink most toward the population mean rate (posterior mean $\alpha/\beta \approx 0.94$; the ratio of point estimates $0.70/0.93 \approx 0.75$). This is the explicit-random-effects form of the Negative Binomial: marginalising $\theta_i$ gives the NB, so R fits it by empirical Bayes via glm.nb — no sampler. Rugby (crossed effects). The Baio & Blangiardo (2010) Six-Nations model: each team carries two crossed random effects (attack & defence) entering every match through different team indices, with a sum-to-zero translation move for identifiability. Recovers the 2014 standings from 15 matches — home advantage $\approx 0.37$ ($\times 1.4$), England strongest attack, Ireland best defence. Each application is cross-checked three ways — from-scratch sampler, PyMC NUTS, and R (lme4::glmer or the NB marginal) — all in agreement.

Downloads

References

• Breslow, N. E. & Clayton, D. G. (1993). Approximate inference in generalized linear mixed models. Journal of the American Statistical Association 88(421), 9–25. — the foundational Poisson-GLMM paper; the Epil random-intercept analysis reproduced here
• Thall, P. F. & Vail, S. C. (1990). Some covariance models for longitudinal count data with overdispersion. Biometrics 46(3), 657–671. — the progabide epilepsy trial (the epil data)
• Baio, G. & Blangiardo, M. (2010). Bayesian hierarchical model for the prediction of football results. Journal of Applied Statistics 37(2), 253–264. — the crossed attack/defence random-effects model applied to the Rugby data
• Spiegelhalter, D. J., Thomas, A., Best, N. & Gilks, W. R. (1996). BUGS Examples, Volume 1. MRC Biostatistics Unit. — the conjugate gamma–Poisson Pumps example
• Gelman, A., Carlin, J. B., Stern, H. S., Dunson, D. B., Vehtari, A. & Rubin, D. B. (2013). Bayesian Data Analysis (3rd ed.). CRC Press. — hierarchical models, the gamma–Poisson conjugacy, and shrinkage
• Roberts, G. O., Gelman, A. & Gilks, W. R. (1997). Weak convergence and optimal scaling of random walk Metropolis algorithms. Annals of Applied Probability 7(1), 110–120. — the RW-Metropolis tuning behind the Metropolis-within-Gibbs blocks

GLMM Sampler — Source Code

"""
Hierarchical Poisson regression (Poisson GLMM with group random intercepts) -- from-scratch.
Metropolis-within-Gibbs. The hierarchical extension of count_reg_mcmc.py; the model behind the
Breslow & Clayton (1993) / BUGS "Epil" analysis.

Model
-----
  y_i ~ Poisson(mu_i),   log mu_i = x_i' beta + b_{g_i},   b_j ~ N(0, sigma_b^2),  j = 1..J groups
  g_i = group (e.g. patient) of observation i.

The group random intercept b_j absorbs within-group correlation and the residual (extra-Poisson)
heterogeneity -- a *log-normal* mixing of the Poisson rate (cf. the Negative Binomial's *gamma* mixing).
It is the right model when counts are clustered (repeated measures per unit).

Sampler (Metropolis-within-Gibbs)
---------------------------------
  1. beta     -- block RW-Metropolis; proposal N(beta, s^2 (I_pool + prior_prec)^{-1}) scaled by the
                 pooled Poisson Fisher information (fixed), s = 2.38/sqrt(k).
  2. b_j      -- per-group RW-Metropolis, vectorised; proposal SD = s_b / sqrt(curvature_j),
                 curvature_j = sum_{i in j} mu_i + 1/sigma_b^2 (recomputed each sweep).
  3. sigma_b^2 -- Inverse-Gamma conjugate from {b_j}.

Prior: beta ~ N(0, prior_sd^2 I);  sigma_b^2 ~ IG(a0, b0).
"""
import numpy as np


def simulate_hpois(J=60, ni=5, beta=(0.5, 0.8, -0.4), sigma_b=0.7, seed=0):
    """Simulate clustered Poisson counts. x1 obs-level, x2 group-level (constant within group)."""
    rng = np.random.default_rng(seed)
    beta = np.asarray(beta, float); k = len(beta)
    group = np.repeat(np.arange(J), ni); n = J * ni
    x1 = rng.normal(size=n)                                  # observation-level covariates
    x2 = rng.normal(size=n)
    X = np.column_stack([np.ones(n), x1, x2])[:, :k]
    b = rng.normal(0, sigma_b, J)
    y = rng.poisson(np.exp(X @ beta + b[group]))
    return y, X, group, b


def _pois_mle(X, y, maxit=100, tol=1e-10):
    beta = np.zeros(X.shape[1])
    for _ in range(maxit):
        mu = np.exp(X @ beta)
        step = np.linalg.solve((X * mu[:, None]).T @ X, X.T @ (y - mu)); beta += step
        if np.max(np.abs(step)) < tol:
            break
    return beta


def hpois_gibbs(y, X, group, R=8000, burn=2000, prior_sd=10.0, a0=2.0, b0=1.0,
                seed=0, s_beta=None, s_b=1.6):
    rng = np.random.default_rng(seed)
    y = np.asarray(y, float); n, k = X.shape
    group = np.asarray(group, int); J = group.max() + 1
    sy = np.bincount(group, weights=y, minlength=J)         # per-group sum of y (fixed)
    s_beta = (2.38 / np.sqrt(k)) if s_beta is None else s_beta

    bpool = _pois_mle(X, y)
    Ipool = (X * np.exp(X @ bpool)[:, None]).T @ X
    Lb = np.linalg.cholesky(np.linalg.inv(Ipool + np.eye(k) / prior_sd ** 2))

    beta = bpool.copy(); b = np.zeros(J); sigma_b = 1.0
    eta0 = X @ beta
    keep = R - burn
    B = np.zeros((keep, k)); SIG = np.zeros(keep); Bm = np.zeros(J); accb = 0; acc_beta = 0
    for r in range(R):
        # 1. beta | b  (block RW-Metropolis)
        cur = X @ beta + b[group]
        ll = np.sum(y * cur - np.exp(cur)) - 0.5 * np.sum(beta ** 2) / prior_sd ** 2
        prop = beta + s_beta * (Lb @ rng.standard_normal(k))
        cp = X @ prop + b[group]
        llp = np.sum(y * cp - np.exp(cp)) - 0.5 * np.sum(prop ** 2) / prior_sd ** 2
        if np.log(rng.random()) < llp - ll:
            beta = prop; acc_beta += 1
        eta0 = X @ beta
        # 2. b_j | beta, sigma_b  (per-group RW-Metropolis, vectorised)
        mu = np.exp(eta0 + b[group])
        curv = np.bincount(group, weights=mu, minlength=J) + 1.0 / sigma_b ** 2
        prop_b = b + s_b * rng.standard_normal(J) / np.sqrt(curv)
        mu_p = np.exp(eta0 + prop_b[group])
        d = sy * (prop_b - b) - np.bincount(group, weights=(mu_p - mu), minlength=J) \
            - 0.5 * (prop_b ** 2 - b ** 2) / sigma_b ** 2
        a_ok = np.log(rng.random(J)) < d
        b = np.where(a_ok, prop_b, b); accb += a_ok.mean()
        # 2b. location move: the likelihood depends only on beta_0 + b_j, so the level mixes slowly
        #     between them. Draw the shared level delta ~ N(mean(b), sigma_b^2/J) and move it into
        #     beta_0 (vague prior), keeping b mean-zero. Exact Gibbs step; fixes the confounding.
        delta = rng.normal(b.mean(), sigma_b / np.sqrt(J))
        beta[0] += delta; b -= delta; eta0 = X @ beta
        # 3. sigma_b^2 | b  (Inverse-Gamma)
        sigma_b = np.sqrt(1.0 / rng.gamma(a0 + J / 2.0, 1.0 / (b0 + 0.5 * np.sum(b ** 2))))
        if r >= burn:
            B[r - burn] = beta; SIG[r - burn] = sigma_b; Bm += b
    return dict(beta=B, sigma_b=SIG, b=Bm / keep,
                accept_beta=acc_beta / R, accept_b=accb / R)


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