# Read in the dataset
library(rms)
cancer <- stata.get("http://biostat.mc.vanderbilt.edu/wiki/pub/Main/CourseBios312/prostaticcancer.dta")

# Kaplan-Meier Plot
plot(survfit(Surv(time, status) ~ treatment, data=cancer), lty=2:3, xlab="Time", ylab="Survival") 
legend("bottomleft", c("Treatment = 1", "Treatment = 2"), lty=2:3, inset=0.05)


# List the observations
cancer[cancer$treatment==2,]
15/16

# Log rank test
survdiff(Surv(time, status) ~ treatment, data=cancer)

# Cox proportional hazard regression
coxph(Surv(time, status) ~ treatment, data=cancer, robust=FALSE)
coxph(Surv(time, status) ~ treatment, data=cancer, robust=TRUE)


# Index
coxph(Surv(time, status) ~ index, data=cancer, robust=FALSE)
coxph(Surv(time, status) ~ index, data=cancer, robust=TRUE)