Seeds — binomial GLMM (logistic regression with over-dispersion)

From-scratch Metropolis-within-Gibbs; the binomial twin of Epil

Module: seeds_gibbs.py. BUGS Vol I "Seeds" (Crowder 1978; Breslow & Clayton 1993). Grouped binomial counts — seeds germinated per plate — that vary more than a single logistic curve allows (over-dispersion). A plate random effect absorbs that extra variation. This completes the hierarchical-GLM grid: hierarchical Gaussian (Rats), hierarchical Poisson (Epil/Pumps/Rugby), and now hierarchical binomial.

The data — Orobanche seed germination (Crowder 1978)

21 plates in a 2×2 factorial: two Orobanche seed types (O. aegyptiaca 75 vs 73) crossed with two root extracts (bean vs cucumber). Each plate reports $r_i$ germinated out of $n_i$ planted.

field	meaning
r	seeds germinated on the plate
n	seeds planted (the binomial denominator; 4–81)
x1	seed type (0 = Oa75, 1 = Oa73)
x2	root extract (0 = bean, 1 = cucumber)

Germination rates range 0–83%, and plates in the same treatment cell still disagree more than binomial sampling predicts — that extra-binomial spread is the over-dispersion the random effect models.

Model & algorithm

$$r_i\sim\text{Binomial}(n_i,p_i),\quad \text{logit}\,p_i=\beta_0+\beta_1 x_{1i}+\beta_2 x_{2i}+\beta_{12}x_{1i}x_{2i}+b_i,\quad b_i\sim N(0,\sigma^2).$$ $b_i$ is the plate random effect — the logistic-scale analogue of the Negative Binomial's gamma mixing. Metropolis-within-Gibbs: (1) $\beta$ block RW-Metropolis (binomial-logistic likelihood, proposal from the pooled logistic Fisher info); (2) $b_i$ per-plate RW-Metropolis (curvature-scaled), re-centred into $\beta_0$; (3) $\sigma^2$ Inverse-Gamma. Self-contained — does not reuse the older hierarchical-Poisson code.

import numpy as np, pandas as pd, matplotlib.pyplot as plt
from scipy.special import expit
from seeds_gibbs import seeds_data, seeds_gibbs, simulate_seeds, summary, _logit_mle

# ---- Section 1: synthetic recovery ----
r, n, X, b = simulate_seeds(I=400, beta=(-0.5, 0.1, 1.3, -0.8), sigma=0.35, seed=1)
o = seeds_gibbs(r, n, X, seed=1)
print('SYNTHETIC (truth beta -0.5/0.1/1.3/-0.8, sigma 0.35):')
for nm, m, s, lo, hi in summary(o['beta'], ['b0','b1','b2','b12']): print('  %-4s %.3f (%.3f)' % (nm, m, s))
print('  sigma %.3f   (acc beta %.2f, b %.2f)' % (o['sigma'].mean(), o['accept_beta'], o['accept_b']))

SYNTHETIC (truth beta -0.5/0.1/1.3/-0.8, sigma 0.35):
  b0   -0.533 (0.046)
  b1   0.092 (0.070)
  b2   1.347 (0.074)
  b12  -0.842 (0.103)
  sigma 0.393   (acc beta 0.31, b 0.57)

# ---- Section 2: the Seeds data; plain logistic (over-dispersed) vs the GLMM ----
r, n, X = seeds_data()
# plain logistic (no random effect) + Pearson over-dispersion
bp = _logit_mle(X, r, n); ph = expit(X @ bp)
pearson = np.sum((r - n*ph)**2 / (n*ph*(1-ph))) / (len(r) - X.shape[1])
W = (X * (n*ph*(1-ph))[:, None]).T @ X; se_glm = np.sqrt(np.diag(np.linalg.inv(W)))
print('plain logistic Pearson dispersion = %.2f  (>1 => over-dispersed)\n' % pearson)

o = seeds_gibbs(r, n, X, R=30000, burn=8000, seed=1)
nm = ['intercept','seed type','root extract','interaction']
print('%-13s %18s %16s' % ('', 'GLMM (Bayes)', 'plain logistic'))
for j, (lab, m, s, lo, hi) in enumerate(summary(o['beta'], nm)):
    print('%-13s %7.3f [%.2f,%.2f]   %7.3f (%.3f)' % (lab, m, lo, hi, bp[j], se_glm[j]))
print('\nrandom-effect sigma = %.3f  95%% CI [%.2f, %.2f]' % (o['sigma'].mean(), np.percentile(o['sigma'],2.5), np.percentile(o['sigma'],97.5)))
print('(canonical: intercept -0.55, seed-type ~0.1, root-extract +1.34, interaction -0.81)')

plain logistic Pearson dispersion = 1.86  (>1 => over-dispersed)

                    GLMM (Bayes)   plain logistic
intercept      -0.560 [-0.85,-0.28]    -0.558 (0.126)
seed type       0.115 [-0.40,0.60]     0.146 (0.223)
root extract    1.344 [0.93,1.78]     1.318 (0.177)
interaction    -0.802 [-1.53,-0.13]    -0.778 (0.306)

random-effect sigma = 0.134  95% CI [0.02, 0.45]
(canonical: intercept -0.55, seed-type ~0.1, root-extract +1.34, interaction -0.81)

fig, ax = plt.subplots(1, 2, figsize=(12, 4.6))
# (A) observed germination by plate within each 2x2 cell + GLMM-fitted cell probability
cell = (X[:,1]*2 + X[:,2]).astype(int)                       # 0..3 = (type,extract) cells
labs = ['Oa75+bean','Oa75+cuc','Oa73+bean','Oa73+cuc']; bm = o['beta'].mean(0)
rng = np.random.default_rng(0)
for c in range(4):
    m = cell == c; jit = c + 0.12*rng.standard_normal(m.sum())
    ax[0].scatter(jit, (r[m]/n[m]), s=20+n[m], color='steelblue', alpha=.6, edgecolor='k', lw=.3)
    x1c, x2c = c//2, c%2; pcell = expit(bm[0]+bm[1]*x1c+bm[2]*x2c+bm[3]*x1c*x2c)
    ax[0].plot([c-0.3, c+0.3], [pcell, pcell], 'firebrick', lw=2)
ax[0].set_xticks(range(4)); ax[0].set_xticklabels(labs, fontsize=8, rotation=15)
ax[0].set_ylabel('germination rate $r/n$'); ax[0].set_title('Plate rates by cell (size ∝ n); red = GLMM fit')
# (B) coefficient forest: GLMM (wider, honest) vs plain logistic (too narrow under over-dispersion)
yy = np.arange(4)
gm = o['beta'].mean(0); glo = np.percentile(o['beta'],2.5,0); ghi = np.percentile(o['beta'],97.5,0)
ax[1].errorbar(gm[1:], yy[1:]+0.12, xerr=[gm[1:]-glo[1:], ghi[1:]-gm[1:]], fmt='o', color='steelblue', capsize=3, label='GLMM (Bayes)')
ax[1].errorbar(bp[1:], yy[1:]-0.12, xerr=1.96*se_glm[1:], fmt='s', color='darkorange', capsize=3, label='plain logistic')
ax[1].axvline(0, color='firebrick', ls='--', lw=1); ax[1].set_yticks(yy[1:]); ax[1].set_yticklabels(nm[1:], fontsize=8)
ax[1].set_xlabel('coefficient (log-odds)'); ax[1].set_title('GLMM vs plain logistic'); ax[1].legend(fontsize=8)
plt.tight_layout(); plt.savefig('seeds.png', dpi=120, bbox_inches='tight'); plt.show()

Further graphs: justify, interpret, and check the model

• Over-dispersion diagnostic — Pearson residuals from the plain logistic; plates beyond ±2 are the extra-binomial spread the random effect is there to absorb.
• The interaction on the probability scale — germination probability vs root extract for each seed type; the non-parallel lines are the −0.81 interaction (cucumber helps Oa75 much more than Oa73) — far more readable than log-odds.
• Posterior-predictive check — the observed dispersion sits in the tail of the plain-logistic null (inadequate) but is reproduced by the GLMM.
• Plate random effects — the per-plate $b_i$ spread, the hierarchical signature.

# --- why a random effect (over-dispersion diagnostic) + the interaction on the probability scale ---
fig, ax = plt.subplots(1, 2, figsize=(12, 4.6))
# (A) Pearson residuals from the plain logistic vs fitted prob; |resid|>2 = extra-binomial spread
ph = expit(X @ bp); pres = (r - n*ph) / np.sqrt(n*ph*(1-ph))
cellc = (X[:,1]*2 + X[:,2]).astype(int); cols = ['steelblue','seagreen','goldenrod','firebrick']
ax[0].scatter(ph, pres, c=[cols[c] for c in cellc], s=25+n, edgecolor='k', lw=.3, zorder=3)
ax[0].axhline(0, color='grey', lw=.7); ax[0].axhline(2, color='firebrick', ls='--', lw=1); ax[0].axhline(-2, color='firebrick', ls='--', lw=1)
ax[0].set_xlabel('plain-logistic fitted probability'); ax[0].set_ylabel('Pearson residual')
ax[0].set_title('Over-dispersion: %d/%d plates beyond ±2  (dispersion %.2f)' % (int((np.abs(pres)>2).sum()), len(r), pearson))
# (B) interaction on the probability scale, with 95% CrI
B = o['beta']
def cellp(x1, x2): return expit(B[:,0] + B[:,1]*x1 + B[:,2]*x2 + B[:,3]*x1*x2)
for x1, lab, col in [(0,'Oa75','steelblue'), (1,'Oa73','firebrick')]:
    ms = [cellp(x1,x2).mean() for x2 in (0,1)]
    lo = [np.percentile(cellp(x1,x2),2.5) for x2 in (0,1)]; hi = [np.percentile(cellp(x1,x2),97.5) for x2 in (0,1)]
    ax[1].errorbar([0,1], ms, yerr=[np.array(ms)-lo, np.array(hi)-ms], marker='o', color=col, capsize=4, lw=2, label=lab)
ax[1].set_xticks([0,1]); ax[1].set_xticklabels(['bean','cucumber']); ax[1].set_xlim(-0.3, 1.3)
ax[1].set_xlabel('root extract'); ax[1].set_ylabel('germination probability')
ax[1].set_title('Interaction: extract effect depends on seed type'); ax[1].legend(title='seed type', fontsize=8)
plt.tight_layout(); plt.savefig('seeds_effects.png', dpi=120, bbox_inches='tight'); plt.show()

# --- posterior-predictive check (does the random effect fix over-dispersion?) + plate random effects ---
fig, ax = plt.subplots(1, 2, figsize=(12, 4.6))
rng = np.random.default_rng(1); I = len(r); kdim = X.shape[1]; ph = expit(X @ bp)
def disp(rr): return np.sum((rr - n*ph)**2 / (n*ph*(1-ph))) / (I - kdim)   # Pearson dispersion vs plain-logistic fit
T_null = np.array([disp(rng.binomial(n.astype(int), ph).astype(float)) for _ in range(3000)])      # plain-logistic replicates
B = o['beta']; SG = o['sigma']; idx = rng.integers(0, len(B), 3000)
T_glmm = np.array([disp(rng.binomial(n.astype(int), expit(X @ B[j] + rng.normal(0, SG[j], I))).astype(float)) for j in idx])
ax[0].hist(T_null, bins=40, density=True, alpha=.6, color='darkorange', label='plain logistic')
ax[0].hist(T_glmm, bins=40, density=True, alpha=.6, color='steelblue', label='GLMM')
ax[0].axvline(pearson, color='k', ls='--', lw=1.5, label=f'observed {pearson:.2f}')
ax[0].set_xlabel('Pearson dispersion (replicated)'); ax[0].set_yticks([])
ax[0].set_title(f'PPC: P(disp>=obs) null={ (T_null>=pearson).mean():.3f}, GLMM={(T_glmm>=pearson).mean():.2f}'); ax[0].legend(fontsize=8)
# (D) plate random-effect caterpillar
BRE = o['bdraws']; bmn = BRE.mean(0); blo = np.percentile(BRE,2.5,0); bhi = np.percentile(BRE,97.5,0); order = np.argsort(bmn)
ax[1].errorbar(bmn[order], range(I), xerr=[bmn[order]-blo[order], bhi[order]-bmn[order]], fmt='o', ms=3, lw=.6, color='steelblue')
ax[1].axvline(0, color='firebrick', ls='--', lw=1)
ax[1].set_xlabel('plate random effect $b_i$ (log-odds)'); ax[1].set_ylabel('plate (sorted)')
ax[1].set_title('Plate effects  (sd $\\approx$ %.2f)' % o['sigma'].mean())
plt.tight_layout(); plt.savefig('seeds_ppc.png', dpi=120, bbox_inches='tight'); plt.show()

Reading the posterior-predictive check

disp is the Pearson dispersion (≈ 1 if the data are truly binomial; observed = 1.86, i.e. over-dispersed). For each model we simulate many replicate datasets from that model and ask what fraction are at least as dispersed as the real data — the posterior-predictive p-value P(disp ≥ obs):

• Plain logistic — P = 0.067: the observed 1.86 lands in the model’s upper tail; a pure binomial simply cannot scatter this much → misfit (small p).
• GLMM — P = 0.36: the plate random effect $b_i$ injects the extra spread, so 1.86 sits centrally in the replicate distribution → adequate (non-extreme p).

A well-calibrated model gives a non-extreme value (roughly 0.1–0.9); p near 0 (or 1) flags systematic misfit. It is a one-sided upper tail (over-dispersion is a one-sided concern), and posterior-predictive p-values are mildly conservative (the data inform both the fit and the check), so the plain-logistic failure is, if anything, understated. Passing this check is what justifies trusting the GLMM’s wider coefficient CIs — ignoring the over-dispersion would overstate significance (the Breslow–Clayton point).

Results & interpretation

• The counts are over-dispersed: the plain logistic Pearson dispersion is >1 — plates within the same cell scatter more than binomial sampling allows (left panel shows the spread).
• Fixed effects (GLMM, matching the canonical Seeds analysis): root extract strongly raises germination (β ≈ +1.34, cucumber > bean), the interaction is clearly negative (β ≈ −0.81 — the extract effect is weaker for seed type Oa73), and seed type alone is ~0. These match glmer / BUGS.
• The random effect matters for inference, not the point estimates: modelling $\sigma$ widens the coefficient intervals vs the plain logistic (right panel) — ignoring over-dispersion would overstate significance (the Breslow–Clayton point, exactly as in Epil).
• σ is small and weakly identified (≈ 0.1–0.3, wide CI): the over-dispersion here is mild but real. (The synthetic case, where σ is well-identified, recovers it cleanly — confirming the sampler.)

Takeaways

• Binomial over-dispersion = a random effect on the logit, the exact analogue of the Poisson GLMM (Epil). Same Metropolis-within-Gibbs machinery, binomial likelihood swapped in.
• Cross-checks: PyMC (pm.Binomial + non-centred random effect) and R (lme4::glmer, with plain glm for the over-dispersion contrast) follow.