Frailty survival — Weibull with per-patient random effects (BUGS "Kidney")

The hierarchical extension of the survival project

Module: frailty_surv.py. BUGS Vol I "Kidney" (McGilchrist & Aisbett 1991): times to infection at the site of a dialysis catheter, 2 recurrences per patient for 38 patients, right-censored. Because each patient contributes two correlated times, we add a patient frailty — a random effect on the hazard. This is the random-effects (hierarchical) extension of survival_python.ipynb.

The implementation is deliberately modular: it is the survival likelihood (survival_mcmc.py) × the random-intercept sampler (hpois_gibbs.py) — the two pieces you've already built, wired together.

The hazard, in brief

Survival models are written in terms of the hazard $h(t)$ — the instantaneous event rate at time $t$ among those still event-free: $$h(t)=\lim_{dt\to0}\frac{P(t\le T<t+dt\mid T\ge t)}{dt}=\frac{f(t)}{S(t)}=-\frac{d}{dt}\log S(t).$$ It is a rate, not a probability (it can exceed 1; units = events per unit time), and it fully determines the survival curve via $S(t)=\exp\!\big(-\int_0^t h(u)\,du\big)$. For the Weibull used here, $h(t)=\lambda\,k\,t^{k-1}$ — rising over time if the shape $k>1$, flat if $k=1$ (the exponential/constant-hazard case).

Everything in this model acts multiplicatively on the hazard:

• a covariate scales it by the hazard ratio $e^{\beta_j}$ — a proportional-hazards effect, the same multiple at every $t$ (so $e^{\beta_{\text{female}}}=0.17$ means a female's hazard is ~17% of a male's, always);
• the frailty scales it by $u_j=e^{w_j}$, the patient's personal multiplier.

(The companion survival_python.ipynb has the fuller treatment, including a fitted-hazard and cumulative-hazard figure.)

What is a "frailty"?

A frailty is an unobserved random multiplier on a subject's hazard — the survival analogue of a random intercept in a mixed model: $$h_{ij}(t)=\underbrace{h_0(t)\,e^{x_{ij}'\beta}}_{\text{observed part}}\times\underbrace{u_j}_{\text{frailty}},\qquad u_j=e^{w_j},\quad w_j\sim N(0,\sigma_w^2).$$ $u_j>1$ marks a frail patient — one whose hazard runs above what their covariates predict (they get infected sooner); $u_j<1$ a robust one. It does two jobs at once:

1. Absorbs unmeasured heterogeneity. Patients differ in immune strength, hygiene, catheter care, comorbidities — none of which are in the data. The frailty is a per-patient catch-all for all of it.
2. Models the within-patient correlation. A patient's two recurrence times share the same $u_j$ ("shared frailty"), so they're correlated. Treating the 76 times as independent would be pseudo-replication — it would make the covariate confidence intervals falsely narrow. (Exactly the lesson from the Epil hierarchical-Poisson model, now on the hazard scale.)

The single number $\sigma_w$ summarises the heterogeneity: $\sigma_w=0$ means patients are exchangeable (no frailty needed); larger $\sigma_w$ means wider spread. Our $\sigma_w\approx0.66$ implies the per-patient hazard multiplier $e^{w}$ ranges roughly [0.3, 3.6] across 95% of patients — a ~12-fold gap between the most robust and the most frail, on top of age/sex/disease. (We use a log-normal frailty, $w\sim N(0,\sigma_w^2)$; a Gamma prior on $u$ is the classic alternative.)

Model

$$h(t_{ij}\mid x,w)=\lambda_{ij}\,k\,t^{k-1},\qquad \lambda_{ij}=\exp(x_{ij}'\beta+w_{p(i)}),\qquad w_j\sim N(0,\sigma_w^2),$$ where $w_j$ is patient $j$'s frailty (shared by both of their observations); $e^{w_j}>1$ means a higher-than-average infection hazard. Covariates $x$: age, sex (female), and disease type (GN/AN/PKD vs Other). Events contribute $h(t)S(t)$, censored cases only $S(t)$.

Algorithm (Metropolis-within-Gibbs) — likelihood × random-intercept sampler

1. $(\beta,\log k)$ — block RW-Metropolis on the Weibull censored log-likelihood (frailty enters as an offset).
2. $w_j$ — per-patient RW-Metropolis, vectorised. The frailty conditional is exactly the hpois random-intercept form: with $A_j=\sum_{i\in j}\delta_{ij}$ ("event count") and $S_j=\sum_{i\in j}e^{x_{ij}'\beta}t_{ij}^{k}$ ("rate"), $\ \log p(w_j)= A_j w_j - S_j e^{w_j} - w_j^2/(2\sigma_w^2)$ — curvature-scaled step.
3. centre $w\to\beta_0$; $\sigma_w^2$ Inverse-Gamma.

The data — BUGS "Kidney" (McGilchrist & Aisbett 1991)

38 patients on portable (home) kidney dialysis, monitored for infection at the catheter insertion site. When an infection occurs the catheter is removed, the infection is cleared, and a new catheter is inserted — producing a second time-to-infection. So each patient contributes two recurrence times (76 rows). A time is censored (status = 0) when the catheter is removed for a reason other than infection (e.g. it stopped functioning) — 18 of the 76.

field	meaning
id	patient (1–38), 2 rows each
time	days to infection (or censoring); range 2–562, median 40
status	1 = infection (event, 58), 0 = censored (18)
age	years (10–69, mean 44); may differ slightly between a patient's two rows
sex	1 = male, 2 = female → we use female. 10 male / 28 female patients
disease	kidney disease: GN (glomerulonephritis), AN (acute nephritis), PKD (polycystic), or Other (baseline)

The sex effect is already stark in the raw data: male records get infected in 90% of cases (median 16 days) vs 71% for females (median 62 days). The model's task is to confirm that gap survives adjustment for disease and age — and for the two-observations-per-patient correlation, which is what the frailty handles.

import numpy as np, pandas as pd, matplotlib.pyplot as plt
from frailty_surv import frailty_weibull_gibbs, summary

d = pd.read_csv('kidney.csv')
age = d['age'].to_numpy(float) - d['age'].mean(); female = (d['sex']==2).astype(float).to_numpy()
GN = (d['disease']=='GN').astype(float); AN = (d['disease']=='AN').astype(float); PKD = (d['disease']=='PKD').astype(float)
X = np.column_stack([np.ones(len(d)), age, female, GN, AN, PKD])
t = d['time'].to_numpy(float); delta = d['status'].to_numpy(float); pat = (d['id']-1).to_numpy(int)
print('n=%d obs, %d patients (2 each), %d events, %d censored' % (len(d), pat.max()+1, int(delta.sum()), int((1-delta).sum())))

o = frailty_weibull_gibbs(t, X, delta, pat, R=25000, burn=5000, seed=1)
print('acceptance: beta %.2f  w %.2f' % (o['accept_beta'], o['accept_w']))
nm = ['intercept','age','female','GN','AN','PKD']
print('\ncovariate effects (hazard ratios):')
for r in summary(o['beta'], nm):
    if r[0] != 'intercept':
        print('  %-8s coef %+.3f [%+.2f, %+.2f]   HR %.2f' % (r[0], r[1], r[3], r[4], np.exp(r[1])))
print('\nshape k     = %.3f [%.2f, %.2f]' % (o['k'].mean(), np.percentile(o['k'],2.5), np.percentile(o['k'],97.5)))
print('frailty sd  = %.3f [%.2f, %.2f]  (>0 => real patient heterogeneity)' % (o['sigma_w'].mean(), np.percentile(o['sigma_w'],2.5), np.percentile(o['sigma_w'],97.5)))

n=76 obs, 38 patients (2 each), 58 events, 18 censored

acceptance: beta 0.36  w 0.63

covariate effects (hazard ratios):
  age      coef +0.002 [-0.03, +0.03]   HR 1.00
  female   coef -1.791 [-2.72, -0.88]   HR 0.17
  GN       coef -0.003 [-1.12, +1.07]   HR 1.00
  AN       coef +0.535 [-0.44, +1.61]   HR 1.71
  PKD      coef -0.889 [-2.40, +0.74]   HR 0.41

shape k     = 1.108 [0.90, 1.33]
frailty sd  = 0.655 [0.39, 1.05]  (>0 => real patient heterogeneity)

fig, ax = plt.subplots(1, 2, figsize=(12, 4.8))
# (1) per-patient frailty (caterpillar): which patients run hot/cold
W = o['W']; wm = W.mean(0); lo = np.percentile(W,2.5,axis=0); hi = np.percentile(W,97.5,axis=0)
order = np.argsort(wm); J = len(wm)
ax[0].errorbar(wm[order], range(J), xerr=[wm[order]-lo[order], hi[order]-wm[order]], fmt='o', ms=3, lw=.6, color='steelblue')
ax[0].axvline(0, color='firebrick', ls='--', lw=1)
ax[0].set_xlabel('patient frailty  $w_j$  (log-hazard)'); ax[0].set_ylabel('patient (sorted)')
ax[0].set_title(f'Per-patient frailty  (sd $\\approx$ {o["sigma_w"].mean():.2f})')
# (2) covariate hazard-ratio forest
labs = ['age','female','GN','AN','PKD']; idx = [1,2,3,4,5]
hrm = [np.exp(o['beta'][:,j]).mean() for j in idx]
hl = [np.exp(np.percentile(o['beta'][:,j],2.5)) for j in idx]; hh = [np.exp(np.percentile(o['beta'][:,j],97.5)) for j in idx]
yy = range(len(labs))
ax[1].errorbar(hrm, yy, xerr=[np.array(hrm)-np.array(hl), np.array(hh)-np.array(hrm)], fmt='o', color='black', capsize=3)
ax[1].axvline(1, color='firebrick', ls='--', lw=1, label='HR = 1')
ax[1].set_yticks(list(yy)); ax[1].set_yticklabels(labs); ax[1].set_xscale('log')
ax[1].set_xlabel('hazard ratio (log scale)'); ax[1].set_title('Covariate effects'); ax[1].legend(fontsize=8)
plt.tight_layout(); plt.savefig('frailty_kidney.png', dpi=120, bbox_inches='tight'); plt.show()

Results & interpretation

• Sex dominates: women have a far lower infection hazard — HR ≈ 0.15–0.17 (coef ≈ −1.8), CI well below 1. PKD patients are also lower-risk (HR ≈ 0.3–0.4); age has no effect.
• Shape $k\approx1.1$ — hazard is close to constant (near-exponential), i.e. the infection risk doesn't change much over follow-up.
• The frailty is real: $\sigma_w\approx0.65$ (CI excludes small values) — substantial patient-to-patient variation beyond the covariates. The caterpillar plot shows individual patients running hot or cold; ignoring this clustering (and the two-per-patient correlation) would understate the uncertainty in the covariate effects.

Takeaways

• Frailty = survival likelihood × random-intercept sampler. The per-patient update is identical in form to the hpois Poisson random intercept — only the likelihood term (Weibull-censored vs Poisson) changed. The project's modular design pays off directly here.
• Bayes keeps the frailty that ML can drop. A Cox-frailty maximum-likelihood fit (coxph(... + frailty(id))) drives $\sigma_w\to0$ (boundary), whereas the Bayesian posterior — and the parametric parfm ML — keep it near 0.6. With only 38 clusters of 2, the variance is weakly identified, so this prior/method sensitivity is the honest headline. Cross-checks (PyMC, R parfm/coxme) follow.

Additional graphs

Three more views that put the fit on the survival scale, stress-test the frailty, and check the model:

1. Survival views — Kaplan–Meier by sex vs the fitted Weibull, and the survival curves for a typical / frail / robust patient (what the frailty spread means in days).
2. Is the frailty real? — the $\sigma_w$ prior vs posterior (the posterior pulls clearly off 0, where the coxph ML estimate collapses), and the covariate hazard ratios with vs without the frailty (modelling the within-patient correlation widens the intervals).
3. Diagnostics — the within-patient time-1/time-2 correlation that motivates the frailty, the shrinkage of unpooled patient estimates toward the population, and a posterior-predictive check (observed KM inside the model's predictive band).

# --- Survival views: (1) KM by sex + fitted, (2) typical vs frail vs robust patient ---
from survival_mcmc import km_estimator
bm = o['beta'].mean(0); kbar = o['k'].mean(); grid = np.linspace(1, 562, 250)
fig, ax = plt.subplots(1, 2, figsize=(12, 4.6))

# (1) KM by sex (empirical) vs population-averaged fitted Weibull (w=0)
for g, lab, col in [(0, 'male', 'steelblue'), (1, 'female', 'firebrick')]:
    msk = female == g
    tk, Sk = km_estimator(t[msk], delta[msk])
    ax[0].step(tk, Sk, where='post', color=col, lw=1.8, label='KM ' + lab)
    Sfit = np.array([np.mean(np.exp(-np.exp(X[msk] @ bm) * tt ** kbar)) for tt in grid])
    ax[0].plot(grid, Sfit, '--', color=col, lw=1.1)
ax[0].set_xlabel('days'); ax[0].set_ylabel('S(t)'); ax[0].set_ylim(0, 1.02)
ax[0].set_title('KM by sex (step) vs fitted Weibull (dashed)'); ax[0].legend(fontsize=8)

# (2) frailty spread: a typical female patient (Other disease, mean age) at u = e^{-sd}, 1, e^{+sd}
xref = np.array([1.0, 0.0, 1.0, 0.0, 0.0, 0.0])             # mean age, female, Other disease
lam0 = np.exp(xref @ bm); sw = o['sigma_w'].mean()
for mult, lab, col in [(np.exp(-sw), 'robust  $u=e^{-\\sigma}$', 'seagreen'),
                       (1.0, 'typical  $u=1$', 'black'),
                       (np.exp(sw), 'frail  $u=e^{+\\sigma}$', 'firebrick')]:
    ax[1].plot(grid, np.exp(-lam0 * mult * grid ** kbar), color=col, lw=1.7, label=lab)
ax[1].set_xlabel('days'); ax[1].set_ylabel('S(t)'); ax[1].set_ylim(0, 1.02)
ax[1].set_title(f'What the frailty means (sd $\\approx$ {sw:.2f})'); ax[1].legend(fontsize=8)
plt.tight_layout(); plt.savefig('frailty_survival.png', dpi=120, bbox_inches='tight'); plt.show()

# --- Is the frailty real? (3) sigma_w posterior vs prior, (4) frailty vs no-frailty HRs ---
from survival_mcmc import weibull_rwm
fig, ax = plt.subplots(1, 2, figsize=(12, 4.6))

# (3) sigma_w: prior (IG(2,1) on variance) vs posterior; coxph ML collapses to ~0
rng2 = np.random.default_rng(0)
sw_prior = np.sqrt(1.0 / rng2.gamma(2.0, 1.0, size=40000))    # IG(a0=2,b0=1) on sigma_w^2
bins = np.linspace(0, 2.5, 60)
ax[0].hist(sw_prior, bins=bins, density=True, alpha=.45, color='grey', label='prior')
ax[0].hist(o['sigma_w'], bins=bins, density=True, alpha=.75, color='steelblue', label='posterior')
ax[0].axvline(o['sigma_w'].mean(), color='k', ls='--', lw=1, label=f"post. mean {o['sigma_w'].mean():.2f}")
ax[0].axvline(0.0, color='firebrick', ls=':', lw=2, label='coxph ML $\\to$ 0')
ax[0].set_xlabel('frailty sd  $\\sigma_w$'); ax[0].set_yticks([]); ax[0].set_xlim(0, 2.5)
ax[0].set_title('Is the frailty real?  $\\sigma_w$ prior vs posterior'); ax[0].legend(fontsize=8)

# (4) covariate HRs with vs without the frailty
nf = weibull_rwm(X, t, delta, R=20000, burn=5000, seed=1)['draws']   # no-frailty Weibull; cols: beta..,logk
labs = ['female', 'GN', 'AN', 'PKD']; idx = [2, 3, 4, 5]; yy = np.arange(len(labs))
for dr_, off, col, lab in [(nf, -0.13, 'darkorange', 'no frailty'), (o['beta'], 0.13, 'steelblue', 'with frailty')]:
    hrm = np.array([np.exp(dr_[:, j]).mean() for j in idx])
    hl = np.array([np.exp(np.percentile(dr_[:, j], 2.5)) for j in idx])
    hh = np.array([np.exp(np.percentile(dr_[:, j], 97.5)) for j in idx])
    ax[1].errorbar(hrm, yy + off, xerr=[hrm - hl, hh - hrm], fmt='o', color=col, capsize=3, label=lab)
ax[1].axvline(1, color='firebrick', ls='--', lw=1)
ax[1].set_yticks(yy); ax[1].set_yticklabels(labs); ax[1].set_xscale('log')
ax[1].set_xlabel('hazard ratio (log scale)'); ax[1].set_title('Effect of modelling the frailty'); ax[1].legend(fontsize=8)
plt.tight_layout(); plt.savefig('frailty_check.png', dpi=120, bbox_inches='tight'); plt.show()

# --- Diagnostics: (5) within-patient correlation, (6) shrinkage, (7) posterior-predictive check ---
J = int(pat.max() + 1); bm = o['beta'].mean(0); kbar = o['k'].mean()
fig, ax = plt.subplots(1, 3, figsize=(15, 4.4))

# (5) within-patient correlation: each patient's two times, coloured by frailty
t1 = np.array([t[pat == j][0] for j in range(J)]); t2 = np.array([t[pat == j][1] for j in range(J)])
sc = ax[0].scatter(np.log(t1), np.log(t2), c=o['w'], cmap='coolwarm', s=45, edgecolor='k', lw=.3)
lims = [np.log(t.min()), np.log(t.max())]; ax[0].plot(lims, lims, 'k:', lw=1)
plt.colorbar(sc, ax=ax[0], label='frailty $w_j$')
ax[0].set_xlabel('log(time 1)'); ax[0].set_ylabel('log(time 2)')
ax[0].set_title('Within-patient correlation\n(colour = estimated frailty)')

# (6) shrinkage: unpooled per-patient frailty (MLE) vs pooled (hierarchical) estimate
base = np.exp(X @ bm) * t ** kbar
A = np.bincount(pat, weights=delta, minlength=J); S = np.bincount(pat, weights=base, minlength=J)
w_unp = np.where(A > 0, np.log(np.maximum(A, 1e-9) / S), np.nan)
xmin = np.nanmin(w_unp) - 0.6
ax[1].scatter(w_unp, o['w'], s=35, color='steelblue', edgecolor='k', lw=.3)
nz = ~(A > 0)
if nz.any():
    ax[1].scatter(np.full(nz.sum(), xmin), o['w'][nz], marker='<', color='grey', label='0 events (unpooled $-\\infty$)')
ll = [xmin, np.nanmax(w_unp) + 0.3]; ax[1].plot(ll, ll, 'k:', lw=1, label='no shrinkage')
ax[1].axhline(0, color='firebrick', ls='--', lw=.8)
ax[1].set_xlabel('unpooled frailty (per-patient MLE)'); ax[1].set_ylabel('pooled $\\hat w_j$')
ax[1].set_title('Shrinkage toward the population'); ax[1].legend(fontsize=7)

# (7) posterior-predictive check: marginal survival band (fresh frailties) vs observed KM
from survival_mcmc import km_estimator
grid = np.linspace(1, 562, 250); tk, Sk = km_estimator(t, delta)
B, K, SG = o['beta'], o['k'], o['sigma_w']; rng = np.random.default_rng(1)
sel = np.arange(0, len(K), max(1, len(K) // 250)); SS = np.zeros((len(sel), len(grid)))
for r, ii in enumerate(sel):
    wnew = rng.normal(0, SG[ii], size=J)                     # fresh patient frailties
    lam = np.exp(X @ B[ii] + wnew[pat])
    SS[r] = np.array([np.mean(np.exp(-lam * tt ** K[ii])) for tt in grid])
ax[2].fill_between(grid, np.percentile(SS, 2.5, 0), np.percentile(SS, 97.5, 0), color='seagreen', alpha=.3, label='posterior predictive 95%')
ax[2].step(tk, Sk, where='post', color='black', lw=1.6, label='observed KM')
ax[2].set_xlabel('days'); ax[2].set_ylabel('S(t)'); ax[2].set_ylim(0, 1.02)
ax[2].set_title('Posterior-predictive check'); ax[2].legend(fontsize=8)
plt.tight_layout(); plt.savefig('frailty_diagnostics.png', dpi=120, bbox_inches='tight'); plt.show()