Surgical — hierarchical binomial for institutional comparison

Shrinkage and ranking of hospital mortality; reuses the Seeds binomial-GLMM engine

Module: seeds_gibbs.py (same engine as the Seeds project — see it for the algorithm; this is an intercept-only binomial GLMM, so no new code). BUGS Vol I "Surgical": deaths in 12 hospitals performing cardiac surgery on infants. The goal is the classic "hospital league table" problem — which hospitals have unusually high mortality? — done right: shrinking unreliable small-hospital rates and showing the uncertainty in any ranking.

The data — infant cardiac-surgery mortality (BUGS Surgical)

12 hospitals; each reports deaths $r_i$ out of $n_i$ operations.

field	meaning
r	deaths
n	operations (volume: 47 to 810 — very uneven)

Raw mortality ranges 0%–14%, but the volumes differ enormously, so the raw rates are not comparable: hospital 1's 0/47 = 0% is based on tiny numbers, while hospital 4's 5.7% rests on 810 operations. The hierarchical model pools information across hospitals to give each a stabilised rate.

Model

$$r_i\sim\text{Binomial}(n_i,p_i),\qquad \text{logit}\,p_i=\mu+b_i,\qquad b_i\sim N(0,\sigma^2).$$ A pure intercept-only binomial GLMM: $\mu$ is the population log-odds, $b_i$ the hospital deviation. This is exactly seeds_gibbs with the design matrix $X$ = a column of ones (one random effect per hospital = per row). The posterior hospital rate $p_i$ is a partial-pooling estimate — a precision-weighted blend of the hospital's own $r_i/n_i$ and the population rate, shrinking small-volume hospitals most.

import numpy as np, pandas as pd, matplotlib.pyplot as plt
from scipy.special import expit
from seeds_gibbs import seeds_gibbs

d = pd.read_csv('surgical.csv'); r = d['r'].to_numpy(float); n = d['n'].to_numpy(float)
X = np.ones((len(d), 1))                                       # intercept-only -> random-intercept model
o = seeds_gibbs(r, n, X, R=30000, burn=8000, seed=1)
mu = o['beta'][:, 0]; P = expit(mu[:, None] + o['bdraws'])    # per-hospital posterior rate draws
raw = r / n; shrunk = P.mean(0); pop = expit(mu).mean()
print('overall mortality %.1f%%   population rate (model) %.1f%%   sigma %.3f   (acc %.2f/%.2f)'
      % (100*r.sum()/n.sum(), 100*pop, o['sigma'].mean(), o['accept_beta'], o['accept_b']))
print('\n%-5s %4s %6s %8s %16s' % ('hosp','n','raw%','shrunk%','95% CrI'))
for i in range(len(d)):
    print('%-5d %4d %6.1f %8.1f   [%.1f, %.1f]' % (i+1, int(n[i]), 100*raw[i], 100*shrunk[i],
          100*np.percentile(P[:,i],2.5), 100*np.percentile(P[:,i],97.5)))

overall mortality 7.4%   population rate (model) 7.4%   sigma 0.321   (acc 0.45/0.57)

hosp     n   raw%  shrunk%          95% CrI
1       47    0.0      5.8   [2.4, 9.2]
2      148   12.2      9.9   [6.7, 14.9]
3      119    6.7      7.2   [4.3, 10.9]
4      810    5.7      6.1   [4.6, 7.8]
5      211    3.8      5.6   [3.1, 8.2]
6      196    6.6      7.1   [4.6, 9.9]
7      148    6.1      6.9   [4.1, 9.9]
8      215   14.4     11.6   [7.3, 16.5]
9      207    6.8      7.1   [4.6, 10.1]
10      97    8.2      7.9   [4.8, 11.9]
11     256   11.3      9.8   [6.9, 13.8]
12     360    6.7      7.0   [4.9, 9.3]

fig, ax = plt.subplots(1, 2, figsize=(12, 4.8))
# (A) shrinkage: raw rate -> partial-pooled rate, point size proportional to volume n
order = np.argsort(raw); xpos = np.arange(len(d))
for k, i in enumerate(order):
    ax[0].plot([k, k], [100*raw[i], 100*shrunk[i]], color='grey', lw=.8, zorder=1)
ax[0].scatter(xpos, 100*raw[order], s=35, color='darkorange', zorder=2, label='raw $r/n$')
ax[0].scatter(xpos, 100*shrunk[order], s=15+0.12*n[order], color='steelblue', edgecolor='k', lw=.3, zorder=3, label='shrunk (partial pool)')
ax[0].axhline(100*pop, color='firebrick', ls='--', lw=1, label='population rate')
ax[0].set_xticks(xpos); ax[0].set_xticklabels([str(i+1) for i in order], fontsize=7)
ax[0].set_xlabel('hospital (sorted by raw rate)'); ax[0].set_ylabel('mortality %')
ax[0].set_title('Shrinkage toward the population (point size ∝ volume)'); ax[0].legend(fontsize=8)
# (B) ranking with uncertainty: posterior rate + 95% CrI, sorted -> intervals overlap heavily
order2 = np.argsort(shrunk); yy = np.arange(len(d))
ax[1].errorbar(100*shrunk[order2], yy, xerr=[100*(shrunk[order2]-np.percentile(P,2.5,0)[order2]),
               100*(np.percentile(P,97.5,0)[order2]-shrunk[order2])], fmt='o', color='steelblue', capsize=3)
ax[1].axvline(100*pop, color='firebrick', ls='--', lw=1, label='population rate')
ax[1].set_yticks(yy); ax[1].set_yticklabels(['H%d'%(i+1) for i in order2], fontsize=7)
ax[1].set_xlabel('posterior mortality % (95% CrI)'); ax[1].set_title('Ranking is uncertain (intervals overlap)'); ax[1].legend(fontsize=8)
plt.tight_layout(); plt.savefig('surgical.png', dpi=120, bbox_inches='tight'); plt.show()

# companion graphs: how much each hospital shrinks (vs volume), and P(worse than average)
fig, ax = plt.subplots(1, 2, figsize=(12, 4.4))
ax[0].scatter(n, 100*np.abs(raw - shrunk), s=40, color='purple', edgecolor='k', lw=.3)
for i in range(len(d)): ax[0].annotate(str(i+1), (n[i], 100*abs(raw[i]-shrunk[i])), fontsize=7, xytext=(3,2), textcoords='offset points')
ax[0].set_xlabel('hospital volume $n$'); ax[0].set_ylabel('shrinkage |raw − shrunk| (pct pts)')
ax[0].set_title('Small hospitals shrink most'); ax[0].set_xscale('log')
exc = (P > pop).mean(0)                                        # P(hospital rate > population rate)
o2 = np.argsort(exc); ax[1].barh(np.arange(len(d)), exc[o2], color=['firebrick' if e>0.5 else 'steelblue' for e in exc[o2]])
ax[1].axvline(0.5, color='k', ls=':', lw=1); ax[1].set_yticks(np.arange(len(d))); ax[1].set_yticklabels(['H%d'%(i+1) for i in o2], fontsize=7)
ax[1].set_xlabel('P(mortality > population average)'); ax[1].set_title('Probability each hospital is above average')
plt.tight_layout(); plt.savefig('surgical_extra.png', dpi=120, bbox_inches='tight'); plt.show()

Results & interpretation

• Partial pooling stabilises the rates. Hospital 1's raw 0% (only 47 ops) is pulled up to ~5.8% — a far more believable estimate than "zero risk"; the high outliers (H8 14.4% → 11.6%, H2 12.2% → 9.9%) are pulled down. The large-volume hospital 4 (810 ops) barely moves — its raw rate is already reliable. Small hospitals shrink most (left companion panel).
• Ranking is genuinely uncertain. The 95% credible intervals overlap heavily (right panel of the main figure) — the data do not support a confident "league table." Only H8 (and to a lesser extent H2/H11) are probably above the population average (exceedance probabilities), while most hospitals are indistinguishable.
• Population mortality ≈ 7.4%, between-hospital SD σ ≈ 0.32 (on the logit scale) — real but moderate variation.

Takeaways

• The institutional-comparison problem solved by shrinkage: raw league tables over-trust small units and manufacture spurious extremes; the hierarchical binomial gives stabilised rates and honest uncertainty about who is really worse.
• Same engine as Seeds — an intercept-only binomial GLMM (X = ones). It's the binomial analogue of the conjugate Pumps Poisson shrinkage, on the logit scale.
• Cross-checks: PyMC (pm.Binomial) and R (lme4::glmer(... ~ 1 + (1|hospital))) follow.