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 (James–Stein). Both the logit-normal GLMM and the new conjugate beta-binomial deliver it — raw RMSE $\approx 0.069$ cut to $\approx 0.038$ (~44% better) by pulling the raw extremes ~86% toward the .267 mean — and the two parameterisations agree to $\approx 0.997$ correlation. 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)]