Survival Data Analysis

library(survival)
library(ISwR)
data(melanom)
str(melanom)

### The Kaplan-Meier estimates
# Surv creats a survival object
mfit <- survfit(Surv(days, status==1), data=melanom)
summary(mfit)
plot(mfit)

## By sex
mfit.bysex <- survfit(Surv(days, status==1)~sex, data=melanom)
plot(mfit.bysex, col=c("black", "grey"), legend.text=c("Female", "Male"))

# log-transformation
plot(mfit.bysex, fun="log", col=c("black", "grey"), legend.text=c("Female", "Male"))

# cloglog-transformation
plot(mfit.bysex, fun="cloglog", col=c("black", "grey"), legend.text=c("Female", "Male"))

### The log-rank test
survdiff(Surv(days, status==1)~sex, data=melanom)

### Cox Regression Model
data(ovarian)
ovcph <- coxph(Surv(futime, fustat) ~ age + factor(resid.ds) + factor(rx), x=T, data=ovarian)
summary(ovcph)

# Kaplan-Meier
ovarian.surv <- survfit(Surv(futime, fustat) ~ rx, data=ovarian)
plot(ovarian.surv, lty=2:3)
legend(600, .3, c("rx=1", "rx=2"), lty=2:3)

# Examine outliers and check proportional hazards assumption
temp <- cox.zph(ovcph)
print(temp)
plot(temp[1]) # Schoenfeld plot for age
split.screen(c(2,2))
plot(temp[3])
screen(2)
plot(temp[2])
screen(3)
plot(temp[1])
screen(4)
plot(ovarian$futime, ovcph$residuals, type="n", xlab="Time", ylab="Mart Res")
text(ovarian$futime, ovcph$residuals, labels=ovarian$rx, cox=.6)