Mice — Weibull survival across 4 groups (BUGS "Mice")¶

Reuses survival_mcmc.py; the new angle is a strongly increasing hazard¶

Same parametric machinery as survival_python.ipynb (Weibull proportional hazards with right-censoring) — see that notebook for the full censored-likelihood derivation and the MCMC algorithm. Here the application is the BUGS Vol I "Mice" experiment: 80 mice in 4 groups (20 each), survival in weeks, with administrative censoring at the end of follow-up. Two things differ from Leuk: the covariate is a 4-level group factor, and the hazard turns out to be strongly increasing (Weibull shape $k\gg1$ — ageing / wear-out), where Leuk's was near-flat.

Model¶

$$T_{i}\sim\text{Weibull}:\ \ h(t\mid g)=\lambda_g\,k\,t^{k-1},\quad \lambda_g=\exp(\beta_0+\beta_{g}),\quad S(t)=e^{-\lambda_g t^{k}},$$ with group 1 as baseline ($\beta_1\equiv0$); $e^{\beta_g}$ is the hazard ratio of group $g$ vs group 1, and the median survival of group $g$ is $(\log 2/\lambda_g)^{1/k}$. Censored mice (still alive at the cutoff) contribute only $S(t)$.

In [1]:
import numpy as np, pandas as pd, matplotlib.pyplot as plt
from survival_mcmc import weibull_mle, weibull_rwm, km_estimator, summary

d = pd.read_csv('mice.csv'); G = d['group'].to_numpy()
t = d['time'].to_numpy(float); delta = d['cens'].to_numpy(float)
X = np.column_stack([np.ones(len(d))] + [(G==g).astype(float) for g in (2,3,4)])  # group 1 baseline
for g in (1,2,3,4):
    m = G==g; print('  group %d: n=%d  events=%d  censored=%d  median.obs=%.0f' % (g,m.sum(),int(delta[m].sum()),int((1-delta[m]).sum()),np.median(t[m])))

mle, se, _ = weibull_mle(X, t, delta)
o = weibull_rwm(X, t, delta, seed=1); dr = o['draws']; k = np.exp(dr[:,-1])
print('\nWeibull shape k = %.3f  95%% CI [%.2f, %.2f]   (k>>1 = strongly increasing hazard / wear-out)' % (k.mean(), np.percentile(k,2.5), np.percentile(k,97.5)))
print('hazard ratios vs group 1:')
for j,g in [(1,2),(2,3),(3,4)]:
    hr = np.exp(dr[:,j]); print('  group %d: HR %.2f [%.2f, %.2f]  P(HR<1)=%.2f' % (g, hr.mean(), np.percentile(hr,2.5), np.percentile(hr,97.5), (hr<1).mean()))
# posterior median survival per group
beta = dr[:,:4]; meds = {}
for g in (1,2,3,4):
    xg = np.array([1, g==2, g==3, g==4], float)
    lam = np.exp(beta @ xg); meds[g] = (np.log(2)/lam)**(1/k)
print('median survival (weeks):  ' + '  '.join('g%d=%.1f[%.0f,%.0f]'%(g,meds[g].mean(),np.percentile(meds[g],2.5),np.percentile(meds[g],97.5)) for g in (1,2,3,4)))

  group 1: n=20  events=19  censored=1  median.obs=24
  group 2: n=20  events=13  censored=7  median.obs=32
  group 3: n=20  events=16  censored=4  median.obs=24
  group 4: n=20  events=17  censored=3  median.obs=21

Weibull shape k = 3.179  95% CI [2.56, 3.82]   (k>>1 = strongly increasing hazard / wear-out)
hazard ratios vs group 1:
  group 2: HR 0.34 [0.15, 0.66]  P(HR<1)=1.00
  group 3: HR 0.74 [0.36, 1.38]  P(HR<1)=0.85
  group 4: HR 1.57 [0.75, 3.00]  P(HR<1)=0.14

median survival (weeks):  g1=24.3[21,28]  g2=34.9[30,42]  g3=27.2[23,32]  g4=21.5[18,25]
In [2]:
fig, ax = plt.subplots(1, 2, figsize=(12, 4.6))
cols = ['steelblue','seagreen','goldenrod','firebrick']; kk = np.exp(mle[-1]); grid = np.linspace(0,42,250)
for g in (1,2,3,4):
    m = G==g; tk, Sk = km_estimator(t[m], delta[m])
    ax[0].step(tk, Sk, where='post', color=cols[g-1], lw=1.7, label=f'group {g}')
    xg = np.array([1, g==2, g==3, g==4], float); lam = np.exp(mle[:4] @ xg)
    ax[0].plot(grid, np.exp(-lam*grid**kk), '--', color=cols[g-1], lw=1.0)
ax[0].set_xlabel('weeks'); ax[0].set_ylabel('S(t)'); ax[0].set_ylim(0,1.02)
ax[0].set_title('KM (step) vs fitted Weibull (dashed) by group'); ax[0].legend(fontsize=8)
mm = [meds[g] for g in (1,2,3,4)]
ax[1].errorbar(range(1,5), [m.mean() for m in mm],
               yerr=[[m.mean()-np.percentile(m,2.5) for m in mm],[np.percentile(m,97.5)-m.mean() for m in mm]],
               fmt='o', color='black', capsize=4)
for g in (1,2,3,4): ax[1].scatter(g, mm[g-1].mean(), s=60, color=cols[g-1], zorder=3)
ax[1].set_xticks(range(1,5)); ax[1].set_xlabel('group'); ax[1].set_ylabel('median survival (weeks)')
ax[1].set_title('Posterior median survival by group (95% CrI)')
plt.tight_layout(); plt.savefig('mice_surv.png', dpi=120, bbox_inches='tight'); plt.show()
No description has been provided for this image

Results¶

  • Shape $k\approx3.2$ (CI ~2.6–3.8), far above 1 — a steeply increasing hazard: mice die faster the longer the experiment runs (ageing/wear-out). This is the qualitative contrast with Leuk ($k\approx1.35$), and it rules out the exponential decisively.
  • Group differences: group 2 survives longest (median ≈ 35 wk; HR ≈ 0.34 vs group 1 — a real protective effect, P(HR<1) high), group 4 shortest (median ≈ 21 wk; HR ≈ 1.6). Groups 1 and 3 sit in between.
  • The censoring is heaviest in group 2 (7 of 20 still alive at cutoff), which is exactly why its raw median understates the benefit — the censored likelihood credits those survivors, and the fitted curve stays high.

Takeaway¶

Methodologically identical to Leuk (same survival_mcmc.py), this example highlights the Weibull shape as a finding in its own right ($k\gg1$ = wear-out) and shows the model handling a multi-level factor. Cross-checks: PyMC and R (survreg) follow.