# =================================================================
# Supplementary Andersen-Gill Analyses
# =================================================================


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$fevc    <- rhDNase$fev - mean(rhDNase$fev[rhDNase$enum == 1])
rhDNase$gtime   <- rhDNase$time2 - rhDNase$time1
rhDNase$status1 <- ifelse(rhDNase$etype == 1, rhDNase$status, 0)

rhDNase[1:10, c("id","enum","etype","time1","time2","gtime",
                "status","status1","trt","fev","fevc")]


# Data includes patients in risk set when they have exacerbations

fit1a <- coxph(Surv(gtime, status1) ~ trt + fevc,
               data=rhDNase, method="breslow")
summary(fit1a)


fit1b <- coxph(Surv(gtime, status1) ~ trt + fevc + frailty(id),
               data=rhDNase, method="breslow")
summary(fit1b)


# Data excludes patients from risk set when they have exacerbations

fit2a <- coxph(Surv(gtime, status1) ~ trt + fevc,
               data=rhDNase, subset=(etype == 1), method="breslow")
summary(fit2a)


fit2b <- coxph(Surv(gtime, status1) ~ trt + fevc + frailty(id),
               data=rhDNase, subset=(etype == 1), method="breslow")
summary(fit2b)