#install and load the survival package
library(survival)

#import data
x=read.table("R://Teaching/Training materials/R course/survival analysis/PBC.dat", header=TRUE)
#show first 10 observations
x[1:10,]

#restrict to first 312 patients
x=x[1:312,]

#(a)
x$delta=ifelse ((x$status==0|x$status==1), 0, 1) 
  
#(b)
x$trt=ifelse (x$drug==2, 0, 1) 

#(c)
x$age_yr=x$age/365

#(d)
x$logbili=log(x$bili)

attach(x)

#(e) compute Kaplan-Meier estimate and confidence intervals
surv=survfit(Surv(futime, delta) ~ trt, conf.type = "log-log")
summary(surv)
#Kaplan-Meier curves 
plot(surv, lty = 1:2, xlab = "Time (days)", ylab = "Survival probability")
#plot(surv, conf.int = T, lty = 1:2, xlab = "Time of remission (weeks)", ylab = "Survival probability")
legend('topright', c("placebo", "drug"), lty = 1:2, lwd = 2)

#(f) median survival time and confidence interval
survfit(Surv(futime, delta) ~ trt, conf.type = "log-log")

#(g) log-rank test for differnces between the two groups
survdiff(Surv(futime, delta) ~ trt)


#(h) stratified log-rank test
survdiff(Surv(futime, delta) ~ trt + strata(stage))

#(i) fit Cox model
fit=coxph(Surv(futime, delta) ~ logbili + albumin + age_yr + edema)
#obtain CI of HRs
#exp(confint(fit))


#(j) deviance residual plot to detect any outliers
dev.resid <- residuals(fit, type="deviance")
plot(x$age_yr, dev.resid, xlab="Age",ylab="Deviance Residuals")


#(k) statistical test of PH assumption
cox.zph(fit)