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

* Generate the relapse variable
psa$relapse <- as.numeric(psa$inrem=="no")
with(psa, table(relapse, inrem))

* Kaplan-Meier Plot
plot(survfit(Surv(obstime, relapse) ~ bss, data=psa), lty=2:4, xlab="Time", ylab="Survival") 
legend("topright", c("BSS = 1", "BSS = 2", "BSS = 3"), lty=2:4, inset=0.05)


* List the observations
psa[psa$bss==1,]

* Continuous bone scan score
m.cont <- coxph(Surv(obstime, relapse) ~ bss, data=psa, robust=TRUE)
m.cont
summary(m.cont)

* Indicator variable for BSS==3 and regression
psa$bss3 <- psa$bss==3
m.3 <- coxph(Surv(obstime, relapse) ~ bss3, data=psa, robust=TRUE)
summary(m.3)


* Indicator variable for BSS=2 and regression
psa$bss2 <- psa$bss==2
m.23 <- coxph(Surv(obstime, relapse) ~ bss2 + bss3, data=psa, robust=TRUE)
summary(m.23)

* A formal method for comparing two models
anova(m.3, m.cont)
anova(m.23, m.cont)