# R-commands for Cox regression
# =============================

# For illustration we will use the melanoma data 
# (treating deaths from other causes as censorings).

# We will use the survival library, so this has to be loaded

# We first read the data:
melanoma=read.table("http://www.uio.no/studier/emner/matnat/math/STK4080/h12/melanoma.txt",header=T)



# ESTIMATION OF REGESION COEFICIENTS
# ----------------------------------

# We first consider the model with sex as the only covariate:
fit.s=coxph(Surv(lifetime,status==1)~sex,data=melanoma)
summary(fit.s)

# Note that since sex is a binary covariate (coded 1 and 2), we get the same estimates if we treat the covariate 
# as a numeric covariate or as a categorical covariate [by using factor(sex) in the coxph-command] 


# Then we consider the model with sex and tumor thickness as covariates:
fit.st=coxph(Surv(lifetime,status==1)~sex+thickn,data=melanoma)
summary(fit.stu)


# Finally we consider the model with sex, tumor thickness, and ulceration:
fit.stu=coxph(Surv(lifetime,status==1)~sex+thickn+ulcer,data=melanoma)
summary(fit.stu)


# The three models may be compared using the likelihood ratio test:
anova(fit.s,fit.st,fit.stu,test="Chisq")

# Note that the computations above is not an example of model building by forward selection. 
# For in that case one would have included the covariates in the order of their significance 
# (first ulcetation, then tumor thickness, and finally sex).



# ESTIMATION OF CUMULATIVE HAZARDS AND SURVIVAL FUNCTIONS
# -------------------------------------------------------

# We first consider sex as the only covariate

# We start by making Nelson-Aalen plots for both genders:
fit.f=coxph(Surv(lifetime,status==1)~strata(sex),data=melanoma)
surv.f=survfit(fit.f)
plot(surv.f,fun="cumhaz", mark.time=F,xlim=c(0,10),ylim=c(0,0.70),
     xlab="Years since operation",ylab="Cumulative hazard",lty=1:2)
legend("topleft",c("Females","Males"),lty=1:2)

# We then fit a Cox model with sex as the only covariate and plot 
# the model based estmates of the cumulative hazards in the same plot:
fit.s=coxph(Surv(lifetime,status==1)~sex,data=melanoma)
surv.s=survfit(fit.s,newdata=data.frame(sex=c(1,2)))
lines(surv.s,fun="cumhaz", mark.time=F,conf.int=F, lty=1:2,col="red")


# We then consider the model with sex, tumor thickness, and ulceration
# and plot the cumulative hazards for the four covariate combinations
#    1) sex=1, thickn=2, ulcer=1
#    2) sex=1, thickn=2, ulcer=2
#    3) sex=2, thickn=5, ulcer=1
fit.stu=coxph(Surv(lifetime,status==1)~sex+thickn+ulcer,data=melanoma)
new.covariates=data.frame(sex=c(1,1,2),thickn=c(2,2,5),ulcer=c(1,2,1))
surv.stu=survfit(fit.stu,newdata= new.covariates)
plot(surv.stu,fun="cumhaz", mark.time=F,xlim=c(0,10), 
     xlab="Years since operation",ylab="Cumulative hazard",lty=1:3)
legend("topleft",c("1","2","3"),lty=1:3)


# To plot the survival functions for the same combinations of the 
# covariates we just omit the "cumhaz" option:
plot(surv.stu, mark.time=F,xlim=c(0,10), xlab="Years since operation",ylab="Survival probability",lty=1:3)
legend("bottomleft",c("1","2","3"),lty=1:3)