# ==========================================================  
# Table 3.1 Estimates of treatment coefficient from several    
#           multiplicative rate function models for rat        
#           tumorgenicity data. (p.105)      
# ==========================================================  


rats <- read.table("c:/WORKSHOP/rats.dat", header=F)
dimnames(rats)[[2]] <- c("id","start","stop","status","enum","trt")
rats$rtrunc  <- rep(NA, nrow(rats))
rats$tstatus <- ifelse(rats$start == 0, 1, 2)

rats[1:12, c("id","start","stop","status","rtrunc","tstatus","enum","trt")]


# Model: Poisson, Baseline Rate: Weibull

wfit <- censorReg(censor(stop, status) ~ trt, data=rats,
                  truncation=censor(start,rtrunc,tstatus),
                  distribution="weibull")
summary(wfit)


# Model: Poisson, Baseline Rate: Piecewise

gd.pw.f <- function(indata) { 
  pid <- sort(unique(indata$id))
   
  data <- matrix(0, nrow=(length(pid)*4), ncol=5)
  for (i in 1:length(pid)) {  
    tmp <- indata[indata$id == pid[i], ]
    etime <- floor( tmp$stop[tmp$status == 1] )
  
    startpos <- 4*(i-1) + 1   
    stoppos  <- 4*i
    data[startpos:stoppos,1] <- rep(pid[i],4)
    data[startpos:stoppos,2] <- c(1,2,3,4)
    data[startpos:stoppos,3] <- c(sum((etime >  0) & (etime <= 30)),
                                  sum((etime > 30) & (etime <= 60)),
                                  sum((etime > 60) & (etime <= 90)),
                                  sum((etime > 90) & (etime <= 122)))
    data[startpos:stoppos,4] <- c(30,30,30,32)
    data[startpos:stoppos,5] <- rep(unique(tmp$trt),4)
  }
  data <- data.frame(data)
  dimnames(data)[[2]] <- c("id","interval","count","len","trt")
  return(data)
}

rats.pois <- gd.pw.f(indata=rats)

pfit <- glm(count ~ offset(log(len)) + factor(interval) + trt,
            family=poisson(link=log), data=rats.pois)
summary(pfit, corr=F)


# Model: Poisson, Baseline Rate: Nonparametric

fit <- coxph(Surv(start,stop,status) ~ trt,
             data=rats, method="breslow")
summary(fit)

cox.zph(fit, transform="identity")
cox.zph(fit, transform="log")


# Model: Mixed Poisson, Baseline Rate: Weibull

library(MASS)
n <- tapply(rats$status, list(rats$id), sum)
n <- data.frame(id=names(n), n=n)
tau <- tapply(rats$stop, list(rats$id), max)
tau <- data.frame(id=names(tau), tau=floor(tau))
rats0 <- rats[rats$enum == 1, c("id","trt")]
rats0 <- merge(rats0, tau, by="id", all.x=T)
rats0 <- merge(rats0, n, by="id", all.x=T)

fit <- glm.nb(n ~ offset(log(tau)) + trt, data=rats0, link="log")
summary(fit)


# Model: Mixed Poisson, Baseline Rate: Nonparametric

fit <- coxph(Surv(start,stop,status) ~ trt + frailty(id), 
             data=rats, method="breslow")
summary(fit)


# Model: Robust, Baseline Rate: Nonparametric

fit <- coxph(Surv(start,stop,status) ~ trt + cluster(id),
             data=rats, method="breslow")
summary(fit)