### Name: censboot
### Title: Bootstrap for Censored Data
### Aliases: censboot cens.return
### Keywords: survival

### ** Examples

data(aml, package="boot")
library(survival)
# Example 3.9 of Davison and Hinkley (1997) does a bootstrap on some
# remission times for patients with a type of leukaemia.  The patients
# were divided into those who received maintenance chemotherapy and 
# those who did not.  Here we are interested in the median remission 
# time for the two groups.
aml.fun <- function(data) {
     surv <- survfit(Surv(time, cens)~group, data=data)
     out <- NULL
     st <- 1
     for (s in 1:length(surv$strata)) {
          inds <- st:(st+surv$strata[s]-1)
          md <- min(surv$time[inds[1-surv$surv[inds]>=0.5]])
          st <- st+surv$strata[s]
          out <- c(out,md)
     }
     out
}
aml.case <- censboot(aml,aml.fun,R=499,strata=aml$group)

# Now we will look at the same statistic using the conditional 
# bootstrap and the weird bootstrap.  For the conditional bootstrap 
# the survival distribution is stratified but the censoring 
# distribution is not. 

aml.s1 <- survfit(Surv(time,cens)~group, data=aml)
aml.s2 <- survfit(Surv(time-0.001*cens,1-cens)~1, data=aml)
aml.cond <- censboot(aml,aml.fun,R=499,strata=aml$group,
     F.surv=aml.s1,G.surv=aml.s2,sim="cond")

# For the weird bootstrap we must redefine our function slightly since
# the data will not contain the group number.
aml.fun1 <- function(data,str) {
     surv <- survfit(Surv(data[,1],data[,2])~str)
     out <- NULL
     st <- 1
     for (s in 1:length(surv$strata)) {
          inds <- st:(st+surv$strata[s]-1)
          md <- min(surv$time[inds[1-surv$surv[inds]>=0.5]])
          st <- st+surv$strata[s]
          out <- c(out,md)
     }
     out
}
aml.wei <- censboot(cbind(aml$time,aml$cens),aml.fun1,R=499,
     strata=aml$group, F.surv=aml.s1,sim="weird")

# Now for an example where a cox regression model has been fitted
# the data we will look at the melanoma data of Example 7.6 from 
# Davison and Hinkley (1997).  The fitted model assumes that there
# is a different survival distribution for the ulcerated and 
# non-ulcerated groups but that the thickness of the tumour has a
# common effect.  We will also assume that the censoring distribution
# is different in different age groups.  The statistic of interest
# is the linear predictor.  This is returned as the values at a
# number of equally spaced points in the range of interest.
data(melanoma, package="boot")
library(splines)# for ns
mel.cox <- coxph(Surv(time,status==1)~ns(thickness,df=4)+strata(ulcer),
     data=melanoma)
mel.surv <- survfit(mel.cox)
agec <- cut(melanoma$age,c(0,39,49,59,69,100))
mel.cens <- survfit(Surv(time-0.001*(status==1),status!=1)~
     strata(agec),data=melanoma)
mel.fun <- function(d) { 
     t1 <- ns(d$thickness,df=4)
     cox <- coxph(Surv(d$time,d$status==1) ~ t1+strata(d$ulcer))
     ind <- !duplicated(d$thickness)
     u <- d$thickness[!ind]
     eta <- cox$linear.predictors[!ind]
     sp <- smooth.spline(u,eta,df=20)
     th <- seq(from=0.25,to=10,by=0.25)
     predict(sp,th)$y
}
mel.str<-cbind(melanoma$ulcer,agec)
# this is slow!
system.time(mel.mod <- censboot(melanoma,mel.fun,R=999,F.surv=mel.surv,
     G.surv=mel.cens,cox=mel.cox,strata=mel.str,sim="model"))
# To plot the original predictor and a 95% pointwise envelope for it
mel.env <- envelope(mel.mod)$point
plot(seq(0.25,10,by=0.25),mel.env[1,], ylim=c(-2,2),
     xlab="thickness (mm)", ylab="linear predictor",type="n")
lines(seq(0.25,10,by=0.25),mel.env[1,],lty=2)
lines(seq(0.25,10,by=0.25),mel.env[2,],lty=2)
lines(seq(0.25,10,by=0.25),mel.mod$t0,lty=1)