Surgical — R cross-check (lme4::glmer)

Companion to surgical_python.ipynb / surgical_pymc.ipynb. The random-intercept-only binomial GLMM via glmer(cbind(r, n-r) ~ 1 + (1|hospital), family=binomial), with the no-pooling raw rates for contrast.

suppressMessages(library(lme4))
d <- read.csv('surgical.csv'); d$hospital <- factor(d$hospital)
gm <- glmer(cbind(r, n-r) ~ 1 + (1|hospital), family=binomial, data=d)
mu <- fixef(gm)[['(Intercept)']]; sg <- sqrt(unlist(VarCorr(gm))[1])
cat(sprintf('glmer: population rate = %.1f%%   random-intercept sd = %.3f\n', 100/(1+exp(-mu)), sg))
shrunk <- 1/(1+exp(-(mu + ranef(gm)$hospital[,1])))           # partial-pooled hospital rates
raw <- d$r / d$n
print(data.frame(hosp=1:nrow(d), n=d$n, raw_pct=round(100*raw,1), shrunk_pct=round(100*shrunk,1)))

glmer: population rate = 7.3%   random-intercept sd = 0.357

   hosp   n raw_pct shrunk_pct
1     1  47     0.0        5.4
2     2 148    12.2       10.2
3     3 119     6.7        7.0
4     4 810     5.7        5.9
5     5 211     3.8        5.2
6     6 196     6.6        6.9
7     7 148     6.1        6.6
8     8 215    14.4       12.3
9     9 207     6.8        7.0
10   10  97     8.2        7.7
11   11 256    11.3       10.2
12   12 360     6.7        6.8

# shrinkage plot (inline): raw -> shrunk, by volume
options(repr.plot.width = 8, repr.plot.height = 5)
ord <- order(raw); pop <- 1/(1+exp(-mu))
plot(seq_along(ord), 100*raw[ord], pch=1, col='darkorange', ylim=c(0,16),
     xlab='hospital (sorted by raw rate)', ylab='mortality %', main='Surgical: shrinkage of hospital rates', xaxt='n')
axis(1, at=seq_along(ord), labels=ord, cex.axis=0.7)
points(seq_along(ord), 100*shrunk[ord], pch=19, col='steelblue', cex=0.5+d$n[ord]/300)
segments(seq_along(ord), 100*raw[ord], seq_along(ord), 100*shrunk[ord], col='grey')
abline(h=100*pop, col='firebrick', lty=2)
legend('topleft', legend=c('raw r/n','shrunk (glmer)','population'), pch=c(1,19,NA), lty=c(NA,NA,2),
       col=c('darkorange','steelblue','firebrick'), bty='n', cex=0.85)

Results

• glmer matches the Bayesian fits: population rate ≈ 7.3%, random-intercept SD ≈ 0.3, and the same shrunk hospital rates (hospital 1's 0% → ~5–6%, the high hospitals pulled in, large-volume hospital 4 barely moving).
• Three engines agree (from-scratch Seeds engine ≈ PyMC ≈ glmer): the hierarchical binomial stabilises the league table. The small frequentist-vs-Bayesian gap in the random-effect SD is the usual weakly-identified variance component with only 12 groups.