msSurv, an R Package for Nonparametric Estimation of Multistate Models

We present an R package, msSurv, to calculate the marginal (that is, not conditional on any covariates) state occupation probabilities, the state entry and exit time distributions, and the marginal integrated transition hazard for a general, possibly non-Markov, multistate system under left-truncation and right censoring. For a Markov model, msSurv also calculates and returns the transition probability matrix between any two states. Dependent censoring is handled via modeling the censoring hazard through observable covariates. Pointwise confidence intervals for the above mentioned quantities are obtained and returned for independent censoring from closed-form variance estimators and for dependent censoring using the bootstrap. Note that this vignette is adapted from the publication on msSurv in the Journal of Statistical Software (Ferguson et al. 2012).


Introduction
Multistate models are systems of multivariate survival data where individuals transition through a series of distinct states following certain paths of possible transitions.These systems are illustrated by a directed graph, where distinct states are treated as nodes and possible transitions are considered directed edges.Transitions between states may be reversible or irreversible while states can be either absorbing or transient.Multistate models have a wide range of application including epidemiology, dentistry, clinical trials, reliability studies in engineering, and medicine where individuals progress through the different states of a disease such as cancer and AIDS.Data in these applications are often subject to right censoring and possibly left truncation.Aalen (1978, see also Nelson 1972) proposed an estimator for the integrated hazard under a broad class of counting process models.Aalen and Johansen (1978) obtained an estimator for the transition probability matrix and subsequently state occupation probabilities through product limit integration of the Nelson-Aalen estimator.Datta and Satten (2001) established that the resulting estimators of state occupation probabilities remained valid even when the process is non-Markovian.Datta and Satten (2002) also proposed an estimator for state occupation probabilities that can handle state dependent censoring and other flexible models through a weighting function based on the censoring scheme.Estimation of state entry and exit distribution functions are also of interest (Pepe 1991;Datta and Ferguson 2011), and can be calculated through normalized sums of state occupation probabilities.
Several R (R Development Core Team 2012) packages are available on the Comprehensive R Archive Network (CRAN, http://www.r-project.org) for use with multistate models.Recently, The Journal of Statistical Software published a special volume on Competing Risks and Multi-State Models featuring papers on the msm (Jackson 2011), mstate (de Wreede et al. 2011), etm (Allignol et al. 2011), and p3state.msm(Meira-Machado and Roca-Pardinas 2011) packages.Other packages currently available include changeLOS (Wrangler et al. 2006) and mvna (Allignol et al. 2008).The msm package provides functions for fitting multistate Markov models to panel count data and offers extensions to hidden Markov multistate models and possibly inhomogeneous Markov models.The mstate package can be applied to right censored and left truncated data in semiparametric or nonparamertric multistate models with or without covariates and it may also be applied to competing risk models.The package offers functions which calculate transition probabilities and standard errors.It also uses Cox regression models to estimate different types of covariate effects.The packages changeLOS, mvna, and etm all provide methods for nonparametric estimation in multistate models.The most specialized package available is changeLOS, which is based on methods described in Schulgen and Schumacher (1996) and computes changes in length of hospital stay.It does offer a function to compute the Aalen-Johansen estimator, but it does not provide variance estimates and cannot be applied to left truncated data.The mvna package computes the Nelson-Aalen estimator of the cumulative transition hazard for any multistate model with right censored, left truncated data, but does not compute transition probability matrices.The etm package calculates the transition probability matrices and corresponding variance estimates for any finite-state multistate model with data subject to right censoring and left truncation.State occupation probabilities may be indirectly found using an initial time of zero in the etm package.
No packages currently available estimate state entry and exit time distributions, nor do they estimate desired quantities like transition probabilities for dependent censoring.These missing methods are addressed in the msSurv package (Ferguson et al. 2012), which we introduce in this paper.msSurv calculates nonparametric estimates of marginal quantities (that is, not conditional on any covariates) for time to event multistate data subject to right censoring and possibly left truncation.We assume left truncation occurs with respect to the total time of an individual in the multistate system, so that individuals are not considered for the estimation of transitions prior to their left truncation time.The main function msSurv() calculates and returns the marginal state occupation probabilities for a general, possibly non-Markov, multistate system, and provides features not currently found in other R packages.Additionally, state entry and exit time distributions are calculated and returned for states where it is appropriate.The function also calculates and returns the marginal integrated transition hazards and the hazard rate functions, and for a Markov model, the transition probability matrix between any two times.Users specify whether the censoring is independent or state dependent and then appropriate calculations of variance and confidence intervals are performed.Pointwise confidence intervals for the above mentioned quantities are obtained and returned for independent censoring from closed-form variance estimators and from bootstrapping for dependent censoring.The function returns an object of S4 class msSurv which includes state and transition information, event times, estimates, variance estimates, confidence intervals, and counting process information.The msSurv object can then be summarized or plotted using methods available in the package.
The rest of this paper is organized as follows.Section 2 describes nonparametric estimation methods which are used in the msSurv package.Section 3 describes the implementation of msSurv and illustrates available functions through examples.Possible future extensions of our software package are discussed in Section 4.

The estimators
Consider a multistate model with a finite state space S = {1, . . ., M } and set of possible transitions between the states in S. The quantities of interest (like state occupation probabilities and state entry and exit distributions) can be calculated for complete data, when available, and also using estimators obtained from censored data when complete data is not known.It is useful to keep track of all the transitions an individual makes before ending in an absorbing state.Let T * ik represent the time of the kth transition for individual i = 1, . . ., n, where T * i0 = 0 and T * ik = ∞ if the ith individual enters the absorbing state before the kth transition is made.Let C i be the right censoring time for the ith individual, L i be the left truncation time for the ith individual, and s ik be the state occupied by the ith individual between times T * ik−1 and T * ik .Let T * i = sup k {T * ik : T * ik < ∞} be the time for the last transition for individual i.The collection of all transition times and states occupied by individual i can be denoted as T * i = (T * ik : k ≥ 1) and s * i = (s ik : k ≥ 1), respectively.Define an indicator of whether the ith individual was never censored, e.g., δ i = I(C i > T * i ), and let T i = min(T * i , C i ).Usually, one assumes that censoring is independent of the multistate system.However, there are times when censoring and hazards of future transitions are affected by time varying covariates Z = Z (t).In these cases, the covariate variables Z explain the dependence and thus future transitions and censoring events behave conditionally independent given the covariate process Z (Datta andSatten 2001, 2002).Modeling of this censoring hazard using the covariate process is discussed in Section 2.3.Andersen et al. (1995) presented formulas for the Nelson-Aalen estimator for the integrated hazard matrix A and the Aalen-Johansen estimator of the state occupation probability matrix of a Markov system.The counting process and the number at risk for data subject to left truncation and right censoring are estimated as

Nelson-Aalen and Aalen-Johansen estimators
and The Nelson-Aalen estimator of the cumulative hazard is given by The Aalen-Johansen estimator of the transition probability matrix of a Markov multistate system is obtained by product integration of A jj ′ (t), i.e., where A = A jj ′ , which reduces to simple empirical proportions for the complete data.
A recursive formula for computing the variance of transition probability can be found in Andersen et al. (1995)(see formula 4.4.19 on p. 295 for details).The resulting estimator is of the Greenwood-type.

State occupation probabilities
The state occupation probability answers the marginal question: "What is the probability that an individual is in state j at time t?" Let p j (t) = P r {s (t) = j} denote the state occupation probabilities where s (t) is the state occupied by an individual at time t, j ∈ {1, ..., M }.For incomplete data, the state occupation probabilities are estimated as where p kj (0, t) is the k, jth element of the matrix P (0, t) = (0,t) I + d A (u) and p k (0) is the initial state occupation proportions for state k.This estimator is in fact the Aalen-Johansen estimator of state occupation and holds its validity regardless of Markovian assumptions (Datta and Satten 2001).Estimation of state occupation probabilities can be extended to dependent censoring and other flexible models by explicitly modeling the censoring process (Datta andSatten 2000, 2002).

Datta-Satten estimators
Datta and Satten (2002) use the principle of inverse probability of censoring weights (Robins and Rotnitzky 1992) to extend the Nelson-Aalen and Aalen-Johansen estimators to data subject to dependent censoring.The resulting Datta-Satten estimators use a weighting function, denoted by K, which is based on a fitted model of the censoring hazards using Aalen's linear hazards model with possibly time-dependent covariates Z (Aalen 1980).
Reweighted estimators for the counting process and the size of the at risk sets of complete data are defined as and where See Datta and Satten (2002) for a formal argument.
For state dependent right censoring (i.e., when Z(t) = s(t)), this is equivalent to estimating the censoring hazard by a state specific Nelson-Aalen estimator of censoring (Datta and Satten 2002).Currently, the msSurv package only implements this type of dependent censoring.
Substituting the formulas for the estimated counting process and the number at risk into Equations 3, 4, and 5 yields the Datta-Satten estimators of integrated transition hazards, transition probabilities for Markov Systems and state occupation probabilities for non-Markov systems under dependent censoring.Variance estimates for the Datta-Satten estimators of transition probabilities are obtained using the bootstrap.

State entry and exit distributions
An important application of state occupation probabilities is computing the state entry and exit time distribution functions.We assume that the only path between the states before and after the state of interest (e.g., state j) is through state j, and that the system is acyclic with respect to state j (i.e., it is not possible to return to state j after leaving it).Suppose X j denotes the possibly unobserved indicator of an individual ever entering state j.Suppose F j and G j denote the state entry and exit time distribution functions, respectively, for the individuals who ever enter state j (i.e., X j = 1).Let S j denote the collection of all states which come after state j in the progressive model.The entry time distribution to state j is estimated by taking the normalized sum of the estimated state occupation probabilities of state j and all the other states that come after j in the system, i.e., where p k (∞) = lim t→∞ p k (t).The denominator in Equation 8is the probability of an individual ever entering state j, while the numerator gives the unconditional probability of having entered state j by time t.The unconditional probability (e.g., the subdistribution function) may be of interest in itself, and the msSurv package has the option to return either the normalized or unnormalized versions.
The exit time distribution from a transient state j is estimated by taking the normalized sum of estimated state occupation probabilities of all states that come after state j in the progressive system, i.e., The denominator here is the same as in Equation 8, and msSurv provides the option of returning either the normalized or unnormalized version of the function.
The variance for the state entry and exit time distributions are obtained using the bootstrap.

Confidence intervals
Pointwise confidence intervals for estimators of transition probabilities and state occupation probabilities follow methods described in Andersen et al. (1995).Let p (s, t) be the transition probability between two states in the system (the subscripts are omitted to simplify the notation) between times s and t and let σ (s, t) be the corresponding variance estimate.Then the linear confidence interval for p (s, t) is defined as where c (α/2) is the upper α/2 percentile of the standard normal distribution.It may be beneficial to consider transformations to improve estimation especially in the case of small sample sizes (Bie et al. 1987;Thomas and Grunkemeier 1975).Borgan and Liestol (1990) suggested a log transformation to improve small sample properties.The resulting formula for the confidence interval is Other transformations to improve small sample properties include the log-log transformation, proposed in the first edition of Kalbfleisch and Prentice (2002) and defined as p (s, t) , and the complementary log-log transformation .
Confidence intervals for state occupation probabilities are calculated using s = 0 in the formulas above.
Confidence intervals for state entry and exit time distributions for state i at time t are calculated using F i (t) and G i (t) instead of p (s, t), with corresponding variance estimate σ(t) obtained using the bootstrap.

msSurv package implementation
The msSurv package may be applied to any general multistate model with data subject to right censoring and possibly left truncation.msSurv is written entirely in the R programming language using S4 classes and methods, and is available for download from CRAN.It can be installed on all operating systems for which the R software is installed.Other package dependencies include the graph package (Gentleman et al. 2010) and the lattice package (Sarkar 2008).
The main function of the package, msSurv, calculates the counting processes and risk sets according to both Andersen et al. (1995) and Datta and Satten (2001), as well as the state occupation probabilities and transition probabilities described in Section 2.5.The msSurv package contains a function to calculate the transition probabilities between two specific times (Pst), a function to display the state occupation probabilities at a specific time t (SOPt), a function to display the state entry and exit time distributions at a specific time t (EntryExit), as well as print, plot,and summary methods for msSurv objects.
We will illustrate the application of the msSurv package through three examples.The first example uses simulated data with independent right censoring, the second example uses a simulated data set with left truncation and independent right censoring, and the third example uses a bone marrow transplant data set from Klein and Moeschberger (2010) with state dependent right censoring.

A 5 state example
We consider a five-state progressive model with the tree structure illustrated in Figure 1.We simulated a data set of 1000 individuals subject to independent right censoring with 60% of individuals starting in state 1 at time 0 and 40% starting in state 2. Those in state 1 remain there until they transition to the transient state 2 or the terminal state 3. Individuals in state 2 remain there until they transition to either terminal state 4 or 5.
A right censoring time is generated for each of the 1000 individuals using the log normal distribution with log mean -0.5 and log standard deviation 2. For individuals starting in state 1, two times are generated using the Weibull distribution using a sample size of 600 and shape parameter of 2 to reflect transition times between states 1 and 2 and states 1 and 3. Times are compared and the minimum time is kept as the event time and the corresponding state is recorded.Then, two additional times are generated to reflect transitions between states 2 and 4 and states 2 and 5.These times are generated using the formula where T 1 is the first transition time, R 2 is a random number generated from U (0, 1) independent of T 1 , D (•) denotes the distribution function for the Weibull distribution with shape parameter 2, D −1 (•) denotes the corresponding quantile function.The minimum of these two times and the censoring times are compared and the minimum time is taken as the event time and the corresponding state is recorded.For individuals starting in state 2, two times are generated using the Weibull distribution with a sample size of 400 and shape parameter of 2 to reflect transitions between states 2 and 4 and states 2 and 5.These times are then compared with the corresponding censoring times and the minimum time is taken as the event time and the state information is recorded.All times were rounded to the fourth decimal place for clarity of presentation.The simulated data is available as RCdata in package msSurv.
We begin by loading the package.

R> library("msSurv")
Now we load the data.

R> data("RCdata")
Data should be in a data frame with column names "id", "stop", "start.stage",and "end.stage",where "id" is the individual's identification number, "stop" is the transition time from state j to j ′ , "start.stage" is the state the individual is transitioning from (i.e., j), and "end.stage" is the state the individual is transitioning to (i.e., j ′ ) and equals 0 if right censored.Now we specify the tree structure for this multistate system.First, input the states in the multistate system as a list or character vector of the state names.Then, store the transition information as a list of possible states with allowed transitions being lists of edges.For terminal states, the lists of edges will be NULL.Nodes correspond to the states in the model and edges refer to the allowed transitions.
Exit distributions calculated for states 1 2 .
Using bootstrap to calculate variances ...
Results of the analysis can be viewed using the print, plot, and summary methods, as well as the Pst, SOPt and EntryExit functions, available for the msSurv object ex1.The print method is discussed in this section while the summary and plot methods are discussed in Sections 3.2 and 3.3, respectively.
The print method gives an overview of the model, specifying the number of transient and absorbing states and identifying the states and possible transitions in the system.
The msSurv package contains two functions for the user to easily access estimators of transition and state occupation probabilities at specific time points.The package also contains a function for the user to access state entry and exit time distribution estimators at specific time points.
The transition probabilities between any two times s and t are computed using the Pst function.The Pst function takes an msSurv object, a starting time s, and an ending time t and prints the transition probability matrix P (s, t).For example, the following code produces the transition probability matrix P (1, 3.1): msSurv, an R Package for Nonparametric Estimation of Multistate Models R> Pst(ex1, s = 1, t = 3.1) Estimate of P(1, 3.1) cols 1 2 3 4 5 1 0 0 0.4805 0.307 0.2124 2 0 0 0.0000 0.586 0.4140 3 0 0 1.0000 0.000 0.0000 4 0 0 0.0000 1.000 0.0000 5 0 0 0.0000 0.000 1.0000 The corresponding covariance matrix for this transition probability is printed by adding the argument covar = TRUE.

R> Pst(ex1, s = 1, t = 3.1, covar = TRUE)
State occupation probabilities at a specific time t are given using the SOPt function.This function takes a msSurv object as well as time t as arguments.The default time t is the maximum event time in the data set.Individuals may start in any state in the system at time 0. The function prints the state occupation probabilities for all states in the system at time t.For example, to find the state occupation probabilities at time t = 0.85.The corresponding variance estimates are provided using the argument covar = TRUE.Variance estimates are found by evaluating the formula in Andersen et al. (1995) at P (0, t) for data where all the individuals start in a single initial state at time 0. The bootstrap is used to find variance estimates of state occupation for data where individuals start in different states of the system.

R> SOPt(ex1, t = 0.85, covar = TRUE)
The EntryExit function in msSurv displays the state entry and exit time distributions at a specific time t.This function takes an msSurv object and time t as arguments and displays the either normalized (if norm = TRUE, the default) or non-normalized (subdistribution) state entry and exit distributions.State entry distributions are only displayed for non-initial states while state exit distributions are only displayed for non-terminal states.Estimates are rounded to four decimal places by default, but the user may specify a different number through the deci argument.For example, the normalized state entry and exit distributions at time 1 are displayed by

R> EntryExit(ex1, t = 1)
The The bootstrap (via the bs = TRUE argument in the call to msSurv) is needed to obtain variance estimates for the state entry and exit time distributions.If this was done, then the var = TRUE argument in the call to the EntryExit function can be used to display the variance estimates.Note that all of the function Pst, SOPt, and EntryExit return their results invisibly, so that users can save and access the estimates for later use if desired.For example, in the following code chunk the results of the EntryExit call are stored in res, which consists of a list with items entry.norm,entry.var.norm,exit.norm, and exit.var.norm.If instead subdistributions were requested, (through the argument norm = FALSE), then the elements in the list would be the same except with 'sub' in place of 'norm'.

Left truncation and right censoring example
We consider an irreversible three-state illness-death model with data subject to independent right censoring and left truncation (Andersen et al. 1995).We simulated a data set of 1000 individuals starting in state 1 at time 0. Individuals remain in state 1 until they transition to the transient state 2 (ill) or the terminal state 3 (death).Individuals in state 2 remain there until they transition to the terminal state 3 (death).
Two times are generated using the Weibull distribution with a sample size of 1000 and shape parameter of 2 to reflect transition times for either illness or death.Right censoring and left truncation times are generated using the log normal distribution with mean -0.5 and standard deviation 2 on the log scale.We assume 20% of individuals have a left truncation time of 0. Only individuals whose left truncation times were less than the terminal event times or censored event times were kept.The left truncation time was taken to be the time the individual entered the study.Individuals were not included in the at risk set before their left truncation time.Times for these individuals were compared and the minimum time was kept as the event time and the corresponding state was recorded.Then, another time was generated reflecting the transition between states 2 and 3 using the formula where T 1 is the first transition time, R 2 is a random number generated from U (0, 1) independent of T 1 , D (•) denotes the distribution function for the Weibull distribution with shape parameter 2, D −1 (•) denotes the corresponding quantile function.These "death" times are then compared to the censoring times and the minimum time is kept as the event time and the corresponding state information is recorded.All times were rounded to the fourth decimal place for clarity of presentation.The simulated data is available as LTRCdata in package msSurv.
Begin by loading the data

R> data("LTRCdata")
Data should be in a data frame with column names "id", "start", "stop", "start.stage",and "end.stage"where "id" is the individual's identification number, "start" is the start time for the period of observation after the individual enters state j (and is the left truncation time for the first observed transition), "stop" is the transition time from state j to j ′ , "start.stage" is the state the individual is transitioning from (i.e., j), and "end.stage" is the state the individual is transitioning to (and equals 0 if right censored).
Exit distributions calculated for states 1 .
When left truncation is present, starting proportions in each state are calculated empirically based on the initial states of those individuals whose left truncation time is prior to the first observed transition time.
Results of the analysis will be illustrated through the summary method available for the msSurv object ex2.The summary method displays information for the state occupation probabilities, state entry and exit distributions, and transition probabilities.Estimates of state occupation probability are displayed with corresponding variance estimates and cofidence intervals (denoted lower.ciand upper.ci in the output).The default settings give these estimates to three decimal places for key percentile event times (minimum, maximum, 25th percentile, median, and 75th percentile).The summary method also displays estimates of the normalized state entry and exit time distirbutions (entry.normand exit.norm,respectively), as well as the corresponding subdistributions (entry.suband exit.sub).
The summary method also provides summary information for each allowed transition in the system.Estimates of the transition probabilities are given along with estimates of the variance and corresponding confidence intervals (lower.ciand upper.ci).In addition, the risk sets (n.risk), number of transitions (n.event, for transitions from one state to a different state), and number still remaining (n.remain, for transitions into the same state) are all displayed based on the counting process calculations in Andersen et al. (1995).Additionally (using DS = TRUE), users can display the weighted counting processes based on the formulas in Datta and Satten (2001)  Confidence intervals provided by default are 95% linear confidence intervals, but the user may change the confidence level to either 90% or 99% by changing the ci.level argument or apply a transformation of "log", "cloglog" or "log-log" by changing the ci.trans argument.
The user may change the number of significant digits through the digits argument.
Summary information for all event times in the data set are displayed using the all = TRUE argument, while summary information at specific timepoints can be requested using the times argument.

R> summary(ex2, all = TRUE)
Users can customize the display of information by using the logical arguments stateocc, trans.pr,and dist to select whether the state occupation probabilities, transition probabilities, and state entry and exit time distributions are displayed.For example, to display only the state occupation probabilities use the following code: R> summary(ex2, trans.pr= FALSE, dist = FALSE)

An example of state dependent censoring
State dependent censoring in msSurv will be illustrated using data on 136 cancer patients who received bone marrow transplants found in Klein and Moeschberger (2010).Following  Datta and Satten (2002), we defined five states based on platelets returning to a normal level, presence of acute graft-versus-host disease (GVHD), and onset of chronic GVHD.
State 2 is entered when acute GVHD develops before the patients platelet level returns to normal.State 3 is entered if the patient's platelet level returns to normal before acute GVHD develops.State 4 is for patients who have both normal platelet levels and acute GVHD while state 5 is for those patients who have chronic GVHD.Patients transition to either state 2 or state 3 and then either state 4 or state 5, as depicted in Figure 2. Patients do not necessarily progress to state 5 as they may remain in any state for any amount of time.Those patients who died or experienced relapse were considered censored for this example.
Our analysis allows for the censoring hazards to vary by state, since relapse or death can be predicted by the different (immunologic) states defined in this system.The cumulative censoring hazard for each state is estimated as the Nelson-Aalen estimator for censoring at each state and used to weigh the counting processes.
The default number of bootstraps is 200 and that is used for this example.
Exit distributions calculated for states 1 .
Using bootstrap to calculate variances ...
The results are displayed graphically using the plot method.The plot method takes msSurv objects and produces plots of estimated quantities using functions available in the lattice package (Sarkar 2008).The plots of the state occupation probabilities for every state in the system are produced by default with their corresponding 95% linear confidence intervals.For state dependent censoring, these confidence intervals are based on variances obtained through bootstrapping.Each state is plotted separately in a single panel.The results are given in Figure 3.

R> if(exists("DepEx")) plot(DepEx)
Users may also specify specific states to be plotted using the states argument.For example, to plot the state occupation probabilities for state 2, the user would type R> if(exists("DepEx")) plot(DepEx, state = "2") To plot the state occupation probabilities for states 2 and 3, use

Plot of State Occupation Probabilites
Event Times State Occupation Probabilities R> if(exists("DepEx")) plot(DepEx, state = c("2","3")) Confidence intervals may be altered using the ci.level and ci.trans arguments which allow the user to specify a different confidence level (0.9 or 0.99) or a different transformation ("log", "log-log", and "cloglog", as described in Section 2.5).Confidence intervals are omitted from the plots using the CI = FALSE argument.
The plot.type argument is used to change the plotted estimators.As previously mentioned, the default is "stateocc" for the state occupation probabilities, but the user may create plots for the transition probabilities ("transprob"), normalized state entry and exit distributions ("entry.norm"and "exit.norm",respectively), and state entry and exit subdistribution functions ("entry.sub"and "exit.sub",respectively).The state entry time distributions are plotted for all states except the initial state by default, but users may specify specific states using the "states" argument.The state exit time distributions are plotted for all nonterminal states by default, but specific non-terminal states may be requested by the user.We included the argument bs = TRUE in the call of msSurv, so that confidence interval estimates of the state entry and exit distributions will be plotted.
The default plot produced for transition probabilities plots estimates and confidence intervals for all possible transitions in the system.
Users may specify specific transitions using the trans argument.For example, if the user wants to plot the transition probabilities for the transitions out of state 1, say the "1 2" and "1 3" transitions, they would type R> if(exists("DepEx")) + plot(DepEx, plot.type= "transprob", trans = c("1 2", "1 3")) For state entry / exit distributions, in the current example the conditions given for ensuring proper calculation of these distributions (see Section 2.4) are only satisfied for states 1 and 5.
To calculate state entry / exit distributions for the transient states, say states 2 and 3, we need to expand the model to differentiate between paths to states 4 and 5 that go through state 2 and paths that go through state 3.The expanded system is plotted in Figure 5, where states 4 and 5 have been relabeled as 4-2 and 5-2 for paths that go through state 2 and 4-3 and 5-3 for paths that go through state 3.The data needs to be modified correspondingly to incorporate these additional states.Note that individuals who start in state 4 are dropped from this model, as these individuals have no impact on the entry / exit distributions from states 2 and 3.The updated tree and data are in objects deptree2 and data.sdc2,respectively.
The updated model is fitted below, and the normalized state entry distribution plots for states 2 and 3 over the time range of 0 to 100 days are given in Figure 6.The corresponding normalized exit distributions are plotted in Figure 7.By default (i.e., with the states="ALL" argument) all states where calculation of entry / exit distributions is valid are plotted.
Using bootstrap to calculate variances ... days.About 40% of these patients never develop GVHD (remain in state 3, Figure 7), with a small percentage subsequently developing acute GHVD (3 → 4 transition in Figure 4) and 60% eventually going on to develop chronic GVHD (3 → 5 transition in Figure 4).A much smaller percentage of patients first develop acute GVHD (1 → 2 transition in Figure 4), with about 40% of these having their platelet levels returning to normal by 70-80 days (peak of 2 → 4 transition in Figure 4).All of the patients who develop acute GVHD first eventually go on to develop chronic GVHD (2 → 5 transition in Figure 4).

An illustration of state-dependent censoring using simulated data
In the previous example based on real data the dependency of the censoring distribution on the state in the system was not too heavy, and the difference between estimates assuming state-dependent and state-independent censoring was relatively minor (data not shown).In this example we illustrate the magnitude of the difference between the two estimates that can occur using simulated data in the same spirit as Datta and Satten (2002).Data were simulated from a three state irreversible illness-death model, where all individuals begin in state 1 (well) at time 0. Individuals in state 1 may subsequently transition to state 2 (illness) or state 3 (death).Individuals with illness remain in state 2 until death (entry into state 3).For each individual, a vector W = (W l , W 2 , W 3 ) was generated from a trivariate normal distribution with mean vector (1.81, 1.7, 1.7) and variance-covariance matrix 0.393 × (0.01I + 0.9111 T ), where I is the identity matrix and 1 is a 3 × 1 vector of ones.The values of W are (log) latent transition times with the following possible transition paths:  2 and expanded in Figure 5.
The hazard for censoring was state dependent with censoring times following a Weibull distribution.In state 1, the censoring hazard was t 1/4 /60, while in state 2 the censoring hazard was t/50.For each observation, a state 1 censoring time C 1 was generated and if C 1 < min{exp(W l ), exp(W 2 )} the observation was censored in state 1.If exp(W 2 ) < C 1 and if state 2 was entered, then a state 2 censoring time was generated from the Weibull distribution conditional on C 2 > exp(W 2 ).Simulations were based on a single data set of 1000 subjects.Figure 8 illustrates the state 2 occupation probability estimates based on the complete (non-censored) data, the Aalen-Johansen estimates, and the Datta-Satten estimates assuming state-dependent censoring.The Datta-Satten estimates closely follow the empirical estimates based on non-censored data, while the Aalen-Johansen estimates deviate appreciably.

Plot of Normalized State Exit Time Distributions
Event Times Normalized State Exit Time Distributions

Discussion
We present a comprehensive R package for nonparametric estimation of general multistate models subject to independent right censoring and possibly left truncation.This package computes the marginal state occupation probabilities and state entry and exit time distributions for possibly non-Markov models.For a Markov model, the R package, msSurv (Ferguson et al. 2012), also calculates and returns the marginal integrated transition hazard and the hazard rate functions, as well as the transition probability matrix between any two states.Motivated by Datta and Satten (2001), msSurv also performs nonparametric estimation for state dependent right censored and possibly left truncated data.Currently no other packages available on CRAN calculate the state entry and exit time distributions for censored multistate data or calculate nonparametric estimates for data subject to state dependent censoring.Pointwise confidence intervals for state occupation probability and transition probability matrices are obtained and returned for independent censoring from closed-form variance estimators and for state dependent censoring using the bootstrap.Pointwise confidence intervals for state entry and exit time distributions are obtained using the bootstrap.The bootstrap is also used to find variance estimators for state occupation probabilities of data subject to state dependent censoring.The msSurv package provides functions to find the state occupation probabilities at a specific time t and to find the transition probabilities between any two times s and t (s ≤ t).Package msSurv is written using S4 classes and methods.Methods are available to print, plot, and summarize the msSurv objects.
The msSurv package has a number of limitations that will be improved upon in future releases.
Currently state dependent censoring is the only censoring scheme incorporated in msSurv, future expansion could include incorporating general covariates into the (dependent) censoring mechanism.The package could also be extended to include estimation for current status data and eventually interval censored data (Datta and Sundaram 2006;Lan and Datta 2010).Other extensions include calculations of the state waiting time distributions (Satten and Datta 2002) and allowing estimation for recurrent event models.

Figure 1 :
Figure 1: The five state model for the simulated multistate data discussed in Section 3.1.

Figure 2 :
Figure 2: The five state model for cancer patients who received bone marrow transplants, discussed in Section 3.3.

Figure 3 :
Figure 3: Plot of state occupation probability estimates and corresponding confidence intervals for each state in the system illustrated in Figure 2.

Figure 7 :
Figure 7: Plot of normalized state exit time distribution estimates and corresponding confidence intervals for states 2 and 3 in the system initially illustrated in Figure 2 and expanded in Figure 5.

Figure 8 :
Figure 8: Plot of state 2 occupation probabilities for the three-state irreversible illness-death model in Section 3.4.The Datta-Satten (DS) estimator closely follows the empirical estimator based on complete (non-censored) data, while the Aalen-Johansen (AJ) estimator deviates noticeably.

of Normalized State Entry Time Distributions
Death occurs before illness at time exp(W 1 ) W 2 < W 1 : Illness occurs at time exp(W 2 ) and death at time exp(W 2 ) + exp(W 3 ) Plot