#################################################################################### # R function to generate conditional survival plots # # Function arguments: # .basekm - an object of class survfit # .at - a vector of times on which to condition # .main - main title for the plot # .xlab - label for the x-axis # .ylab - label for the y-axis, defaults to "Survival probability" # .lwd - line width, defaults to 1 # .mark - controls whether censoring times are marked on the curves, defaults to F #################################################################################### cond.surv=function(.basekm,.at,.main=NULL,.xlab=NULL,.ylab="Survival probability", .lwd=1,.mark=F){ if(class(.basekm) != "survfit"){stop("Argument to .basekm must be of class survfit")} if(max(.at)>max(.basekm\$time)){stop(paste( "Argument to .at specifies value(s) outside the range of observed times;", "the maximum observed time is", round(max(.basekm\$time),2)))} plot(.basekm,conf.int=F,xlab=.xlab,ylab=.ylab,main=.main,lwd=.lwd,mark.time=.mark) nt=length(.at) fitkm=list() for(i in 1:nt){ fitkm[[i]]=survfit(formula=as.formula(.basekm\$call\$formula), data=eval(.basekm\$call\$data), start.time=.at[i]) lines(fitkm[[i]],conf.int=F,col=i+1,lwd=.lwd,mark.time=.mark) abline(v=i,lty=3) } } #################################################################################### # R function to calculate conditional survival with a 95% confidence interval # # Function arguments: # .basekm - an object of class survfit # .t1 - the time on which to condition # .t2 - the survival time to estimate # # e.g. if .t1=2 and .t2=5, this will give you the probability of surviving to # year 5 conditional on having already survived to year 2 #################################################################################### cs.ci=function(.basekm, .t1, .t2){ if(class(.basekm) != "survfit"){stop("Argument to .basekm must be of class survfit")} if(max(.t1)>max(.basekm\$time)){stop(paste( "Argument to .t1 specifies a value outside the range of observed times;", "the maximum observed time is", round(max(.basekm\$time),2)))} if(max(.t2)>max(.basekm\$time)){stop(paste( "Argument to .t2 specifies a value outside the range of observed times;", "the maximum observed time is", round(max(.basekm\$time),2)))} cs<-summary(.basekm, times=c(.t1, .t2))\$surv[2]/summary(.basekm, times=c(.t1, .t2))\$surv[1] cs.sq<-cs^2 d<-.basekm\$n.event[.basekm\$time>=.t1 & .basekm\$time<=.t2 & .basekm\$n.event>0] r<-.basekm\$n.risk[.basekm\$time>=.t1 & .basekm\$time<=.t2 & .basekm\$n.event>0] dr<-d/(r*(r-d)) var.cs<-cs.sq*sum(dr) ci<-cs+c(-1,1)*(qnorm(0.975)*sqrt(var.cs)) if(ci[1]<0){warning("Lower bound of CI has been truncated to 0")} if(ci[2]>1){warning("Upper bound of CI has been truncated to 1")} ci.cs<-round(ci,2) if(ci.cs[1]<0){ci.cs[1]<-0} if(ci.cs[2]>1){ci.cs[2]<-1} return(list(cs=round(cs,2), cilow=ci.cs[1], cihigh=ci.cs[2])) }