library(survival)
#**********************************************************
# Data excerpted from NIH Multi Center AIDS Cohort Study. *
# Cohort consists of gay and bisexual men.  Variables are *
# happy: self-reported level of happiness at enrollment   *
# time: time to death or loss to followup,                *
# entry: age at first follow-up, status: 1 for death, 0 ow*
# life: finer self-reported happiness measure             *
# aaids: age at development of aids.                      *
#**********************************************************
mcac<-as.data.frame(scan("mcac.dat",
   what=list(borny=0,life=0,happy=0,status=0,entry=0,time=0,aaids=0)))
# 202 of the individuals in this study are lost to followup immediately
# after enrollment.  Delete them.
mcac<-mcac[mcac$entry<mcac$time,]
# Put the correct and incorrect survival fits on the same
# plot.  Both require a plot area of the same size, so it
# doesn't matter which we do first.  The calls to plot and
# lines below actually call plot.survfit and lines.survfit,
# and in plot.survfit, conf.int and mark.time defaults
# are true.  They are false in lines.survfit.  To get both
# lines in the same format, turn these features off in plot.
plot(survfit(Surv(time,status),data=mcac),conf.int=FALSE,mark.time=FALSE)
# Denote delayed entry by a first argument to Surv.
lines(survfit(Surv(entry,time,status),data=mcac))
coxph(Surv(entry,time,status)~happy,data=mcac)
#**********************************************************
# Data from http://lib.stat.cmu.edu/datasets/csb, Case    *
# Studies in Biometry, ch. 14, Interpretation of a        *
# Leukemia Trial Stopped Early, Emerson, et al.           *
# Variables are age, treatment arm (0 for D or 1 for I)   *
# indicator for remission and time of remission,          *
# indicator of transplant and time of transplant, time to *
# end of study and indicator of death, fab (a measure of  *
# disease status), gender(1=M,0=F), white blood count,    *
# kscore (a measure of disease status), and number of     *
# chemotherapy rounds before remission (-1 => missing)    *
#**********************************************************
emerson<-as.data.frame(scan("emerson.dat",what=list(
   age=0,ttt=0,crind=0,crtime=0,transs=0,transp=0,time=0,
   status=0,fab=0,gender=0,wbc=0,kscore=0,chem=0)))
emerson$fab[emerson$fab==-1]<-NA
emerson$wbc[emerson$wbc==-1]<-NA
emerson$wbcl<-1*(emerson$wbc<13.5)
cat('Fixed-covariate Cox regression\n')
coxph(Surv(time,status)~ttt,data=emerson)
cat('Fixed-covariate Cox regression\nIncorrect Use of Transpant Status\n')
coxph(Surv(time,status)~ttt+transs,data=emerson)
emerson2<-emerson1<-emerson
emerson1$entry<-rep(0,length(emerson$time))
emerson2$entry<-emerson$transp
emerson1$time<-emerson$transp
# Non-transplant status is 0 if patient got transplant or
# final status was 0, and 1 otherwise.
emerson1$status<-(1-emerson$transs)*emerson$status
emerson1$transi<-rep(0,length(emerson$time))
emerson2$transi<-rep(1,length(emerson$time))
emersonbig<-rbind(emerson1,emerson2[emerson$transs==1,])
coxph(Surv(entry,time,status)~ttt+transi,data=emersonbig,method='breslow')
#
# Remission indicator is treated the same way as transplant, and I won't
# repeat it.  I'll present the model with remission indicator and time
# since remission.
# 
emerson$id<-seq(length(emerson$time))
alltimes<-sort(unique(emerson$time))
start<-0
for(tt in alltimes){
   newemerson<-emerson[emerson$time>=tt,]
   newemerson$entry<-rep(start,length(newemerson$time))
   newemerson$cri<-1*(newemerson$crtime<tt)
   newemerson$crt<-newemerson$cri*(tt-newemerson$crtime)
   newemerson$status<-newemerson$status*(newemerson$time==tt)
   newemerson$time<-rep(tt,length(newemerson$time))
   if(start==0){
      verybigemerson<-newemerson
   }else{
      verybigemerson<-rbind(verybigemerson,newemerson)
   }
   start<-tt
}
coxph(Surv(entry,time,status)~ttt+cri+crt,data=verybigemerson,method='breslow')