myeloid.Rnw

\documentclass{article}[11pt]
\usepackage{amsmath}
\addtolength{\textwidth}{1in}
\addtolength{\oddsidemargin}{-.5in}
\setlength{\evensidemargin}{\oddsidemargin}

\newcommand{\code}[1]{\texttt{#1}}
\title{Myeloid data}
\author{Terry Therneau}
\date{19 Sept 2023}

\begin{document}
\maketitle

<<setup, echo=FALSE>>=
library(knitr)
opts_chunk$set(comment=NA, tidy=FALSE, highlight=FALSE, echo=FALSE,
               fig.width=4.5, fig.height=3, fig.path="figures/",
               device="pdf", dev.args=list(pointsize=8),
               cache=FALSE,   background="#ffffff",
               warning=FALSE, error=FALSE, message=FALSE, prompt=TRUE,
               strip.white=FALSE, mymar=TRUE)
options(contrasts= c("contr.treatment", "contr.poly"),
        show.signif.stars = FALSE, continue=" ", width=65)

# because "mymar" is set to TRUE above, this hook will be run on every chunk.
# it sets the margins, font, etc to values that work
knit_hooks$set(mymar= function(before, options, envir) {
    if (before) {
        look.for <- c("mar", "font", "cex")  # options we want
        plist <- options[match(look.for, names(options), nomatch=0)]
        if (is.null(plist$mar)) plist$mar <- c(4, 4, .5, .5)
        if (is.null(plist$cex)) plist$cex <- 1
        do.call("par", plist)
    } else NULL
})

table2 <- function(...) table(..., useNA = "ifany")
library(survival)
library(splines)
library(xtable)

crisk <- function(what, horizontal = TRUE, ...) {
    nstate <- length(what)
    connect <- matrix(0, nstate, nstate,
                      dimnames=list(what, what))
    connect[1,-1] <- 1  # an arrow from state 1 to each of the others
    if (horizontal) statefig(c(1, nstate-1),  connect, ...)
    else statefig(matrix(c(1, nstate-1), ncol=1), connect, ...)
}
@

\section{Myeloid data}
The \code{myeloid} data set contains data from a clinical trial
in subjects with  acute myeloid leukemia.  To protect patient confidentiality
the data set in the survival package has been slightly perturbed, but 
results are essentially unchanged.
In this comparison of two conditioning regimens, the
canonical path for a subject is initial therapy $\rightarrow$
complete response (CR) $\rightarrow$
hematologic stem cell transplant (SCT) $\rightarrow$
sustained remission, followed by relapse or death.
Not everyone follows this ideal path, of course.

<<overall>>=
myeloid[1:5,]
@ 
The first few rows of data are shown above.
The data set contains the follow-up time and status at last follow-up 
for each subject, along with the time to transplant
(txtime),
complete response (crtime) or relapse after CR (rltime).
Subject 1 did not receive a transplant, as shown by the NA value,
and subject 2 did not achieve CR.

\begin{figure}
<<sfit0, echo=FALSE>>=
sfit0 <- survfit(Surv(futime,death) ~ trt, myeloid)
sfit0b <- survfit(Surv(futime, death) ~ sex, myeloid)
sfit0c <- survfit(Surv(futime, death) ~ flt3, myeloid)
oldpar <- par(mfrow=c(2,2), mar=c(5,5,1,1))
plot(sfit0, xscale=365.25, col=1:2, lwd=2,  xmax=4*365.25, fun="event",
     xlab="Years from randomization", ylab="Death")
legend(730, .25, c("Trt A", "Trt B"), col=1:2, lwd=2, lty=1, bty='n')
plot(sfit0b, xscale=365.25, col=1:2, lwd=2,  xmax=4*365.25, fun="event",
     xlab="Years from randomization", ylab="Death")
legend(730, .25, c("Female", "Male"), col=1:2, lwd=2, lty=1, bty='n')
plot(sfit0c, xscale=365.25, col=1:3, lwd=2,  xmax=4*365.25, fun="event",
     xlab="Years from randomization", ylab="Death")
legend(730, .25, c("TKD only", "ITD low", "ITD high"), col=1:3, lwd=2, 
       lty=1, bty='n')
par(oldpar)
@ 
  \caption{Overall survival curves for the myeloid study, by treatment
    arm, sex, and FLT3 mutation group.}
  \label{sfit0}
\end{figure}

Overall survival curves for the data are shown in figure \ref{sfit0}.
The difference between the treatment arms A and B
is substantial.
A goal of this analysis is to better understand this difference. 
The presence of mutations in the FLT3 gene was an eligibility criteria
for the study; subtypes of the mutation were a stratification factor and
are also significant predictors of the outcome.
This is reflected in the simple Cox model below.

<<cfit0>>=
cfit0 <- coxph(Surv(futime, death) ~ trt + strata(flt3), myeloid)
cfit1 <- coxph(Surv(futime, death) ~ trt + sex + flt3, myeloid)
print(cfit1, digits=2)
@ 

\section{Multi-state curves}

The full multi-state data set can be created with the 
\code{tmerge} routine.
<<sfit0a, echo=TRUE>>=
mdata <- tmerge(myeloid[,1:4], myeloid,  id=id,  death= event(futime, death),
                sct = event(txtime), cr = event(crtime), 
                relapse = event(rltime))
temp <- with(mdata, cr + 2*sct  + 4*relapse + 8*death)
table(temp) 
@ 

Our check shows that there is one subject who had CR and stem cell transplant
on the same day (temp=3).
In the multi-state data set the software does not allow the subject to have
two transitions at the same exact time, e.g., entry:cr and an immediate 
cr:transplant transition.
To avoid length 0 intervals, we break the tie so that complete response (CR)
happens first.  The amount of time shift added below is completely arbitrary,
you could for instance use only .1 day; the choice will have an impact on summary
statistics such as mean time in state of course, but the effect is small.
(Students may be surprised to see anomalies like this, since they never appear
in textbook data sets.  In real data such issues always appear.)

<<sfit0b, echo=TRUE>>=
tdata <- myeloid  # temporary working copy
tied <- with(tdata, (!is.na(crtime) & !is.na(txtime) & crtime==txtime))
tdata$crtime[tied] <- tdata$crtime[tied] -1
mdata <- tmerge(tdata[,1:4], tdata,  id=id,  death= event(futime, death),
                sct = event(txtime), cr = event(crtime), 
                relapse = event(rltime),
                priorcr = tdc(crtime), priortx = tdc(txtime),
                priorrel = tdc(rltime))
temp <- with(mdata, cr + 2*sct  + 4*relapse + 8*death)
table(temp)
mdata$event <- factor(temp, c(0,1,2,4,8),
                       c("none", "CR", "SCT", "relapse", "death"))

mdata[1:7, c("id", "trt", "tstart", "tstop", "event", "priorcr", "priortx")]
@ 

Subject 1 has a CR on day 44, relapse on day 113, death on day 235 and
did not receive a stem cell transplant. 
The data for the first three subjects looks good.  
Check it out a little more thoroughly using survcheck.

<<check, echo=TRUE>>=
survcheck(Surv(tstart, tstop, event) ~1, mdata, id=id)
@ 

The second table shows that no subject had more than one CR, SCT,
relapse, or death; the intention of the study was to count only the first
of each of these, so this is as expected.
Several subjects visited all four intermediate states.  
The transitions table shows 11 subjects who achieved CR \emph{after} stem
cell transplant and another 106 who received a transplant before
achieving CR, both of which are deviations from the ``ideal'' pathway.
No subjects went from death to another state or had further follow-up time
after death (which is good).

For investigating the data we would like to add a set of alternate endpoints.
\begin{enumerate}
  \item The competing risk of CR and death, ignoring other states.  This
    is used to estimate the fraction who ever achieved a complete response.
  \item The competing risk of SCT and death, ignoring other states.
  \item An endpoint that distinguishes death after SCT from death 
    before SCT.
\end{enumerate}
Each of these can be accomplished by adding further outcome variables to 
the data set, we do not need to change the time intervals.

<<newevent>>=
temp1        <- with(mdata, ifelse(priorcr, 0, c(0,1,0,0,2)[event]))
mdata$crstat <- factor(temp1, 0:2, c("none", "CR", "death"))

temp2        <- with(mdata, ifelse(priortx, 0, c(0,0,1,0,2)[event]))
mdata$txstat <- factor(temp2, 0:2, c("censor", "SCT", "death"))

temp3     <- with(mdata, c(0,0,1,0,2)[event] + priortx)
mdata$tx2 <- factor(temp3, 0:3,
                    c("censor", "SCT", "death w/o SCT", "death after SCT"))
@ 

\begin{figure}
<<curve1>>=
# I want to have the plots in months, it is simpler to fix time
#  once rather than repeat xscale many times
tdata$month <-  tdata$futime * 12 /365.25
mdata$mstart <- mdata$tstart * 12 /365.25
mdata$mstop  <- mdata$tstop * 12 /365.25


sfit1 <- survfit(Surv(month, death)  ~ trt, tdata) # survival
sfit2 <- survfit(Surv(mstart, mstop, crstat) ~ trt, 
                 data= mdata, id = id) # CR
sfit3 <- survfit(Surv(mstart, mstop, txstat) ~ trt, 
                 data= mdata, id =id) # SCT

layout(matrix(c(1,1,1,2,3,4), 3,2), widths=2:1)
oldpar <- par(mar=c(5.1, 4.1, 1.1, .1))

mlim   <- c(0, 48) # and only show the first 4 years
plot(sfit2[,"CR"], xlim=mlim, 
         lty=3, lwd=2, col=1:2, xaxt='n',
     xlab="Months post enrollment", ylab="Fraction with the endpoint")
lines(sfit1, mark.time=FALSE, xlim=mlim,
      fun='event', col=1:2, lwd=2)

lines(sfit3[,"SCT"], xlim=mlim, col=1:2, 
          lty=2, lwd=2)

xtime <- c(0, 6, 12, 24, 36, 48)
axis(1, xtime, xtime) #axis marks every year rather than 10 months
temp <- outer(c("A", "B"), c("CR", "transplant", "death"),  paste)
temp[7] <- ""
legend(25, .3, temp[c(1,2,7,3,4,7,5,6,7)], lty=c(3,3,3, 2,2,2 ,1,1,1),
       col=c(1,2,0), bty='n', lwd=2)
abline(v=2, lty=2, col=3)

# add the state space diagrams
par(mar=c(4,.1,1,1))
crisk(c("Entry", "CR", "Death"), alty=3)
crisk(c("Entry", "Tx", "Death"), alty=2)
crisk(c("Entry","Death"))
par(oldpar)
layout(1)
@
  \caption{Overall survival curves: time to death, to transplant (Tx), 
    and to complete response (CR).
    Each shows the estimated fraction of subjects who have ever reached the
    given state.  The vertical line at 2 months is for reference.
    The curves were limited to the first 48 months to more clearly show 
    early events. The right hand panel shows the state-space model for each
  pair of curves.}
  \label{curve1}
\end{figure}

This data set is the basis for our first set of curves, which are shown in
figure \ref{curve1}.
The plot overlays three separate \code{survfit} calls: standard survival
until death, complete response with death as a competing risk, and
transplant with death as a competing risk.
For each fit we have shown only one state: the fraction
who have died, fraction ever in CR, and fraction ever to receive transplant,
respectively.
Most of the CR events happen before 2 months (the green
vertical line) and most of the additional CRs
conferred by treatment B appear to occur between months 2 and 8.
This visual impression is however contradicted by a simple tabulation of the
CR events shown below.

<<crtable>>=
ctemp <- subset(mdata, event=="CR")
with(ctemp, table(trt, "CR month"= cut(mstop, c(0:3, 6,20))))
@  
 
Most transplants happen after 2 months, which is consistent with the
clinical guide of transplant after CR.
The survival advantage for treatment B begins between 4 and 6 months, 
which argues that it could be at least partially a consequence of the
additional CR events.

The association between a particular curve and its corresponding state space
diagram is critical.  As we will see below, many different models are
possible and it is easy to get confused.
Attachment of a diagram directly to each curve, as was done above,
will not necessarily be day-to-day practice, but the state space should
always be foremost.  If nothing else,  draw it on a scrap of paper and tape it 
to the side of the terminal when creating a data set and plots.

\begin{figure}
<<badfit>>=
badfit <- survfit(Surv(mstart, mstop, event=="SCT") ~ trt, 
                       id=id, mdata, subset=(priortx==0))

layout(matrix(c(1,1,1,2,3,4), 3,2), widths=2:1)
oldpar <- par(mar=c(5.1, 4.1, 1.1, .1))
plot(badfit, fun="event", xmax=48, xaxt='n', col=1:2, lty=2, lwd=2,
     xlab="Months from enrollment", ylab="P(state)")
axis(1, xtime, xtime)
lines(sfit3[,2], xmax=48, col=1:2, lwd=2)
legend(24, .3, c("Arm A", "Arm B"), lty=1, lwd=2,
       col=1:2, bty='n', cex=1.2)

par(mar=c(4,.1,1,1))
crisk(c("Entry", "transplant"), alty=2, cex=1.2)
crisk(c("Entry","transplant", "Death"), cex=1.2)
par(oldpar)
layout(1)
@ 
 \caption{Correct (solid) and invalid (dashed) estimates of the number
  of subjects transplanted.}
 \label{badfit}
\end{figure}
 
Figure \ref{badfit} shows the transplant curves overlaid with the naive KM that
censors subjects at death.  There is no difference in the initial portion as
no deaths have yet intervened, but the final portion of the curve overstates the 
transplant outcome by more than 10\%. 
\begin{enumerate}
  \item The key problem with the naive estimate is that subjects who die can 
    never have a transplant.  The result of censoring them 
    is an estimate of the ``fraction who would
    be transplanted, if death before transplant were abolished''.  This is not
    a real world quantity.
  \item In order to estimate this fictional quantity one needs to assume that
    death is uninformative with respect to future disease progression.  The
    early deaths in months 0--2, before transplant begins, are however a very
    different class of patient.  Non-informative censoring is untenable.
\end{enumerate}
We are left with an unreliable estimate of an uninteresting quantity.
Mislabeling any true state as censoring is always a mistake, one that 
will not be repeated here.
(There are cases where the ``what if'' curve would be of interest for policy
or other reasons, e.g., expected survival if lung cancer were eliminated.
The necessary assumptions to estimate such a curve are however very strong,
and often untenable.)

\begin{figure}
<<cr2>>=
state3 <- function(what, horizontal=TRUE, ...) {
    if (length(what) != 3) stop("Should be 3 states")
    connect <- matrix(c(0,0,0, 1,0,0, 1,1,0), 3,3,
                      dimnames=list(what, what))
    if (horizontal) statefig(1:2, connect, ...)
    else statefig(matrix(1:2, ncol=1), connect, ...)
}

temp <- c(1,2,1,3,3)[as.numeric(mdata$event)]
temp[mdata$priorrel==1] <- 1  # don't treat death as a new event after relapse
mdata$cr2 <- factor(temp, 1:3, c("none", "CR", "death/relapse"))

crsurv <- survfit(Surv(mstart, mstop, cr2) ~ trt,
                  data= mdata, id=id, influence=TRUE)
# the influence=TRUE allows for a standard error of sojourn time, later

layout(matrix(c(1,1,2,3), 2,2), widths=2:1)
oldpar <- par(mar=c(5.1, 4.1, 1.1, .1))
plot(sfit2[,2], lty=3, lwd=2, col=1:2, xmax=12, 
     xlab="Months", ylab="CR")
lines(crsurv[,2], lty=1, lwd=2, col=1:2)
par(mar=c(4, .1, 1, 1))
crisk( c("Entry","CR", "Death"), alty=3)
state3(c("Entry", "CR", "Death/Relapse"))

par(oldpar)
layout(1)
@ 
  \caption{Models for `ever in CR' and `currently in CR';
    the only difference is an additional transition.
    Both models ignore transplant.}
  \label{cr2}
\end{figure}

Complete response is a goal of the initial therapy; figure \ref{cr2}
looks more closely at this.
As was noted before arm B has an increased number of responses.
The duration of response is also increased: 
the solid curves show the number of subjects still in response, and
we see that they spread apart as much as the dotted ``ever in response''
curves.
The figure shows only the first year in order to better visualize
the details, but continuing the curves out to 48 months reveals a similar
pattern  
 ($.77- .68= .11$ vs $.53 - .40= .13$ at 12 months, .11 vs .09 at 48).

<<cr2b, echo=TRUE>>=
print(crsurv, rmean=48, digits=2)
@ 

<<cr2c, echo=FALSE>>=
temp <- summary(crsurv, rmean=48)$table
delta <- round(temp[4,3] - temp[3,3], 2)
@ 

@ 
The restricted mean time in the CR state is extended by 
\Sexpr{round(temp[4,3], 1)} - \Sexpr{round(temp[3,3], 1)} =
\Sexpr{round(delta,1)} months.  
A question which immediately gets asked is whether this difference
is ``significant'', to which there are two answers. 
The first and more important is to ask whether 5 months is an important gain
from either a clinical or patient perspective. 
The overall restricted mean survival for the study is approximately
30 of the first 48 months post entry;
on this backdrop an extra 5 months in CR might or might not be viewed as an
meaningful advantage.
The less important answer is to test whether the apparent gain is sufficiently
rare from a mathematical point of view, i.e., ``statistical'' significance.
The standard errors of the two values are
\Sexpr{round(temp[3,4],1)} and \Sexpr{round(temp[4,4],1)},
and since they are based
on disjoint subjects the values are independent, leading to a standard error
for the difference of $\sqrt{1.1^2 + 1.1^2} = 1.6$.
The 5 month difference is more than 3 standard errors, so highly significant.


The code for the figure created yet another event
variable so as to ignore transitions to the transplant state.
For the new \code{cr2} variable `transplant' becomes 'none', and relapse and
death are combined.  Transplant is essentially ignored.
Below we show a small part of the data set. 
The \code{event} variable is the overall code for the subject's progression
through states, while \code{crstat} and \code{txstat} deal with the competing
risk for CR/death and SCT/death, respectively, which are shown in
figure \ref{curve1}.
The \code{crstat} variable pick up the first of CR or death, and codes all
other lines as `none'.  
The \code{cr2} event is used for the duration of CR curves in figure \ref{cr2},
it captures CR and death/relapse, with all other rows as 'none'.

<<showstate>>=
subset(mdata, id<6, c(id, tstart, tstop, event, crstat, cr2, txstat))
@ 

For the estimates themselves, do we need to remove redundant lines. That is
should the CR/death curve be one of the two lines below and not the other?
It turns out that the results are the same.  None of the ``extra'' rows after
CR contain a transition and hence they all get ignored, since they are not
in the risk set for any possible transition.

<<codec, echo=TRUE, eval=FALSE>>=
curve1 <- survfit(Surv(tstart, tstop, crstat) ~ trt, mdata, id=id)
curve2 <- survfit(Surv(tstart, tstop, crstat) ~ trt, mdata, id=id,
                  subset = (priorcr ==0))
@ 


\begin{figure}
<<txsurv>>=
event2 <- with(mdata, ifelse(event=="SCT" & priorcr==1, 6,
               as.numeric(event)))
event2 <- factor(event2, 1:6, c(levels(mdata$event), "SCT after CR"))
txsurv <- survfit(Surv(mstart, mstop, event2) ~ trt, mdata, id=id,
                  subset=(priortx ==0))
dim(txsurv) # number of strata by number of states
txsurv$states # Names of states

layout(matrix(c(1,1,1,2,2,0),3,2), widths=2:1)
oldpar <- par(mar=c(5.1, 4.1, 1,.1))
plot(txsurv[,c(3,6)], col=1:2, lty=c(1,1,2,2), lwd=2, xmax=48,
     xaxt='n', xlab="Months", ylab="Transplanted")
axis(1, xtime, xtime)
legend(15, .13, c("A, transplant without CR", "B, transplant without CR",
                 "A, transplant after CR", "B, transplant after CR"),
       col=1:2, lty=c(1,1,2,2), lwd=2, bty='n')

# Add the state figure
sname <- c("Entry", "CR", "Transplant", "Transplant")
xypos <- cbind(c(1/2, 3/4, 1/4, 3/4),
                    c(5/6, 1/2, 1/2, 1/6))
connect <- matrix(0,4,4, dimnames=list(sname, sname))
connect[1, 2:3] <- 1
connect[2,4] <- 1
statefig(xypos, connect)

par(oldpar)
layout(1)
@ 
   \caption{Transplant status of the subjects, broken down by whether it
    occurred before or after CR.}
  \label{txsurv}
\end{figure}

In summary
  \begin{itemize}
    \item Arm B adds further complete responses (about 10\%);
      206/317 = 65\% achieve CR in arm A vs. 248/329 = 75\% in arm B.
    \item The difference in 4 year survival is about 6\%.
    \item There is approximately 2 months longer average duration of CR (of 48).
  \end{itemize}

CR $\rightarrow$ transplant is the target treatment path for a
patient; given the improvements listed above
why does figure \ref{curve1} show no change in the number transplanted?
Figure \ref{txsurv} shows the transplants broken down by whether this
happened before or after complete response.
Most of the non-CR transplants happen by 10 months.
One possible explanation is that once it is apparent to the
patient/physician pair that CR is not going to occur, they proceed forward with
other treatment options.  
The extra CR events on arm B lead to 
a consequent increase in transplant as well, but at a later time of 12--24 
months: for a subject in CR we can perhaps afford to defer the transplant date.

Computation is again based on a manipulation of the event variable: in this
case dividing the transplant state into two sub-states based on the presence
of a prior CR.  The code makes use of the time-dependent covariate 
\code{priorcr}.
(Because of scheduling constraints within a hospital it is unlikely that
a CR that is within a few days prior to transplant could have affected the
decision to schedule a transplant, however.  An alternate breakdown that
might be useful would be ``transplant without CR or within 7 days after CR''
versus those that are more than a week later. 
There are many sensible questions that can be asked.)


\begin{figure}
<<sfit4>>=
sfit4 <- survfit(Surv(mstart, mstop, event) ~ trt, mdata, id=id)
sfit4$transitions
layout(matrix(1:2,1,2), widths=2:1)
oldpar <- par(mar=c(5.1, 4.1, 1,.1))
plot(sfit4, col=rep(1:4,each=2), lwd=2, lty=1:2, xmax=48, xaxt='n',
     xlab="Months", ylab="Current state")
axis(1, xtime, xtime)
text(c(40, 40, 40, 40), c(.51, .13, .32, .01),
     c("Death", "CR", "Transplant", "Recurrence"), col=c(4,1,2,3))

par(mar=c(5.1, .1, 1, .1))
sname <- c("Entry", "CR", "SCT", "Relapse", "Death")
connect <- matrix(0, 5, 5, dimnames=list(sname, sname))
connect[1, -1] <- c(1,1,1, 1.4)
connect[2, 3:5] <- c(1, 1.4, 1)
connect[3, c(2,4,5)] <- 1
connect[4, c(3,5)]  <- 1
statefig(matrix(c(1,3,1)), connect, cex=.8)
par(oldpar)
layout(1)
@ 
  \caption{The full multi-state curves for the two treatment arms.}
  \label{sfit4}
\end{figure}

Figure \ref{sfit4} shows the full set of state occupancy probabilities for the
cohort over the first 4 years.  At each point in time the curves
estimate the fraction of subjects currently in that state.
 The total who are in the transplant state peaks at
about 9 months and then decreases as subjects relapse or die; 
the curve rises
whenever someone receives a transplant and goes down whenever someone
leaves the state.
At 36 months treatment arm B (dashed) has a lower fraction who have died,
the survivors are about evenly split between those who have received a
transplant and those whose last state is a complete response
(only a few of the latter are post transplant). 
The fraction currently in relapse -- a transient state -- is about 5\% for
each arm.
The figure omits the curve for ``still in the entry state''.
The reason is that 
at any point in time the sum of the 5 possible states is 1 ---
everyone has to be somewhere.  Thus one of the curves
is redundant, and the fraction still in the entry state is the least 
interesting of them.

\section{Multi-state models}
We now repeat some of the above using MSH models.  
An advantage in principle is that we can jointly adjust for treament and
FLT3 group, though the fact that the study is randomized and thus balanced
for FLT3 suggests that any changes in results will be modest.
The MLT also assumes proportional hazards, i.e., a constant coefficient or
hazard ratio over time;
figure \ref{sfit0} suggests that this may be questionable, for the first half 
year at least, for both treatment and FLT3.
A formal test of PH is not rejected, however.

<<crmod>>=
crfit1 <- coxph(Surv(mstart,mstop, crstat) ~ trt + flt3, mdata, id=id)
crfit2 <- coxph(Surv(mstart,mstop, cr2) ~ trt + flt3, mdata, id=id)

print(crfit1, digits=2)

print(crfit2, digits=2)
@ 

The above are fits for the two models in the right margin of figure \ref{cr2}.
It is interesting that FLT3 group has minimal effect on the probability
of reaching CR, or on the rate of death without CR, 
but a major effect on the transition from CR to death.
We can see this clearly in figure \ref{cr3b}, and in the sojourn time table
below.  
The range between the three FLT3 levels is smaller for time in the entry
state (1.6 months) than for time in the CR state (12.1), while the difference
between treatments A and B is fairly constant across FLT3.  

<<cr3a>>=
cdummy <- expand.grid(flt3 = levels(mdata$flt3), trt=c("A","B"))
crsurv2 <- survfit(crfit2, newdata= cdummy)
temp <- summary(crsurv2, rmean=48)$table
temp <- array(temp[1:12, 3], dim=c(3,4),
              dimnames=list(levels(mdata$flt3), NULL)) 
temp2 <- cbind(temp[,1:2], temp[,2]-temp[,1], NA, temp[,3:4], temp[,4]-temp[,3])
colnames(temp2) <- c("Entry A", "Entry B", "diff", "    ", 
                     "CR A", "CR B", "diff")
print(round(temp2, 1), na.print="")
@ 

\begin{figure}
<<cr3b>>=
plot(crsurv2[,2], col=1:3, lty=c(1,1,1,2,2,2), lwd=2, xaxt='n', xmax=48,
     xlab="Months from entry", ylab="Currently in CR")
axis(1, xtime, xtime)
legend(5, .2, levels(mdata$flt3), col=1:3, lty=1, bty='n', lwd=2)
@   
  \caption{Estimated probability of being in the CR state.  Treatment A =
    solid line, B = dashed.}
 \label{cr3b}
\end{figure}

I was puzzled at first why the coefficients for the entry:CR transition in
the first fit are different than the the second.  The answer is the 13 subjects
who go directly from s0 to relapse without passing throug CR: in the first fit
they remain at risk for a CR event for just a bit longer.  
The resulting small increase in a few risk sets
moves the coefficients just a smidge.  
(Per the study design, subjects were supposed to only be counted as a relapse 
if they had experienced CR, so this was not anticipated.)

For AJ modeling transplant we had divided transplant into two states: before 
and after CR.  There are then 6 states, the Aalen-Johansen estimate 
with treatment and FLC3 as
covariates creates 6 x 6= 36 curves, we ignored all but the two SCT states
when plotting and printing.
An MPH model is focused instead on the transitions: for this splitting the
SCT state only adds complexity, replacing the CR:SCT transition with
CR:(SCT after CR), a 1:1 substitution, but also splitting the 3 transitions
from SCT to CR, relapse, and death into 2.

Instead start by focusing our question.  
We already know that more non-CR transplants come from arm A than B, 
that more CR transplant come from B than A,  and that FLT3 appears to
have little effect on CR rates.  What further effect might FLT3 have
on transplant rates?
One possible effect, both before and after CR, would be acceleration of
transplant, for instance if the procedure is being delayed until failure seems
immenent.  For non CR patients the delay may be in hope of CR, for those
in CR there may be a desire to let the patient enjoy their treatment
holiday for a bit longer. 
For this question, simplify the states by ignoring relapse. The table below
shows coefficients.

<<txfit1>>=
temp <- c(1,2,3,1,4)[mdata$event]
mdata$event3 <- factor(temp, 1:4, c('none', "CR","SCT", "death"))
survcheck(Surv(mstart, mstop, event3) ~1, mdata, id=id)$transitions

txsurv <- coxph(Surv(mstart, mstop, event3) ~trt + flt3, mdata, id=id)
# Our interest lies in the 1:3 and 2:3 transitions

temp <- summary(txsurv)$coefficients
temp2 <- temp[grepl("1:3", rownames(temp)) | grepl("2:3", rownames(temp)),]
temp3 <- rbind(temp2[1:3,], NA, temp2[4:6,])
print(signif(temp3[,c(1,2,5,6)], 2), na.print='')
@ 

The FLT3 level appears to have little or no effect on the rate of entry:SCT
transitions, but a strong effect on CR:SCT transitions.
This may simply be an aspect of shorter duration of CR for the more severe
mutations. 
Treatment B has a somewhat lower rate of entry:SCT transitions.  We can only
speculate, but perhaps this may be an anticipatory effect of the higher CR
rate, i.e., subjects on B not yet declared to be CR are nevertheless 
less likely to be seen as a clear failure and promoted to early transplant.

\end{document}