Toenail — R cross-check (lme4::glmer, and the quadrature-point caveat)¶

Companion to toenail_python.ipynb / toenail_pymc.ipynb. The mixed logistic via glmer(y ~ treat*month + (1|patient), family=binomial). This dataset is Lesaffre & Spiessens' (2001) example that the estimated σ depends on the number of Gauss-Hermite quadrature points (nAGQ).

In [1]:
suppressMessages(library(lme4))
d <- read.csv('toenail.csv'); d$patient <- factor(d$patient)
cat('sigma vs number of quadrature points (nAGQ) -- the Lesaffre-Spiessens point:\n')
for (q in c(1, 5, 20, 50)) {
  gm <- glmer(y ~ treat*month + (1|patient), family=binomial, data=d, nAGQ=q)
  cat(sprintf('  nAGQ=%2d:  sigma = %.2f   month = %.3f   treat:month = %.3f\n',
              q, sqrt(unlist(VarCorr(gm))[1]), fixef(gm)[['month']], fixef(gm)[['treat:month']]))
}

sigma vs number of quadrature points (nAGQ) -- the Lesaffre-Spiessens point:

  nAGQ= 1:  sigma = 4.56   month = -0.400   treat:month = -0.137
  nAGQ= 5:  sigma = 3.69   month = -0.382   treat:month = -0.134
  nAGQ=20:  sigma = 4.00   month = -0.391   treat:month = -0.137
  nAGQ=50:  sigma = 4.01   month = -0.391   treat:month = -0.137
In [2]:
gm <- glmer(y ~ treat*month + (1|patient), family=binomial, data=d, nAGQ=50)
print(round(summary(gm)$coefficients, 3))
cat(sprintf('\nrandom-intercept sd (nAGQ=50) = %.2f\n', sqrt(unlist(VarCorr(gm))[1])))

            Estimate Std. Error z value Pr(>|z|)
(Intercept)   -1.618      0.434  -3.726    0.000
treat         -0.161      0.584  -0.275    0.783
month         -0.391      0.044  -8.810    0.000
treat:month   -0.137      0.068  -2.011    0.044

random-intercept sd (nAGQ=50) = 4.01

Results¶

  • glmer matches the Bayesian fits (at enough quadrature points): month ≈ −0.39, treat × month ≈ −0.14, σ ≈ 4 — agreeing with the from-scratch PG-Gibbs and PyMC.
  • The Lesaffre–Spiessens caveat is visible: the estimate of σ is numerically unstable at low quadrature and only settles near 4.0 once nAGQ ≥ 20. Here it swings from 4.56 at nAGQ=1 (the Laplace approximation) to 3.69 at nAGQ=5 before converging to 4.00–4.01 at nAGQ=20–50 — over-shooting first, then under-shooting, not a clean one-directional bias. The direction is specific to this optimiser and dataset (Lesaffre & Spiessens' original software under-estimated σ monotonically with too few points); what is general, and what matters, is that a large random-effect variance is genuinely hard to integrate with few nodes, so the low-nAGQ numbers cannot be trusted. This is why the exact MCMC samplers (PG-Gibbs, NUTS), which never approximate that integral, are reassuring here — both land at σ ≈ 4.1, consistent with the well-resolved glmer fit.
  • All engines agree on the science: onycholysis declines strongly over time, terbinafine somewhat faster, with large between-patient heterogeneity.