#install and load the survival package
library(survival)


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


#install and load the survival package
library(survival)

#compute Kaplan-Meier estimate and confidence intervals
surv=survfit(Surv(time, delta) ~ group, conf.type = "log-log")
summary(surv)

#Kaplan-Meier curves 
plot(surv, lty = 1:2, xlab = "Time of remission (weeks)", ylab = "Survival probability")
#plot(surv, conf.int = T, lty = 1:2, xlab = "Time of remission (weeks)", ylab = "Survival probability")
legend('topright', c("control", "6-MP"), lty = 1:2, lwd = 2)

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

#log-rank test for differnces between the two groups
survdiff(Surv(time, delta) ~ group)



###############################################################

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

#stratified log-rank test
survdiff(Surv(time, status) ~ treatment + strata(age), data=x)


###############################################################

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

fit=coxph(Surv(time, status) ~ age + factor(sex) + bun + ca + hb + pcells + factor(protein), data = x)

#obtain CI of beta coefficients
confint(fit)

#obtain CI of HRs
exp(confint(fit))


##ties in Cox model
#sort by time
sorted=x[order(x[,2],decreasing=FALSE),]
sorted[1:10,]

#Efron method for handling ties (default)
coxph(Surv(time, status) ~ age + factor(sex) + bun + ca + hb + pcells + factor(protein), data = x, ties=c("efron"))

#exact method for handling ties (default)
coxph(Surv(time, status) ~ age + factor(sex) + bun + ca + hb + pcells + factor(protein), data = x, ties=c("exact"))




##deviance residual plot to detect any outliers

dev.resid <- residuals(fit, type="deviance")
plot(x$age, dev.resid, xlab="Age",ylab="Deviance Residuals")
#lines(lowess(x$age, dev.resid))


#schoenfeld residual plot to check PH assumption, x axis is observed EVENT time
sch.resid <- residuals(fit, type="schoenfeld")
plot(sorted[sorted$status==1,]$time, sch.resid[,1], xlab="Time",ylab="Schoenfeld residuals for age")
abline(h=0)


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