# R-commands for competing risks (and Markov chains)
# =================================================

# For illustration we will use the melanoma data and consider 
# the two causes of death 1: dead from malignant melanoma and 
# 2: death from other causes (for men and women combined).

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

# We will use the R package "mstate", so this has to be installed and loaded.


# The Cuminc-command may be used find the Kaplan-Meier estimate 
# for overall survival and the empirical cumulative incidence functions:
fit.ci=Cuminc(time="lifetime",status="status",data=melanoma,failcodes=c(1,4))
fit.ci

# There are no specific plotting commands that support the Cuminc-command, so we have to do the plotting "from skratch": 
plot(fit.ci$time,fit.ci$Surv,type="s",lty=1,xlim=c(0,12),ylim=c(0,1),xlab="Years since operation",ylab="Probability")
lines(fit.ci$time,fit.ci$CI.1,type="s",lty=2)
lines(fit.ci$time,fit.ci$CI.4,type="s",lty=3)



# We will also show how we may get the estimates by applying 
# commands that may be used (with minor modifications) for general Markov chains

# We first define the possible transitions between the states
#(for general Markov chains this is done by the command "transMat"):
tmat=trans.comprisk(2,c("alive","melanoma","other"))

# We then extend the melanoma data-frame with indicators for the two causes of death:
melanoma$mel.status=as.numeric(melanoma$status==1)
melanoma$other.status= as.numeric(melanoma$status==4)

# We convert the data-frame to long format, where there are two lines
# for each person (one for each cause of death):
melanomalong=msprep(time=c(NA,"lifetime","lifetime"), status=c(NA,"mel.status","other.status"), data=melanoma,trans=tmat)

# We then fit a stratified cox-model with no covariates, extract the 
# cumulative cause-specific hazard and make a plot of these
cox.mel=coxph(Surv(Tstart,Tstop,status)~strata(trans),data=melanomalong, method="breslow")
haz.mel=msfit(cox.mel,trans=tmat)
plot(haz.mel,xlab="Years since operation",xlim=c(0,12))

# We the obtain the estimated transition probabilities
prob.mel=probtrans(haz.mel,predt=0)

# We may list the estimated probabilities by:
prob.mel[[1]]

# The transition probabilities may be plotted in a stacked manner: 
plot(prob.mel,xlim=c(0,12),type="stacked")

# Or they may be plotted as separate estimates:
plot(prob.mel,xlim=c(0,12),type="single")