Survival Variable Selection with Censored Data

Python · R · Part 10 of the variable-selection arc — time-to-event outcomes

Overview

Part 10 of the variable-selection arc: time-to-event outcomes. The last outcome type brings a complication no earlier notebook faced — censoring. When a study ends, many patients are still alive: we know only that their survival time exceeds their follow-up, not what it is. A selection method must use that partial information correctly or it will be badly biased. The classic testbed is the Mayo Clinic PBC data (primary biliary cirrhosis): 276 patients, 16 candidate covariates, only 40% observed deaths. The famous Mayo "natural history" model, built by painstaking stepwise selection, kept exactly five — bilirubin, albumin, age, edema, prothrombin time. Can Bayesian variable selection recover them? We fit a Weibull proportional-hazards model with a spike-and-slab prior (Part 1's method extended to censored data), cross-checked against PyMC and R's BMA::bic.surv (model averaging over Cox models).

\log L = \sum_i\Big[\,\delta_i\big(\log\alpha + (\alpha-1)\log t_i + \eta_i\big) - t_i^{\alpha}e^{\eta_i}\,\Big], \qquad \eta_i = \beta_0 + x_i'\beta

Weibull proportional hazards, and how censoring enters

A Weibull proportional-hazards model with hazard $h(t\mid x)=\alpha t^{\alpha-1}\exp(\beta_0+x'\beta)$ , an event indicator $\delta_i$ (1 = death observed, 0 = censored) and follow-up $t_i$ . In the log-likelihood the cumulative-hazard term $t_i^\alpha e^{\eta_i}$ is subtracted for every subject, while only an observed death adds the log-hazard — which is exactly how a censored patient contributes "survived at least this long" ( $\log S(t_i)=-t_i^\alpha e^{\eta_i}$ ) rather than "died at $t_i$ ". A spike-and-slab prior sits on the coefficients as before. The Weibull likelihood is not conjugate, so $\beta$ (and the baseline $\beta_0,\alpha$ ) are updated by adaptive random-walk Metropolis with the spike/slab prior variance set by the current $\gamma$ — but the selection step is unchanged from Part 1: $\gamma_j\mid\beta_j$ is a Bernoulli from the spike/slab density ratio, because that conditional never sees the likelihood. So the whole apparatus is Part-1 SSVS with a Metropolis coefficient move and a survival log-likelihood. Strongly right-skewed labs are log-transformed and all covariates standardized. Engine: a self-contained survsel.py.

Results

The Mayo model, recovered. The from-scratch Weibull-PH search selects all five factors of the hand-built Mayo natural-history model — bilirubin (inclusion 1.00), age, edema, albumin and prothrombin time (the last just over the line at 0.52) — alongside histologic stage and urinary copper, both established PBC markers. The search found the accepted prognostic set from a Bayesian posterior, censoring and all, and the Kaplan–Meier curves confirm it: split by bilirubin tertile or by edema, the observed survival curves separate dramatically. A famous null result: treatment (D-penicillamine) has inclusion ≈ 0.10 — not selected — matching the trial's conclusion that the drug did not improve survival; selection correctly declines an ineffective covariate. Three engines, one prognostic set: the from-scratch Weibull-PH SSVS, PyMC's version, and R's bic.surv (which averages over Cox models) select the same core six — bilirubin, age, edema, albumin, copper and histologic stage — with prothrombin time right on the selection boundary (in for the Weibull at 0.52, just shy for PyMC at ≈0.47 and Cox BMA at 0.41). That a fully parametric Weibull and a semiparametric Cox agree on which covariates matter confirms the selection is a property of the proportional-hazards structure, not the baseline-hazard form. Censoring, handled cleanly: with 60% of times censored, the cumulative-hazard term does the work — a censored patient contributes the probability of surviving past their follow-up, and the selection is unbiased for it (treating censored times as deaths would corrupt every inclusion probability). This completes the arc's outcome types — continuous, counts and tables, multivariate time series, binary, and now time-to-event — the same spike-and-slab idea traveling to each by swapping in the appropriate likelihood, with the inclusion step untouched.

Notebooks

Python Survival SSVS from scratch — a Weibull proportional-hazards spike-and-slab with the censored likelihood and Metropolis coefficient moves, plus a PyMC cross-check, on the Mayo PBC data (Python) Open → R R cross-check — BMA::bic.surv, Bayesian model averaging over Cox proportional-hazards models (R) Open →

Downloads

survsel.py The from-scratch survival-selection engine — Weibull-PH censored log-likelihood, Metropolis coefficient / baseline moves, and the Part-1 spike-and-slab indicator step (new self-contained module) pbc_data.csv Mayo Clinic PBC trial — 276 patients, 16 covariates, time to death with right-censoring (survival::pbc)

References

George, E. I. & McCulloch, R. E. (1993). Variable selection via Gibbs sampling. Journal of the American Statistical Association 88(423), 881–889. — the spike-and-slab SSVS carried over to the censored likelihood
Volinsky, C. T., Madigan, D., Raftery, A. E. & Kronmal, R. A. (1997). Bayesian model averaging in proportional hazard models. Journal of the Royal Statistical Society C 46(4), 433–448. — the R BMA::bic.surv Cox model averaging
Cox, D. R. (1972). Regression models and life-tables. Journal of the Royal Statistical Society B 34(2), 187–220. — proportional-hazards regression
Kaplan, E. L. & Meier, P. (1958). Nonparametric estimation from incomplete observations. Journal of the American Statistical Association 53(282), 457–481. — the survival-curve estimator
Dickson, E. R., Grambsch, P. M., Fleming, T. R., Fisher, L. D. & Langworthy, A. (1989). Prognosis in primary biliary cirrhosis: model for decision making. Hepatology 10(1), 1–7. — the Mayo natural-history model and its five prognostic factors
Therneau, T. M. & Grambsch, P. M. (2000). Modeling Survival Data: Extending the Cox Model. Springer. — the PBC dataset (survival::pbc) and proportional-hazards modelling