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 — a per-group adjustment that is itself treated as a draw from a common distribution rather than as a free parameter, which is what makes the groups inform each other instead of being fitted in isolation — 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 — the Poisson forces the variance to equal the mean, and real counts are almost always more variable than that, so the excess has to be modelled by something). 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 the Negative Binomial's 0.10 to 0.15–0.18, and the progabide rate-ratio interval moves onto the significance boundary. (The comparison is against the NB rather than the plain Poisson deliberately: the single-level Poisson SE on this data is 0.048, and its own overdispersion has to be fixed before the correlation question can even be posed — see count regression.) All three engines put the rate ratio at $0.73\text{–}0.74$ and all three land within a whisker of the boundary, but on different sides of it: glmer just inside ($p=0.036$, SE 0.150), PyMC just inside (95% CrI $[0.52,\ 0.99]$, P(reduction) $0.978$), and the from-scratch Gibbs just outside ($[0.53,\ 1.07]$, P(reduction) $0.943$, SE 0.180). That spread across engines is the finding — an effect whose significance depends on which of three near-identical fits you happen to run is not an effect the data has settled. Modelling correlation changed the conclusion, not just the standard error — the classic Breslow–Clayton lesson; the single-level significant effect was pseudo-replication — counting 4 visits from one patient as 4 independent pieces of evidence when they are largely 1, which inflates the apparent sample size and shrinks the standard errors to match. 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, an attack strength and a defence strength, which enter every match through different team indices. Crossed rather than nested: each effect does not sit inside one group, because every team meets every other, so a single match draws on two teams' effects at once. A sum-to-zero translation move keeps them identified. 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)]