Bayesian Weibull Proportional-Hazards Model

Model

Time-to-event data bring a genuinely new ingredient over every earlier example — censoring. A subject still event-free at the end of follow-up contributes only the survival probability $S(t)$, not a density. With event indicator $\delta_i$ (1 = event observed, 0 = censored), the likelihood splits in two: events contribute the hazard × survival, censored cases only the survival. The Weibull proportional-hazards form gives $h(t\mid x) = \lambda_i k\, t^{k-1}$ with $\lambda_i = e^{x_i'\beta}$ — so $e^{\beta_j}$ is the hazard ratio for covariate $j$, and the shape $k$ controls whether hazard rises ($k>1$), falls ($k<1$), or is constant ($k=1$, the exponential model).

• $\delta_i$ — event indicator (1 = event observed, 0 = right-censored)
• $e^{\beta_j}$ — hazard ratio for covariate $j$
• $k$ — Weibull shape: $k>1$ rising, $k<1$ falling, $k=1$ exponential (constant hazard)

Survival is a count model underneath

A nice structural fact: the Weibull log-likelihood is a Poisson GLM kernel in disguise ($\delta$ as 0/1 "counts", offset $\log t^k$) — the piecewise-exponential = Poisson trick behind the BUGS Leuk example, linking survival analysis directly back to the count-model arc. Everything else is familiar GLM machinery: MLE by direct maximisation, a Bayesian RW-Metropolis on $(\beta, \log k)$ with a full-covariance proposal from the MLE Hessian, and a Kaplan–Meier curve for nonparametric comparison — the model-free survival estimate, a staircase that drops by the right fraction at each observed event and simply steps over each censored subject, removing them from the pool without recording a death. Its companion, the Nelson–Aalen estimator, does the same job on the cumulative hazard scale.

Notebooks

Three applications share the same survival_mcmc.py machinery, each surfacing a different facet of the Weibull. Leuk — treatment effect (BUGS Leuk). 42 leukemia patients randomised to 6-MP (mercaptopurine) vs. placebo, time to relapse in weeks. The placebo arm relapsed fast and fully (all 21 events, median 8 weeks); the 6-MP arm did far better, with 12 of 21 still in remission at study end (censored). The control relapse hazard is several times higher than 6-MP's, with $P(\text{control HR}>1)=1.00$ and a mildly increasing baseline ($k\approx 1.35$). The exact multiple depends on a convention worth stating once, because it recurs throughout this section. The fitted coefficient is $\beta_{\text{control}} = 1.73$, so exponentiating the estimate gives a hazard ratio of $e^{1.73} = 5.65$ — equivalently a 6-MP HR of $0.18$ — and that is the figure R's survreg and coxph report and the one the Cox-via-Poisson page compares against. Averaging the hazard ratio itself over the posterior instead gives $E[e^\beta] = 6.31$ $[2.51,\ 14.05]$, a 6-MP HR of $0.158$, which is what the notebook's summary line prints. The second is larger because $e^\beta$ is a convex function of a parameter that is uncertain, so the posterior of the ratio is right-skewed and its mean exceeds the ratio at the mean. Neither is wrong; they answer slightly different questions, and the gap widens as the coefficient gets less certain. Ignoring the censored patients would understate the benefit — the censored likelihood credits them with surviving at least that long. Mice — a multi-level factor & wear-out shape (BUGS Mice). 80 mice in 4 groups, survival in weeks. The covariate is a 4-level group factor, and the shape comes out strongly increasing ($k\approx 3.2$) — an ageing/wear-out hazard that decisively rules out the exponential, the qualitative contrast with Leuk's near-flat $k$. Group 2 survives longest (median $\approx 35$ wk, HR $\approx 0.3$), group 4 shortest. Reliability — engineering life metrics (Meeker & Escobar shock absorber). A single sample (no covariates) of 38 vehicles' distance-to-failure, with 23 of 38 suspended (61% censored). The shape $k\approx 3.3$ signals a clear wear-out failure mode, and the engineering quantities read straight off the posterior — characteristic life $\eta\approx 25{,}550$ km (the distance by which 63% have failed — the Weibull's natural scale), $B_{10}\approx 12{,}700$ km (the distance by which the first 10% have failed, the standard warranty quantity), mean time to failure $\approx 22{,}900$ km, and a reliability at 10,000 km of $0.948$ $[0.880,\ 0.989]$ — with a Weibull probability plot confirming adequacy. Worth noting against the 61% censoring: all three life metrics fall inside the observed range, below the 28,100 km the longest-surviving unit reached, so what the censoring costs here is precision rather than reach — B10's interval runs $9{,}200\text{–}16{,}100$ km because only 15 units actually failed. It is the curve past 28,100 km that is carried entirely by the fitted shape, which is why $k$ rather than the fit to the plotted points is the assumption worth arguing about. Every application is cross-checked three ways — from-scratch (MLE + Bayes), PyMC (censored Weibull via pm.Potential, NUTS), and R (survival::survreg, which fits the identical Weibull but reports it in accelerated-failure-time form — coefficients describing how a covariate stretches or compresses the time axis rather than how it scales the hazard, so survreg's numbers are $-\beta/k$ against the ones here and agree once converted) — with R's semiparametric Cox coxph as an extra benchmark on Leuk that recovers the same HR without assuming a Weibull baseline.

Downloads

References

• Weibull, W. (1951). A statistical distribution function of wide applicability. Journal of Applied Mechanics 18(3), 293–297. — the Weibull distribution and its shape parameter
• Cox, D. R. (1972). Regression models and life-tables. Journal of the Royal Statistical Society: Series B 34(2), 187–220. — proportional hazards and the semiparametric Cox benchmark
• Kaplan, E. L. & Meier, P. (1958). Nonparametric estimation from incomplete observations. Journal of the American Statistical Association 53(282), 457–481. — the Kaplan–Meier survival curve used for comparison
• Freireich, E. J. et al. (1963). The effect of 6-mercaptopurine on the duration of steroid-induced remissions in acute leukemia. Blood 21(6), 699–716. — the leukemia remission trial (the BUGS Leuk data)
• Meeker, W. Q. & Escobar, L. A. (1998). Statistical Methods for Reliability Data. Wiley. — the shock-absorber life data and the η / B₁₀ / MTTF reliability metrics
• Kalbfleisch, J. D. & Prentice, R. L. (2002). The Statistical Analysis of Failure Time Data (2nd ed.). Wiley. — parametric survival, censoring, and PH-vs-AFT parameterisations

Sampler Module — Source Code

"""
Parametric survival regression (Weibull / exponential proportional hazards) with RIGHT-CENSORING --
from-scratch MLE + Bayesian RW-Metropolis. The BUGS Vol I "Leuk" example (Gehan leukemia, 6-MP vs
placebo) and the PyMC Bayesian-survival case study.

The new ingredient vs every earlier project is CENSORING: a subject still in remission at the end of
follow-up contributes only the survival probability S(t), not a density. With event indicator delta_i
(1 = relapse observed, 0 = censored):
    L = prod_i  h(t_i)^{delta_i} * S(t_i)
i.e. events contribute the hazard x survival, censored cases only the survival.

Weibull proportional-hazards parameterisation
---------------------------------------------
  h(t | x) = lambda_i * k * t^{k-1},   lambda_i = exp(x_i' beta),   S(t) = exp(-lambda_i t^k)
  log L = sum_i delta_i [ x_i'beta + log k + (k-1) log t_i ]  -  sum_i exp(x_i'beta) t_i^k
k = Weibull shape: k>1 increasing hazard, k<1 decreasing, k=1 = EXPONENTIAL (constant hazard).
exp(beta_j) is the HAZARD RATIO for covariate j. Note the log-lik is a Poisson GLM kernel in disguise
(delta as 0/1 "counts", offset log t^k) -- the piecewise-exponential/Poisson trick behind BUGS Leuk.

Priors (Bayes): beta ~ N(0, prior_sd^2),  log k ~ N(0, logk_sd^2)  [logk_sd defaults to 1].
The default N(0,1) on log k is mildly informative and fine for hazards near constant, but it pulls
noticeably on strong wear-out shapes (k ~ 3 sits 1.1 SD out), so the reliability example relaxes it.
"""
import numpy as np
from scipy.optimize import minimize


def simulate_weibull(n=600, beta=(-1.0, -1.3), k=1.4, cens_scale=3.0, seed=0):
    """X = [1, treat]; right-censor by an independent exponential follow-up time."""
    rng = np.random.default_rng(seed)
    beta = np.asarray(beta, float)
    treat = rng.integers(0, 2, n).astype(float)
    X = np.column_stack([np.ones(n), treat])
    lam = np.exp(X @ beta)
    T = (-np.log(rng.random(n)) / lam) ** (1.0 / k)          # Weibull event times
    C = rng.exponential(cens_scale, n)                        # random right-censoring
    t = np.minimum(T, C); delta = (T <= C).astype(float)
    return t, X, delta


def _negloglik(params, X, t, delta, fix_k):
    k = 1.0 if fix_k else np.exp(params[-1])
    beta = params if fix_k else params[:-1]
    xb = X @ beta
    ll = np.sum(delta * (xb + np.log(k) + (k - 1) * np.log(t))) - np.sum(np.exp(xb) * t ** k)
    return -ll


def _num_hess(f, x, eps=1e-4):
    n = len(x); H = np.zeros((n, n)); fx = f(x)
    for i in range(n):
        for j in range(i, n):
            xpp = x.copy(); xpp[i] += eps; xpp[j] += eps
            xpm = x.copy(); xpm[i] += eps; xpm[j] -= eps
            xmp = x.copy(); xmp[i] -= eps; xmp[j] += eps
            xmm = x.copy(); xmm[i] -= eps; xmm[j] -= eps
            H[i, j] = H[j, i] = (f(xpp) - f(xpm) - f(xmp) + f(xmm)) / (4 * eps ** 2)
    return H


def weibull_mle(X, t, delta, fix_k=False):
    k0 = X.shape[1] + (0 if fix_k else 1)
    x0 = np.zeros(k0)
    res = minimize(_negloglik, x0, args=(X, t, delta, fix_k), method='Nelder-Mead',
                   options=dict(xatol=1e-8, fatol=1e-8, maxiter=20000))
    H = _num_hess(lambda p: _negloglik(p, X, t, delta, fix_k), res.x)
    se = np.sqrt(np.diag(np.linalg.inv(H)))
    return res.x, se, -res.fun


def weibull_rwm(X, t, delta, R=12000, burn=3000, prior_sd=10.0, fix_k=False, seed=0, s=1.0,
                logk_sd=1.0):
    rng = np.random.default_rng(seed)
    mle, _, _ = weibull_mle(X, t, delta, fix_k)
    H = _num_hess(lambda p: _negloglik(p, X, t, delta, fix_k), mle)
    L = np.linalg.cholesky(np.linalg.inv(H)) * s
    p = mle.copy(); d = len(p)

    def logpost(par):
        beta = par if fix_k else par[:-1]
        lp = -0.5 * np.sum(beta ** 2) / prior_sd ** 2
        if not fix_k:
            lp += -0.5 * par[-1] ** 2 / logk_sd ** 2   # log k ~ N(0, logk_sd^2)
        return -_negloglik(par, X, t, delta, fix_k) + lp

    keep = R - burn; P = np.zeros((keep, d)); acc = 0; cur = logpost(p)
    for it in range(R):
        prop = p + L @ rng.standard_normal(d)
        lp = logpost(prop)
        if np.log(rng.random()) < lp - cur:
            p, cur = prop, lp; acc += 1
        if it >= burn:
            P[it - burn] = p
    return dict(draws=P, accept=acc / R, fix_k=fix_k)


def km_estimator(t, delta):
    """Kaplan-Meier product-limit estimate. Returns (times, survival)."""
    order = np.argsort(t); t, delta = t[order], delta[order]
    uniq = np.unique(t[delta == 1])
    n = len(t); S = 1.0; ts = [0.0]; Ss = [1.0]
    for ti in uniq:
        at_risk = np.sum(t >= ti)
        d = np.sum((t == ti) & (delta == 1))
        S *= (1 - d / at_risk); ts.append(ti); Ss.append(S)
    return np.array(ts), np.array(Ss)


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)]