mstate : An R Package for the Analysis of Competing Risks and Multi-State Models

Multi-state models are a very useful tool to answer a wide range of questions in survival analysis that cannot, or only in a more complicated way, be answered by classical models. They are suitable for both biomedical and other applications in which time-to-event variables are analyzed. However, they are still not frequently applied. So far, an important reason for this has been the lack of available software. To overcome this problem, we have developed the mstate package in R for the analysis of multi-state models. The package covers all steps of the analysis of multi-state models, from model building and data preparation to estimation and graphical representation of the results. It can be applied to non-and semi-parametric (Cox) models. The package is also suitable for competing risks models, as they are a special category of multi-state models. This article oﬀers guidelines for the actual use of the software by means of an elaborate multi-state analysis of data describing post-transplant events of patients with blood cancer. The data have been provided by the EBMT (the European Group for Blood and Marrow Transplantation). Special attention will be paid to the modeling of diﬀerent covariate eﬀects (the same for all transitions or transition-speciﬁc) and diﬀerent baseline hazard assumptions (diﬀerent for all transitions or equal for some).


Introduction
Recently, multi-state and competing risks models have gained considerable popularity in survival analysis.In the first place, this popularity is due to the fact that in comparison to classical models, these models describe the disease/recovery process of patients in more detail, thus yielding more insight.In the second place, these models are useful in a clinical setting for the prediction of survival duration for specific patients.These predictions can eas-ily be updated if information about events later occurring to such a patient becomes available.The influence of covariates can be taken into account.
Although the relevance of these models is clearly recognized, their application by non-statisticians has so far been limited.An important reason for this is the lack of available software.For this reason, we have developed a software package in R (R Development Core Team 2010), called mstate, which can be used for the different phases of the description and analysis of competing risks and multi-state models.It is primarily designed for non-and semi-parametric Cox models.However, we have paid special attention to the flexibility of the functions: they can be used independently of each other to enable users to combine their own functions with those of mstate when considering models different from ours.mstate is available from the Comprehensive R Archive Network at http://CRAN.R-project.org/package=mstate.
The mathematics underlying the package, its philosophy and features and an example have been discussed in our previous paper de Wreede et al. (2010).A more general introduction into competing risks and multi-state models can be found in Putter et al. (2007).The current article builds on these previous articles.It describes in detail how the functions work by means of a more elaborate example.The mathematical concepts related to the functions will only be mentioned briefly in the relevant places.The use of the package is illustrated by the analysis of a model describing the disease/recovery process of leukemia patients, of which model several variants will be discussed.The code for all of them will be given, and it will be explained how this should be adapted for other models.Data, clinical questions and model are introduced in Section 2. We consider estimation and prediction both for a non-parametric model (Section 3) and a semi-parametric model including several relevant clinical covariates (Section 4.1).The use of transition-specific covariates will be explained.As a special semi-parametric model, we consider a proportional baseline hazards model, which is useful for reducing the number of parameters in the model (Section 4.2).In Section 4.3 we discuss reduced rank models and two functions for simulation and bootstrapping in the context of multi-state models.Finally, the analysis of competing risks models by means of mstate is discussed briefly in Section 4.4.

Data, questions, and model
We consider survival after a transplant treatment of patients suffering from a blood cancer.The data have been provided by the EBMT (the European Group for Blood and Marrow Transplantation); they are available in mstate as ebmt4.The present data set has been compiled to illustrate the models and the software.To facilitate this illustration, only patients with complete covariate information and a reasonable amount of information about intermediate events have been included.Although the current data set mimics a real-life situation and although the order of magnitude of the outcomes is correct, the clinical meaning of the results of the analyses is restricted.To avoid misinterpretation, we have abstracted from the actual disease, covariate values and intermediate events.
Three intermediate events are included in the model: Recovery (Rec), an Adverse Event (AE) and a combination of the two (AE and Rec).It is to be expected that recovery improves the prognosis and an adverse event deteriorates it.The model is suitable to show the size of these effects, and to capture the influence of their timing and of the covariates on the prognosis.

Death
Figure 1: A six-states model for leukemia patients after bone marrow transplantation Moreover, it shows what happens when both the positive and negative event take place, as compared to one or none of them.The models under study here can easily be made more specific in the case of real applications.For instance, instead of Recovery, Engraftment can be included, and instead of Adverse Event, Acute Graft-versus-Host Disease.
We consider 2279 patients who were treated between 1985 and 1998.Four prognostic factors are known at baseline for all patients (see Table 1).They are: donor-recipient match (where gender mismatch is defined as female donor, male recipient), prophylaxis, year of transplant and age at transplant in years.All these covariates are treated as time-fixed categorical covariates.The distribution of the values of the covariates over the patients in the data set is shown in Table 1.
A multi-state approach is particularly appropriate for these data, since it can help to model both the disease-related and the treatment-related morbidity and mortality.These are modeled here by including the intermediate events recovery and adverse event.Information about the occurrence of these events is used to update the prognosis of the patients.
We consider the following six-states model (see Figure 1): 1. Alive and in remission, no recovery or adverse event; 2. Alive in remission, recovered from the treatment; 3. Alive in remission, occurrence of the adverse event; 4. Alive, both recovered and adverse event occurred; 5. Alive, in relapse (treatment failure); 6. Dead (treatment failure).
All patients start in state 1.States 5 and 6 are called absorbing: once the patient has entered one of them, she/he stays there.This leaves us with a model with 12 transitions.
The data have been made suitable for a multi-state analysis by some small adjustments.Since the model does not allow patients to enter two states at the same time, we have set the death indicator of patients with simultaneous relapse and death to 0, because although these events were reported at the same time, in reality the patients must have experienced the relapse before their death.For those with equal time to the adverse event and time to Rec we have lowered the time to AE by half a day to avoid a transition with only very few events.Two new variables have been created to express the time of entry in state 4 (AE and Rec) and the accompanying status indicator: recae and recae.srespectively.
The data are available in mstate in wide format.This means that each row in the data corresponds to a single subject.
R> library("mstate") R> data("ebmt4  6 1995-1998 20-40 no no gender mismatch The columns rec, ae, rel and srv are time variables, indicating the time measured in days post-transplant to recovery, AE, relapse and death respectively in case of an event, or last follow-up otherwise.The .s-variables are the corresponding status variables (1 for an event, 0 for censoring).For instance, patient 1 had recovered after 22 days (transition from state 1 to state 2) and was censored after 995 days without a further event.Patient 2 experienced the adverse event after 12 days (transition from state 1 to state 3), then recovery after 29 days (transition from state 3 to state 4) and a relapse after 422 days (transition from state 4 to state 5).Finally, he/she died after 579 days, but this last event is not relevant to the model, because the patient had already reached an absorbing state.
Diverse clinical questions can be answered by fitting this model to the data, such as: how does age influence the different phases of the disease/recovery process?How does the occurrence of the adverse event after 2 months change the prognosis of 10-year survival for a patient?Which risks should be monitored most carefully for a patient who has recovered after one month, taking into account certain covariate values?
The model can be described by means of a transition matrix, which is in this case a 6-by-6 matrix.A number at entry (g, h) of the matrix represents a possible transition from state g to state h.These numbers range here from 1 to 12, because the model has 12 transitions.These transition numbers are used in the various stages of the analysis.If a transition between two states is not allowed, the entry becomes NA.In the present format, the data are not yet suitable for a multi-state analysis.First they have to be recoded into 'long format'.In this format, each subject has as many rows as transitions for which he/she is at risk.The function msprep() transforms a data frame in wide format into one in long format.Arguments for msprep() are a time and status vector indicating the time of entry in every state and the accompanying status indicator.The keep argument contains the names of the covariates that will be used in the analysis.The output is an object of class 'msdata': a data frame in long format, which has the transition matrix as an attribute.This object has its own print() method.
The call to msprep() and output for the first patient are as follows (covariate match is not shown in the output): R> msebmt <-msprep(data = ebmt, trans = tmat, time = c(NA, "rec", "ae", + "recae", "rel", "srv"), status = c(NA, "rec.s","ae.s", "recae.s",+ "rel.s","srv.s"),keep = c("match", "proph", "year", "agecl")) R> msebmt[msebmt$id == 1, c(1:8, 10:12)] In Section 4.1, a model will be introduced in which it is assumed that the covariates have different effects on each transition.This can be achieved by creating transition-specific covariates.They are derived from the covariates at baseline as follows: each covariate Z is split up into as many covariates Z gh as there are transitions in the model.For the transition from state g to state h, Z gh is equal to Z; for all other transitions, Z gh = 0.
The function expand.covs()expands the covariates specified by the user on the basis of an object of class 'msdata'.Factors are expanded into dummy variables.The option longnames controls whether the names of the dummy variables are based on the levels of the categorical covariates (TRUE) or on numbers (FALSE).In our example, we obtain six new covariates for each transition: one dummy variable each for donor-recipient match and prophylaxis, and two each for year of transplant and age.For the coding of the dummy variables, consider agecl, which has three values.These are coded by agecl1 and agecl2.For the reference category (1985)(1986)(1987)(1988)(1989) agecl1 and agecl2 are both 0; for 1990-1994, agecl1 = 1 and agecl2 = 0 and for 1995-1998, agecl1 = 0 and agecl2 = 1.This results in 72 new covariates altogether.
We illustrate this process by showing a selection of the long format data for the first patient, leaving out the values at baseline and the transition-specific counterparts of match, proph and agecl.We also omit year1.1 to year1.12, because they are all 0.
Finally, for future modeling we convert time into years.

A non-parametric model
First we consider a non-parametric model, in which we ignore the influence of covariates.The basic quantities of interest are the transition intensities or hazard rates.Denote by g → h a transition from state g to state h, by X(t) the state occupied at time t and by α gh (t) the corresponding transition intensity.The transition intensity expresses the instantaneous risk of a transition from state g into state h at time t.It is defined as (1) This definition makes an implicit assumption that the multi-state model is Markovian, which implies that the probability of going to a future state depends only on the present state and not on the history.The cumulative transition hazard for transition g → h is defined as A gh (t) = t 0 α gh (u)du.Let A(t) be a matrix with dimensions S × S, in which S denotes the number of states.This matrix has elements A gh (t) (g = h), and diagonal elements For the estimation of the cumulative hazard, we assume we have data with independent censoring.The Nelson-Aalen estimator A gh (t) of the cumulative hazard for transition g → h is now given by A gh (t) = where t i indicates the event times, dN gh (t i ) is the observed number of transitions from state g to state h at time t i , and Y g (t i ) is the number of subjects at risk for a transition from state g at t i .The Nelson-Aalen estimator makes jumps of size ∆ A gh (t i ) = dN gh (t i )/Y g (t i ) at the event times t i .The estimates of the diagonal elements of A(t) are given by A gg (t) = − h =g A gh (t).
We use the function coxph() from the survival package (Therneau and Lumley 2010) to estimate the cumulative hazards.
In principle, the transition intensities could also be estimated separately, but the combined use of long format data and a single stratified copxh object makes further calculations easier.
The output of coxph() is the input for mstate's function msfit().It estimates transition hazards and their associated (co)variances.The output of msfit() is an object of class 'msfit'.This is a list with elements Haz (with the estimated cumulative hazard values at all event times), varHaz (with the covariances of each pair of estimated cumulative hazards at each event time point, i.e., cov( A gh (t), A kl (t))), and trans, in which the transition matrix is stored for further use.The (co)variances of the estimated cumulative hazards may be computed in two different ways: by means of the Aalen estimator or by means of the Greenwood estimator.
An advantage of the Greenwood estimator is the fact that it yields exact multinomial standard errors for the transition probabilities when there is no censoring.The two estimators give almost equal results in all practical applications.For details see Section IV R> msf0 <-msfit(object = c0, vartype = "greenwood", trans = tmat) The 'msfit' class has its own summary() and plot() methods.By default, summary() prints head and tail of the cumulative hazards of all transitions.Since there are twelve transitions in the model, its output is not shown here in the interest of space.The head and tail of the Haz and varHaz items look as follows (time is now given in years): R> head(msf0$Haz) time Haz trans 1 0.002737851 0.000000000 1 2 0.008213552 0.000000000 1 3 0.010951403 0.000000000 1 4 0.013689254 0.000000000 1 5 0.016427105 0.000443066 1 6 0.019164956 0.001333142 1 The last time point in the list indicates the last time point in the data, either of an event or of a censoring.The plot() function produces a graph of all estimated cumulative hazards in different colors: R> plot(msf0, las = 1, lty = rep(1:2, c(8, 4)), + xlab = "Years since transplantation") When we choose the "aalen" option for the vartype argument in msfit(), the estimated cumulative hazards are the same as before, but the standard errors are slightly different.
The function probtrans() calculates the estimated transition probabilities, and optionally the standard errors and/or the covariances of the transition probabilities.In the case of non-parametric models, the user can choose between Greenwood or Aalen standard errors in the method argument, in accordance with the choice in msfit().Just as in the case of the estimates of the hazards, the estimates of the transition probabilities themselves do not depend on this choice.
The argument predt gives the starting time for prediction, that is, the starting time for the calculation of the transition probabilities.Two directions of prediction have been implemented, which can be specified by the direction argument: "forward" (the default) and "fixedhorizon".Any string starting with "fo" or "fi" is sufficient to distinguish between the two options.The "forward" option means that the prediction is made from predt; in P(s, t), time s remains fixed at the value predt, while time t varies from s to the last (possibly censored) time point in the data.The "fixedhorizon" option means that the prediction is made for predt; in P(s, t), time t remains fixed at the value predt and time s varies from 0 to predt.The use of the fixed horizon option will be illustrated in Section 4.1.
The default output of probtrans() is an object of class 'probtrans', which is a list of (S + 1) elements, where S is the number of states.Element g of the list contains all transition probabilities starting from state g and going to all states, and their associated standard errors (g = 1, . . ., S).The last element of the list contains again the transition matrix.If the option covariance=TRUE is chosen, an additional element is added to the list containing all estimated covariances of the transition probabilities.The 'probtrans' class has two methods, summary() and plot().By default the summary() method prints the head and tail of the estimated transition probabilities for each of the states g.If the transition probabilities from a particular starting state are required, the argument from must be added.These results show how the prognosis of a patient depends on his/her starting state and on the moment that is taken as the starting point for prediction.
The plot() method enables the user to show the transition probabilities in several ways.The output of probtrans() is perhaps most conveniently interpretable when plotted in a figure with stacked transition probabilities.The argument ord is used to specify an informative ordering of the transition probabilities to be stacked.The numbers in ord correspond to the states specified in the transition matrix.The argument type specifies the type of plot; we ask here for the space between adjacent curves to be filled with suitable colors (darker is more serious).For this we used the colorspace package (Ihaka et al. 2009;Zeileis et al. 2009).
R> library("colorspace") R> statecols <-heat_hcl(6, c = c(80, 30), l = c(30, 90), + power = c(1/5, 2))[c(6, 5, 3, 4, 2, 1)] R> ord <-c(2, 4, 3, 5, 6, 1) R> plot(pt0, ord = ord, xlab = "Years since transplantation", + las = 1, type = "filled", col = statecols[ord]) Figure 2 shows the result.The distance between two adjacent curves represents the probability of being in the corresponding state.The particular order chosen makes it possible to combine A comparison of Figure 2 and Figure 3 clearly shows that the fact that patient 1 has not had any adverse event in the first 100 days post-transplant has improved his/her prognosis considerably; notably, his/her probability of long-term relapse-free survival has increased significantly.On the other hand, the long-term relapse-free survival of patient 2 is unchanged by the fact that he/she has experienced the adverse event, its negative impact being balanced by the good news of still being alive and relapse-free at 100 days.The relapse probability for patient 2 is somewhat smaller than that for patient 1 and the most likely scenario for 2 is that he/she will have no further events.
When reading the figures, one must keep in mind that we consider no further information about the patients who have had a relapse as this is considered an absorbing state; many of them may die as well.Moreover, being in an intermediary state, such as AE, must be interpreted as being alive after having experienced the entering event; this does not necessarily mean that the patient continues to suffer from the event.
For the analysis of non-parametric multi-state models, several other R packages are also available, but they all have some limitations.Among them are mvna (Allignol et al. 2008) to calculate the Nelson-Aalen estimator of the cumulative hazard and etm (Allignol et al. 2011) to calculate the Aalen-Johansen estimator of the transition probability matrix and its variance-covariance matrix.Contrary to mstate, mvna and etm cannot be applied to semiparametric models.

Semi-parametric models 4.1. A model with transition-specific covariates
Next, we will show how the prediction of the transition probabilities can be improved by taking covariates at baseline into account.The mstate package supports the analysis of typespecific Cox models.Events are of the same 'type' or 'stratum' if they share a baseline hazard.
In this section, we consider a model in which 'type' is equivalent to transition: each transition has its own baseline hazard.In Section 4.2, we consider a so-called 'proportional baseline hazards model'.In both models, covariates can have the same effect for all transitions or different effects for different transitions; in the latter case, transition-specific covariates are needed.For details see de Wreede et al. (2010).
The model that we consider in this section is a transition-specific Cox model: where gh indicates a transition from state g to state h, α gh,0 (t) is the baseline hazard for this transition, Z is the vector of covariates at baseline and Z gh is the vector of transition-specific covariates (see page 7 for covariates expansion).This model specifies different covariate effects for the different transitions, as well as separate baseline transition hazards for each transition.
For the estimators of the regression coefficients and baseline hazards, see de Wreede et al.
In our leading example, we have tested for each covariate separately whether its effect is the same for all transitions or different across transitions by means of a likelihood ratio test.The results suggest a model in which all covariates have transition-specific effects.We refer to this model as the full model.In the call to coxph(), each of the expanded covariates has to be included.To specify that each transition has its own baseline transition intensity, + strata(trans) has to be added to the covariates.agecl2.12+ strata(trans), data = msebmt, method = "breslow")

R>
The estimated regression coefficients of the covariates and their standard errors for each of the transitions are shown in Table 2.For each covariate, the estimated effects are positive for some transition hazards and negative for others.The use of transition-specific covariates is very convenient to observe such effects.A model without transition-specific covariates could be estimated by substituting these by the basic covariates in the function call above.
We will now illustrate the use of probtrans() to do prediction for two example patients (see Table 3).Patient A has a (relatively) good prognosis and patient B has a bad prognosis.The prediction will be done in two steps.First, msfit() is used to estimate the transition intensities specific to these two patients.Similar to survfit() in the survival package, a newdata argument can be defined, specifying the values of the covariates of the patient.This newdata data frame looks somewhat different from newdata in survfit().In msfit(), newdata needs to have as many rows as the number of transitions and each row needs to contain the values of all the (transition-specific) covariates used in the coxph object.An additional column strata specifies to which stratum in the coxph object each transition corresponds.In this case, the values of strata are equal to those of trans (from the 'msdata' object) because every transition has its own baseline hazard.The fastest way to create data frames for these example patients is by selecting a patient from the data set with the correct characteristics of interest.Alternatively one could specify the basic covariates and apply expand.covs() to obtain transition-specific covariates.
For the current model, msfit() can only calculate Aalen-type standard errors, because the Greenwood estimator is not defined for models with covariates.The code for patient A is shown below; the code for patient B is similar.The first commands select the covariates and inherit the factor levels of the first patient in the data set with the required characteristics.
From a clinical point of view, it is worthwhile to update this prediction if more information becomes known.Assume that both patients are in state 3 (AE, no Rec or relapse) at 100 days post-transplant (compare Figure 3, right-hand side).
R> pt100A <-probtrans(msfA, predt = 100/365.25)For patient A the prognosis of relapse-free survival (i.e., being in state 1, 2, 3 or 4) is about the same as the prognosis just after transplant, but the distribution of the probabilities is different from the previous situation.By surviving the first 100 days, the prospects of patient B regarding relapse-free survival have somewhat improved.
The software also enables the user to make a different kind of predictions for individual patients.Suppose we are interested in dynamic prediction of 10-year relapse-free survival (RFS) probabilities.The question is how these prediction probabilities change as more information about intermediate events becomes known in the course of time.We consider again patient A and we want to study how 10-year RFS probabilities change when the patient experiences the adverse event 60 days (0.164 years) post-transplant and is recovered from the treatment 80 days (0.219 years) post-transplant.The 10-year survival probabilities can all be calculated by setting the direction in probtrans() to "fixedhorizon" (or in fact any string starting with "fi").
R> ptA10yrs <-probtrans(msfA, predt = 10, direction = "fixedhorizon") R> head(ptA10yrs We can extract the dynamic predictions of 10-year RFS probabilities for this specific patient directly.The 'probtrans' object ptA10yrs contains estimates of P gh (s, 10).For the 10-year RFS probabilities the quantity 1 − (P g5 (s, 10) + P g6 (s, 10)) is needed, in which g indicates the state in which the patient is at time s.For 0 ≤ s < 0.164 (time of the AE) patient A is in state 1 and we extract 1 − (P 15 (s, 10) + P 16 (s, 10)) from the first list element ptA10yrs [[1]].For 0.164 ≤ s < 0.219 (time of Recovery), patient A is in state 3 and now we need 1−(P 35 (s, 10)+ P 36 (s, 10)) from the third list element ptA10yrs [[3]].Finally, for 0.219 ≤ s < 10, we extract the relevant information from ptA10yrs [[4]].The results of this fixed horizon prediction are plotted in Figure 6 (s ≤ 2 years).They show the 10-year RFS probabilities for patient A based on his/her history as described above.The dashed lines show what the survival probabilities would have been if the patient had not had a transition.We see that her/his prognosis becomes worse the moment when the adverse event occurs and improves when recovery takes place.When time progresses, the prognosis improves irrespective of the current state simply because he/she has already survived a potentially dangerous period.Note that these curves, although tending to increase in general, do not have to be monotonely increasing, as the shape of the black curve shows (see also van Houwelingen and Putter 2008).
The function probtrans() only needs to be called once, since all the information from different starting states is present in a single 'probtrans' object.Without further computations we could study changes in 10-year RFS probabilities for patients with the same covariate values but different event histories.

A proportional baseline hazards model
In the full model discussed in Section 4.1, 12 baseline hazards and 72 covariate effects have to be estimated.These numbers can be reduced in several ways.One of them is to consider a model in which some of the baseline hazards are assumed to be proportional.Alternatively, the number of regression parameters to be estimated can be reduced by applying reduced rank techniques; these techniques, well known in regression theory, have been adapted to the multi-state context.This will be explored in Section 4.3.
In the current section, we consider a model in which we assume that all transitions into the Relapse state have a common baseline hazard and that all transitions into the Death state have another common baseline hazard.These are reasonable assumptions from a clinical point of view and they can be checked using standard methods.In formula, these are expressed as α gh,0 (t) = δgh α 1h,0 (t), g = 2, 3, 4; h = 5, 6.The transition intensities from Tx to Relapse and to Death serve as baseline transition intensities.To estimate the δgh s, we need a specific kind of time-dependent covariates Zg (t) in the regression model, g = 2, 3, 4. Within a stratum or type, this covariate distinguishes between different transitions into the same state: Zg (t) is 1 after the patient has moved to state g and 0 otherwise.These covariates may be expanded into transition-specific covariates as before.The proportionality is then expressed by the coefficient βgh of Zgh (t), in which Zgh (t) indicates the transition-specific covariate for the transition g → h expanded from Zg (t) : exp( βgh ) = δgh (for details see de Wreede et al. (2010)).Each of the other transitions also determines one stratum.
Other covariates in the new model can again be analyzed both as generic covariates and as transition-specific covariates.For clinical reasons, the last option will again be further pursued here.The new model has 6 baseline hazards instead of the 12 of the model of Section 4.1.The hazard rates, transitions probabilities and their asymptotic standard errors can be calculated as before.
The data now have to be extended with a strata column (the column name is irrelevant), indicating which transitions are together in one type or stratum, and with the time-dependent covariates indicating for which transition within a stratum the patient is at risk.Once the data have been adjusted, the analysis is very similar to that shown before.
The call to coxph() now includes the same 72 transition-specific covariates as in the full model (cfull), plus the six covariates measuring the effects of occurrences of intermediate events on relapse and death.Instead of stratifying by transition, we stratify by strata, which means that six baseline hazards are estimated instead of twelve.The first part of the function call with all transition specific covariates is equal to that of cfull, the last part becomes +Z2.6+Z2.7+Z3.9+Z3.10+Z4.11+Z4.12+strata(strata).The coefficients of this Cox model are very close to those of the full model of Section 4.1.Of the new covariates Z gh , only the coefficient of Z3.10 is significant (p = 0.0024); the hazard ratio of Z3.10 equals 2.60, which means that, adjusted for the other covariates, the occurrence of the adverse event increases the rate of dying with a factor of 2.60.The main advantage of proportional baseline hazards models is precisely that we obtain such a measure for the impact of intermediate events on the final outcome.
The prediction for Patient A is very similar to the procedure discussed before.Only the values of the strata column need to be adjusted and the time-dependent covariates Z2.6, . . ., Z4.12 need to be added.

Reduced rank models and simulation
In this section two more specialized functions of mstate are illustrated.The first of these is useful to obtain a lower dimensional representation of the regression coefficients of the full model of Section 4.1.In our example, 72 coefficients were estimated.Table 2 does not give a clear overview of the structure of the covariate effects on the transitions.By reducing the rank R of the matrix B of regression coefficients we reduce the number of parameters to be estimated.The matrix B can be factorized as B = AΓ .This implies that the number of free parameters that need to be estimated is reduced from p × K to R(p + K − R), where p and K denote respectively the number of covariates and the number of transitions in the model.For more details see Fiocco et al. (2005) and Fiocco et al. (2008).
We will illustrate the rank 1 model, which is formulated as for transitions numbered k = 1, . . ., K. In this model all covariates have the same effect (given by the parameter vector α) on each transition apart from the proportionality coefficients γ k .The factor α Z can be seen as a prognostic score for a patient with a vector of covariates Z; this prognostic score determines how likely a patient is to experience an event.The parameter γ k determines the size of the effect of the prognostic score on transition k.In this example, we use a model in which each transition has its own baseline hazard, but this is not necessary.The prognostic index in the Alpha item of rr1 is higher for later years of transplantation and for younger age.The coefficients for this risk score in Gamma are negative and of substantial size for all transitions into death (and for many into relapse).This means that in this case lower values of Alpha (for instance higher age) correspond to higher death rates, showing that a higher prognostic index does not necessarily imply a higher risk of an adverse event.
In contrast, the effects of Alpha on the "good" event of recovery (Tx → Rec and AE → Rec+AE) are positive.The matrix of regression coefficients B obtained by B = AΓ is given in rr1$Beta.

R> rr1$Beta
Tx If the rank 1 model is too simplistic, reduced rank models of higher rank may be investigated.
The standard errors of the estimates can be obtained by means of a bootstrap procedure.The bootstrap is useful in cases where asymptotic variances of estimators are not available in closed form or may be very complicated to compute.In mstate, a non-parametric bootstrap procedure has been implemented through the function msboot().Its main argument theta is a function that should return a real-valued scalar or vector.The function theta itself should have an object of class 'msdata' as its first argument.msboot() is modeled after boot() in the boot package, but one important difference is that in the multi-state context of msboot(), the dataframe contains multiple rows of data for the same individual.
As an example, we calculate the standard errors of the B matrix of the reduced rank regression of rank 1 illustrated above.For this purpose, we define theta as the stacked vector representation of B. First, a function is defined that turns B into a vector: R> rr1beta <-function(data) { + rr1 <-redrank(Surv(Tstart, Tstop, status) ~match + + proph + year + agecl, data = data, R = 1, print.level= 0) + return(as.vector(rr1$Beta))+ } R> th <-rr1beta(msebmt) Now msboot() is called with rr1beta() as theta argument.The other arguments include data (an object of class 'msdata'), id (needed to identify individuals in the data), and B, which specifies the number of bootstrap simulations.
R> set.seed(1234)R> msb <-msboot(theta = rr1beta, data = msebmt, id = "id", + B = 500, verbose = 0) The result is a matrix with B columns, each of which contains the result of the theta function applied to a bootstrap data set (see Fiocco et al. (2008) for details).It can be used to assess the bias and compute the standard errors of B. For instance, the standard errors are calculated as R> sqrt(apply(msb, 1, var)) Finally, it is possible to use simulation to calculate transition probabilities.In Markov models this is not necessary, because calculation of transition probabilities is implemented in probtrans(), but especially in models where the Markov assumption is not fulfilled, simulation is very useful.In mstate, a function mssample() is provided to do this.We will illustrate this function by recomputing transition probabilities for patient A in Section 4. The code below uses simulation to approximate the probabilities already calculated in ptA (see p. 17).
The argument tvec indicates a vector of time points at which the probabilities are to be calculated.

Competing risks models
The mstate package has been designed for general multi-state models.Since competing risks models are special cases of multi-state models, all the functionality for general multi-state models also applies to competing risks models.In particular, cumulative incidences for causespecific proportional hazards models may be obtained using msfit() and probtrans().Two functions in mstate are designed specifically for competing risks models: trans.comprisk(),which defines a transition matrix for competing risks models, and Cuminc(), which calculates non-parametric cumulative incidence functions and associated standard errors, possibly for subgroups defined by a categorical covariate.The vignette of mstate contains example code for the analysis of competing risks data using mstate; type vignette("Tutorial", package= "mstate") to access the vignette.

Discussion
Although multi-state models are a very useful tool to answer a wide range of questions in survival analysis, they are not frequently applied.So far, an important reason for this has been the lack of available software.For this reason we have developed a package in R that offers the user the opportunity to explore different kinds of multi-state models and estimate their parameters of interest on the basis of a regular data set containing the times to event of the events of interest and optionally covariate values.The functions in the package are flexible, which means that they can easily be combined with user-written software in cases when models not covered by mstate are studied.
In mstate, we restrict ourselves to non-and semi-parametric models.This means that for parametric models other software is needed, such as the R package msm developed by Christopher Jackson (see Jackson 2010Jackson , 2011)).In this article, we have explored the role of the mstate functions in the different phases of a multi-state analysis: model building, data preparation, exploration of different covariate effects and baseline assumptions, estimation of hazards, transition probabilities and associated standard errors.In particular, it has been explained how predictions can be updated if more information becomes known (dynamic prediction).
This possibility is an important extra feature of multi-state models compared to classical survival models.Moreover, several ways of presenting the outcomes in figures have been shown.Finally, we have presented some functions for reduced rank modeling and for bootstrap and simulation procedures.
In future releases of mstate we plan to implement prediction in Markov renewal models, with and without covariates.Formulas similar to the Aalen-Johansen estimator currently implemented in probtrans() also exist for this type of models, but are not yet available in any software package.Other methods are also planned to be included in mstate: dynamic prediction using landmarking of van Houwelingen (2007), vertical modeling in competing risks models of Nicolaie et al. (2010), the Fine-Gray regression method for competing risks with left truncation developed by Geskus (2010) and methods for quality-of-life adjusted survival in the context of multi-state models.

Figure
Figure 2: Stacked transition probabilities

Figure 3 :
Figure 3: Non-parametric estimates of stacked transition probabilities at 100 days posttransplant.On the left starting state is 1 (transplant), on the right starting state is 3 (AE).

Figure 6 :
Figure 6: Dynamic prediction of 10-year relapse-free survival for patient A, semi-parametric model

Figure 7 :
Figure 7: Standard errors for the transition probabilities AE → Relapse and AE → Death (starting time: 100 days), both for the full model and the proportional baseline hazards model An object of class 'msdata'

]
An object of class 'msdata' This matrix B may be compared with that of the full model shown in Table2.Differences between the reduced rank matrix and Table 2 may be explained both by parameter estimate uncertainty of both models and by a possible lack of fit of the reduced rank model of rank 1.