Practical exercise 5
In this exercise 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. See here for a description of the coding of the data.
To read the data into R you may give the command:
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"))
To perform the estimation, we need to convert the data-frame to a long format, where there are 2-3 lines for each patient:
bonemarrow.long=msprep(time=c(NA,"TP","T2"), status=c(NA,"P","DF"),data=bonemarrow,trans=tmat)
Then we fit a stratified Cox-model with no covariates (stratifying on the type of transition), extract the Nelson-Aalen estimates for the cumulative transition intensities and make a plot of these
cox.bm=coxph(Surv(Tstart,Tstop,status)~strata(trans),data=bonemarrow.long, method="breslow")
haz.bm=msfit(cox.bm,trans=tmat)
plot(haz.bm,xlab="Days post-transplant",xlim=c(0,1000))
To find the Aalen-Johansen estimates of the transition probabilities, we give the command:
prob.bm=probtrans(haz.bm,predt=0)
We may extract the empirical transition probabilities from state "transplanted" (state 1) by the command prob.bm[[1]], and similarly we obtain the empirical transition probabilities form state "recovered" (state 2) by prob.bm[[2]].
The transition probabilities may be plotted in a stacked manner or as separate estimates by the commands:
plot(prob.bm,type="stacked",ord=c(2,1,3))
plot(prob.bm,type="single",xlim=c(0,1000))
a) Perform the commands given above. Make sure that you understand what each of the commands and plots give you, and compare the results with those of example 3.16 in the ABG-book.
b) Above we used the probtrans command with the option predt=0. This gives transition probabilities from time s = 0 to time t. Also compute and plot the transition probabilities for s = 14, 28 and 42. What can you learn from these plots?