R-TUTORIAL ON THE ILLNESS-DEATH MODEL

This tutorial is based on Practical exercise 5, STK4080, 2014.

In this tutorial we will see how we may find the Nelson-Aalen estimates for the cumulative transition intensities and the Aalen-Johansen estimates for the transition probabilities in a Markov illness-death model; cf. section 3.4.2 in the ABG-book. To this end we will use the bone marrow transplantation data described in example 1.13 and used for illustration in example 3.16. 

A description of the coding of the data is given here:

https://www.uio.no/studier/emner/matnat/math/STK4080/h14/bone-marrow-info.html

To read the data into R you may give the commands: 

path="http://www.uio.no/studier/emner/matnat/math/STK4080/h14/bone-marrow.txt"
bonemarrow=read.table(path,header=T)

The data contain information for 137 patients with acute leukemia who had a bone marrow transplantation. The patients were followed for a maximum of seven years, and times to relapse and death were recorded. It was also recorded if and when the platelet count of a patient returned to a self-sustaining level. The possible events for a patient may be described by an illness-death model without recovery with the three states "transplanted", "platelet recovered", and "relapsed or dead". A patient starts out in state "transplanted" at time zero when he/she gets the bone marrow transplant. If the platelets recover, the patient moves to state "platelet recovered", and if he/she then relapses or dies, the patient moves on to state "relapsed or dead". If the patient relapses or dies without the platelets returning to a normal level, he moves directly from state "transplanted" to state "relapsed or dead". 

To find the Nelson-Aalen estimates and the empirical transition probabilities (i.e. Aalen-Johansen estimates) we will use the 'mstate' package so this has to be installed and loaded.

We start out by defining the states and the possible transitions between them. This is achieved by the command: 

tmat=transMat(list(c(2,3),c(3),c()), names=c("transplanted","recovered","relapsed.dead"))

giving the output:
> tmat
               to
from            transplanted recovered relapsed.dead
  transplanted            NA         1             2
  recovered               NA        NA             3
  relapsed.dead           NA        NA            NA


To perform the estimation, we need to convert the data-frame to a long format, where there are 2-3 lines for each patient: 

bonemarrowlong=msprep(time=c(NA,"TP","T2"), status=c(NA,"P","DF"), data=bonemarrow,trans=tmat)

Output:
> head(bonemarrowlong)
An object of class 'msdata'

Data:
  id from to trans Tstart Tstop time status
1  1    1  2     1      0    13   13      1
2  1    1  3     2      0    13   13      0
3  1    2  3     3     13  2081 2068      0
4  2    1  2     1      0    18   18      1
5  2    1  3     2      0    18   18      0
6  2    2  3     3     18  1602 1584      0

Thus, for each individual there is one line for each possible transition (these are numbered 1,2,3, as seen from 'tmat' above). Consider individual 1: This individual was at risk for both transitions from 1 to 2 and 1 to 3 for 13 time units. The revealed transition was to state 2 (status = 1) and not to state 3 (status = 0). Then from time 13 to time 2081, individual 1 was at risk for transition from 2 to 3. There was instead a censoring, so status = 0 for this transition.


Now we fit a stratified Cox-model with no covariates (stratifying on the type of transition), and then we extract the Nelson-Aalen estimates for the cumulative transition intensities and make a plot of these cumulative transition probabiities- These correspond to Fig. 3.20 p. 120 in ABG.

cox.bm=coxph(Surv(Tstart,Tstop,status)~strata(trans),data=bonemarrowlong, method="breslow")
haz.bm=msfit(cox.bm,trans=tmat)
plot(haz.bm,xlab="Days post-transplant",xlim=c(0,1000),legend.pos="bottomright")

To find the Aalen-Johansen estimates of the transition probabilities, we give the command: 

prob.bm=probtrans(haz.bm,predt=0) 

We may then list the estimated transition probabilities from state 1 and state 2, respectively, by:

prob.bm[[1]]
prob.bm[[2]]

The transition probabilities may be plotted in a stacked manner: 

plot(prob.bm,type="stacked",ord=c(2,1,3))

Or they may be plotted as separate estimates:

plot(prob.bm,type="single",xlim=c(0,1000))

Above we used the probtrans command with the option predt=0. This gives transition probabilities from time s = 0 to time t. One may also compute and plot the transition probabilities for other values of s, for example s = 14, 28 and 42, see below. What can you learn from these plots?

s=14:
prob.bm=probtrans(haz.bm,predt=14)
plot(prob.bm,type="single",xlim=c(0,1000))


REFERENCE:

de Wreede, Liesbeth C., Marta Fiocco, and Hein Putter. "mstate: an R package for the analysis of competing risks and multi-state models." Journal of statistical software 38.7 (2011): 1-30.