# ==========================================================  
# Table 5.4 Modulated Markov models for pulmonary 
#           exacerbations (p.192)
# ==========================================================  

library(survival)


rhDNase <- read.table("rhDNase.dat", header=F)
dimnames(rhDNase)[[2]] <- c("id","trt","fev","fev2","time1","time2","status",
                            "etype","enum","enum1","enum2")
rhDNase <- data.frame(rhDNase[rhDNase$etype == 1,], row.names=NULL)
rhDNase$fevc <- rhDNase$fev - mean(rhDNase$fev[rhDNase$enum1==1])

rhDNase[1:10,]


# Model 1
fit1 <- coxph(Surv(time1,time2,status) ~ trt + fevc,
              data=rhDNase, method="breslow")
summary(fit1)


# Model 2
fit2 <- coxph(Surv(time1,time2,status) ~ trt + fevc + frailty(id),
              data=rhDNase, method="breslow")
summary(fit2)
fit2$history


# Function to create data frame with time-dependent covariate
# Description of variables in data frame, rhDNase80 :
#   id      patient ID
#   time1   start time
#   time2   stop time
#   status  1 if transition time; 0 if censoring time
#   enum    cumulative number of lines for each id
#   trt     1 if rhDNase; 0 if placebo
#   fevc    FEV (centered)
#   It      0 if time > 80; 1 if time <= 80
#   Nt      1 if N(t-) >= 1; 0 if N(t-) = 0

gettddata.f <- function(indata) {
  indata <- indata[,c("id","time1","time2","status","enum1","trt","fevc")]

  pid <- sort(unique(indata$id))

  data <- NULL
  for (i in 1:length(pid)) {
    tmp <- indata[indata$id == pid[i],]

    dataj <- NULL
    for (k in 1:length(tmp$id)) {
      if ( (tmp$time1[k] < 80) & (tmp$time2[k] >= 80) ) {
        if (tmp$time2[k] == 80) {
          time1  <- tmp$time1[k]
          time2  <- tmp$time2[k]
          status <- tmp$status[k]
          It <- 1
        }
        else {  
          time1  <- c(tmp$time1[k], 80)
          time2  <- c(80, tmp$time2[k])
          status <- c(0, tmp$status[k])
          It     <- c(1,0)        
        }
      }
      else {
        time1  <- tmp$time1[k]
        time2  <- tmp$time2[k]
        status <- tmp$status[k]
        It     <- ifelse(tmp$time1[k] > 80, 0, 1)
      }
      nlen <- length(time1)
      id   <- rep(tmp$id[k], nlen)
      trt  <- rep(tmp$trt[k], nlen)
      fevc <- rep(tmp$fevc[k], nlen)

      Nt <- ifelse(tmp$enum1[k] == 1, 0, 1) 
      Nt <- rep(Nt, nlen)

      dataj <- rbind(dataj, data.frame(id,time1,time2,status,trt,fevc,It,Nt))
    }
    dataj$enum <- 1:length(dataj$id)

    data <- rbind(data, dataj[,c("id","time1","time2","status","enum","trt","fevc","It","Nt")])
  }
  return(data)
}

rhDNase80 <- gettddata.f(indata=rhDNase)
rhDNase80$trtIt0 <- rhDNase80$trt*(rhDNase80$It == 1)
rhDNase80$trtIt1 <- rhDNase80$trt*(rhDNase80$It == 0)  
rhDNase80[1:18,]


# Model 3
fit3 <- coxph(Surv(time1,time2,status) ~ trtIt0 + trtIt1 + fevc + frailty(id),
              data=rhDNase80, method="breslow")
summary(fit3)


# Model 4
fit4 <- coxph(Surv(time1,time2,status) ~ trtIt0 + trtIt1 + fevc + Nt,                
              data=rhDNase80, method="breslow")
summary(fit4)