# ============================================================
# A.4.1 Timeline Plots
# R code to create Figure A.2
#
# December 11, 2018
# =============================================================
# Note: rhDNase.dat dataset is available on the book website
# =============================================================


# -------------------------------------------------------------
# Create an input dataframe for event.chart function
# -------------------------------------------------------------

createdata.f <- function(indata) {
  pid <- sort(unique(indata$id))
  nlen <- length(pid)

  etype1 <- etype2 <- matrix(NA, nrow=nlen, ncol=5)
  dstart <- rep(NA, nlen)
  dstop  <- rep(NA, nlen)
  etime  <- matrix(NA, nrow=nlen, ncol=20)
  trt    <- rep(NA, nlen)
  for (i in 1:nlen) {
    datai <- indata[indata$id == pid[i],]

    trt[i]    <- datai$trt[1] 
    dstart[i] <- min(datai$time1)
    dstop[i]  <- max(datai$time2)

    tt <- NULL
    for (j in 1:length(datai$id)) {
      tt <- c(tt, c(datai$time1[j], datai$time2[j]))
    }

    etime[i,1:length(tt)] <- tt

    if ( sum((datai$etype == 1) & (datai$status == 1)) > 0 ) {
       times <- datai$time2[(datai$etype == 1) & (datai$status == 1)]
       etype1[i,1:length(times)] <- times
    }

    if ( sum((datai$etype == 2) & (datai$status == 1)) > 0 ) {
      times <- datai$time2[(datai$etype == 2) & (datai$status == 1)]
      etype2[i,1:length(times)] <- times
    }
  }

  outdata <- data.frame(id=pid, trt, etype1, etype2, dstart, dstop, etime)
  dimnames(outdata)[[2]] <- c("id", "trt", paste("e",1:5,sep=""), paste("ef",1:5,sep=""),
                              "dstart","dstop", paste("etime",1:20,sep=""))
  return(outdata)
}


rhDNase <- read.table("rhDNase.dat", header=T)


copdEC <- createdata.f(indata=rhDNase)
copdEC <- copdEC[copdEC$id %in% c(493313, 589303, 601320, 766303,
                                  493303, 589310, 765303, 768309),]
copdEC <- copdEC[,c("id",paste("e",1:3,sep=""),paste("ef",1:3,sep=""),
                    "dstart","dstop",paste("etime",1:12,sep=""))]
copdEC <- copdEC[order(copdEC$id),]
copdEC$id <- c(1:nrow(copdEC))


# -------------------------------------------------------------
# Use event.chart function to create Figure A.2
# -------------------------------------------------------------

library(Hmisc)

pdf("figureA2.pdf", width=10, height=6)
par(mai=c(0.8,0.8,0.4,0.4))

event.chart(copdEC, subset.c=c(paste("e",1:3,sep=""),paste("ef",1:3,sep=""),
                               "dstart","dstop",paste("etime",1:12,sep="")),
            sort.by="id", 
            y.idlabels=1, x.julian=TRUE, y.lim.extend =c(0,0.5),
            point.cex=c(rep(1.2,6), 1.5, 1,5, rep(0.1,12)),
            point.pch=c(rep(19,3), rep(1,3), 3, 3, rep(NA,12)),
            x.lab="DAYS  SINCE  RANDOMIZATION", y.lab="PATIENT ID",
            titl ="EVENT  CHART  FOR  8  SELECTED  INDIVIDUALS",
            line.lty=2, line.lwd=1, 
            cex=1, legend.plot=TRUE, legend.location="i",
            line.add = matrix(paste("etime",1:12,sep=""),2,6),
            line.add.lty=rep(c(1,2),3),
            line.add.col=rep("black",6), 
            line.add.lwd=rep(c(3,1),3),
            legend.point.pch=c(19,1),
            legend.point.text=c("EXACERBATION","EXACERBATION FREE"),
            legend.point.at=list(x=c(0,50), y=c(9,9.5)),
            legend.cex=1, legend.bty="n")

dev.off()