# EXERCISE 3.13
# ====================

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


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


# We first define the possible transitions between the states
tmat=trans.comprisk(2,c("in place","expelled","removed"))

# We then extend the iud data-frame with indicators for the two causes for ending iud use:
iud$expelled=as.numeric(iud$status==1)
iud$removed= as.numeric(iud$status==2)

# We convert the data-frame to long format, where there are two lines
# for each person (one for each cause of ending iud use):
iudlong=msprep(time=c(NA,"time","time"), status=c(NA,"expelled","removed"), data=iud,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.iud=coxph(Surv(Tstart,Tstop,status)~strata(trans),data=iudlong, method="breslow")
haz.iud=msfit(cox.iud,trans=tmat)
plot(haz.iud,xlab="Days since insertion")

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

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

# The transition probabilities may be plotted in a stacked manner: 
plot(prob.iud,type="stacked")

# Or they may be plotted as separate estimates:
plot(prob.iud,type="single")