\documentclass{article}[11pt] \usepackage{Sweave} \usepackage{amsmath} \addtolength{\textwidth}{1in} \addtolength{\oddsidemargin}{-.5in} \setlength{\evensidemargin}{\oddsidemargin} %\VignetteIndexEntry{Multi-state models and competing risks} \SweaveOpts{keep.source=TRUE, fig=FALSE} % Ross Ihaka suggestions \DefineVerbatimEnvironment{Sinput}{Verbatim} {xleftmargin=2em} \DefineVerbatimEnvironment{Soutput}{Verbatim}{xleftmargin=2em} \DefineVerbatimEnvironment{Scode}{Verbatim}{xleftmargin=2em} \fvset{listparameters={\setlength{\topsep}{0pt}}} \renewenvironment{Schunk}{\vspace{\topsep}}{\vspace{\topsep}} % I had been putting figures in the figures/ directory, but the standard % R build script does not copy it and then R CMD check fails \SweaveOpts{prefix.string=compete,width=6,height=4} \newcommand{\myfig}[1]{\includegraphics[height=!, width=\textwidth] {compete-#1.pdf}} \setkeys{Gin}{width=\textwidth} <>= options(continue=" ", width=60) options(SweaveHooks=list(fig=function() par(mar=c(4.1, 4.1, .3, 1.1)))) pdf.options(pointsize=10) #text in graph about the same as regular text options(contrasts=c("contr.treatment", "contr.poly")) #ensure default library("survival") @ \title{Multi-state models and competing risks} \author{Terry Therneau \and Cynthia Crowson \and Elizabeth Atkinson} \newcommand{\code}[1]{\texttt{#1}} \begin{document} \maketitle \section{Multi-state models} \begin{figure} <>= oldpar <- par(mar=c(.1, .1, .1, .1), mfrow=c(2,2)) cmat1 <- matrix(c(0,0,1,0), nrow=2, dimnames=list(c("Alive", "Dead"), c("Alive", "Dead"))) statefig(c(1,1), cmat1, cex=1.5) state4 <- c("0", "1", "2", "...") cmat4 <- matrix(0, 4,4, dimnames= list(state4, state4)) cmat4[1,2] <- cmat4[2,3] <- cmat4[3,4] <- 1 statefig(c(1,1,1,1), cmat4, bcol=c(1,1,1,0), cex=c(1.5, 1.5, 1.5, 3)) states <- c("A", "D1", "D2", "D3") cmat2 <- matrix(0, 4,4, dimnames=list(states, states)) cmat2[1, -1] <- 1 statefig(c(1,3), cmat2, cex=1.5) state3 <- c("Health", "Illness", "Death") cmat3 <- matrix(0, 3, 3, dimnames = list(state3, state3)) cmat3[1,2] <- cmat3[2,1] <- cmat3[-3, 3] <- 1 statefig(c(1,2), cmat3, offset=.03, cex=1.5) state4 <- c("0", "1", "2", "...") cmat4 <- matrix(0, 4,4, dimnames= list(state4, state4)) cmat4[1,2] <- cmat4[2,3] <- cmat4[3,4] <- 1 statefig(c(1,1,1,1), cmat4, bcol=c(1,1,1,0), cex=c(1.5, 1.5, 1.5, 3)) par(oldpar) @ \caption{Four multi-state models. The upper left panel depicts simple survival, the upper right depicts sequential events, the lower left is an example of competing risks, and the lower right panel is a multi-state illness-death model.} \label{mfig1} \end{figure} A multi-state model is used to model a process where subjects transition from one state to the next. For instance, a standard survival curve can be thought of as a simple multi-state model with two states (alive and dead) and one transition between those two states. A diagram illustrating this process is shown in the top left corner of figure \ref{mfig1}. In these types of diagrams, each box is a state and each arrow is a possible transition. The lower left diagram depicts a classic competing risk analysis, where all subjects start on the left and each subject can make a single transition to one of 3 terminal states. The bottom right diagram shows a common multi-state situation known as the illness-death model with recovery. Finally, the upper right diagram represents sequential events, such as repeated infections in the CGD study. In that case one subject had 8 events so there are 9 states corresponding to entry into the study (0 infections) and the first, second, \ldots, eighth events. As will be shown below, there are often multiple choices for the state and transition diagram, and for some data sets it is revealing to look at a problem from multiple views. In addition to deciding the diagram that best matches the research questions, the two other primary decisions are the choice of time scale for the fits, e.g., time from entry in the study vs. time from entry in the state, and what covariates will be used. \section{Multi-state curves} \subsection{Aalen-Johansen estimate} As a starting point for the analysis, it is important to compute and plot estimates of $p(t)$, which is a vector containing the probability of being in each of the states at time $t$. If there is no censoring then $p$ becomes a simple tabulation at time $t$ of all the states. For the general case, we compute this using the Aalen-Johansen estimate via the \code{survfit} function. Mathematically the estimate is simple. For each unique time where an event occurs, form the transition matrix $T(t)$ with elements or rates of $\lambda_{ij}(t) =$ the fraction of subjects who transition from state $i$ to $j$ at time $t$, among those in state $i$ just prior to time $t$. ($T$ is equal to the identity matrix at any time point without an observed transition.) Then \begin{equation} p(t) = p(0) \prod_{s \le t} T(s) \label{AJ} \end{equation} where $p(0)$ is the initial distribution of subjects. Let's work this out for the simple two-state alive $\rightarrow$ {death} model. Let $n(t)$ be the number of subjects still at risk at time $t$ and $d(t)$ the number of deaths at that time point. All subjects start in the alive state and thus $p(0) = (1,0)$ and the transition matrix is \begin{equation*} T(s) = \left( \begin{array}{cc} \frac{n(s)- d(s)}{n(s)} & \frac{d(s)}{n(s)} \\\\ 0 & 1 \end{array} \right) \end{equation*} The two rows are ``start in state 1 (alive)'' and ``start in state 2'' and the two colums are ``finish in state 1'', ``finish in state 2 (death)''. The second row corresponds to the fact that death is an absorbing state: the probability of a death $\rightarrow$ alive transition (lower left element) is 0. Writing out the matrices for the first few transitions and multiplying them leads to \begin{equation} p_1(t) = \prod_{s \le t} \left[n(s) - d(s)\right] /n(s) \label{km} \end{equation} which we recognize as the Kaplan-Meier estimate of survival. For the two state alive-dead model the Aalen-Johansen (AJ) estimate has reprised the KM. In the competing risks case $p(t)$ has an alternate form known as the \emph{cumulative incidence} (CI) function \begin{equation} CI_k(t) = \int_0^t \lambda_k(u) S(u) du \label{cuminc} \end{equation} where $\lambda_k$ is the incidence function for outcome $k$, and $S$ is the overall survival curve for ``time to any endpoint''. Repeating the same matrix exercise for the competing risks, i.e. writing out the Aalen-Johansen computation, exactly recovers the CI formula. The CI is also a special case of the Aalen-Johansen. (The label ``cumulative incidence'' is one of the more unfortunate ones in the survival lexicon, since we normally use `incidence' and `hazard' as interchangeable synonyms but the CI is \emph{not} a cumulative hazard.) The AJ estimate is very flexible; subjects can visit multiple states during the course of a study, subjects can start after time 0 (delayed entry), and they can start in any of the states. The \code{survfit} function implements the AJ estimate and will handle all these cases. The standard error of the estimates is computed using an infinitesimal jackknife. Let $D(t)$ be a matrix with one row per subject and one column per state. Each row contains the \emph{change} in $p(t)$ corresponding to subject $i$, i.e., the derivative of $p$ with respect to the $i$th subject's case weight $dp(t)/dw_i$. Then $V(t) = D'WD$ is the estimated variance-covariance matrix of the estimates at time $t$, where $W$ is a diagonal matrix of observation weights. If a single subject is represented by multiple rows in the data set, then $D$ is first collapsed to have one row per subject, the new row for subject $i$ is the sum of the rows for the observations that represented the subject. This is essentially the same algorithm as the robust variance for a Cox model. For simple two state alive -> dead model without case weights, the IJ estimate of variance is identical to the traditional Greenwood estimate for the variance of the survival curve $S$. (This was a surprise when we first observed it; proving the equivalence is not straightforward.) The $p(t)$ vector obeys the obvious constraint that its sum at any time is equal to one; i.e., each person has to be somewhere. We will use the phrase \emph{probability in state} or simply $p$ when referring to the vector. In the simple two state model Pr(alive) is the usual KM survival estimate, and we have $p_1(t) = 1- p_2(t)$, Pr(alive) = 1 - Pr(dead). Plots for the 2 state case sometimes choose to show Pr(alive) and sometimes Pr(dead). Which one is used often depends on a historical whim of the disease specialty; cardiology journals for instance quite often use Pr(event) resulting in curves that rise starting from zero, while oncology journals invariably use Pr(alive) giving curves that fall downhill from 1. The \code{survfit} routine's historical default for the 2 state case is to print and plot Pr(alive)= $p_1(t)$, which reflects that the author of the routine was working primarily in cancer trials at the time said default was chosen. For simple survival we have gotten used to the idea of using Pr(dead) and Pr(alive) interchangeably, but that habit needs to be left behind for multi-state models, as curves of $1-p_k(t)$ = probability(any other state than $k$) are not useful. In the multi-state case, some curves will rise and then fall. For competing risks the curve for the initial state (leftmost in the diagram) is rarely included in the final plot. Since the curves sum to 1, the full set is redundant. Pr(nothing yet) is usually the least interesting of the set and so it is left off to make the plot less busy. The remaining curves in the competing risks case rise from 0. (This bothers some researchers as it `just looks wrong' to them.) \subsection{Examples} \begin{figure} <>= par(mar=c(.1, .1, .1, .1)) frame() oldpar <- par(usr=c(0,100,0,100)) # first figure xx <- c(0, 10, 10, 0) yy <- c(0, 0, 10, 10) polygon(xx +10, yy+70) temp <- c(60, 80) for (j in 1:2) { polygon(xx + 30, yy+ temp[j]) arrows(22, 70 + 3*j, 28, temp[j] +5, length=.1) } text(c(15, 35, 35), c(75, 65, 85),c("Entry", "Death", "PCM")) text(25, 55, "Competing Risk") # Second figure polygon(xx +60, yy+70) for (j in 1:2) { polygon(xx + 80, yy+ temp[j]) arrows(72, 70+ 3*j, 78, temp[j] +5, length=.1) } text(50+ c(15, 35, 35), c(75, 65, 85),c("Entry", "Death", "PCM")) arrows(85, 78, 85, 72, length=.1) text(75, 55, "Multi-state 1") # third figure polygon(xx+10, yy+25) temp <- c(15, 35) for (j in 1:2) { polygon(2*xx +30, yy + temp[j]) arrows(22, 25 + 3*j, 28, temp[j] +5, length=.1) } text(c(15, 40, 40), c(30, 20, 40),c("Entry", "Death w/o PCM", "PCM")) polygon(2*xx + 60, yy+temp[2]) arrows(52, 40, 58, 40, length=.1) text(70, 40, "Death after PCM") text(40, 10, "Multi-state 2") @ \caption{Three models for the MGUS data.} \label{mfig2} \end{figure} Start with a simple competing risks problem as illustrated in the first diagram of figure \ref{mfig2}. The \code{mgus2} data set contains the time to plasma cell malignancy (PCM) and/or death for 1384 subjects diagnosed with monoclonal gammopathy of undetermined significance (MGUS). Survival and progression time are in months. The code below creates an ordinary Kaplan-Meier curve of post-diagnosis survival for these subjects, along with a histogram of age at diagnosis. The mean age at diagnosis is just over 70 years. <>= oldpar <- par(mfrow=c(1,2)) hist(mgus2$age, nclass=30, main='', xlab="Age") with(mgus2, tapply(age, sex, mean)) mfit1 <- survfit(Surv(futime, death) ~ sex, data=mgus2) mfit1 plot(mfit1, col=c(1,2), xscale=12, mark.time=FALSE, lwd=2, xlab="Years post diagnosis", ylab="Survival") legend("topright", c("female", "male"), col=1:2, lwd=2, bty='n') par(oldpar) @ The xscale and yscale arguments to \code{plot.survfit} affect only the axis labels, not the data. Further additions to the plot region such as \code{legend}, \code{lines}, or \code{text} remain in the original scale. This simplifies programmatic additions such as adding another curve to the plot, while making interactive additions such as a legend somewhat less simple. As a second model for these subjects we will use competing risks (CR), where PCM and death without malignancy are the two terminal states, as shown in the upper left of figure \ref{mfig2}. In a CR model we are only interested in the first event for each subject. Formally we are treating progression to a PCM as an \emph{absorbing state}, i.e., one that subjects never exit. We create a variable \code{etime} containing the time to the first of progression, death, or last follow-up along with an event variable that contains the outcome. The starting data set \code{mgus2} has two pairs of variables \code{(ptime, pstat)} that contain the time to progression and \code{(futime, status)} that contain the time to death or last known alive. The code below creates the necessary \code{etime} and \code{event} variables, then computes and plots the competing risks estimate. <>= mgus2$etime <- with(mgus2, ifelse(pstat==0, futime, ptime)) event <- with(mgus2, ifelse(pstat==0, 2*death, 1)) mgus2$event <- factor(event, 0:2, labels=c("censor", "pcm", "death")) table(mgus2$event) mfit2 <- survfit(Surv(etime, event) ~ sex, data=mgus2) print(mfit2, rmean=240, scale=12) mfit2$transitions plot(mfit2, col=c(1,2,1,2), lty=c(2,2,1,1), mark.time=FALSE, lwd=2, xscale=12, xlab="Years post diagnosis", ylab="Probability in State") legend(240, .6, c("death:female", "death:male", "pcm:female", "pcm:male"), col=c(1,2,1,2), lty=c(1,1,2,2), lwd=2, bty='n') @ The \code{mfit2} call is nearly identical to that for an ordinary Kaplan-Meier, with the exception of the \code{event} variable. \begin{enumerate} \item The event variable was created as a \emph{factor}, whereas for ordinary single state survival the status is either 0/1 or TRUE/FALSE. The first level of the factor must be censoring, which is the status code for those whose follow-up terminated without reaching a new endpoint. Codes for the remaining states can be in any order. The labels for the states are unrestricted, e.g., the first one does not have to be ``censor''. (It will be treated as `no event at this time', whatever the label.) \item A simple print of the \code{mfit2} object shows the order in which the curves will be displayed. This information was used to choose the line types and colors for the curves. \item The \code{mfit2} object contains curves for all the states, but by default the entry state will not be plotted. The remaining curves all start at 0. \item The transitions component of the result is useful as a data check, e.g., if it showed a transition from death to PCM. \item Each subject's initial state can be specified by the \code{istate} argument. When this is omitted all subjects are assumed to start from an entry state named '(s0)'. This default name is an amalgam of ``state 0'', a label sometimes used in textbooks for the leftmost box, and an R convention of using () for constructed names, e.g., `(Intercept)' from an lm() call. \end{enumerate} The printout shows that a male subject will spend, on average, 8.7 of his first 20 years post diagnosis in the entry state, 1.1 years in the PCM state and 10.3 of those 20 in the death state. If a cutoff time is not given the default is to use the maximum observed time for all curves, which is 424 months in this case. The result of a multi-state \code{survfit} is a matrix of probabilities with one row per time and one column per state. These are contained in the probability-in-state (\code{pstate}) component of the resulting survfit object. Since the three MGUS states of entry/pcm/death must sum to 1 at any given time (everyone has to be somewhere), one of the three curves is redundant and the ``fraction still in the entry state'' curve is normally the least interesting. By default, any state with the label `(s0)' is left off of the plot; which is the default value of the \code{noplot} option of \code{plot.survfit}. One can easily plot all of the states by setting the option to NULL. A common mistake with competing risks is to use the Kaplan-Meier separately on each event type while treating other event types as censored. The next plot is an example of this for the PCM endpoint. <>= pcmbad <- survfit(Surv(etime, pstat) ~ sex, data=mgus2) plot(pcmbad[2], lwd=2, fun="event", conf.int=FALSE, xscale=12, xlab="Years post diagnosis", ylab="Fraction with PCM") lines(mfit2[2,"pcm"], lty=2, lwd=2, mark.time=FALSE, conf.int=FALSE) legend(0, .25, c("Males, PCM, incorrect curve", "Males, PCM, competing risk"), col=1, lwd=2, lty=c(1,2), bty='n') @ There are two problems with the \code{pcmbad} fit. The first is that it attempts to estimate the expected occurrence of plasma cell malignancy (PCM) if all other causes of death were to be disallowed. In this hypothetical world it is indeed true that many more subjects would progress to PCM (the incorrect curve is higher), but it is also not a world that any of us will ever inhabit. This author views the result in much the same light as discussions of survival after the zombie apocalypse. The second problem is that the computation for this hypothetical case is only correct if all of the competing endpoints are independent, a situation which is almost never true. We thus have an unreliable estimate of an uninteresting quantity. The competing risks curve, on the other hand, estimates the fraction of MGUS subjects who \emph{will experience} PCM, a quantity sometimes known as the lifetime risk, and one which is actually observable. The last example chose to plot only a subset of the curves, something that is often desirable in competing risks problems to avoid a ``tangle of yarn'' plot that simply has too many elements. This is done by subscripting the \code{survfit} object. For subscripting, multi-state curves behave as a matrix with the outcomes as the second subscript. The columns are in order of the levels of \code{event}, i.e., as displayed by our earlier call to \code{table(event)}. The first subscript indexes the groups formed by the right hand side of the model formula, and will be in the same order as simple survival curves. Thus \code{mfit2[2,1]} corresponds to males (2) and the PCM endpoint (1). Curves are listed and plotted in the usual matrix order of R. A third example using the MGUS data treats it as a multi-state model and it shown in the upper right of figure \ref{mfig2}. In this version a subject can have multiple transitions and thus multiple rows in the data set. In this case it is necessary to identify which data rows go with which subject via the \code{id} argument of \code{survfit}; valid estimates of the curves and their standard errors both depend on this. Our model looks like the illness-death model of figure \ref{mfig1} but with ``PCM'' as the upper state and no arrow for a return from that state to health. The necessary data set will have two rows for any subject who has further follow-up after a PCM and a single row for all others. The data set is created below using the \code{tmerge} function, which is discussed in detail in another vignette. We need to decide what to do with the 9 subjects who have PCM and death declared at the same month. (Some of these were cancer cases discovered at autopsy.) They slipped through without comment in the earlier competing risks analysis; only when setting up this second data set did we notice the ties. Looking back at the code, the prior example counted these subjects as a progression. In retrospect this is a defensible choice: even though undetected before death, the disease must have been present for some amount of time previous and so progression did occur first. For the multi-state model we need to be explicit in how this is coded since a sojourn time of 0 within a state is not allowed. Below we push the progression time back by .1 month when there is a tie, but that amount is entirely arbitrary. Since there are 3 possible transitions we will call the data set \code{data3}. <>= ptemp <- with(mgus2, ifelse(ptime==futime & pstat==1, ptime-.1, ptime)) data3 <- tmerge(mgus2, mgus2, id=id, death=event(futime, death), pcm = event(ptemp, pstat)) data3 <- tmerge(data3, data3, id, enum=cumtdc(tstart)) with(data3, table(death, pcm)) @ The table above shows that there are no observations in \code{data3} that have both a PCM and death, i.e., the ties have been resolved. The last \code{tmerge} line above creates a variable \code{enum} which simply counts rows for each person, which is often useful. <>= temp <- with(data3, ifelse(death==1, 2, pcm)) data3$event <- factor(temp, 0:2, labels=c("censor", "pcm", "death")) mfit3 <- survfit(Surv(tstart, tstop, event) ~ sex, data=data3, id=id) print(mfit3, rmean=240, digits=2) mfit3$transitions plot(mfit3[,"pcm"], mark.time=FALSE, col=1:2, lty=1:2, lwd=2, xscale=12, xlab="Years post MGUS diagnosis", ylab="Fraction in the PCM state") legend(40, .4, c("female", "male"), lty=1:2, col=1:2, lwd=2, bty='n') @ This plot is quite different in that it shows the fraction of subjects \emph{currently} in the PCM state. Looking at our multi-state diagram this is the fraction of subjects in the upper right PCM box. The curve goes up whenever someone enters the box (progression) and down when they leave (death). Myeloma survival was quite short during the era of this study --- most subjects died within 1 year of PCM --- and the proportion currently in the PCM state rarely rises above 2 percent. The result of \code{print(mfit3)} reveals, as expected, less time spent in the PCM state. In the prior \code{mfit2} model, subjects who enter that state remain there for the duration; in this one they quickly pass through. It is worthwhile to check the \code{transitions} table in the output, simply as a data check, since an error in creating the input data can lead to surprising counts and an even more surprising curve. In this case it shows subjects going from the entry state to PCM and death along with transitions from PCM to death. This is as expected. We have often found the three curve display shown below useful in the case of a transient state. It combines the results from competing risk model used above along with a second fit that treats death after PCM as a separate state from death before progression, the \emph{multi-state 2} model of figure \ref{mfig2}. In this plot the fraction of subjects currently in the PCM state is shown by the distance between the two curves. Only males are shown in the plot to minimize overlap. <>= # Death after PCM will correspond to data rows with # enum = 2 and event = death d2 <- with(data3, ifelse(enum==2 & event=='death', 4, as.numeric(event))) e2 <- factor(d2, labels=c("censor", "pcm", "death w/o pcm", "death after pcm")) mfit4 <- survfit(Surv(tstart, tstop, e2) ~ sex, data=data3, id=id) plot(mfit2[2,], lty=c(1,2), xscale=12, mark.time=FALSE, lwd=2, xlab="Years post diagnosis", ylab="Probability in State") lines(mfit4[2,4], mark.time=FALSE, col=2, lty=1, lwd=2, conf.int=FALSE) legend(200, .5, c("Death w/o PCM", "ever PCM", "Death after PCM"), col=c(1,1,2), lty=c(2,1,1), lwd=2, bty='n', cex=.82) @ \subsection{Further notes} The Aalen-Johansen method used by \code{survfit} does not account for interval censoring, also known as panel data, where a subject's current state is recorded at some fixed time such as a medical center visit but the actual times of transitions are unknown. Such data requires further assumptions about the transition process in order to model the outcomes and has a more complex likelihood. The \code{msm} package, for instance, deals with data of this type. If subjects reliably come in at regular intervals then the difference between the two results can be small, e.g., the \code{msm} routine estimates time until progression \emph{occurred} whereas \code{survfit} estimates time until progression was \emph{observed}. \begin{itemize} \item When using multi-state data to create Aalen-Johansen estimates, individuals are not allowed to have gaps in the middle of their time line. An example of this would be a data set with (0, 30, pcm] and (50,70, death] as the two observations for a subject where the time from 30-70 is not accounted for. \item Subjects must stay in the same group over their entire observation time, i.e., variables on the right hand side of the equation cannot be time-dependent. \item A transition to the same state is allowed, e.g., observations of (0,50, 1], (50, 75, 3], (75, 89, 4], (89, 93, 4] and (93, 100, 4] for a subject who goes from entry to state 1, then to state 3, and finally to state 4. However, a warning message is issued for the data set in this case, since stuttering may instead be the result of a coding mistake. The same result is obtained if the last three observations were collapsed to a single row of (75, 100, 4]. \end{itemize} \section{Rate models} For simple two-state survival, the Cox model leads to three relationships \begin{align} \lambda(t) &= \lambda_0(t) e^{X\beta} \label{hazard} \\ \Lambda(t) &= \Lambda_0(t) e^{X\beta} \label{cumhaz}\\ S(t) &= \exp(-\Lambda(t)) \label{surv} \end{align} where $\lambda$, $\Lambda$ and $S$ are the hazard, cumulative hazard and survival functions, respectively. There is a single linear predictor which governs both the rate $\lambda$ (the arrow in figure \ref{mfig1}) and probability of residing in the left hand box of the figure. For multi-state models this simplicity no longer holds; proportional hazards does not lead to proportional $p(t)$ curves. There is a fundamental dichotomy in the analysis namely that \begin{itemize} \item hazards can be computed one at a time, \item probability in state must be done for all states at once. \end{itemize} The analysis of multi-state data has four key steps. In order of importance: \begin{enumerate} \item Draw a box and arrow figure describing the model. \item Think through the rates (arrows). \begin{enumerate} \item Which covariates should be attached to each rate? Sometimes a covariate is important for one transition, but not for another. \item For which transitions should one or more of the covariates be constrained to have the same coefficient? Sometimes there will be a biologic rationale for this. For other studies an equivalence is forced simply because we have too many unknowns and cannot accommodate them all. (This is the same reason that models often contain very few interaction terms). \item Which, if any, of the transitions should share the same baseline hazard? Most of the time the baseline rates are all assumed to be different. \item Should there be random effects, and if so what is an appropriate correlation structure? Do some pairs of transitions have a shared effect, some pairs separate effects and others no random effect? Mixed effects Cox models tend to need larger sample size --- does the data set have enough events? \end{enumerate} \item Build an appropriate data set. \item Fit the data. Examine multiple summaries of the model fit, including the predicted occupancy curves. \end{enumerate} Step 1 is key to the entire endeavor. We saw in figure \ref{mfig2} and the examples above that multiple views of a multi-state process can be useful, and this will hold for modeling as well. Step 3 will often be the one that demands the most attention to detail. \subsection{MGUS example} Start with the simplest model for the MGUS data: a competing risks model (upper left diagram of figure \ref{mfig2}), distinct baseline hazards for the two rates, no shared coefficients, and three covariates. <>= options(show.signif.stars = FALSE) # display intelligence cfit1 <- coxph(Surv(etime, event) ~ age + sex + mspike, mgus2, id=id) cfit1$cmap print(cfit1, digits=2) @ This above call fits both endpoints at once. The resulting R object contains a single vector of coefficients, the \code{cmap} component documents which goes with which endpoint and is used by the print routine to organize the printout. We see that the effect of age and sex on non-PCM mortality is profound, which is not a surprise given the median starting age of \Sexpr{median(mgus2$age)}. Risk rises \Sexpr{round(exp(10*coef(cfit1)[4]),1)} fold per decade of age and the death rate for males is \Sexpr{round(exp(coef(cfit1)[5]),1)} times as great as that for females. The size of the serum monoclonal spike has almost no impact on non-PCM mortality. A 1 unit increase changes mortality by only 6\%. The mspike size has a major impact on progression, however; each 1 gram change increases risk by \Sexpr{round(exp(coef(cfit1)[3]) ,1)} fold. The interquartile range of \code{mspike} is 0.9 gram so this risk increase is clinically important. The effect of age on the progression rate is much less pronounced, with a coefficient only 1/4 that for mortality, while the effect of sex on progression is completely negligible. The effect of sex on the \emph{lifetime} probability of PCM is not zero, however. Because females live longer, a female with MGUS will on average spend more total years at risk for PCM than the average male, and so has a larger lifetime risk of PCM. The average rate of progression is about 1\% per year, as shown below, while the mean post diagnosis lifetime is 19 months longer for females. The overall effect is a 1.6\% increase in lifetime risk. <>= pfit1 <- pyears(Surv(ptime, pstat) ~ sex, mgus2, scale=12) round(100* pfit1$event/pfit1$pyears, 1) # PCM rate per year temp <- summary(mfit1, rmean="common") #print the mean survival time round(temp$table[,1:6], 1) @ The same Cox model coefficients can also be obtained by fitting separate models for the PCM and death endpoint, censoring cases that fail due to the other cause. Hazards can be computed one at a time. When computing $p(t)$, on the other hand, all the rates must be considered at once, it is necessary to use the joint fit \code{cfit1} above. We create predicted curves for four hypothetical subjects below. <>= dummy <- expand.grid(sex=c("F", "M"), age=c(60, 80), mspike=1.2) dummy csurv <- survfit(cfit1, newdata=dummy) dim(csurv) @ The resulting set of Aalen-Johansen estimates can be indexed as though it were an array. There is only one stratum (the Cox model did not have strata), four hypothetical subjects, and 3 states. (When there is only a single stratum, as here, the subscript method allows that index to be omitted, for backwards compatability). <>= plot(csurv[,'pcm'], xmax=25*12, xscale=12, xlab="Years after MGUS diagnosis", ylab="PCM", col=1:2, lty=c(1,1,2,2), lwd=2) legend(10, .14, outer(c("female", "male "), c("diagnosis at age 60", "diagnosis at age 80"), paste, sep=", "), col=1:2, lty=c(1,1,2,2), bty='n', lwd=2) @ The individual survival curves that would result from fitting each endpoint separately are not actually of interest, since each will be a Cox model analog of the pcmbad curve we criticized earlier. The coxph fits showed that sex has nearly no effect on the hazard of PCM, i.e., on any given day the risk of conversion for a male is essentially the same as for a female of the same age. Yet we see above that the fitted Cox models predict a higher lifetime risk for females; this is due to a longer lifespan. The age effect on lifetime risk that is far from proportional: the 80 year curves flatten out while predictions for the younger subjects do not. Very few subjects acquire PCM more than 15 years after a MGUS diagnosis at age 80 for the obvious reason: very few of them will still be alive. <>= # Print out a M/F results at 20 years temp <- summary(csurv, time=20*12)$pstate cbind(dummy, PCM= round(100*temp[,,2], 1)) @ The above table shows that females are predicted to have a higher risk of 20 year progression, even though their hazard at any given moment is nearly identical to males. The female/male difference at 20 years is on the order of our ``back of the napkin'' person-years estimate of 1\% progression per year * 1.7 more years of life for the females, but the progression fraction varies substantially by group. \section{Fine-Gray model} For the competing risk case the Fine-Gray model provides an alternate way of looking at the data. As we saw above, the impact of a particular covariate on the final values $P$ can be complex, even if the models for the hazards are relatively simple. The primary idea of the Fine-Gray approach is to turn the multi-state problem into a collection of two-state ones. In the upper right diagram of figure \ref{mfig2}, draw a circle around all of the states except the chosen endpoint and collapse them into a single meta-state. For the MGUS data these are \begin{itemize} \item Model 1 \begin{itemize} \item left box: All subjects in the entry or ``death first'' state \item right box: PCM \end{itemize} \item Model 2 \begin{itemize} \item left box: All subjects in the entry or ``PCM first'' state \item right box: Death (without PCM) \end{itemize} \end{itemize} An interesting aspect of this is that the fit can be done as a two stage process: the first stage creates a special data set while the second fits a weighted \code{coxph} or \code{survfit} model to the data. The data set can be created using the \code{finegray} command. <>= # (first three lines are identical to an earlier section) etime <- with(mgus2, ifelse(pstat==0, futime, ptime)) event <- with(mgus2, ifelse(pstat==0, 2*death, 1)) event <- factor(event, 0:2, labels=c("censor", "pcm", "death")) pcmdat <- finegray(Surv(etime, event) ~ ., data=mgus2, etype="pcm") pcmdat[1:4, c(1:3, 11:14)] deathdat <- finegray(Surv(etime, event) ~ ., data=mgus2, etype="death") dim(pcmdat) dim(deathdat) dim(mgus2) @ The \code{finegray} command has been used to create two data sets: one for the PCM endpoint and one for the death endpoint. In each, four new variables have been added containing a survival time \code{(fgstart, fgstop, fgstatus)} with an `ordinary' status of 0/1, along with a case weight and a large number of new rows. We can use this new data set as yet another way to compute multi-state survival curves, though there is no good reason to use this rather roundabout approach instead of the simpler \code{survfit(Surv(etime, event) \textasciitilde sex)}. <>= # The PCM curves of the multi-state model pfit2 <- survfit(Surv(fgstart, fgstop, fgstatus) ~ sex, data=pcmdat, weight=fgwt) # The death curves of the multi-state model dfit2 <- survfit(Surv(fgstart, fgstop, fgstatus) ~ sex, data=deathdat, weight=fgwt) @ Inspection shows that the two new curves are almost identical to the prior estimates based on the Aalen-Johansen estimate (\code{mfit2}), and in fact would be identical if we had accounted for the slightly different censoring patterns in males and females (by adding \code{strata(sex)} to the right hand side of the \code{finegray} formulas). A Cox model fit to the constructed data set yields the Fine-Gray models for PCM and for death. <>= fgfit1 <- coxph(Surv(fgstart, fgstop, fgstatus) ~ sex, data=pcmdat, weight= fgwt) summary(fgfit1) fgfit2 <- coxph(Surv(fgstart, fgstop, fgstatus) ~ sex, data=deathdat, weight= fgwt) fgfit2 mfit2 <- survfit(Surv(etime, event) ~ sex, data=mgus2) #reprise the AJ plot(mfit2[,'pcm'], col=1:2, lwd=2, xscale=12, conf.times=c(5, 15, 25)*12, xlab="Years post diagnosis", ylab="Fraction with PCM") ndata <- data.frame(sex=c("F", "M")) fgsurv1 <- survfit(fgfit1, ndata) lines(fgsurv1, fun="event", lty=2, lwd=2, col=1:2) legend("topleft", c("Female, Aalen-Johansen", "Male, Aalen-Johansen", "Female, Fine-Gray", "Male, Fine-Gray"), col=1:2, lty=c(1,1,2,2), bty='n') # rate models with only sex cfitr <- coxph(Surv(etime, event) ~ sex, mgus2, id=id) rcurve <- survfit(cfitr, newdata=ndata) #lines(rcurve[, 'pcm'], col=6:7) # makes the plot too crowsded @ The FG model states that males have a less \emph{observed} PCM, by a factor of \Sexpr{round(exp(coef(fgfit1)), 2)}, and that this hazard ratio is constant over time. An overlaid plot of the non-parametric Aalen-Johansen estimates for the PCM state (from \code{survfit}) along with predicted curves from the Fine-Gray model show that proportional hazards is not unreasonable for this particular fit. The predicted values from the rate model, computed just above but not plotted on the curve, also fit well with the data. When there is only a single categorical 0/1 covariate the Fine-Gray model reduces to Gray's test of the subdistribution function, in the same way that a \code{coxph} fit with a single categorical predictor is equivalent to the log-rank test. The mathematics behind the Fine-Gray estimate starts with the functions $F_k(t) = p_k(t)$, where $p$ is the probability in state function estimated by the AJ estimate. This can be thought of as the distribution function for the improper random variable $T^*= I(\mbox{event type}=k)T + I(\mbox{event type}\ne k)\infty$. Fine and Gray refer to $F_k$ as a subdistribution function. In analogy to the survival probability in the two state model define \begin{equation} \gamma_k(t) = - d \log[1-F_k(t)]/dt \label{FG} \end{equation} and assume that $\gamma_k(t;x) = \gamma_{k0}(t) \exp(X\beta)$. In a 2 state alive $\longrightarrow$ death model, $\gamma$ becomes the usual hazard function $\lambda$. In the same way that our multivariate Cox model \code{cfit1} made the simplifying assumption that the impact of male sex is to increase the hazard for death by a factor of \Sexpr{round(exp(coef(cfit1)[5]), 2)}, independent of the subject's age or serum mspike value, the Fine-Gray model assumes that each covariate's effect on $\log(1-F)$ is a constant, independent of other variables. Both model's assumptions are wonderfully simplifying with respect to understanding a covariate, since we can think about each covariate separately from all the others. This is, of course, under the assumption that the model is correct: that additivity across covariates, linearity, and proportional hazards all hold. In a multi-state model, however, these assumptions cannot hold for both the per-transition and Fine-Gray models formulations at the same time; if true for one, they will not be true for the other. Now consider a multivariate fit on age, sex, and serum m-spike. <>= fgfit2a <- coxph(Surv(fgstart, fgstop, fgstatus) ~ age + sex + mspike, data=pcmdat, weight=fgwt) fgfit2b <- coxph(Surv(fgstart, fgstop, fgstatus) ~ age + sex + mspike, data=deathdat, weight=fgwt) round(rbind(PCM= coef(fgfit2a), death=coef(fgfit2b)), 3) @ The Fine-Gray fits show an effect of all three variables on the subdistribution rates. Males have a lower lifetime risk of PCM before death and a higher risk of death before PCM, while a high serum m-spike works in the opposite direction. The Cox models showed no effect of sex on the instantaneous hazard of PCM and none for serum m-spike on the death rate. However, as shown in the last section, the Cox models do predict a greater lifetime risk for females. We had also seen that older subjects are less likely to experience PCM due to the competing risk of death; this is reflected in the FG model as a negative coefficient for age. Now compute predicted survival curves for the model, and show them alongside the predictions from the multi-state Cox model. <>= oldpar <- par(mfrow=c(1,2)) dummy <- expand.grid(sex= c("F", "M"), age=c(60, 80), mspike=1.2) fsurv1 <- survfit(fgfit2a, dummy) # time to progression curves plot(fsurv1, xscale=12, col=1:2, lty=c(1,1,2,2), lwd=2, fun='event', xlab="Years", ylab="Fine-Gray predicted", xmax=12*25, ylim=c(0, .15)) legend(1, .15, c("Female, 60", "Male, 60","Female: 80", "Male, 80"), col=c(1,2,1,2), lty=c(1,1,2,2), lwd=2, bty='n') plot(csurv[,2], xscale=12, col=1:2, lty=c(1,1,2,2), lwd=2, xlab="Years", ylab="Multi-state predicted", xmax=12*25, ylim=c(0, .15)) legend(1, .15, c("Female, 60", "Male, 60","Female: 80", "Male, 80"), col=c(1,2,1,2), lty=c(1,1,2,2), lwd=2, bty='n') par(oldpar) @ The predictions as a function of age group are quite different for the Fine-Gray model: new PCM cases are predicted 20+ years after diagnosis in both the old and young age groups, while they are predicted to cease in the multi-state fit. The average of all four curves is nearly the same at each age, but the global proportional hazards assumption of the FG model forces the curves to remain parallel. We can check the proportional hazards assumption of the models using the \code{cox.zph} function, linearity of the continuous variables age and mspike by using non-linear terms such as \code{pspline} or \code{ns}, and additivity by exploring interactions. All are obvious and important next steps. For instance, the proportional hazards assumption for age shows clear violations. <>= zph.fgfit2a <- cox.zph(fgfit2a) zph.fgfit2a plot(zph.fgfit2a[1]) abline(h=coef(fgfit2a)[1], lty=2, col=2) @ A further weakness of the Fine-Gray approach is that since the two endpoints are modeled separately, the results do not have to be consistent. Below is a graph of the predicted fraction who have experienced neither endpoint. For subjects diagnosed at age 80 the Fine-Gray models predict that more than 100\% will either progress or die by 30 years. Predictions based on the Aalen-Johansen approach do not have this issue. <>= fsurv2 <- survfit(fgfit2b, dummy) # time to progression curves xtime <- 0:(30*12) #30 years y1a <- 1 - summary(fsurv1, times=xtime)$surv #predicted pcm y1b <- 1 - summary(fsurv2, times=xtime)$surv #predicted deaths before pcm y1 <- (y1a + y1b) #either matplot(xtime/12, y1, col=1:2, lty=c(1,1,2,2), type='l', xlab="Years post diagnosis", ylab="FG: either endpoint") abline(h=1, col=3) legend("bottomright", c("Female, 60", "Male, 60","Female: 80", "Male, 80"), col=c(1,2,1,2), lty=c(1,1,2,2), lwd=2, bty='n') @ The primary strength of the Fine-Gray model with respect to the Cox model approach is that if lifetime risk is a primary question, then the model has given us a simple and digestible answer to that question: ``females have a 1.2 fold higher lifetime risk of PCM, after adjustment for age and serum m-spike''. This simplicity is not without a price, however, and these authors are not proponents of the approach. There are five issues. \begin{enumerate} \item The attempt to capture a complex process as a single value is grasping for a simplicity that does not exist for many (perhaps most) data sets. The necessary assumptions in a multivariate Cox model of proportional hazards, linearity of continuous variables, and no interactions are strong ones. For the FG model these need to hold for a combined process --- the mixture of transition rates to each endpoint --- which turns out to be a more difficult barrier. \item The sum of predictions need not be consistent. \item From the per-transition Cox model one can work forward and compute $p(t)$, the occupancy probabilities for each state over time; both the hazard ratios and $p$ are useful summaries of the data. We don't have tools to work backwards from a Fine-Gray fit to the per transition hazards. \item The approach is viable only for competing risks and not for other multi-state models. \item The risk sets are odd. \end{enumerate} The last of these is perhaps the most frequently listed issue with the Fine-Gray model, but it is actually a minor complaint. The state probabilities $p(t)$ in a multi-state model are implicitly fractions of the total population we started with: someone who dies in month 1 is still a part of the denominator for the fraction of subjects with PCM at 20 years. In the Fine-Gray formulas this subject explicitly appears in risk set denominators at a later time, which looks odd but is more of an artifact. The first issue is substantial, however, and checking the model assumptions of a Fine-Gray fit is mandatory. The second point is alarming, but it normally does not have a practical impact unless there is long follow-up. \section{Shared coefficients} To fit risk models that have shared coefficients or baseline hazards we use an extended formula notation for the \code{coxph} function. For the simple competing risks MGUS fit above, assume that we wanted to add hemoglobin to the fit, with a common coefficient for both the PCM and death endpoint. (Anemia is a feature of both PCM and old age.) <>= hgfit <- coxph(list(Surv(etime, event) ~ age + sex + mspike, 1:2 + 1:3 ~ hgb / common), data = mgus2, id = id) hgfit @ Notice that the hemoglobin coefficient is identical for the PCM and death models. A closer look at the result shows that the coefficient vector is of length 7 and not 8, orchestrated by the \code{cmap} component. The variance matrix likewise is 7 by 7. <>= hgfit$coef hgfit$cmap @ In a \code{coxph} formula list, the first element is the default formula which will be applied to all the transitions. The second, third, etc. elements are pseudo-formulas that have a set of transitions as the ``response'', and usual covariate formulas on the right. For instance, if we wanted creatinine to be a covariate for only the transition to PCM a third line of \code{1:2 \textasciitilde creat} could be added. There is flexibility in how the left hand side is written. The following are all valid forms for the second line of our formula. <>= list(1: c("pcm", "death") ~ hgb / common, 1:0 ~ hgb / common, state("(s0)"): (2:3) ~ hgb / common) @ A '0' is shorthand for ``any state''. Transitions in the formula that do not occur in the data are silently ignored, e.g., the implicit 1:1 transition of our second line above. For the three transition model based on data3, the following fits first a model with 3 baseline hazards and 6 coefficients. The second fit assumes a common baseline hazard for the two transitions to death. <>= cox3a <- coxph(Surv(tstart, tstop, event) ~ age + sex, data3, id=id) cox3b <- coxph(list(Surv(tstart, tstop, event) ~ age + sex, 1:3 + 2:3 ~ 1/ common), data= data3, id= id) cox3b$cmap @ For this data set model \code{cox3b} makes very little medical sense, since death rates after PCM do not even remotely resemble those for patients without a plasma cell malignancy. (The fact that the package allows something does not make it a good idea.) In other data sets a shared baseline for all the transitions to death can be an effective model. Be aware that though shared baselines for two transitions that end in the same state is acceptable, a shared hazard for two transitions that originate in the same state is not compatable with the partial likelihood formulation of a Cox model, and will lead to surprising results (often a large number of NA coefficients). An alternative approach to the above is to fit the joint model using a special data set. To reproduce the \code{hgfit} model above, create a stacked data set with $2n$ observations. The first $n$ rows are the data set we would use for a time to PCM analysis, with a simple 0/1 status variable encoding the PCM outcome. The second $n$ rows are the data set we would have used for the `death before PCM' fits, with status encoding the death-before-PCM endpoint. A last variable, \code{group}, is `pcm' for the first $n$ observations and `death' for the remainder. Then fit a model <>= temp1 <- data.frame(mgus2, time=etime, status=(event=="pcm"), group='pcm') temp2 <- data.frame(mgus2, time=etime, status=(event=="death"), group="death") stacked <- rbind(temp1, temp2) allfit <- coxph(Surv(time, status) ~ hgb + (age + sex+ mspike)*strata(group), data=stacked) all.equal(allfit$loglik, hgfit$loglik) @ This fits a common effect for hemoglobin (hgb) but separate age and sex effects for the two endpoints, along with separate baseline hazards. \section{Other software} \subsection{The mstate package} As the number of states + transitions (arrows + boxes) gets larger then the `by hand' approach used above for creating a stacked data set becomes a challenge. (It is still fairly easy to do, just not as easy to ensure it has been done \emph{correctly}.) The \code{mstate} package starts with a definition of the matrix of possible transitions and uses that to drive tools that build and analyze a stacked data set in a more automated fashion. We find it a little more difficult to use than a \code{coxph} model with a multistate status variable. (The fact that we like our own child best should be no surprise, however). One current disadvantage of the survival package is that the Aalen-Johansen curves from a multi-state coxph model currently do not include a variance estimate, whereas those from mstate do have a variance. A current restriction in R is that packages on the \emph{recommended} list (of which survival is one) should not depend on any packages outside that list. This precludes adding an mstate example to this vignette, even though it is a natural fit. \subsection{The \code{msm} package} There are two broad classes of multi-state data: \begin{itemize} \item Panel data arises when subjects have regular visits, with the current state assessed at each visit. We don't know when the transitions between states occur, or if other states may have been visited in the interim --- only the subject's state at specific times. This is also known as interval censored data. \item Survival data arises when we observe the transition times; death, for example. \end{itemize} The overall model (boxes and arrows), the quantities of interest (transition rates and $p(t)$), and the desired printout and graphs are identical for the two cases. Much of the work in creating a data set is also nearly the same. The underlying likelihood equations and resulting analytical methods for solving the problem are, however, completely different. The \code{msm} package addresses panel data, while \code{survival}, \code{mstate}, and many others are devoted to survival data. \section{Conclusions} When working with acute diseases, such as advanced cancer or end-stage liver disease, there is often a single dominating endpoint. Ordinary single event Kaplan-Meier curves and Cox models are then efficient and sufficient tools for much of the analysis. Such data was the primary use case for survival analysis earlier in the authors' careers. Data with multiple important endpoints is now common, and multi-state methods are an important addition to the statistical toolbox. As shown above, they are now readily available and easy to use. It is sometimes assumed that the presence of competing risks \emph{requires} the use of a Fine-Gray model (we have seen it in referee reports), but this is not correct. The model may often be useful, but is one available option among many. Grasping the big picture for a multi-state data set is always a challenge and we should make use of as many tools as possible. It is not always easy to jump between observed deaths, hazard rates, and lifetime risk. We are often reminded of the story of the gentleman on his 100th birthday who proclaimed that he was looking forward to many more years ahead because ``I read the obituaries every day, and you almost never see someone over 100 listed there''. %\bibliographystyle{plain} %\bibliography{refer} \end{document}