#Chapter 5.

library(survival)

#Table 5.1, p. 166. 

#Creating the dummy variables for hercoc.

print(table(uis$hercoc))

uis$hercoc1 <- 1*(uis$hercoc==1)
uis$hercoc2 <- 1*(uis$hercoc==2)
uis$hercoc3 <- 1*(uis$hercoc==3)
uis$hercoc4 <- 1*(uis$hercoc==4)


# The Log-rank test and Likelihood ratio test for the main effect of hercoc. 
# (Should be done for all the variables with
# more than two categories (ie ivhx, agecat, beckcat and drugcat.)

hercoc.ph<-coxph(formula = Surv(time, censor) ~ hercoc2+hercoc3+hercoc4, data = uis,
  method="breslow", na.action=na.exclude)
print(summary(hercoc.ph))

# The medians by hercoc group (not exactly the same as with STATA).

hercoc.surv <-survfit(formula = Surv(time, censor) ~ hercoc, data = uis, 
  na.action=na.exclude)
print(hercoc.surv)

uis$ivhx1 <- 1*(uis$ivhx==1)
uis$ivhx2 <- 1*(uis$ivhx==2)
uis$ivhx3 <- 1*(uis$ivhx==3)

ivhx.ph<-coxph(formula = Surv(time, censor) ~ ivhx2+ivhx3, data = uis,
  method="breslow", na.action=na.exclude)
print(summary(ivhx.ph))

ivhx.surv<-survfit(formula = Surv(time, censor) ~ ivhx, data = uis, 
  na.action=na.exclude)
print(ivhx.surv)

race.ph<-coxph(formula = Surv(time, censor) ~ race, data = uis,
  method="breslow", na.action=na.exclude)
print(summary(race.ph))

race.surv<-survfit(formula = Surv(time, censor) ~ race, data = uis, 
  na.action=na.exclude)
print(race.surv)

treat.ph<-coxph(formula = Surv(time, censor) ~ treat, data = uis,
  method="breslow", na.action=na.exclude)
print(summary(treat.ph))

treat.surv<-survfit(formula = Surv(time, censor) ~ treat, data = uis, 
  na.action=na.exclude)
print(treat.surv)

site.ph<-coxph(formula = Surv(time, censor) ~ site, data = uis,
  method="breslow", na.action=na.exclude)
print(summary(site.ph))

site.surv<-survfit(formula = Surv(time, censor) ~ site, data = uis, 
  na.action=na.exclude)
print(site.surv)

#Creating categorical variable for age: agecat.

uis$agecat[uis$age < 28] <- 1
uis$agecat[uis$age >= 28 & uis$age < 33] <- 2
uis$agecat[uis$age >= 33 & uis$age < 38] <- 3
uis$agecat[uis$age >= 38] <- 4
print(table(uis$agecat))

# Creating dummy variables for agecat.

uis$agecat1 <- 1*(uis$age < 28)
uis$agecat2 <- 1*(uis$age >= 28 & uis$age < 33)
uis$agecat3 <- 1*(uis$age >= 33 & uis$age < 38)
uis$agecat4 <- 1*(uis$age >= 38)


agecat.ph<-coxph(formula = Surv(time, censor) ~ agecat2+agecat3+agecat4, data = uis, method="breslow",
  na.action=na.exclude)
print(summary(agecat.ph))

agecat.surv<-survfit(formula = Surv(time, censor) ~ agecat, data = uis, 
  na.action=na.exclude)
print(agecat.surv)

#Creating categorical variable for becktota: beckcat.

uis$beckcat[uis$becktota < 10] <- 1
uis$beckcat[uis$becktota >= 10 & uis$becktota < 15] <- 2
uis$beckcat[uis$becktota >= 15 & uis$becktota < 25] <- 3
uis$beckcat[uis$becktota >= 25] <- 4
print(table(uis$beckcat))


#Creating the dummy variables for beckcat.

uis$beckcat1 <- 1*(uis$becktota < 10)
uis$beckcat2 <- 1*(uis$becktota >= 10 & uis$becktota < 15)
uis$beckcat3 <- 1*(uis$becktota >= 15 & uis$becktota < 25)
uis$beckcat4 <- 1*(uis$becktota >= 25)

beckcat.ph<-coxph(formula = Surv(time, censor) ~ beckcat2+beckcat3+beckcat4, 
  data = uis, method="breslow", na.action=na.exclude)
print(summary(beckcat.ph))

beckcat.surv<-survfit(formula = Surv(time, censor) ~ beckcat, data = uis, 
  na.action=na.exclude)
print(beckcat.surv)

#Creating the categorical variable drugcat.

uis$drugcat[uis$ndrugtx < 2] <- 1
uis$drugcat[uis$ndrugtx >= 2 & uis$ndrugtx < 4] <- 2
uis$drugcat[uis$ndrugtx >= 4 & uis$ndrugtx < 7] <- 3
uis$drugcat[uis$ndrugtx >= 7] <- 4
print(table(uis$drugcat))

#Creating the dummy variables for drugcat.

uis$drugcat1 <- 1*(uis$ndrugtx < 2)
uis$drugcat2 <- 1*(uis$ndrugtx >= 2 & uis$ndrugtx < 4)
uis$drugcat3 <- 1*(uis$ndrugtx >= 4 & uis$ndrugtx < 7)
uis$drugcat4 <- 1*(uis$ndrugtx >= 7)

drugcat.ph<-coxph(formula = Surv(time, censor) ~ drugcat2+drugcat3+drugcat4, 
  data = uis, method="breslow", na.action=na.exclude)
print(summary(drugcat.ph))

drugcat.surv<-survfit(formula = Surv(time, censor) ~ drugcat, data = uis, 
  na.action=na.exclude)
print(drugcat.surv)

# Table 5.2, p. 167.

uis$age5 <- uis$age/5
age5.ph<-coxph(formula = Surv(time, censor) ~ age5, data = uis, method="breslow",
  na.action=na.exclude)
print(summary(age5.ph))

uis$beck10 <- uis$becktota/10
beck10.ph<-coxph(formula = Surv(time, censor) ~ beck10, data = uis, method="breslow",
  na.action=na.exclude)
print(summary(beck10.ph))

uis$drug5 <- uis$ndrugtx/5
drug5.ph<-coxph(formula = Surv(time, censor) ~ drug5, data = uis, method="breslow",
  na.action=na.exclude)
print(summary(drug5.ph))

# Table 5.3, p. 168.

full.ph <- coxph(formula = Surv(time, censor) ~ age + becktota + ndrugtx + hercoc2 +
  hercoc3 + hercoc4 + ivhx2 + ivhx3 + race + treat + site, data=uis, method="breslow", 
  na.action=na.exclude)
print(summary(full.ph))

# Table 5.4, p. 168.

reduced1.ph <- coxph(formula = Surv(time, censor) ~ age + becktota + ndrugtx + 
   ivhx2 + ivhx3 + race + treat + site, data=uis, method="breslow", 
  na.action=na.exclude)
print(summary(reduced1.ph))

# Table 5.5, p. 169.

reduced2.ph <- coxph(formula = Surv(time, censor) ~ age + becktota + ndrugtx + 
   ivhx3 + race + treat + site, data=uis, method="breslow", 
  na.action=na.exclude)
print(summary(reduced2.ph) )

# Table 5.6, p. 170.

cat.ph <- coxph(formula=Surv(time, censor) ~ agecat2+agecat3+agecat4+beckcat2+beckcat3+
  beckcat4+drugcat2+drugcat3+drugcat4+ ivhx3+site+race+treat, method="breslow", data=uis, 
  na.action=na.exclude)
print(summary(cat.ph))

# Fig. 5.2a, p. 175.  
# The fitted line appears to be approximately linear thus there is no transformation of age 
# needed.

age.ph <- coxph(formula = Surv(time, censor) ~ becktota + ndrugtx + 
   ivhx3 + race + treat + site, data=uis, method="breslow", 
  na.action=na.exclude)
uis$resid <- residuals.coxph(age.ph, type="martingale")

windows(6,6)
plot(uis$age, uis$resid, xlab="Age", ylab="Martingale Residuals", 
     ylim=c(-4, 1.0))
lines(lowess(uis$age[ !is.na(uis$age) & !is.na(uis$resid)],
       uis$resid[ !is.na(uis$age) & !is.na(uis$resid)]))

# Fig. 5.3a, p. 175.
# The fitted line appears to be approximately linear thus there is no transformation of becktota 
# needed.

becktota.ph <- coxph(formula = Surv(time, censor) ~ age + ndrugtx + 
   ivhx3 + race + treat + site, data=uis, method="breslow", 
  na.action=na.exclude)
uis$resid <- residuals.coxph(becktota.ph, type="martingale")

windows(6, 6)
plot(uis$becktota, uis$resid, xlab="Becktota",ylab="Martingale Residuals", ylim=c(-4, 1.0))
lines(lowess(uis$becktota[ !is.na(uis$becktota) & !is.na(uis$resid)], 
  uis$resid[ !is.na(uis$becktota) & !is.na(uis$resid)]))

# Fig. 5.4a, p. 175.
# The fitted line has a noticable squiggle in the beginning which makes us very nervous and 
# it is an indicator that we can't just include ndrugtx in the model; we should transform it
# before including it in the final model.

ndrugtx.ph <- coxph(formula = Surv(time, censor) ~ age + becktota + 
   ivhx3 + race + treat + site, data=uis, method="breslow", 
  na.action=na.exclude)
uis$resid <- residuals.coxph(ndrugtx.ph, type="martingale")

windows(6,6)
plot(uis$ndrugtx, uis$resid, xlab="Ndrugtx",ylab="Martingale Residuals", ylim=c(-4, 1.0))
lines(lowess(uis$ndrugtx[ !is.na(uis$ndrugtx) & !is.na(uis$resid)], 
  uis$resid[ !is.na(uis$ndrugtx) & !is.na(uis$resid)]))

# Creating ndrugfp1 and ndrugfp2, p. 174. 

uis$ndrugfp1 <- 1/((uis$ndrugtx+1)/10)
uis$ndrugfp2 <- (1/((uis$ndrugtx+1)/10))*log((uis$ndrugtx+1)/10)

# Table 5.8, p. 174. 

reduced3.ph <- coxph(formula=Surv(time, censor)~ age + becktota + ndrugfp1 + ndrugfp2 +
   ivhx3 + race + treat + site, data=uis, method="breslow", na.action=na.exclude)
print(summary(reduced3.ph))

# Creating the interaction variables, p. 178. 

uis$agesite <- uis$age*uis$site
uis$racesite <- uis$race*uis$site
uis$agefp1 <- uis$age*uis$ndrugfp1
uis$agefp2 <- uis$age*uis$ndrugfp2

# Table 5.10, p. 179.

interaction.ph <- coxph(formula=Surv(time, censor) ~ age + becktota + + ndrugfp1 + ndrugfp2 + 
   ivhx3 + race + treat + site + age:site + race:site + age:ndrugfp1 + age:ndrugfp2,
 data=uis, method="breslow", na.action=na.exclude)
print(summary(interaction.ph))

# Table 5.11, p. 179.

reducedinter.ph <- coxph(formula=Surv(time, censor) ~ age + becktota + + ndrugfp1 + ndrugfp2 +
 ivhx3 + race + treat + site + agesite + racesite,
 data=uis, method="breslow", na.action=na.exclude)
print(summary(reducedinter.ph))