# Plots and tests based on cumulative sums of martingale residuals
# ----------------------------------------------------------------

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

# We use the timereg package
library(timereg)

# We first consider a model with ulceration and thickness (not log-transformed) 
fit.ut=cox.aalen(Surv(lifetime,status==1)~prop(ulcer)+prop(thickn),  
     data=melanoma, weighted.test=0, residuals=1,rate.sim=0,n.sim=1000)

#Check of log-linearity
resids.ut=cum.residuals(fit.ut,data=melanoma,cum.resid=1,n.sim=1000)
plot(resids.ut,score=2,xlab="Tumor thickness")
summary(resids.ut)

# We then check log-linearity for a model with log-transformed thickness 
fit.ult=cox.aalen(Surv(lifetime,status==1)~prop(ulcer)+prop(log2(thickn)),  
     data=melanoma, weighted.test=0, residuals=1,rate.sim=0,n.sim=1000)
resids.ult=cum.residuals(fit.ult,data=melanoma,cum.resid=1,n.sim=1000)
plot(resids.ult,score=2,xlab="log tumor thickness")
summary(resids.ult)

# Finally we check proportional hazards:
par(mfrow=c(1,2))
plot(fit.ult,score=T,xlab="Years since operation")
summary(fit.ult)


# Stratified Cox regression
# -------------------------

# We fit a model where we stratify on  ulceration use log-thickness as covariate 
fit.strat=coxph(Surv(lifetime,status==1)~log2(thickn)+strata(ulcer), data=melanoma)
summary(fit.strat)

# We may plot the cumuative baseline hazards for the two ulceration strata:
par(mfrow=c(1,1))
baseline.covar=data.frame(thickn=1)
surv.strat=survfit(fit.strat,newdata=baseline.covar)
plot(surv.strat,fun="cumhaz", mark.time=F,xlim=c(0,10), 
     xlab="Years since operation",ylab="Cumulative hazard",lty=1:2)


# Check of log-linearity using smoothing splines
# ----------------------------------------------

# We will check log-linearity of thickness in a model with tumor thickness and ulceration as covariates.
# We then used a penalized smoothing spline for the effect of thickness:
fit.pstu=coxph(Surv(lifetime,status==1)~factor(ulcer)+pspline(thickn), data=melanoma)
print(fit.pstu)
par(mfrow=c(1,2))
termplot(fit.pstu,se=T)

# Both the formal test (from the print command) and the plot (from the termplot command) 
# indicate that the effect of thickness is NOT log-linear. 


# The model checks above indicate that it may be better to use log2-thickness 
# as covariate so we define this covariate and log-linearity
melanoma$log2thick=log2(melanoma$thickn)
fit.pslogtu=coxph(Surv(lifetime,status==1)~factor(ulcer)+pspline(log2thick), data=melanoma)
print(fit.pslogtu)
par(mfrow=c(1,2))
termplot(fit.pslogtu,se=T)

# Now both the formal test (from the print command) and the plot (from the termplot command) 
# indicate that the effect of log2-thickness is reasonably log-linear. 



# Check of proportional hazards
# -----------------------------

# A simple method is to fit a stratified Cox model, where we stratify according to the levels of 
# a categorical covariate, and plot the stratum-specific cumulative baseline hazard estimates on a log-scale.
# If we have proportional hazards, the curves should be fairly parallel (i.e. have a constant vertical distance).

# We will use two small R-function to obtain the cumulative hazards for the different strata and the
# corresponding times where the cumulative hazard jumps
# (There may exist easier methods, but I have not made a systematic search here.)

cumhaz.cox=function(cox.surv)
{
no.strata=length(cox.surv$strata)
strata.ind=NULL
for (i in 1:no.strata){
strata.ind=c(strata.ind,rep(i,cox.surv$strata[i]))}
cumhaz.cox=vector('list',length=no.strata)
for (g in 1:no.strata) {
cumhaz.cox[[g]]=cox.surv$cumhaz[strata.ind==g]}
return(cumhaz.cox)
}

time.cox=function(cox.surv)
{
no.strata=length(cox.surv$strata)
strata.ind=NULL
for (i in 1:no.strata){
strata.ind=c(strata.ind,rep(i,cox.surv$strata[i]))}
time.cox=vector('list',length=no.strata)
for (g in 1:no.strata) {
time.cox[[g]]=cox.surv$time[strata.ind==g]}
return(time.cox)
}


# Then we check proportional hazards for ulceration:
par(mfrow=c(1,1))
fit.logtstu=coxph(Surv(lifetime,status==1)~strata(ulcer)+log2thick, data=melanoma)
surv.logtstu=survfit(fit.logtstu,newdata=data.frame(log2thick=0))
cumhaz.logtstu=cumhaz.cox(surv.logtstu)
time.logtstu=time.cox(surv.logtstu)
plot(time.logtstu[[1]],log(cumhaz.logtstu[[1]]),type='s',xlim=c(0,10),xlab="Years since operation",
     ylab="log cumulative hazard", main="Ulceration")
lines(time.logtstu[[2]],log(cumhaz.logtstu[[2]]),type='s',lty=2)

# The curves are fairly parallel, except for the first two years


# Then we check for thickness:
fit.sttu=coxph(Surv(lifetime,status==1)~factor(ulcer)+strata(grthick), data=melanoma)
surv.sttu=survfit(fit.sttu,newdata=data.frame(ulcer=1))
cumhaz.sttu=cumhaz.cox(surv.sttu)
time.sttu=time.cox(surv.sttu)
plot(time.sttu[[1]],log(cumhaz.sttu[[1]]),type='s',xlim=c(0,10),ylim=c(-4,0.5), 
     xlab="Years since operation", ylab="log cumulative hazard", main="Tumor thickness")
lines(time.sttu[[2]],log(cumhaz.sttu[[2]]),type='s',lty=2)
lines(time.sttu[[3]],log(cumhaz.sttu[[3]]),type='s',lty=3)

# There is a tendency that the curves are closer together for large values 
# of t than for smaller ones, indicating a non-proportional effect of thickness


# We then make a formal test for proportionality of the covariates. 
# This is done by, for each covariate x, adding time-dependent covariate x*log(t), 
# and testing whether the time-dependent covariates are significant using a score test:
fit.logtu=coxph(Surv(lifetime,status==1)~factor(ulcer)+log2thick, data=melanoma)
cox.zph(fit.logtu,transform='log')

# The test indicates that the effect of tumor-thickness is not proportional.
# The estimate we get for log-thicness in then a weighted average of the time-varying effect


# We also make plots that give nonparametric estimates of the (possible) time dependent effect of the covariates:
par(mfrow=c(1,2))
plot(cox.zph(fit.slogtu))