#Chapter 4

library(survival)

#Table 4.2, p. 119.
#Note1: The coxph uses Efron as default,
#furthermore, the coxph functions only supply the CI for the hazard ratio.

drug.ph <- coxph(formula=Surv(time, censor)~drug, data=hmohiv, method="breslow")
print(summary(drug.ph))

#Creating the dummy variables for categories of age.

hmohiv$age2[hmohiv$age<30 | hmohiv$age>34] <- 0
hmohiv$age2[hmohiv$age>=30 & hmohiv$age<=34] <- 1

hmohiv$age3[hmohiv$age<35 | hmohiv$age>39] <- 0
hmohiv$age3[hmohiv$age>=35 & hmohiv$age<=39] <- 1

hmohiv$age4[hmohiv$age<40] <- 0
hmohiv$age4[hmohiv$age>=40] <- 1

# Table 4.4, p. 122 and table 4.5, p. 123.

dummy.ph <- coxph(formula=Surv(time, censor)~age2+age3+age4, data=hmohiv, method="breslow")
print(summary(dummy.ph))

# Creating deviation from mean coded variables for the categories of age.

hmohiv$age1[ hmohiv$age >= 30] <- 0
hmohiv$age1[ hmohiv$age < 30] <- 1

hmohiv$age2d <- hmohiv$age2
hmohiv$age2d[hmohiv$age1 == 1] <- -1

hmohiv$age3d <- hmohiv$age3
hmohiv$age3d[hmohiv$age1 == 1] <- -1

hmohiv$age4d <- hmohiv$age4
hmohiv$age4d[hmohiv$age1 == 1] <- -1

# Table 4.7, p. 127.

dummyd.ph <- coxph(formula=Surv(time, censor)~age2d+age3d+age4d, data=hmohiv, method="breslow")
print(summary(dummyd.ph))

# Table 4.8, p. 129.

age.ph <- coxph(formula=Surv(time, censor)~age, data=hmohiv, method="breslow")
print(summary(age.ph))

# Creating the variable drug in the uis data set.

uis$drug[uis$ivhx==1] <- 0
uis$drug[uis$ivhx>1] <- 1


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

# Table 4.9, p. 133.
#Note: The estimates for the crude model is slightly different since it only excluded observations based on 
#whether or not they were missing in the variable drug.  So, there are only 18 observations that have been 
#excluded from calculations whereas in the adjusted model there were 23 observations excluded from the calculations.

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

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

# Table 4.10, p. 135.
#Note: again the crude model is slightly different since it only excluded observations that had missing entries
#for treat (and treat did not have any missing values).  The two other models excluded observations that were 
#missing for either age or treat and therefore 5 observations were excluded from the calculations.

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

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

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

# Fig. 4.2, p. 139. 
# Note: First we fit the cox ph model using drug as a predictor.  Then we use the survfit function to produce 
# the baseline survival function.  Using the baseline survival we then can calculate the survival curve for the 
# group where drug=1.  

crude.ph <- coxph(formula=Surv(time, censor)~drug, data=hmohiv, method="breslow", na.action=na.exclude)
fit <- survfit(crude.ph, conf.type="none")
fit$surv2 <- fit$surv^exp(.779)

windows(6,6)
plot( fit$time, fit$surv, xlab="Survival Time (Months)", ylab="Survival Probability", 
     ylim=c(0, 1.0), type="s")
points(fit$time, fit$surv2,  type="s")

#Table 4.12, p. 148. 
hmohiv$agec <- hmohiv$age - 35
full.model <- coxph(formula=Surv(time, censor)~agec+drug, data=hmohiv, 
              na.action=na.exclude, method="breslow")
print(full.model)

#Creating centered age and ndrugtx variables
uis$agec <- uis$age - 30
uis$ndrugtxc <- uis$ndrugtx - 3

#
full.model <- coxph(formula=Surv(time, censor)~treat+agec+drug+ndrugtxc, data=uis, 
              na.action=na.exclude, method="breslow")
print(full.model)
fit <- survfit(full.model, conf.type="none")
fit$surv2 <- fit$surv^exp(-0.2271)

windows(6,6)
plot( fit$time, fit$surv, xlab="Time to Drug Use From Admission (Days)", ylab="Survival Probability", 
     ylim=c(0, 1.0), xlim=c(0, 1200), type="s")
points(fit$time, fit$surv2,  type="s")