######################################################################### # R FUNCTION: Plotconfidenceband.loghazardratio.r # # DESCRIPTION: When sourced in R, this function: # # (1) Inputs a 2-arm clinical trial dataset in the file datafile; # # (2) Calls the function confidenceband.loghazardratio() contained in the file # confidenceband.loghazardratio.r, to compute the estimated log hazard ratio # (treatment 2/treatment 1) and 95% pointwise and 80% and 95% simultaenous confidence # intervals about the log hazard ratio over time; and # # (3) Plots the estimated log hazard ratio over time with 95% pointwise and # 80%, 95% simultaneous confidence intervals; this plot is written to the # file figureloghazardratio.ps # # The method used to compute the simultaneous confidence intervals is described # in Gilbert P, Wei LJ, Kosorok MR, Clemens JD. Simultaneous inference on the # contrast of two hazard functions with censored observations. Biometrics, 2002; # 58:773-780. ######################################################################### # # INPUT variables into the function Plotconfidenceband.loghazardratio.r # # datafile = The name of the text file that contains the dataset, with # n rows (n = number of subjects in the trial) and 3 columns. # First column: the minimum of the censoring time and failure time. # Second column: failure indicator: 0 if censored, 1 if event. # Third column: treatment indicator: 1 = treatment 1 (e.g., control), # 2 = treatment 2 (e.g., active). This variable must be coded as 1 and 2. # # [u1,u2] = the time interval over which the confidence band is computed # [t1,t2] = the time interval over which the confidence band could be computed (at least # as wide as [u1,u2]). t1 should be chosen larger than both of the smallest # infection times for each group. t2 should be chosen smaller than both of # the largest infection times for each group. # In practice it is often fine to set t1=u1; t2=u2 # N = number of simulated Gaussian variates used for computing the # simultaneous confidence bands (N = 1000 is adequate) # band11 = second derivative bandwidth for treatment 1, as on page 249 of Andersen, Borgan, # Gill, and Keiding (1993), used to calculate the optimal bandwidth band1 # for the treatment 1 sample. # If band11 = 0, the optimal band11 is calculated using a bootstrap procedure. # band1 = the bandwidth for the first sample. The default value is # the optimal bandwidth, calculated using the formula in ABGK (1993) # band12 = second derivative bandwidth for treatment 2, as on page 249 of Andersen, Borgan, # Gill, and Keiding (1993), used to calculate the optimal bandwidth band2 # for the treatment 2 sample. # If band12 = 0, the optimal band12 is calculated using a bootstrap procedure. # band2 = the bandwidth for the second sample. The default value is # the optimal bandwidth, calculated using the formula in ABGK (1993) # ngrid = the number of time-points spanning [t1,t2] for computing # the log hazard ratio and the confidence bands. The default value is 50. # biasadjust = T means that the hazard rate estimates are adjusted by the asymptotic # bias correction term (4.2.27) in ABGK. # tailsl = T means that the method described on page 251 of # ABGK is used to calculate the log hazard ratio estimate # and the confidence intervals/bands in the lower tail [t1,u1] # tailsu = T means that the method described on page 251 of # ABGK is used to calculate the log hazard ratio estimate # and the confidence intervals/bands in the upper tail [u2,t2] # nboot = Number of bootstrap samples used for the cross-validation # procedure for estimating the bandwidth parameters used to # estimate the second derivatives band11 of the hazard functions # Nv = 0 implies that an asymptotic variance formula will be used for the # variance function v() # > 0 implies that the simulation procedure will be used to calculate # the variance function v(). In this case, Nv is the number of iterations # used to calculate the simulated estimate of v. # yliml = lower limit of y-axis. If yliml=0, then it is set at the value of # the lower 95% simultaneous confidence band at t1. # ylimu = upper limit of y-axis. If ylimu=0, then it is set at the value of # the upper 95% simultaneous confidence band at t2. ################################################################################### Plotconfidenceband.loghazardratio <- function(datafile,Rx,u1,u2,t1,t2,N=1000,band11=0, band1=0,band12=0,band2=0,ngrid=50,biasadjust=T,tailsl=T,tailsu=T,nboot=500,Nv=0,yliml=0,ylimu=0) { tempmat <- matrix(scan(datafile),ncol=3,byrow=T) time <- tempmat[,1] delta <- tempmat[,2] Rx <- tempmat[,3] # Make sure t1<=u1 and u2<=t2 if (t1 > u1) { u1 <- t1 } if (t2 > u2) { u2 <- t2 } # Make sure t1 is big enough m1 <- sort(time[Rx==1 & delta==1])[2] m2 <- sort(time[Rx==2 & delta==1])[2] t1 <- max(m1,m2) # Make sure t2 is small enough m1 <- sort(time[Rx==1 & delta==1])[length(sort(time[Rx==1 & delta==1]))-1] m2 <- sort(time[Rx==2 & delta==1])[length(sort(time[Rx==2 & delta==1]))-1] t2 <- min(m1,m2) ## Compute the estimates and bands: source('confidenceband.loghazardratio.r') answermat <- confidenceband.loghazardratio(time,delta,Rx,u1,u2,t1,t2,N,band11,band1,band12,band2, ngrid,biasadjust,tailsl,tailsu,nboot,Nv) time <- answermat[,1] loghazratio <- answermat[,2] lowint <- answermat[,3] upint <- answermat[,4] lowband <- answermat[,5] upband <- answermat[,6] lowband80 <- answermat[,9] upband80 <- answermat[,10] yliml <- ifelse(yliml==0,quantile(lowband,.025),yliml) ylimu <- ifelse(ylimu==0,quantile(upband,.975),ylimu) postscript("figureloghazratio.ps",horizontal=T) par(las=1,cex=1.2) plot(time,lowband80,xlab="Time",ylab="Log hazard ratio",ylim=c(yliml,ylimu),type='n') lines(time,loghazratio,lwd=2.3,col="red") lines(time[!is.na(lowint)],lowint[!is.na(lowint)],lty=2,lwd=2.3,col="navy") lines(time[!is.na(upint)],upint[!is.na(upint)],lty=2,lwd=2.3,col="navy") lines(time[!is.na(lowband80)],lowband80[!is.na(lowband80)],lty=3,lwd=2.3,col="purple") lines(time[!is.na(upband80)],upband80[!is.na(upband80)],lty=3,lwd=2.3,col="purple") lines(time[!is.na(lowband)],lowband[!is.na(lowband)],lty=4,lwd=2.3,col="green") lines(time[!is.na(upband)],upband[!is.na(upband)],lty=4,lwd=2.3,col="green") abline(h=0) legend(quantile(time,.05),yliml + 0.2*(ylimu-yliml),legend=c("95% pointwise CIs","80% simultaneous CIs","95% simultaneous CIs"),lty=c(2,3,4),col=c("navy","purple","green"),lwd=2.3) title(main="Estimate of log hazard ratio over time",cex=1.2) dev.off() } #############################################################################