Bayesian Hierarchical Survival — Frailty Models

Model

The hierarchical (random-effects) extension of the Weibull survival model: each subject carries a frailty — an unobserved random multiplier on their hazard, the survival analogue of a random intercept in a mixed model. The implementation is deliberately modular — it is literally the survival likelihood (from the Weibull example) × the random-intercept sampler (from the hierarchical Poisson example), wired together.

h(t_{ij}\mid x, w) = h_0(t)\, e^{x_{ij}'\beta}\, u_j, \qquad \lambda_{ij} = \exp(x_{ij}'\beta + w_{p(i)}), \qquad w_j \sim N(0, \sigma_w^2)

$u_j = e^{w_j}$ — patient $j$ 's frailty (shared by their repeated observations)
$e^{\beta_j}$ — hazard ratio for covariate $j$ ; $k$ — Weibull shape
$\sigma_w$ — frailty spread; $\sigma_w=0$ ⇒ subjects exchangeable (no frailty)

What is a frailty?

A frailty $u_j = e^{w_j}$ ( $w_j \sim N(0,\sigma_w^2)$ ) does two jobs at once. It absorbs unmeasured heterogeneity — patients differ in immune strength, hygiene, comorbidities, none in the data — a per-subject catch-all. And it models within-subject correlation: a subject's repeated event times share the same $u_j$ ("shared frailty"), so treating them as independent would be pseudo-replication that falsely narrows the covariate intervals (exactly the lesson from the Epil hierarchical-Poisson model, now on the hazard scale). $u_j > 1$ marks a frail subject (hazard above what covariates predict), $u_j < 1$ a robust one; $\sigma_w$ summarises the spread.

Sampler — survival likelihood × random-intercept sampler

The key modularity: the per-patient frailty update is identical in form to the hierarchical-Poisson random intercept — only the likelihood term changed (Weibull-censored instead of Poisson). With "event count" $A_j = \sum_{i\in j}\delta_{ij}$ and "rate" $S_j = \sum_{i\in j} e^{x_{ij}'\beta} t_{ij}^k$ , the conditional $\log p(w_j) = A_j w_j - S_j e^{w_j} - w_j^2/(2\sigma_w^2)$ is exactly the hpois form. The block $(\beta, \log k)$ is RW-Metropolis on the censored Weibull likelihood (frailty as an offset); a centring move de-confounds $w$ from $\beta_0$ ; $\sigma_w^2$ is a conjugate Inverse-Gamma draw.

Block	Draw	Update
(1)	$(\beta, \log k) \mid \cdot$	Block RW-Metropolis — censored Weibull likelihood, frailty as offset
(2)	$w_j \mid \cdot$	Vectorised RW-Metropolis — the hpois random-intercept conditional
(3)	$\sigma_w^2 \mid \cdot$	Inverse-Gamma conjugate from the frailties

Re-centre $w$ into $\beta_0$ each sweep (location move; the frailty level is confounded with the intercept).

Notebooks

Applied to the BUGS Kidney data (McGilchrist & Aisbett 1991): times to infection at a dialysis catheter site, 2 recurrences per patient for 38 patients (76 rows, 18 censored). The sex effect dominates — women have a far lower infection hazard (HR $\approx 0.15\text{–}0.2$ across the four engines), and it is the only covariate the data settle: its interval, $[0.07,\ 0.41]$ on the hazard-ratio scale, is nowhere near 1. Nothing else in the model is resolved. Age is flat (HR $1.00$ , interval $[0.97,\ 1.03]$ ) — genuinely no effect rather than an unresolved one. The disease-type indicators all have wide intervals straddling 1 in both directions: PKD's point estimate is low at HR $0.41$ and AN's is high at $1.71$ , but their intervals run $[0.09,\ 2.10]$ and $[0.64,\ 5.00]$ respectively, so with 38 patients neither is a finding. The hazard is near-constant ( $k\approx 1.1$ ). The frailty is real: $\sigma_w\approx 0.65$ implies the per-patient hazard multiplier spans roughly $[0.3, 3.6]$ across 95% of patients — a ~12-fold robust-to-frail gap on top of the covariates. The honest headline is a method sensitivity: a Cox-frailty maximum-likelihood fit drives $\sigma_w\to 0$ (the variance component at the boundary artefact: $\sigma_w$ cannot be negative, so when the data are thin the likelihood is often still rising as it reaches zero and maximum likelihood simply stops there, reporting no heterogeneity rather than a small amount of it — with only 38 clusters of 2 this is exactly the thin case, and a Bayesian posterior, which averages rather than maximises, does not collapse the same way), whereas the Bayesian posterior — and parametric ML via R's parfm — keep it near 0.6. Cross-checked across from-scratch Gibbs, PyMC (non-centered frailty, NUTS), R parfm (parametric Weibull log-normal frailty), and R coxme (semiparametric Cox frailty).

Python Weibull frailty from scratch — Metropolis-within-Gibbs, survival-likelihood × random-intercept, Kidney data (Python) Open → Python PyMC cross-check — non-centered per-patient frailty, censored Weibull via pm.Potential (Python) Open → R R cross-check — parfm (parametric Weibull frailty) + coxme (semiparametric Cox frailty) (R) Open →

Downloads

frailty_surv.py From-scratch Weibull frailty — Metropolis-within-Gibbs (builds on survival_mcmc.py) kidney.csv BUGS Kidney data (McGilchrist & Aisbett 1991) — 38 patients, 2 recurrences each

References

McGilchrist, C. A. & Aisbett, C. W. (1991). Regression with frailty in survival analysis. Biometrics 47(2), 461–466. — the Kidney catheter-infection data and the log-normal frailty model
Vaupel, J. W., Manton, K. G. & Stallard, E. (1979). The impact of heterogeneity in individual frailty on the dynamics of mortality. Demography 16(3), 439–454. — the paper that introduced "frailty" as unobserved heterogeneity
Clayton, D. G. (1978). A model for association in bivariate life tables. Biometrika 65(1), 141–151. — the shared (gamma) frailty for within-cluster association
Hougaard, P. (2000). Analysis of Multivariate Survival Data. Springer. — the standard monograph on frailty and multivariate survival
Therneau, T. M., Grambsch, P. M. & Pankratz, V. S. (2003). Penalized survival models and frailty. Journal of Computational and Graphical Statistics 12(1), 156–175. — the coxme semiparametric Cox-frailty cross-check
Munda, M., Rotolo, F. & Legrand, C. (2012). parfm: parametric frailty models in R. Journal of Statistical Software 51(11), 1–20. — the parfm parametric Weibull-frailty cross-check

Frailty Sampler — Source Code

"""
Frailty survival: Weibull proportional hazards with a per-subject random effect (FRAILTY) and
right-censoring -- from-scratch Metropolis-within-Gibbs. BUGS Vol I "Kidney" (McGilchrist & Aisbett
1991): recurrence times of catheter-insertion infection, 2 per patient, with patient frailties.

This is the hierarchical (random-effects) extension of survival_mcmc.py -- and it is literally the
SURVIVAL LIKELIHOOD (from survival_mcmc.py) x the RANDOM-INTERCEPT SAMPLER (from hpois_gibbs.py).

Model
-----
  h(t_ij | x, w) = lambda_ij * k * t^{k-1},  lambda_ij = exp(x_ij' beta + w_{p(i)}),  w_j ~ N(0, sigma_w^2)
  events contribute h(t)S(t), censored only S(t); w_j is the frailty of patient j (shared by their obs).

Sampler (Metropolis-within-Gibbs)
---------------------------------
  1. (beta, log k)  -- block RW-Metropolis (Weibull censored loglik with the frailty as an offset);
                       proposal = chol(inv(Hessian)) from the no-frailty MLE.
  2. w_j            -- per-patient RW-Metropolis, VECTORISED. The frailty conditional is exactly the
                       hpois random-intercept form: "event count" A_j = sum_{i in j} delta_ij, and
                       "rate" S_j = sum_{i in j} exp(x_ij'beta) t_ij^k, so
                       loglik_j(w) = A_j w - S_j e^{w} - w^2/(2 sigma_w^2);
                       curvature-scaled step s_w / sqrt(S_j e^{w} + 1/sigma_w^2).
  +  centre w -> beta_0 (location move).
  3. sigma_w^2      -- Inverse-Gamma conjugate from {w_j}.
"""
import numpy as np
from survival_mcmc import weibull_mle, _negloglik, _num_hess


def frailty_weibull_gibbs(t, X, delta, patient, R=15000, burn=5000, seed=0,
                          prior_sd=10.0, logk_sd=2.0, a0=2.0, b0=1.0, s_w=1.3, s_beta=0.6):
    rng = np.random.default_rng(seed)
    t = np.asarray(t, float); delta = np.asarray(delta, float); patient = np.asarray(patient, int)
    n, p = X.shape; J = patient.max() + 1; logt = np.log(t)
    A = np.bincount(patient, weights=delta, minlength=J)        # events per patient (fixed)

    mle, _, _ = weibull_mle(X, t, delta)                        # no-frailty init + proposal
    H = _num_hess(lambda q: _negloglik(q, X, t, delta, False), mle)
    Lprop = np.linalg.cholesky(np.linalg.inv(H)) * s_beta
    beta = mle[:p].copy(); logk = mle[p]; w = np.zeros(J); sigma = 0.5

    def wll(b, lk, off):                                        # Weibull censored log-likelihood
        k = np.exp(lk); xb = X @ b + off
        return np.sum(delta * (xb + lk + (k - 1) * logt)) - np.sum(np.exp(xb) * t ** k)

    keep = R - burn
    B = np.zeros((keep, p)); K = np.zeros(keep); SIG = np.zeros(keep); W = np.zeros((keep, J))
    accb = 0; accw = 0.0
    for it in range(R):
        off = w[patient]
        # 1. (beta, log k) block RW-Metropolis
        cur = wll(beta, logk, off) - 0.5*np.sum(beta**2)/prior_sd**2 - 0.5*logk**2/logk_sd**2
        st = Lprop @ rng.standard_normal(p + 1); bp = beta + st[:p]; lkp = logk + st[p]
        pr = wll(bp, lkp, off) - 0.5*np.sum(bp**2)/prior_sd**2 - 0.5*lkp**2/logk_sd**2
        if np.log(rng.random()) < pr - cur:
            beta, logk = bp, lkp; accb += 1
        # 2. w_j vectorised RW-Metropolis (hpois random-intercept form)
        k = np.exp(logk); base = np.exp(X @ beta) * t ** k
        Sj = np.bincount(patient, weights=base, minlength=J)
        curv = Sj * np.exp(w) + 1.0 / sigma ** 2
        wp = w + s_w * rng.standard_normal(J) / np.sqrt(curv)
        d = A * (wp - w) - Sj * (np.exp(wp) - np.exp(w)) - 0.5 * (wp**2 - w**2) / sigma**2
        ok = np.log(rng.random(J)) < d; w = np.where(ok, wp, w); accw += ok.mean()
        m = w.mean(); w -= m; beta[0] += m                     # centre frailties -> intercept
        # 3. sigma_w^2 Inverse-Gamma
        sigma = np.sqrt(1.0 / rng.gamma(a0 + J/2.0, 1.0 / (b0 + 0.5*np.sum(w**2))))
        if it >= burn:
            i = it - burn; B[i] = beta; K[i] = np.exp(logk); SIG[i] = sigma; W[i] = w
    return dict(beta=B, k=K, sigma_w=SIG, w=W.mean(0), W=W,
                accept_beta=accb/R, accept_w=accw/R)


def summary(draws, names):
    return [(nm, draws[:, j].mean(), draws[:, j].std(), *np.percentile(draws[:, j], [2.5, 97.5]))
            for j, nm in enumerate(names)]