Count regression — R benchmark¶
glm(family=poisson) and MASS::glm.nb vs the from-scratch sampler / PyMC¶
Companion to count_reg_python.ipynb (from-scratch MLE + RW-Metropolis) and count_reg_pymc.ipynb. We refit the same data with R's standard GLM tools to confirm the three engines agree.
Parameter mapping: MASS::glm.nb reports theta = the NB size parameter = our r = PyMC's alpha (Var = μ + μ²/θ). Synthetic data are shared via count_synth_pois.csv / count_synth_nb.csv; the real data is epil.csv.
.libPaths(c("C:/Users/user/R/win-library/4.6", .libPaths()))
suppressMessages(library(MASS))
# --- §1 Poisson (equidispersed synthetic) ---
d1 <- read.csv("count_synth_pois.csv")
m1 <- glm(y ~ x1 + x2, family = poisson, data = d1)
cat("Poisson GLM (R) vs from-scratch MLE/Bayes & PyMC:\n")
print(round(summary(m1)$coef[, 1:2], 3))
cat("\n(from-scratch/PyMC: 0.531/0.807/-0.387, SE 0.019/0.014/0.014)\n")
Poisson GLM (R) vs from-scratch MLE/Bayes & PyMC:
Estimate Std. Error (Intercept) 0.531 0.019 x1 0.807 0.014 x2 -0.387 0.014
(from-scratch/PyMC: 0.531/0.807/-0.387, SE 0.019/0.014/0.014)
1. Poisson agrees exactly¶
R's glm(poisson) gives intercept 0.531 (0.019), x1 0.807 (0.014), x2 −0.387 (0.014) — identical to the from-scratch IRLS/MLE, the from-scratch RW-Metropolis, and PyMC NUTS. (Same likelihood; IRLS is the same Newton step.)
# --- §2 Overdispersed synthetic: Poisson dispersion + glm.nb ---
d2 <- read.csv("count_synth_nb.csv")
mp <- glm(y ~ x1 + x2, family = poisson, data = d2)
disp <- sum(residuals(mp, "pearson")^2) / mp$df.residual
mnb <- glm.nb(y ~ x1 + x2, data = d2)
cat(sprintf("Pearson dispersion (Poisson) = %.2f (>>1)\n\n", disp))
cat("glm.nb coefficients (R):\n"); print(round(summary(mnb)$coef[, 1:2], 3))
cat(sprintf("\nglm.nb theta (= r = alpha) = %.2f [from-scratch r 2.10, PyMC 2.09, truth 2.0]\n", mnb$theta))
Pearson dispersion (Poisson) = 2.17 (>>1)
glm.nb coefficients (R):
Estimate Std. Error (Intercept) 0.491 0.025 x1 0.748 0.025 x2 -0.432 0.023
glm.nb theta (= r = alpha) = 2.10 [from-scratch r 2.10, PyMC 2.09, truth 2.0]
2. Overdispersion: glm.nb recovers θ = r¶
Pearson dispersion ≈ 2.17 flags overdispersion. glm.nb recovers θ = 2.1 (= our r = PyMC α, truth 2) and coefficients 0.491 / 0.748 / −0.432 with SEs 0.025 / 0.025 / 0.023 — matching the from-scratch NegBin and PyMC. R's glm.nb SEs also confirm the from-scratch SE fix (the numerical-Hessian SEs were the correct ones).
# --- §3 Epil real data: Poisson vs glm.nb (the treatment effect) ---
e <- read.csv("epil.csv"); e$lBase <- log(e$Base / 4); e$lAge <- log(e$Age)
mep <- glm(seizures ~ lBase + Trt + lAge + V4, family = poisson, data = e)
menb <- glm.nb(seizures ~ lBase + Trt + lAge + V4, data = e)
dispE <- sum(residuals(mep, "pearson")^2) / mep$df.residual
cat(sprintf("EPIL Pearson dispersion (Poisson) = %.2f\n\n", dispE))
cat("Poisson GLM:\n"); print(round(summary(mep)$coef[, 1:2], 3))
cat("\nglm.nb:\n"); print(round(summary(menb)$coef[, 1:2], 3))
cat(sprintf("\nglm.nb theta (= r) = %.2f\n", menb$theta))
ci <- confint.default(menb)["Trt", ]
cat(sprintf("\nTREATMENT (progabide): glm.nb Trt = %.3f SE %.3f rate ratio %.2f 95%% CI [%.2f, %.2f]\n",
coef(menb)["Trt"], summary(menb)$coef["Trt",2], exp(coef(menb)["Trt"]), exp(ci[1]), exp(ci[2])))
EPIL Pearson dispersion (Poisson) = 4.71
Poisson GLM:
Estimate Std. Error (Intercept) -2.340 0.404 lBase 1.224 0.033 Trt -0.017 0.048 lAge 0.579 0.110 V4 -0.160 0.055
glm.nb:
Estimate Std. Error (Intercept) -1.184 0.800 lBase 1.060 0.068 Trt -0.233 0.100 lAge 0.369 0.227 V4 -0.152 0.115
glm.nb theta (= r) = 2.63
TREATMENT (progabide): glm.nb Trt = -0.233 SE 0.100 rate ratio 0.79 95% CI [0.65, 0.96]
Data vs fitted, by arm (R glm.nb). Left: observed mean seizures per 2-week visit (points ± 95% CI) vs the NegBin-fitted means (dashed), placebo vs progabide — the fit separates the arms (placebo ≈ 9.4, progabide ≈ 7.3 at visits 1–3) where the raw points overlap. Right: the covariate-adjusted treatment effect (typical patient, mean baseline & age), placebo ≈ 6.8 vs progabide ≈ 5.4, RR 0.79 with 95% CIs — identical to the from-scratch and PyMC figures.
# Data vs fitted, by treatment arm (mirrors the from-scratch notebook)
e$fit <- fitted(menb)
ag <- aggregate(cbind(seizures, fit) ~ visit + Trt, e, mean)
se <- aggregate(seizures ~ visit + Trt, e, function(z) sd(z)/sqrt(length(z)))
par(mfrow = c(1, 2)); cols <- c('steelblue', 'firebrick')
plot(NA, xlim=c(1,4), ylim=c(5,12), xlab='visit', ylab='mean seizures / 2 weeks', main='Observed vs fitted, by arm')
for (t in 0:1) {
o <- ag[ag$Trt==t, ]; o <- o[order(o$visit), ]; s <- se[se$Trt==t, ]; s <- s[order(s$visit), ]
arrows(o$visit, o$seizures-1.96*s$seizures, o$visit, o$seizures+1.96*s$seizures, angle=90, code=3, length=.04, col=cols[t+1])
points(o$visit, o$seizures, col=cols[t+1], pch=19); lines(o$visit, o$fit, col=cols[t+1], lty=2, lwd=2)
}
legend('bottomleft', c('placebo obs','progabide obs','placebo fit','progabide fit'),
col=c(cols,cols), pch=c(19,19,NA,NA), lty=c(NA,NA,2,2), bty='n', cex=.8)
nd0 <- data.frame(lBase=mean(e$lBase), Trt=0, lAge=mean(e$lAge), V4=0); nd1 <- nd0; nd1$Trt <- 1
p0 <- predict(menb, nd0, type='link', se.fit=TRUE); p1 <- predict(menb, nd1, type='link', se.fit=TRUE)
m <- c(exp(p0$fit), exp(p1$fit))
lo <- c(exp(p0$fit-1.96*p0$se.fit), exp(p1$fit-1.96*p1$se.fit)); hi <- c(exp(p0$fit+1.96*p0$se.fit), exp(p1$fit+1.96*p1$se.fit))
bp <- barplot(m, names.arg=c('placebo','progabide'), col=cols, ylim=c(0, max(hi)*1.1),
ylab='predicted seizures / 2 weeks', main=sprintf('Adjusted effect (RR %.2f)', exp(coef(menb)['Trt'])))
arrows(bp, lo, bp, hi, angle=90, code=3, length=.08); par(mfrow=c(1,1))
3. Epil: R confirms the treatment finding¶
| engine | Poisson Trt (SE) | NegBin Trt (SE) | r / θ |
|---|---|---|---|
| from-scratch | −0.017 (0.048) | −0.233 (0.101) | 2.6 |
| PyMC | −0.017 (0.048) | −0.234 (0.105)* | 2.6 |
R (glm / glm.nb) |
−0.017 (0.048) | −0.233 (0.100) | 2.63 |
*PyMC value is the posterior mean/SD. Pearson dispersion ≈ 4.71; glm.nb θ = 2.63.
All three engines agree exactly. Crucially, R's glm.nb gives the NegBin Trt SE as 0.100, confirming the from-scratch SE fix (the earlier BFGS-hess_inv value of 0.055 was wrong; the numerical-Hessian value 0.101 is right). The treatment effect is rate ratio ≈ 0.79 (a ~21% seizure reduction), but — as the Python notebook stresses — this single-level NegBin ignores the 4 correlated visits per patient, so the proper inference needs patient random effects (the hierarchical Poisson model, Breslow–Clayton 1993 / BUGS Epil).
Conclusion¶
Four engines — Newton/IRLS MLE, from-scratch RW-Metropolis, PyMC NUTS, and R glm/MASS::glm.nb — give the same Poisson and Negative Binomial fits on synthetic and real (Epil) data. glm.nb's θ is exactly our dispersion r. The count-regression core is fully validated; the natural next step is the hierarchical Poisson extension.