#load the survival library

library(survival)


# Table 2.2, p. 32

mini.surv <- survfit(Surv(time, censor)~1, mini, conf.type="none")
print(summary(mini.surv))

# Fig. 2.1, p. 32

windows(6,6)
plot(mini.surv, lty = 1, cex=.7, cex.lab=.7, cex.axis=.7 )

# Table 2.3, p. 35
hmohiv.surv <- survfit(Surv(time, censor)~ 1, hmohiv, conf.type="none")
print(summary(hmohiv.surv))

# Fig. 2.2, p. 34
windows(6,6)
plot(hmohiv.surv, lty=1, xlab="Time", ylab="Survival Probability", cex=.7, cex.lab=.7, cex.axis=.7 )
 
# Create timecat, the categorical variable of time in the hmohiv dataset.
# Note: In order to include the lower end point in the interval you need to specify a cut off
# point that is slightly lower.  Also the lower endpoint of the whole data set needs to be 
# smaller than the smallest number in the whole dataset.  Likewise, the greatest number has 
# to be bigger than the largest number in the whole data set.  
# The object hmohiv$timecat1 is categorical with categories 1,2,...,11 but we would like the 
# categories to represent months so we must multiply them by 6 so now the categories are 6,12, ..., 66.

hmohiv$agecat <- cut(hmohiv$age,  breaks=c(19.9, 29, 34, 39, 54.1) )

hmohiv$timecat1 <- cut(hmohiv$time, 
               breaks = c(-.5,5.5,11.5,17.5,23.5,29.5,35.5,41.5,47.5,53.5,59.5,66.5), labels=F)
hmohiv$timecat <- 6*hmohiv$timecat1

# Table 2.4, p. 38

timecat.surv <- survfit( Surv(timecat, censor)~ 1, hmohiv, conf.type="none")
print(summary(timecat.surv))

# Fig. 2.3, p. 39
windows(6,6)
plot( timecat.surv, xlab="Survival Time (Months)", ylab="Proportion Surviving", cex=.7, cex.lab=.7, cex.axis=.7 ) 

# Fig. 2.5, p. 46
# Note: Could only graph the pointwise confidence bands using the log-log survivorship function.

time.conf <- survfit(Surv(time, censor)~ 1, hmohiv, conf.type="log-log")

windows(6,6)
plot(time.conf, lty=1:2, xlab="Survival Time (Months)", ylab="Survival Probability", cex=.7, cex.lab=.7, cex.axis=.7 )
legend(20, 1.0, c("Kaplan-Meier Curve", "Point-wise 95% Confidence Interval"), lty=1:2, cex=.7 ) 

# Table 2.6, p. 

timeconf.surv <- survfit( Surv(time, censor)~ 1, hmohiv, conf.type="log-log")
print(summary(timeconf.surv))

# The mean of survivorship function, p. 57

print(timeconf.surv)

# Fig. 2.7, p. 58

timestrata.surv <- survfit(Surv(time, censor)~ strata(drug), hmohiv, conf.type="log-log")
windows(6,6)
plot(timestrata.surv, lty=c(1,3), xlab="Time", ylab="Survival Probability", cex=.7, cex.lab=.7, cex.axis=.7)
legend(40, 1.0, c("Drug - No", "Drug - Yes") , lty=c(1,3), cex=.7 ) 

# Using data set minitest in table 2.8, p. 63.

#Table 2.10, testing for differences in survival between the two drug groups.
# Note: When rho=0 then the test is the log-rank test. However, when rho=1 then 
# the weight is similar
# to the Peto and Prentice test.  Rho=.5 appears to be close to the Tarone-Ware test.

print(survdiff(Surv(time, censor) ~ drug, data=minitest, rho=0))
print(survdiff(Surv(time, censor) ~ drug, data=minitest, rho=.5))
print(survdiff(Surv(time, censor) ~ drug, data=minitest, rho=1))

# Table 2.11, p. 65
# Testing for differences between drug group using the G-rho family of tests.  
# Note: When rho=0 then the test is the log-rank test. However, when rho=1 then the weight is similar
# to the Peto and Prentice test.  Rho=.5 appears to be close to the Tarone-Ware test.

print(survdiff(Surv(time, censor) ~ drug, data=hmohiv, rho=0))
print(survdiff(Surv(time, censor) ~ drug, data=hmohiv, rho=.5))
print(survdiff(Surv(time, censor) ~ drug, data=hmohiv, rho=1))

# Creating the categorical age variable, agecat.

hmohiv$agecat <- cut(hmohiv$age,  breaks=c(19.9, 29, 34, 39, 54.1) )

# Table 2.12, p. 65.

age.surv <- survfit(Surv(time, censor)~ strata(agecat), hmohiv, conf.type="log-log")
print(age.surv)

# Fig. 2.8, p. 69

windows(6,6)
plot(age.surv, lty=c(6, 1, 4, 3), xlab="Time", ylab="Survival Probability")
legend(40, 1.0, c("Group 1", "Group 2", "Group 3", "Group 4")  , lty=c(6, 1, 4, 3) ) 

# Table 2.14, p. 70

print(survdiff(Surv(time, censor) ~ agecat, data=hmohiv, rho=0))
print(survdiff(Surv(time, censor) ~ agecat, data=hmohiv, rho=.5))
print(survdiff(Surv(time, censor) ~ agecat, data=hmohiv, rho=1))