Bootstrapping the Proportional Hazards Check

The following R function fits a Cox proportional hazards (PH) model, checks the PH assumption (using Schoenfeld residuals and the R function cox.zph), if PH violated corrects with stratification, and then bootstraps this entire model selection process to allow computation of valid confidence intervals of the predicted survival.

This function is a work in progress, only contains the bare minimums, and may have to be altered for an individual analysis.

Usage

Arguments

formula

a character expression (given in quotations) that is converted to a formula object (to be called by coxph() ), with the response on the left of a '~' operator, and the terms on the right. The response must be a survival object as returned by the 'Surv' function.

number.strata: an integer of at least 2 which will determine how many strata to divide continuous predictors into if they violate the PH assumption. ph.crit.val: if the GLOBAL test of PH in cox.zph is less than this value, then the model stratifies over the predictor which contributed most to this violation (as evidenced by the lowest p-value in cox.zph). surv.time: the time at which to estimate the survival probability surv.x: the levels of predictor variables Nboots: number of bootstrap replications

Details All categorical variables need to be specified as factor variables. If stratification occurs over a factor variable, then the levels of this variable are chosen as factors. If stratification occurs over a continuous variable, then strata are created by dividing observations into strata based on the ordered value of the continuous predictor and the number of strata specified by the user (number.strata). For example, if number.strata=5, then 5 strata are used. The continuous variable which was used for stratification is left in the model to account for any residual information. This function returns surv.est, the estimated survival at a given time (surv.time) for given predictor variables (surv.x), as well as surv.estb, this estimate for all Nboots bootstrap replications. A confidence interval can be computed using quantiles of surv.estb. The function also returns badcovb, which specifies which predictors were used for stratification in each bootstrap replication, 0 if no stratification.

create.strata<-function(a,number){ floor(((rank(a,ties.method="first")-1)/length(a))*number) }

coxphboot<-function(formula,number.strata,ph.crit.val,surv.time,surv.x,Nboots) {

formula<-as.formula(formula) mod<-coxph(formula) express1<-as.character(formula) S<-eval(parse(text=express1[2])) z<-predict(mod, type='terms') f.short<-coxph(S ~ z) phtest<-survival::cox.zph(f.short, transform='identity')$table

badcov<-0 if(phtest[ncol(z)+1,3]<ph.crit.val) { badcov<-which(phtest[,3]==min(phtest[1:ncol(z),3])) }

m<-Design(model.frame(formula)) X<-m[,-1] classes<-attributes(m)$Design$assume.code Xnew<-X Xnew[1,]<-surv.x

if (badcov>0) { if (classes[badcov]==5) { stratnum<-X[,badcov] tmp1<-strsplit(express1, "\\ ") tmp13[1+(badcov-1)*2]<-"strata(stratnum)" express3<-formula(paste(express1[2],express1[1],paste(tmp13,collapse=""))) mod<-coxph(express3) } if (classes[badcov]==1) { stratnum<-create.strata(X[,badcov],number.strata) express2<-"+strata(stratnum)" express3<-as.formula(paste(express1[2],express1[1],express1[3],express2)) mod<-coxph(express3) } }

if (badcov==0) { stuff<-summary(survival::survfit(mod,newdata=Xnew[1,])) surv.est<-stuff$surv[stuff$time==max(stuff$time[stuff$time<surv.time])] } if (badcov>0) { if (classes[badcov]==5){ stuff<-summary(survival::survfit(mod,newdata=Xnew[1,][-badcov])) surv.strat<-surv.x[badcov] } if (classes[badcov]==1){ stuff<-summary(survival::survfit(mod,newdata=Xnew[1,])) surv.strat<-create.strata(c(X[,badcov],surv.x[badcov]),number.strata)[length(X[,badcov])+1] } surv.est<-stuff$surv[stuff$time==max(stuff$time[stuff$time<surv.time& stuff$strata==paste("stratnum=",as.character(surv.strat),sep="")])& stuff$strata==paste("stratnum=",as.character(surv.strat),sep="")] }

surv.estb<-badcovb<-NULL for (j in 1:Nboots) { cat(" Bootstrap: ", j, "\n")

N<-length(X[,1]) bootsamp<-sample(c(1:N),N,replace=TRUE)

yb<-S[,1][bootsamp] db<-S[,2][bootsamp] Sb<-Surv(yb,db) Xb<-X[bootsamp,1:ncol(X)] rownames(Xb)<-1:N Xb2 <- data.frame(cbind(S = Sb, Xb))

formula1<-as.formula(paste("Sb~",express1[3])) modb <- coxph(formula1, data=Xb2) z<-predict(modb, type='terms') f.short<-coxph(Sb ~ z) phtest<-survival::cox.zph(f.short, transform='identity')$table

badcov<-0 if(phtest[ncol(z)+1,3]<ph.crit.val) { badcov<-which(phtest[,3]==min(phtest[1:ncol(z),3])) }

if (badcov>0) { if (classes[badcov]==5) { Xb2$stratnum<-Xb[,badcov] tmp1<-strsplit(express1, "\\ ") tmp13[1+(badcov-1)*2]<-"strata(stratnum)" express3<-formula(paste("Sb~",paste(tmp13,collapse=""))) modb<-coxph(express3, data=Xb2) } if (classes[badcov]==1) { Xb2$stratnum<-create.strata(Xb[,badcov],number.strata) express2<-"+strata(stratnum)" express3<-as.formula(paste("Sb~", express1[3], express2)) modb<-coxph(express3, data=Xb2) } }

if (badcov==0) { stuff<-summary(survival::survfit(modb,newdata=Xnew[1,])) surv.estb[j]<-stuff$surv[stuff$time==max(stuff$time[stuff$time<surv.time])] } if (badcov>0) { if (classes[badcov]==5){ stuff<-summary(survival::survfit(modb,newdata=Xnew[1,][-badcov])) surv.strat<-surv.x[badcov] } if (classes[badcov]==1){ stuff<-summary(survival::survfit(modb,newdata=Xnew[1,])) surv.strat<-create.strata(c(Xb[,badcov],surv.x[badcov]),number.strata)[length(Xb[,badcov])+1] } surv.estb[j]<-stuff$surv[stuff$time==max(stuff$time[stuff$time<surv.time& stuff$strata==paste("stratnum=",as.character(surv.strat),sep="")])& stuff$strata==paste("stratnum=",as.character(surv.strat),sep="")] }

tmp1<-strsplit(express1, "\\ ") badcovb[j]<-badcov if (badcov>0) { badcovb[j]<-tmp13[1+(badcov-1)*2] } } return(surv.est,surv.estb,badcovb) }