# Estiation for competing risks using the mstate package
#-------------------------------------------------------

# We will consider data on causes of death for Norwegian males

# Read data on causes of death and extract the data for males
path="http://folk.uio.no/borgan/abg-2008/data/causes_death.txt"
norw.death=read.table(path, header=T)
males=norw.death[norw.death$sex==1, ]
n.males=dim(males)[1]

#  The mstate package needs to be downloaded 
library(mstate)

# We define the possible transitions between the states
# (for general Markov chains this is done by the command "transMat"):
tmat=trans.comprisk(4,c("Alive","Cancer","Cardiovasc","Medical","Alcohol"))

# We convert the data-frame to long format, where there are four lines
# for each person (one for each competing cause):
males.long=msprep(time=c(NA,"agestop","agestop","agestop","agestop"), 
     status=c(NA,"dead1","dead2","dead3","dead4"),data=males ,trans=tmat,
     start=list(state=rep(1,n.males),time=males$agestart))

# Note that in mstate the states are numbered 1, 2, etc, so that 
# state 1 corresponds to "alive"


# We then fit a stratified Cox-model with no covariates, extract the cumulative 
# cause-specific hazard and make a plot of these (cf. Figure 3.6 in the ABG-book)
cox.males=coxph(Surv(Tstart,Tstop,status)~strata(trans),data=males.long, method="breslow")
haz.males=msfit(cox.males,trans=tmat)
plot(haz.males,xlab="Age",xlim=c(40,70))

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

# We may list the estimated probabilities by:
prob.males[[1]]
 
# The transition probabilities may be plotted in a stacked manner: 
plot(prob.males, type="stacked")
 
# Or they may be plotted as separate estimates:
plot(prob.males,type="single")