# R script for lecture 8 in STK4900/9900

# Slide 13-14: Kaplan-Meier estimator ------------------------------------------

library(survival)

time <- c(2,6,12,54,56,68,89,96,96,125,128,131,140,141,143,145,146,148,162,168,173,181)
cens <- c(1,1,1,1,0,1,1,1,1,0,0,0,0,0,1,0,1,0,0,1,0,0)

fit.surv <- survfit(Surv(time,cens)~1)
summary(fit.surv)

plot(fit.surv, conf.int=FALSE, mark.time=TRUE)



# Slide 16: Kaplan-Meier estimator ---------------------------------------------

fit.surv <- survfit(Surv(time,cens)~1, conf.type="plain")

plot(fit.surv, conf.int=TRUE, mark.time=TRUE)


# Slide 18: Kaplan-Meier estimator ---------------------------------------------

print(fit.surv)



# Slide 20-21: Comparing two groups --------------------------------------------

time <- c(2,3,4,7,10,22,28,29,32,37,40,41,54,61,63,71,127,140,146,158,167,182,2,
       6,12,54,56,68,89,96,96,125,128,131,140,141,143,145,146,148,162,168,173,181)
cens <- c(1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,0,0,0,0,0,0,
       1,1,1,1,0,1,1,1,1,0,0,0,0,0,1,0,1,0,0,1,0,0)
group <- c(rep(1,22),rep(2,22))

fit.surv.both <- survfit(Surv(time,cens)~group, conf.type="plain")

plot(fit.surv.both, lty=1:2, xlab="Time (months)", ylab="Survival", mark.time=TRUE)
legend(5,0.2,c("control","treatment"),lty=1:2)

print(fit.surv.both)



# Slide 25: Comparing two groups (cont) -------------------------------------

survdiff(Surv(time,cens)~group)



# Slide 33-35: Cox regression --------------------------------------------

melanom <- read.table("http://www.uio.no/studier/emner/matnat/math/STK4900/data/melanoma.dat", header=TRUE)
fit.cox <- coxph(Surv(lifetime,status==1)~factor(sex)+thickn, data=melanom)
summary(fit.cox)

fit.cox.sex <- coxph(Surv(lifetime,status==1)~factor(sex), data=melanom)

anova(fit.cox.sex,fit.cox,test="Chisq")