#########################################################################
# This file gives example R code that uses the R function Plotconfidenceband.loghazardratio.r
# to create a plot of an estimated log hazard ratio over time, using the method of
# Gilbert et al. (2002, Biometrics).
#
# The files Plotconfidenceband.loghazardratio.r and confidenceband.loghazardratio.r
# should be placed in the same directory as this code is run.

#########################################################################
library(survival)

source('Plotconfidenceband.loghazardratio.r')

# Make up a clinical trial dataset, written to the file datafile

time <- c(rexp(250,rate=1),rweibull(250,shape=0.5,scale=0.9))
delta <- ifelse(runif(500) < .9,1,0)
Rx <- c(rep(1,250),rep(2,250))

write(t(cbind(time,delta,Rx)),file="datafile",ncolumns=3)

# Set input parameters:
N <- 1000
nboot <- 250

# Set u1=t1 to the maximum of the 2nd smallest failure times for the 2 groups
m1 <- sort(time[Rx==1 & delta==1])[2]
m2 <- sort(time[Rx==2 & delta==1])[2]
t1 <- max(m1,m2)
u1 <- t1

# Set u2=t2 to the minimum of the 2nd largest failure times for the 2 groups
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)
u2 <- t2

band11 <- (u2-u1)/3
band12 <- (u2-u1)/3

# Default settings for the rest of the parameters:
band1 <- 0
band2 <- 0
ngrid <- 50
biasadjust <- T
tailsl <- T
tailsu <- T
Nv <- 0
yliml <- 0
ylimu <- 0

# Create the plot:

answer <- Plotconfidenceband.loghazardratio('datafile',Rx,u1,u2,t1,t2,N,band11,
band1,band12,band2,ngrid,biasadjust,tailsl,tailsu,nboot,Nv,yliml,ylimu)

# The created plot is placed in the file figureloghazratio.ps, in the same directory.