Survival selection — R cross-check

Bayesian variable selection, Part 10 (R)

Validates surv_python.ipynb. BMA::bic.surv (Volinsky et al. 1997) does Bayesian model averaging over Cox proportional-hazards models: it searches the model space, scores each Cox model by BIC (a Laplace approximation to its marginal likelihood), and returns posterior inclusion probabilities and the top models. It is an independent engine and a different survival model (semiparametric Cox vs the from-scratch parametric Weibull), so agreement on the selected covariates shows the result depends on the proportional-hazards structure, not the baseline-hazard form.

Data: pbc_data.csv — Mayo PBC: 276 patients, 16 covariates, time to death with right-censoring.

.libPaths('C:/Users/user/R/win-library/4.6'); suppressMessages({library(BMA); library(survival)})
d <- read.csv('pbc_data.csv')
for (v in c('bili','alk.phos','ast','copper','trig','protime')) d[[v]] <- log(d[[v]])   # skewed labs -> log
cov <- c('age','sex','ascites','hepato','spiders','edema','bili','albumin',
         'copper','alk.phos','ast','trig','platelet','protime','stage','trt')

b <- bic.surv(x=d[,cov], surv.t=d$years, cens=d$event, strict=FALSE, OR=20)   # BMA over Cox models
ip <- b$probne0; names(ip) <- b$namesx
cat('bic.surv posterior inclusion probabilities (%):\n'); print(round(sort(ip, decreasing=TRUE)))
cat('\nmost probable model:', paste(b$namesx[b$which[1, ]], collapse=' + '), '\n')
sel_r <- sort(names(ip)[ip > 50])                                 # bic.surv selected covariates (probne0 > 50%)
py    <- sort(c('bili','age','edema','albumin','copper','stage')) # Python Weibull SSVS consensus set
cat(sprintf('bic.surv selects (probne0 > 50%%): %s\n', paste(sel_r, collapse=', ')))
cat(sprintf('cross-check -> Python Weibull SSVS: %d of %d covariates shared (prothrombin time on the boundary in both engines).\n',
            length(intersect(sel_r, py)), length(py)))

options(repr.plot.width=8, repr.plot.height=5.2)
ips <- sort(ip)/100; sel <- ips > 0.5                             # highlight the selected set, derived from the fit
barplot(ips, horiz=TRUE, las=1, xlim=c(0,1), col=ifelse(sel,'firebrick','grey65'),
        xlab='posterior inclusion probability', main='bic.surv (Cox BMA): PBC prognostic factors')
abline(v=0.5, lty=2, col='grey40')

Warning message:
"package 'BMA' was built under R version 4.6.1"

Warning message:
"package 'leaps' was built under R version 4.6.1"

Warning message:
"package 'robustbase' was built under R version 4.6.1"

Warning message:
"package 'inline' was built under R version 4.6.1"

Warning message:
"package 'rrcov' was built under R version 4.6.1"

bic.surv posterior inclusion probabilities (%):

    bili      age    edema   copper  albumin    stage  protime      ast 
     100       98       86       83       81       61       41       14 
     sex  ascites  spiders     trig platelet alk.phos   hepato      trt 
       4        4        4        3        3        3        3        3 

most probable model: age + edema + bili + albumin + copper + stage 

bic.surv selects (probne0 > 50%): age, albumin, bili, copper, edema, stage

cross-check -> Python Weibull SSVS: 6 of 6 covariates shared (prothrombin time on the boundary in both engines).

Same prognostic factors, a Cox model

• bic.surv reproduces the selection. Bilirubin (100%), age (98%), edema (86%), copper (83%), albumin (81%) and stage (61%) are selected; treatment, sex and the rest sit near zero — matching the from-scratch Weibull-PH SSVS and PyMC. The most probable Cox model is age + edema + bili + albumin + copper + stage.
• Model-form robustness. The from-scratch engine assumes a parametric Weibull baseline hazard; bic.surv uses the semiparametric Cox partial likelihood. They select the same covariates because variable selection under proportional hazards is about the linear predictor $x'\beta$, which both share — the baseline hazard cancels out of the comparison.

References

• Volinsky, C. T., Madigan, D., Raftery, A. E. & Kronmal, R. A. (1997). Bayesian model averaging in proportional hazard models. Applied Statistics 46, 433–448.
• Raftery, A. E. (1995). Bayesian model selection in social research. Sociological Methodology 25, 111–163.
• Dickson, E. R. et al. (1989). Prognosis in primary biliary cirrhosis. Hepatology 10, 1–7.