Bootstrapping coxph after model averaging

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?

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.