library(survival) library(KMsurv) data(bmt) head(bmt) #t2 disease free survival. Time to relapse or death #tp time to platelet recovery #group different patientgroups #convert data to start-stop format. Two lines for each patient where status of platelet recovery changes during disease free survival time bmt.timedep=rbind(bmt,bmt) n=dim(bmt)[1] status.timedep=zp.timedep=bmt.start=bmt.stop=rep(0,2*n) count=0 for (i in 1:n){ if ((bmt$tp[i]0)){#two rows - platelet recovery observed before t2 bmt.timedep[count+1,]=bmt[i,]#repeat fixed covariates bmt.timedep[count+2,]=bmt[i,] bmt.start[count+1]=0 bmt.stop[count+1]=bmt$tp[i] bmt.start[count+2]=bmt$tp[i] bmt.stop[count+2]=bmt$t2[i] zp.timedep[count+1]=0 zp.timedep[count+2]=1 status.timedep[count+1]=0 status.timedep[count+2]=bmt$d3[i] count=count+2 } if (bmt$tp[i]>=bmt$t2[i]){#one row - platelet recovery not observed bmt.timedep[count+1,]=bmt[i,] bmt.start[count+1]=0 bmt.stop[count+1]=bmt$t2[i] zp.timedep[count+1]=0 status.timedep[count+1]=bmt$d3[i] count=count+1 } if (bmt$tp[i]==0){#one row platelet recovery already occurred at time 0 bmt.timedep[count+1,]=bmt[i,] bmt.start[count+1]=0 bmt.stop[count+1]=bmt$t2[i] zp.timedep[count+1]=1 status.timedep[count+1]=bmt$d3[i] count=count+1 } } bmt.timedep=data.frame(cbind(bmt.start,bmt.stop,status.timedep,zp.timedep,bmt.timedep$group)[1:count,]) bmt.timedep[1:10,] names(bmt.timedep)=c("start","stop","status","zp","group") head(bmt.timedep) #fit model with fixed covariate using the two data sets. fitorig=coxph(Surv(t2,d3)~factor(group),data=bmt) summary(fitorig) fit=coxph(Surv(start,stop,status)~factor(group),data=bmt.timedep) summary(fit) #identical results. #forkert analyse hvor vi blot bruger indicator dp for recovery fitwrong=coxph(Surv(t2,d3)~factor(group)+dp,data=bmt) summary(fitwrong) #now model with timedependent platelet fitzp=coxph(Surv(start,stop,status)~factor(group)+zp,data=bmt.timedep) summary(fitzp) factor(group)2 -0.4971 0.6083 0.2892 -1.719 0.085656 . factor(group)3 0.3818 1.4650 0.2676 1.427 0.153626 zp -1.1202 0.3262 0.3292 -3.403 0.000666 *** #slight discrepancy from results in KM Table 9.1 page 299. #models dp overestimates benefit of platelet recovery fitwrong$loglik[2] -356.1788 fitzp$loglik[2] -361.8592 #likelihood biggest for wrong model. It fits better in terms of likelihood but does not make sense. #Why does it not make sense to use indicator for platelets returned to normal ? #Suppose this had no impact of hazard. Then we would still see positive association between platelet recovery and survival time - simply because if a patient dies early, there is less chance of seeing platelet recovery. #Now check for proportional hazards using timedependent, cf. Table 9.5 page 304. z2group=as.numeric(bmt$group==2) z3group=as.numeric(bmt$group==3) fit.group.time=coxph(Surv(t2,d3)~z2group+z3group+tt(z2group)+tt(z3group),data=bmt,tt=function(x,t,...) x*log(t)) summary(fit.group.time) fit.group=coxph(Surv(t2,d3)~z2group+z3group,data=bmt) summary(fit.group) -2*(fit.group$loglik-fit.group.time$loglik)[2]#second component is likelihood at final values. First componentis at initial values. 1.943971 on 2 df.#KM 1.735 Wald statistic. #let us try sex patient z3 and sex donor z4 z3z4=bmt$z3*bmt$z4 fit.sex.time=coxph(Surv(t2,d3)~z3+z4+z3z4+tt(z3)+tt(z4)+tt(z3z4),data=bmt,tt=function(x,t,...) x*log(t)) summary(fit.sex.time) fit.sex=coxph(Surv(t2,d3)~z3+z4+z3z4,data=bmt) -2*(fit.sex$loglik-fit.sex.time$loglik)[2] #0.225. KM Wald: 0.220 #Waiting time to transplant z7 fit.wait=coxph(Surv(t2,d3)~z7,data=bmt) fit.wait.timedep=coxph(Surv(t2,d3)~z7+tt(z7),data=bmt,tt=function(x,t,...) x*log(t)) summary(fit.wait.timedep) #z=-0.071 z^2=Wald=0.005041 -2*(fit.wait$loglik-fit.wait.timedep$loglik)[2] #0.004956894. KM Wald: 0.005 #Waiting time z7*t instead of z7*logt: fit.wait=coxph(Surv(t2,d3)~z7,data=bmt) fit.wait.timedep=coxph(Surv(t2,d3)~z7+tt(z7),data=bmt,tt=function(x,t,...) x*t) summary(fit.wait.timedep) #same conclusion. #let us try age z1 (age in years) fit.age=coxph(Surv(t2,d3)~z1,data=bmt) fit.age.timedep=coxph(Surv(t2,d3)~tt(z1),data=bmt,tt=function(x,t,...) x*365.25+t) #same model except for different scaling of age. coef(fit.age.timedep)*365.25 0.01122251 coef(fit.age) 0.01122251 ################################################################### wgroup Disease Group 1-ALL, 2-AML Low Risk, 3-AML High Risk t1 Time To Death Or On Study Time t2 Disease Free Survival Time (Time To Relapse, Death Or End Of Study) d1 Death Indicator 1-Dead 0-Alive d2 Relapse Indicator 1-Relapsed, 0-Disease Free d3 Disease Free Survival Indicator 1-Dead Or Relapsed, 0-Alive Disease Free) ta Time To Acute Graft-Versus-Host Disease da Acute GVHD Indicator 1-Developed Acute GVHD 0-Never Developed Acute GVHD) tc Time To Chronic Graft-Versus-Host Disease dc Chronic GVHD Indicator 1-Developed Chronic GVHD 0-Never Developed Chronic GVHD tp Time To Chronic Graft-Versus-Host Disease#time to platelet recovery dp Platelet Recovery Indicator 1-Platelets Returned To Normal, 0-Platelets Never Returned to Nor- mal z1 Patient Age In Years z2 Donor Age In Years z3 Patient Sex: 1-Male, 0-Female z4 Donor Sex: 1-Male, 0-Female z5 Patient CMV Status: 1-CMV Positive, 0-CMV Negative z6 Donor CMV Status: 1-CMV Positive, 0-CMV Negative z7 Waiting Time to Transplant In Days z8 FAB: 1-FAB Grade 4 Or 5 and AML, 0-Otherwise z9 Hospital: 1-The Ohio State University, 2-Alferd , 3-St. Vincent, 4-Hahnemann z10 MTX Used as a Graft-Versus-Host- Prophylactic: 1-Yes, 0-No