Hierarchical Poisson — R cross-check

lme4::glmer (Poisson GLMM) vs the from-scratch Gibbs sampler

Companion to hpois_python.ipynb (from-scratch Metropolis-within-Gibbs). We refit the same data with R's standard mixed-model tool, lme4::glmer(family = poisson) with a patient random intercept (1 | patient). Synthetic data is shared via hpois_synth.csv; the real data is epil.csv.

Frequentist vs Bayesian note. glmer uses Laplace-approximate maximum likelihood (point estimate of σ_b, asymptotic SEs); the from-scratch sampler is fully Bayesian (marginalises σ_b uncertainty). They agree closely, with a revealing difference at the significance boundary (see Section 2).

.libPaths(c("C:/Users/user/R/win-library/4.6", .libPaths()))
suppressMessages({library(lme4); library(MASS)})

# --- §1 synthetic: recover beta and sigma_b ---
d <- read.csv("hpois_synth.csv")
m <- glmer(y ~ x1 + x2 + (1 | group), data = d, family = poisson)
cat("glmer fixed effects (truth 0.5 / 0.8 / -0.4):\n"); print(round(fixef(m), 3))
cat(sprintf("sigma_b = %.3f   (truth 0.70;  from-scratch 0.68)\n", sqrt(VarCorr(m)$group[1])))

glmer fixed effects (truth 0.5 / 0.8 / -0.4):

(Intercept)          x1          x2 
      0.512       0.807      -0.415 

sigma_b = 0.675   (truth 0.70;  from-scratch 0.68)

1. synthetic recovery matches

glmer: β = 0.512 / 0.807 / −0.415 and σ_b = 0.675 — essentially identical to the from-scratch Gibbs (0.525 / 0.811 / −0.417, σ_b 0.679) and the truth. The two engines agree on the GLMM fit.

# --- §2 Epil: hierarchical Poisson, and the treatment effect across models ---
e <- read.csv("epil.csv"); e$lBase <- log(e$Base/4); e$lAge <- log(e$Age); e$patient <- factor(e$patient)
mp   <- glm(seizures ~ lBase + Trt + lAge + V4, family = poisson, data = e)            # single-level Poisson
mnb  <- glm.nb(seizures ~ lBase + Trt + lAge + V4, data = e)                            # single-level NegBin
mgl  <- glmer(seizures ~ lBase + Trt + lAge + V4 + (1 | patient), data = e, family = poisson)  # hierarchical

cat("glmer (hierarchical) fixed effects:\n"); print(round(summary(mgl)$coef[, 1:2], 3))
cat(sprintf("sigma_b = %.3f\n\n", sqrt(VarCorr(mgl)$patient[1])))
tab <- rbind(
  Poisson      = summary(mp)$coef['Trt', c(1,2,4)],
  NegBin       = summary(mnb)$coef['Trt', c(1,2,4)],
  glmer_hier   = summary(mgl)$coef['Trt', c(1,2,4)])
colnames(tab) <- c('Trt', 'SE', 'p'); tab <- cbind(tab, RR = exp(tab[,'Trt']))
cat("Treatment effect across models:\n"); print(round(tab, 3))

glmer (hierarchical) fixed effects:

            Estimate Std. Error
(Intercept)   -1.087      1.197
lBase          1.027      0.101
Trt           -0.315      0.150
lAge           0.332      0.343
V4            -0.160      0.054

sigma_b = 0.516

Treatment effect across models:

              Trt    SE     p    RR
Poisson    -0.017 0.048 0.727 0.983
NegBin     -0.233 0.100 0.019 0.792
glmer_hier -0.315 0.150 0.036 0.730

# treatment-effect forest: SE widens single-level -> hierarchical
est <- tab[, 'Trt']; se <- tab[, 'SE']; labs <- c('Poisson (1-level)','NegBin (1-level)','glmer (hierarchical)')
cols <- c('steelblue','steelblue','firebrick')
plot(est, 1:3, xlim=range(c(est-1.96*se, est+1.96*se, 0.1)), ylim=c(0.5,3.5), yaxt='n',
     pch=19, col=cols, xlab='Trt coefficient (95% CI)', ylab='', main='Epil: treatment effect by model')
arrows(est-1.96*se, 1:3, est+1.96*se, 1:3, angle=90, code=3, length=.06, col=cols)
abline(v=0, lty=3, col='gray'); axis(2, at=1:3, labels=labs, las=1, cex.axis=.8)
text(est, 1:3+0.18, sprintf('p=%.3f', tab[,'p']), cex=.8)

2. Epil: R confirms the result, and pins down the boundary

Hierarchical glmer fixed effects (matching the from-scratch Bayesian in parentheses):

coef	glmer (SE)	from-scratch Bayes (SD)
log(Base/4)	1.027 (0.101)	1.13 (0.10)
Trt	−0.315 (0.150)	−0.30 (0.18)
log(Age)	0.332 (0.343)	0.20 (0.45)
V4	−0.160 (0.054)	−0.16 (0.06)
σ_b	0.516	0.576

Treatment effect across the three models:

model	Trt	SE	RR	p / 95% CI
Poisson (single-level)	−0.017	0.048	0.98	p ≈ 0.7
NegBin (single-level)	−0.233	0.100	0.79	p ≈ 0.02 (sig)
glmer (hierarchical)	−0.315	0.150	0.73	p = 0.036 (borderline)

The story is the same as the from-scratch notebook: adding the patient random intercept widens the treatment SE (0.10 → 0.15) and pushes the effect from clearly significant toward the boundary.

The instructive nuance — frequentist vs Bayesian disagree at the margin: glmer's Laplace ML gives p = 0.036 (just significant), while the full-Bayes posterior gives RR 0.74 with 95% CrI [0.53, 1.07] (just not significant). The difference is exactly the extra uncertainty the Bayesian carries — it marginalises over σ_b (and small-cluster, 4-visits-per-patient effects that the Laplace approximation under-states). Robust conclusion: once within-patient correlation is modelled the progabide effect is fragile — sitting right on the 5% boundary, no longer the clear effect the single-level NegBin reported. Whether you call it "significant" depends on the inference framework, which is itself the cautionary lesson.

Conclusion

lme4::glmer and the from-scratch Gibbs agree on the GLMM fit (β, σ_b, the per-visit V4 effect). Both show the hierarchical model inflates the treatment SE ~50% vs single-level. The treatment effect lands on the significance boundary — frequentist marginally in, Bayesian marginally out — a textbook demonstration that modelling correlation changes conclusions, and that the verdict here is genuinely borderline.