Bootstrapping coxph after model averaging

Question

I have selected 7 coxph models using AIC, and did model averaging to obtain model averaged parameter estimates. I wanted to plot the survival curve from the averaged model, and I found someone with my problem, whose answer was given by Terry Therneau:

You can fit a Cox model with fixed coefficients. Assume "fixbeta" are the coefficients from your model averaging, then do

ffit <- coxph(Surv(time, status) ~ x1 + x2 + .... , data=mydata,
init=fixbeta, iter=0)
sfit <- survfit(fit)

The standard errors in sfit are incorrect of course. One could bootstrap the entire model creation process to get accurate values.

My question is: how do I "bootsrap the entire model creation process" using R?

Answer 1

You have made two false assumptions: fitting more than one model is a good idea, and variable selection using AIC is a good idea. Model averaging, when the models being averaged are all from the same family (here, Cox PH), is an effort-intensive way to obtain a penalized single model. A penalized full model will perform as well or better than model averaging, and the process is easier to bootstrap and much easier to interpret. You can use the R survival package's coxph for penalization.

First check whether you need to use penalization (shrinkage) at all. How many events do you have? How many candidate predictors? Are you sure all predictors act linearly? Have you considered much simpler data reduction procedures (masked to survival time)?