Seeds — R cross-check (lme4::glmer, with glm over-dispersion contrast)¶

Companion to seeds_python.ipynb / seeds_pymc.ipynb. glmer(..., family=binomial) fits the binomial GLMM; plain glm shows the over-dispersion the random effect addresses.

In [1]:
suppressMessages(library(lme4))
d <- read.csv('seeds.csv'); d$plate <- factor(d$plate)
gl <- glm(cbind(r, n-r) ~ x1*x2, family=binomial, data=d)              # plain logistic
cat(sprintf('plain logistic: residual deviance %.1f on %d df  -> dispersion %.2f\n',
            deviance(gl), df.residual(gl), sum(residuals(gl,'pearson')^2)/df.residual(gl)))
gm <- glmer(cbind(r, n-r) ~ x1*x2 + (1|plate), family=binomial, data=d)  # binomial GLMM
print(summary(gm)$coefficients)
cat(sprintf('\nrandom-effect sd (plate) = %.3f\n', sqrt(unlist(VarCorr(gm))[1])))

plain logistic: residual deviance 33.3 on 17 df  -> dispersion 1.86

               Estimate Std. Error    z value     Pr(>|z|)
(Intercept) -0.54848370  0.1660825 -3.3024769 9.583497e-04
x1           0.09742568  0.2773573  0.3512642 7.253902e-01
x2           1.33681076  0.2361777  5.6601911 1.512045e-08
x1:x2       -0.81004130  0.3841745 -2.1085247 3.498563e-02

random-effect sd (plate) = 0.235
In [2]:
# SEs: plain glm (too narrow) vs glmer (accounts for over-dispersion)
cat('coefficient SEs:  glm vs glmer\n')
se_glm <- summary(gl)$coefficients[,2]; se_glmer <- summary(gm)$coefficients[,2]
print(round(data.frame(glm=se_glm, glmer=se_glmer, ratio=se_glmer/se_glm), 3))

coefficient SEs:  glm vs glmer

              glm glmer ratio
(Intercept) 0.126 0.166 1.318
x1          0.223 0.277 1.243
x2          0.177 0.236 1.331
x1:x2       0.306 0.384 1.254

Results¶

  • Over-dispersion is present (plain-logistic Pearson dispersion > 1).
  • glmer matches the from-scratch and PyMC fixed effects: root extract ≈ +1.34, interaction ≈ −0.81, seed type ≈ 0; plate random-effect SD ≈ 0.2–0.3 (ML; small, as expected — and a touch larger/smaller than the Bayesian posterior mean depending on the variance prior, since σ is weakly identified with 21 plates).
  • The random effect widens the SEs vs plain glm (ratio > 1) — the inferential payoff of modelling over-dispersion.
  • Three engines agree (from-scratch ≈ PyMC ≈ glmer) on the binomial GLMM; the hierarchical-binomial completes the GLM family alongside the hierarchical-Poisson set.