Bayesian Hierarchical Binomial

Model

A binomial GLMM — logistic regression on grouped counts with a group random effect that absorbs over-dispersion. Grouped binomial counts (e.g. seeds germinated per plate) routinely vary more than a single logistic curve allows; a per-group random effect $b_i$ on the logit soaks up that extra-binomial spread. It is the exact binomial twin of the hierarchical-Poisson Epil model — the logistic-scale analogue of the Negative Binomial's gamma mixing — and completes the hierarchical-GLM grid: hierarchical Gaussian (Rats), hierarchical Poisson (Epil/Pumps/Rugby), and now hierarchical binomial. This is the BUGS Vol I Seeds example (Crowder 1978; Breslow & Clayton 1993).

• $r_i, n_i$ — successes and trials on group (plate) $i$
• $b_i$ — group random effect on the logit; absorbs extra-binomial over-dispersion
• $x_i$ — 2×2 factorial design (seed type, root extract, interaction); $\sigma$ — over-dispersion scale

Sampler — Metropolis-within-Gibbs

The same Metropolis-within-Gibbs structure as the hierarchical Poisson, with the binomial-logistic likelihood swapped in: a block RW-Metropolis for the fixed effects $\beta$ (proposal from the pooled logistic Fisher information), a vectorised per-plate RW-Metropolis for the random effects $b_i$ (curvature-scaled, re-centred into $\beta_0$), and a conjugate Inverse-Gamma draw for $\sigma^2$.

A conjugate alternative — beta-binomial

The logit-normal GLMM above is the non-conjugate route (a normal on the logit, sampled by Metropolis). The beta-binomial is its conjugate twin: put a Beta prior directly on the success probability and the per-unit posterior is exact, so only the two hyperparameters — population mean $\mu = \alpha/(\alpha+\beta)$ and concentration $\kappa = \alpha+\beta$ — need Metropolis (the "almost-pure Gibbs" structure of Pumps). Marginally $r_i \sim$ Beta-Binomial — the binomial cousin of the Negative Binomial. This is the exact binomial mirror of the count-model conjugacy contrast (gamma-Poisson/Pumps vs Poisson-GLMM/Epil); the baseball example fits both and shows they give the same shrinkage. Module: betabin_gibbs.py.

Three applications

The same hierarchical-binomial idea applied three ways — a covariate model, an institutional-comparison model, and a two-parameterisation shrinkage study:

Seeds

The Orobanche germination experiment (Crowder, 1978): 21 plates in a 2×2 factorial — two seed types crossed with two root extracts — each reporting $r_i$ germinated of $n_i$ planted. Plates in the same cell disagree more than binomial sampling allows, so a per-plate random effect on the logit absorbs the extra-binomial spread. This is the covariate application: the random effect changes inference on the fixed effects — modelling $\sigma$ widens the coefficient intervals, so ignoring the over-dispersion would overstate significance (the Breslow–Clayton point).

Surgical

Infant cardiac-surgery mortality across 12 hospitals with wildly uneven volumes (47–810 operations). An intercept-only version of the same model (design matrix = a column of ones) for the classic "hospital league-table" problem. Partial pooling stabilises each rate — a raw 0/47 = 0% is pulled up to a believable ~6%, and high-volume hospitals barely move — and, crucially, the credible intervals overlap heavily, so the data do not support a confident ranking. This is the institutional-comparison application.

Baseball

Efron & Morris (1975): 18 players' batting averages after only 45 at-bats, used to predict their rest-of-season average. Shrinking the raw averages toward the group mean beats them decisively — the textbook James–Stein / Efron–Morris result. It fits the model two ways — the conjugate beta-binomial and the logit-normal GLMM — and shows the two parameterisations give the same shrinkage. This is the two-parameterisation study that ties the conjugate and non-conjugate routes together.

Notebooks

Three applications share the same hierarchical-binomial idea — a covariate model, an institutional-comparison model, and a two-parameterisation shrinkage study. Seeds — over-dispersion in a designed experiment (BUGS Seeds). The Orobanche germination data (Crowder 1978): 21 plates in a 2×2 factorial — two seed types crossed with two root extracts — each reporting $r_i$ germinated of $n_i$ planted, with plates in the same cell disagreeing more than binomial sampling predicts. The fixed effects match the canonical analysis: root extract strongly raises germination ($\beta\approx +1.34$), the interaction is clearly negative ($\beta\approx -0.81$), and seed type alone is $\approx 0$. The Breslow–Clayton point: modelling $\sigma$ widens the coefficient intervals versus a plain logistic — ignoring over-dispersion would overstate significance. Surgical — shrinkage & ranking (BUGS Surgical). An intercept-only version of the same model (design matrix = a column of ones) for the classic "hospital league table" problem: infant cardiac-surgery deaths in 12 hospitals with very uneven volumes (47–810 operations). Partial pooling stabilises the rates — hospital 1's raw 0/47 = 0% is pulled to a believable ~5.8%, the high-volume hospital barely moves — and, crucially, the credible intervals overlap heavily: the data do not support a confident ranking. It is the binomial, logit-scale analogue of the conjugate Pumps Poisson shrinkage. Baseball — two parameterisations of shrinkage (Efron–Morris 1975). The famous result: shrinking 18 players' 45-at-bat batting averages toward the group mean predicts their rest-of-season average far better than the raw averages. This is the James–Stein result, and it was genuinely shocking when it appeared: for three or more quantities estimated at once, pulling every estimate toward the common mean beats using each unit's own data, always, no matter what the true values are. It is the frequentist face of exactly the shrinkage a hierarchical model produces. Both the logit-normal GLMM and the new conjugate beta-binomial deliver it — raw RMSE $0.069$ cut to $0.038$ (~44% better) by pulling the raw extremes ~86% toward the .267 mean — and the two parameterisations agree to $\approx 0.997$ correlation.

The notebook also reports the number that makes the result honest, and it deserves saying out loud. Alongside the raw averages and the shrunken ones it scores complete pooling — simply predicting that every player will hit the group mean, discarding each player's own record entirely. That scores $0.0397$, against partial pooling's $0.0384$. So of the whole $0.031$ improvement over the raw averages, $0.029$ comes from declining to trust 45 at-bats at all, and only $0.0013$ from the refinement of trusting them partially. This does not undo the James–Stein point, it sharpens it: the estimator wins because 45 at-bats carry very little signal about a batting average, and the sophisticated part of the machinery is claiming a small share of the credit. On data where the units were better measured, the balance would shift the other way — which is exactly what the shrinkage factor of 86% is reporting. Every fit is cross-checked three ways — from-scratch Gibbs, PyMC (pm.Binomial / pm.BetaBinomial), and R (lme4::glmer, plus glm/VGAM contrasts).

Downloads

References

• Efron, B. & Morris, C. (1975). Data analysis using Stein's estimator and its generalizations. Journal of the American Statistical Association 70(350), 311–319. — the baseball batting-average shrinkage (James–Stein) example
• Crowder, M. J. (1978). Beta-binomial ANOVA for proportions. Journal of the Royal Statistical Society: Series C 27(1), 34–37. — the Orobanche Seeds germination data analyzed here
• Skene, A. M. & Wakefield, J. C. (1990). Hierarchical models for multicentre binary response studies. Statistics in Medicine 9(8), 919–929. — hierarchical binomial for the Surgical multi-hospital comparison
• 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 binomial GLMM and why over-dispersion widens coefficient intervals
• Spiegelhalter, D. J., Thomas, A., Best, N. & Gilks, W. R. (1996). BUGS Examples, Volume 1. MRC Biostatistics Unit. — the canonical Seeds and Surgical reference fits
• 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 binomial models, partial pooling, and the beta-binomial hyperprior

GLMM Sampler — Source Code

"""
Seeds: binomial GLMM (logistic regression with group random effects + over-dispersion) -- from-scratch
Metropolis-within-Gibbs. BUGS Vol I "Seeds" (Crowder 1978; Breslow & Clayton 1993). The BINOMIAL twin of
the hierarchical-Poisson Epil model -- a plate random effect absorbs extra-binomial variation.

Model
-----
  r_i ~ Binomial(n_i, p_i),   logit(p_i) = x_i' beta + b_i,   b_i ~ N(0, sigma^2),  i = 1..I plates
  x_i = [1, seed-type, root-extract, seed-type*extract]  (a 2x2 factorial with interaction)
The plate random effect b_i captures OVER-DISPERSION: grouped binomial counts vary more than a single
logistic curve allows; b_i is the logistic-scale analogue of the NegBin's gamma or the Poisson-GLMM's
log-normal mixing. Self-contained (does not reuse/modify the older hierarchical-Poisson module).

Sampler (Metropolis-within-Gibbs)
---------------------------------
  1. beta   -- block RW-Metropolis (binomial-logistic loglik); proposal N(beta, s^2 (I_pool+prior_prec)^-1)
               from the pooled logistic Fisher info X' diag(n p(1-p)) X (computed once). s = 2.38/sqrt(k).
  2. b_i    -- per-plate RW-Metropolis, vectorised; step SD = s_b / sqrt(curvature_i),
               curvature_i = n_i p_i(1-p_i) + 1/sigma^2.  + centre b -> beta_0 (location move).
  3. sigma^2 -- Inverse-Gamma conjugate from {b_i}.
"""
import numpy as np
import pandas as pd
from scipy.special import expit


def seeds_data(csv='seeds.csv'):
    d = pd.read_csv(csv)
    X = np.column_stack([np.ones(len(d)), d.x1, d.x2, d.x1 * d.x2])   # 1, type, extract, interaction
    return d['r'].to_numpy(float), d['n'].to_numpy(float), X


def simulate_seeds(I=300, beta=(-0.5, 0.1, 1.3, -0.8), sigma=0.35, nbar=40, seed=0):
    rng = np.random.default_rng(seed); beta = np.asarray(beta, float)
    x1 = rng.integers(0, 2, I).astype(float); x2 = rng.integers(0, 2, I).astype(float)
    X = np.column_stack([np.ones(I), x1, x2, x1 * x2])
    b = rng.normal(0, sigma, I); p = expit(X @ beta + b)
    n = rng.integers(nbar // 2, nbar * 2, I).astype(float)
    r = rng.binomial(n.astype(int), p).astype(float)
    return r, n, X, b


def _logit_mle(X, r, n, maxit=100):
    beta = np.zeros(X.shape[1])
    for _ in range(maxit):
        p = expit(X @ beta); w = n * p * (1 - p)
        step = np.linalg.solve((X * w[:, None]).T @ X + 1e-8 * np.eye(X.shape[1]), X.T @ (r - n * p))
        beta += step
        if np.max(np.abs(step)) < 1e-10:
            break
    return beta


def seeds_gibbs(r, n, X, R=20000, burn=5000, prior_sd=10.0, a0=0.001, b0=0.001,
                seed=0, s_beta=None, s_b=1.6):
    rng = np.random.default_rng(seed); I, k = X.shape
    s_beta = (2.38 / np.sqrt(k)) if s_beta is None else s_beta
    bp = _logit_mle(X, r, n); p0 = expit(X @ bp)
    Ipool = (X * (n * p0 * (1 - p0))[:, None]).T @ X
    Lb = np.linalg.cholesky(np.linalg.inv(Ipool + np.eye(k) / prior_sd ** 2))

    beta = bp.copy(); b = np.zeros(I); sigma = 0.3
    keep = R - burn
    B = np.zeros((keep, k)); SIG = np.zeros(keep); Bm = np.zeros(I); BRE = np.zeros((keep, I))
    accb = 0.0; acc_beta = 0
    for it in range(R):
        # 1. beta  (block RW-Metropolis)
        eta = X @ beta + b
        ll = np.sum(r * eta - n * np.logaddexp(0, eta)) - 0.5 * np.sum(beta ** 2) / prior_sd ** 2
        prop = beta + s_beta * (Lb @ rng.standard_normal(k)); ep = X @ prop + b
        llp = np.sum(r * ep - n * np.logaddexp(0, ep)) - 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_i  (per-plate RW-Metropolis, vectorised, curvature-scaled)
        eta = eta0 + b; p = expit(eta)
        curv = n * p * (1 - p) + 1.0 / sigma ** 2
        prop_b = b + s_b * rng.standard_normal(I) / np.sqrt(curv)
        ep = eta0 + prop_b
        d = r * (prop_b - b) - n * (np.logaddexp(0, ep) - np.logaddexp(0, eta)) \
            - 0.5 * (prop_b ** 2 - b ** 2) / sigma ** 2
        ok = np.log(rng.random(I)) < d; b = np.where(ok, prop_b, b); accb += ok.mean()
        m = b.mean(); b -= m; beta[0] += m                       # centre random effects -> intercept
        # 3. sigma^2  (Inverse-Gamma)
        sigma = np.sqrt(1.0 / rng.gamma(a0 + I / 2.0, 1.0 / (b0 + 0.5 * np.sum(b ** 2))))
        if it >= burn:
            B[it - burn] = beta; SIG[it - burn] = sigma; Bm += b; BRE[it - burn] = b
    return dict(beta=B, sigma=SIG, b=Bm / keep, bdraws=BRE,
                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)]