TPmsm : Estimation of the Transition Probabilities in 3-State Models

One major goal in clinical applications of multi-state models is the estimation of transition probabilities. The usual nonparametric estimator of the transition matrix for non-homogeneous Markov processes is the Aalen-Johansen estimator (Aalen and Johansen 1978). However, two problems may arise from using this estimator: ﬁrst, its standard error may be large in heavy censored scenarios; second, the estimator may be inconsistent if the process is non-Markovian. The development of the R package TPmsm has been motivated by several recent contributions that account for these estimation problems. Estimation and statistical inference for transition probabilities can be performed using TPmsm . The TPmsm package provides seven diﬀerent approaches to three-state illness-death modeling. In two of these approaches the transition probabilities are estimated conditionally on current or past covariate measures. Two real data examples are included for illustration of software usage.


Introduction
In many longitudinal studies it is often of interest to investigate time to a certain event.In medicine the event is an ultimate outcome, such as diagnosis of "death" of the patient or "relapse of the disease".In addition to this primary event of interest one may observe also a number of intermediate ("transient") states, such as "local recurrence" and "distant metastasis" in cancer studies.Analysis of such studies where individuals may experience several events can be successfully performed using a multi-state model (MSM).An MSM is a stochastic process (X(t), t ∈ T ) with a finite state space, where X(t) represents the state occupied by the process at time t ≥ 0. Graphically, these models are represented by diagrams with rectangular boxes and arrows between them indicating respectively possible states and possible transitions.In general, the future state transitions of an MSM may depend on past events.However, for the special case of a Markov model the past and future are independent given its present state.There exists an extensive literature on MSMs.Main contributions include books by Andersen, Borgan, Gill, and Keiding (1993) and Hougaard (2000).Recent reviews on this topic may be found in the papers by Hougaard (1999), Commenges (1999), Andersen and Keiding (2002), Putter, Fiocco, and Geskus (2007) and Meira-Machado, de Uña-Álvarez, Cadarso-Suárez, and Andersen (2009).
The simplest form of an MSM is the mortality model for survival analysis with only two states "alive" and "dead" with a single transition.Other common models include the progressive three-state model, the illness-death model and the competing risks model.The illness-death model is probably the most used model in the literature, in particular for studying progression of many diseases.This model describes the dynamics of healthy subjects who may move to an intermediate "diseased" state before entering into a terminal absorbing state.Many longitudinal medical data with multiple endpoints can be reduced to this structure.Thus, methods developed for the illness-death model have a wide range of applications.There are several issues of interest in an illness-death multi-state model: study of the relationship between covariates and disease evolution; estimation of transition probabilities; state occupation probabilities; distributions of time spent in each state, among other topics.In this paper we will focus on the inference for the transition probabilities.These quantities provide estimates of the clinical prognosis of a patient at a given point in disease progression, allowing long-term predictions of the process.Aalen and Johansen (1978) introduced a nonparametric estimator for these quantities for nonhomogeneous Markov models.Their estimation method extends the Kaplan-Meier estimator (Kaplan and Meier 1958) to Markov chains.As for the Kaplan-Meier estimator, the standard error of the Aalen-Johansen estimator may be large when the censoring is heavy, particularly with a small sample size.To overcome this problem, Moreira, de Uña-Álvarez, and Meira-Machado (2013) propose a modification of Aalen-Johansen estimator based on presmoothing (Dikta 1998), which allows for a variance reduction in the presence of censoring.In a recent paper, Meira-Machado, de Uña-Álvarez, and Cadarso-Suárez ( 2006) introduce a substitute for the Aalen-Johansen estimator in the case of a non-Markov illness-death model.They show that when the Markov assumption does not hold, their estimator may behave much better than the Aalen-Johansen which may be systematically biased.The idea behind their estimator is to weight the data by the Kaplan-Meier weights pertaining to the distribution of the total survival time of the process.However, by removing the Markov condition, the proposed substitute for the Aalen-Johansen estimator provides undesirably large standard errors.This problem becomes worse when there is a large proportion of censored data.In order to overcome this problem, Amorim, de Uña-Álvarez, and Meira-Machado (2011) propose a modification of the Meira-Machado estimator based on presmoothing.Other estimators were proposed to estimate the transition probabilities.A valid estimator was provided by Keilegom, de Uña-Álvarez, and Meira-Machado (2011) for a progressive three-state model.This methodology assumes that the vector of gap times (time in State 1 and time in State 2) satisfies the nonparametric location-scale regression model, allowing for the transfer of tail information from lightly censored areas to heavily ones.All these approaches assume independent censoring and do not account for the influence of covariates.To this regard in a recent work, in a regression setup, Meira-Machado, de Uña-Álvarez, and Somnath (2012) introduce feasible estimation methods for the transition probabilities in an illness-death model conditionally on current or past covariate measures.
Software for multi-state survival analysis has been developed recently.A comprehensive list of the available packages at the Comprehensive R Archive Network (CRAN) can be seen in the CRAN task view "Survival Analysis" (Allignol and Latouche 2014).An issue of the Journal of Statistical Software, entirely devoted to these models, was published in 2011 (Putter 2011).In R (R Core Team 2014) several packages provide functions for estimating the transition probabilities.The timereg package (Scheike and Martinussen 2006;Scheike and Zhang 2011) can be used to obtain the cumulative incidence probability of a specific cause of failure in competing risks data.It also provides an estimate of its variance at each fixed time point, and constructs (1 − ³)100% simultaneous confidence bands over a given time interval.The package cmprsk (Gray 2014) can also be used to obtain the same quantities.The package etm (Allignol, Schumacher, and Beyersmann 2011) computes and displays the transition probabilities for the Aalen-Johansen estimator.This package also features a Greenwoodtype estimator of the covariance matrix.The package msm (Jackson 2011) can be used to obtain estimates for the transition probabilities in time-homogeneous Markov models.The package p3state.msm(Meira-Machado and Roca-Pardiñas 2011) enables the user to perform inference in the illness-death model.The main feature of the package is its ability for obtaining non-Markov estimates for the transition probabilities.Finally, the msSurv package (Ferguson, Datta, and Brock 2012) can be used to estimate the state occupation probabilities along with the corresponding variance estimates, and lower and upper confidence intervals.All of the existing software presents, however, some limitations in practice.Most software assumes the process to be Markovian and assumes independent censoring.Furthermore such software does not account for the influence of covariates.In addition, a comparison between different packages is rather difficult because each of the current programs requests its own data structure.This paper describes the R package TPmsm (Araújo, Roca-Pardiñas, and Meira-Machado 2014) which is available from CRAN at https://CRAN.R-project.org/package=TPmsm.The package aims at implementing nonparametric and semiparametric estimators for the transition probabilities in 3-state models.The package provides the so-called Aalen-Johansen estimator typically assumed in Markov processes but it also covers alternative methods which have been proved to be consistent even without the Markov assumption.Inverse censoring probability reweighting is used in some methods to deal with right censoring.These approaches lead to consistent estimators even in the presence of dependent censoring.Finally, two different estimators are implemented that account for the influence of covariates.Bootstrap confidence bands are provided for all methods.In this article we explain and illustrate how numerical and graphical output for all methods can be obtained using the TPmsm package.
In Section 2 we introduce the notation for the illness-death stochastic model and describe in detail the proposed estimation methods.In Section 3 we describe the implementation of the methods in package TPmsm.Some of the methods are illustrated using generated data in Section 4. Finally, Section 5 illustrates the package's capabilities using two real data examples, and Section 6 gives some concluding remarks and proposals for future work.

Methodological background
In this paper we consider the progressive illness-death model depicted in Figure 1.We assume that all subjects are in State 1 at time t = 0, and that they may either visit State 2 at some time point; or not, going directly to the absorbing state (State 3).The stochastic behavior of the process is represented by a random vector (T 12 , T 13 , T 23 ), where T hj is the potential transition from State h to State j, 1 ≤ h < j ≤ 3, in which T 23 is the sojourn time in State 2. In this model we have two competing transitions 1 → 2 and 1 → 3. Let the sojourn time in State 1 be denoted by Z = min(T 12 , T 13 ).The survival time of the process is given by T = I(T 12 ≤ T 13 )(T 12 + T 23 ) + I(T 12 > T 13 )T 13 .However, censoring may occur due to follow-up limitations, lost cases and so on.Because of censoring, one observes ( Z, T , ∆ 1 , ∆) where Here C denotes the potential censoring time, which we assume to be independent of the process (that is, C and (Z, T ) are assumed to be independent).
Given two time points s < t, define the transition probabilities as p hj (s, t) = P(X(t) = j|X(s) = h).The transition between the three stochastic states is illustrated in Figure 1 Details about the inference for the transition intensities can be seen in Andersen and Perme (2008).
The transition probabilities can also be estimated nonparametrically or semiparametrically using the notation shown in the top of this section.The expressions for the transition probabilities are given by

Aalen-Johansen estimator
The transition probabilities may be estimated via the nonparametric (Aalen-Johansen estimator) model.This can be thought as the generalization of the Kaplan-Meier approach (Kaplan and Meier 1958) for the simple mortality model and was proposed by Aalen and Johansen (1978) for general non-homogeneous Markov models with a finite number of states.
Explicit formulae of the Aalen-Johansen estimator for the illness-death model are available (Borgan 1998).Here we give alternative expressions for this estimator using the notation introduced above.The Aalen-Johansen (AJ) estimate of the transition probability p 11 (s, t) is the Kaplan-Meier estimator where Similarly, the estimate of the transition probability p 22 (s, t) is the Kaplan-Meier estimator where Finally, the estimator for p 12 (s, t) is given by where

Presmoothed Aalen-Johansen estimator
The standard error of the Aalen-Johansen estimator may be large when the censoring is heavy, particularly with a small sample size.Interestingly, the variance of the Aalen-Johansen estimator may be reduced by presmoothing (Dikta 1998).Presmoothing the Aalen-Johansen estimator (Moreira et al. 2013) involves replacing the censoring indicators (in the transition probabilities p 11 (s, t) and p 22 (s, t)) by a smooth fit (e.g., using logistic regression).Then, the corresponding presmoothed Aalen-Johansen (PAJ) estimator of p 11 (s, t) is given by where m 0n ( Z) stands for an estimator of the conditional probability of the event ∆ 1 = 1 given Z; which can be estimated using logistic regression.The presmoothed version of (2) given by where m 1n ( Z, T ) stands for an estimator of the conditional probability of the event ∆ = 1 given ( Z, T ) and given that transition 1 → 2 is observed.Finally the transition probability p 12 (s, t) can be estimated by plugging ( 4) and ( 5) into Equation 3. In the limit case of no presmoothing, the presmoothed Aalen-Johansen estimator reduces to the time-honored Aalen-Johansen estimator, which has become the standard tool for estimating the transition probabilities in Markovian processes.Moreira et al. (2013) derive the consistency of the PAJ estimator which may be much more efficient than the original AJ estimator.
The original and the presmoothed AJ estimators are consistent in Markov models.If the Markov property assumption is violated, then the consistency of the time-honored Aalen-Johansen estimator and of its presmoothed version cannot be ensured in general.Alternative methods that do not rely on the Markov assumption are presented below.

Kaplan-Meier weighted estimator
Recently Meira-Machado et al. (2006) verified that in non-Markov situations, the use of Aalen-Johansen estimators to empirically estimate the transition probabilities may be inappropriate.These authors propose, in the scope of the illness-death model, alternative "Markov-free" estimators for the transition probabilities, which do not rely on the Markov assumption.The idea behind estimation is to use the Kaplan-Meier estimator pertaining to the distribution of the total time to weight the bivariate data.The proposed estimator (Kaplan-Meier weighted estimator, KMW) is given by where W i (and W 1i ) are Kaplan-Meier weights attached to T i (respectively, Z i ) when estimating the marginal distribution of T (respectively, Z) from ( T i , ∆ i )'s (respectively, ( Z i , ∆ 1i )).The expression for the Kaplan-Meier weights, W i , is given by et al. (2006) derive large sample properties of these estimators which may be generalized to more complicated non-Markov processes.

Kaplan-Meier presmooth weighted estimator
Recently, Amorim et al. (2011) propose a modification of estimator ( 6)-( 8) based on presmoothing, which allows for a variance reduction in the presence of censoring.Basically, this method is obtained by replacing the censoring indicator variables in the expression of the Kaplan-Meier weights by a smooth fit of a binary regression.In this estimator (Kaplan-Meier presmooth weighted estimator, KMPW) the presmoothed Kaplan-Meier weights are given by Here, m(x, y) ) belongs to a parametric (smooth) family of binary regression curves, e.g., logistic.Our package provides the results assuming that m denotes a logistic regression model (KMPW).In practice, we assume that m(x, y) = m(x, y; ´) where ´is a vector of parameters which typically will be computed by maximizing the conditional likelihood of the ∆ 2 's given ( T 1 , T ) for those with ∆ 1 = 1.In the limit case of no presmoothing, the KMPW estimator reduces to the KMW estimator.Conditions under which both estimators are consistent is fully discussed in papers by Meira-Machado et al. (2006) and Amorim et al. (2011).In the latter paper the authors compare the performance of the presmoothed (semiparametric) estimator with the purely nonparametric estimator (without presmoothing) and concluded that the presmoothed estimator gains efficiency.The advantages of presmoothing are more clearly seen with an increasing censoring degree and at the distribution's right tail.In general, presmoothing introduces some bias in estimation, while reducing the variance.This bias component is larger when there is some misspecification in the chosen parametric model.Importantly, the validity of a given model for presmoothing can be checked graphically or formally, by applying a goodness-of-fit test.This implies that the risk of introducing a large bias through a misspecified model can be controlled in practice.

Inverse probability of censoring weighted estimator
To account for the influence of covariates, Meira-Machado et al. (2012) introduce estimation methods for the transition probabilities conditionally on current or past measures which we denote by X.The authors provide two competing nonparametric regression estimators for the conditional transition probabilities, p hj (s, t|X), both valid under certain regularity conditions even when the system is non-Markovian.The two estimators use different schemes of inverse censoring probability reweighting to deal with right censoring.In both estimators, local smoothing is done by introducing regression weights that are either based on a local constant (i.e., Nadaraya-Watson) or a local linear regression.To introduce these estimators, we need to introduce first the distribution function of C given X, G X .Let G X i denote the conditional distribution function of C | X = X i and let ĜX i stand for its estimator.This can be done using the estimator introduced by Beran (1981), with where N W 0i (x, a n ) are the Nadaraya-Watson (NW) weights, K 0 is a known probability density function and a n is a sequence of bandwidths.This estimator reduces to the so-known Kaplan-Meier (Kaplan and Meier 1958) estimator when all weights are equal.Then, the inverse probability censoring weighted estimators (IPCW) are given by where N W 1i (x, b n ) are NW weights as introduced above and ĜX i ( Alternatively local linear weights can also be introduced. An alternative approach that also accounts for the influence of covariates is based on the Lin, Sun, and Ying (1999) approach for the bivariate distribution function.Then, a different set of estimators (LIN) are given by where ĤX stands for the Kaplan-Meier estimator of the conditional distribution of C given X based on the ( Z i , 1 − ∆ 1i )'s.This estimator is defined in the same way as Ĝx .Here we assume that C is independent of (Z, T ) given X; this assumption does not exclude the possibility of dependent censoring.The performance of the two estimators has been evaluated through simulations, showing that they are valid even when the system is non-Markov or conditionally non-Markov.Simulation results show that the general performance difference between the two methods is quite small, and both methods perform quite well.However, one of the two approaches (the LIN-based one) has the drawback of occasionally providing nonmonotone curves for transition probabilities which are indeed monotone and, therefore, its practical use is less recommendable.

Location-scale estimator
Other estimators were proposed to estimate the transition probabilities.A valid estimator was provided by Keilegom et al. (2011).This methodology assumes that the vector of gap times (Z, Y ), where Y = T − Z, satisfies the nonparametric location-scale regression model, allowing for the transfer of tail information from lightly censored areas to heavily ones.An automatic bandwidth procedure was proposed by Meira-Machado, Roca-Pardiñas, Keilegom, and Cadarso-Suárez (2013) for this methodology.
Consider the nonparametric location-scale regression model (LS) Y = m(Z) + Ã(Z)ϵ, where the functions m and Ã are 'smooth', and ϵ is independent of Z.Under this model, the authors propose a nonparametric estimator of the distribution of the error variable, F ϵ , to introduce nonparametric estimators for the transition probabilities.They use a Kaplan-Meier estimator of F ϵ based on the ( Êi , ∆ i )'s (where Êi = ( which is the key for the construction of an estimator for the conditional distribution of the second gap time, ).The location and scale functionals are estimated using an extension of the Beran (1981) estimator, which copes with censoring in the first gap time.Then, estimators for the transition probabilities can be obtained from the following expressions: , where F 1 (⋅) is the marginal distribution of the first gap time, which we may estimate by the Kaplan-Meier estimator based on the ( Zi , ∆ 1i )'s.Meira-Machado et al. (2013) suggest that the transfer of tail information may improve the estimation of the transition probabilities especially in points where the uncensored information is scarce.The authors compared the location-scale method with the estimator by Meira-Machado et al. (2006) in several scenarios.It was found that when the deviation from the location-scale model was only minor, the location-scale method outperforms the Kaplan-Meier weighted estimator (Meira-Machado et al. 2006).However, when the model deviates a lot from a location-scale model, the later method becomes better.Another drawback of the location-scale model is that this method can only be used in the progressive three-state model.

Occupation probabilities
Another important target in multi-state modeling is the estimation of the state occupation probabilities.For the illness-death model there are in essence three state occupation proba-bilities to calculate, p 11 (0, t), p 12 (0, t) and p 13 (0, t).Datta and Satten (2001) show that these quantities can be estimated using Aalen-Johansen estimators even when the process is not Markov.Though all methods introduced in the previous sections provide valid estimators for these quantities, the Markovian approaches (AJ and PAJ) are recommended.

Package description
The TPmsm software package contains functions that calculate estimates for the transition probabilities.As mentioned in Section 2, this software package can be used to implement seven methods (AJ, PAJ, KMW, KMPW, IPCW, LIN and LS).This software package is intended to be used within the statistical software program R (R Core Team 2014).TPmsm is composed of several functions that allow users to obtain estimates and plots of the transition probabilities.Table 1 provides a summary of some of the functions in the package.Details on the usage of these functions can be obtained with the corresponding help pages.setThreadsTP Specifies the number of threads used by default in parallel sections.
setPackageSeedTP Set the initial package seed.
Table 1: Summary of functions in the package.
It should be noted that to apply the methods described in Section 2 one needs the following variables: time1, event1, Stime and event.A single covariate can also be included (they are necessary only for IPCW and LIN methods).The variable time1 represents the observed time in State 1 ("healthy"), and event1 the corresponding status/censoring indicator (if the survival time is a censored observation, the value is 0 and otherwise the value is 1).The variable Stime represents the total survival time (time to the absorbing state).If event1 = 0, then the total survival time is equal to the observed time in State 1.The variable event is the final status of the individual (takes the value 1 if the final event of interest is observed and 0 otherwise).

Data generation
Users may use the function dgpTP to generate data from the illness-death model.We assume that all individuals are in the "healthy" state at time t = 0. Therefore, the patient's history (or course) may be divided into two groups according to whether the disease occurred (that is, passing through State 2) (1 → 2 → 3) or not (1 → 3).We separately consider these two possible subgroups of individuals.For the first subgroup of individuals, the successive gap times (Z, T − Z) can be simulated from two of the most known copula functions: the Farlie-Gumbel-Morgenstern copula with exponential marginals and the bivariate Weibull distribution.
In the following, using the dgpTP function we will simulate data from the illness-death model using Gumbel's bivariate exponential distribution (dist = "exponential") with unit exponential margins (dist.par= c(1, 1)).The parameter ¹ controls for the amount of dependency between the gap times (Z, T − Z).Theoretical correlation between the gap times can be obtained using the corrTP function.For the second subgroup of individuals (those that go directly from State 1 to State 3), the corresponding survival time is simulated according to an exponential with rate parameter 1.The computation and the implementation of the proposed estimator involves the construction of pointwise confidence intervals by means of a bootstrap approach and in some cases the choice of an appropriate bandwidth.Thus, some of the methods implemented in package TPmsm can be computationally demanding.To obtain the point estimation and the pointwise confidence intervals, efficient algorithms were developed and implemented in the C programming language.The most computationally demanding parts of the code, namely those that involve the bootstrap and cross-validation techniques, were parallelized by means of the OpenMP API.This should considerably increase performance on multi-core/multithreading machines.To ensure the reproducibility of the results reported in the paper, two threads were considered (setThreadsTP(2)).The random number generator with multiple independent streams ((L'Ecuyer 1999),(L'Ecuyer, Simard, Chen, and Kelton 2002) and (Karl, Eubank, Milovanovic, Reiser, and Young 2014)) was implemented for parallel computation of uniform pseudorandom numbers.Package TPmsm own implemementation of a random number generator makes it independent of R, requiring a different function for defining a seed.The function setPackageSeedTP requires a vector of six integers.library("TPmsm"); setThreadsTP(2); seed <-c(2718, 3141,5436,6282,8154,9423); setPackageSeedTP(seed); sim_data_exp <-dgpTP(n = 1000, corr = 0, dist = "exponential", dist.par= c(1, 1), model.cens= "uniform", cens.par= 3, state2.prob= 0.5); This input command will simulate 1000 observations (n = 1000) assuming no correlation in Gumbel's bivariate exponential distribution (corr = 0), using an independent uniform censoring time (model.cens= "uniform"), according to model U (0, 3) (cens.par= 3).The use of corr = 0 in Gumbel's bivariate exponential distribution leads to independent gap times and therefore to Markov data.The proportion of transitions into State 2 is given by the argument state2.prob(a value of 1 corresponds to the progressive three-state model).
To obtain the estimates for the methods proposed in Section 2 we can use the functions shown in Table 1.As in the simulation by Amorim et al. (2011)  P(Z≤s,T >s) = 0.666.The following two input commands provide the estimates for the AJ and PAJ methods.Since the process is Markovian these are the recommended approaches.With these input commands we obtain the estimates for the transition matrix together with 95% (conf.level= 0.95) pointwise confidence intervals (conf = TRUE) using 1000 bootstrap replicates (n.boot = 1000).The construction of the pointwise confidence intervals is obtained by randomly sampling the n items from the original data set with replacement.This can be achieved using the percentile bootstrap interval (method.boot= "percentile") or using the basic bootstrap interval (method.boot= "basic").By default all functions use the percentile bootstrap method (Davison and Hinkley 1997).
In addition to the numerical results graphical output can also be obtained.This will be shown in the next section using two data sets: the widely used and well-known colon cancer data and data from a bladder cancer study.Details about these data sets are given below.

Examples of application
To illustrate our estimators we consider two real data sets.One of these data sets comes from the well-known colon cancer study which is freely available as part of the R survival package (Therneau and Grambsch 2000;Therneau 2014).In addition to this data set we also use data from a bladder cancer study (Byar 1980) conducted by the Veterans Administration Cooperative Urological Research Group.

Colon cancer data
For illustration, we apply some of the proposed methods of Section 2 to data from a large clinical trial on Duke's stage III patients, affected by colon cancer, that underwent a curative surgery for colorectal cancer (Moertel, Fleming, McDonald et al. 1990).In this study, some of these patients have residual cancer, which leads to disease recurrence and death (in some cases).From the total of 929 patients, 468 developed a recurrence and among these 414 died.38 patients have died of causes unrelated to their disease and without evidence of recurrence.
The remaining 423 patients contributed with censored survival times.For each individual, an indicator of his/her final vital status (censored or not), the survival times (time to recurrence, time to death) from the entry of the patient in the study (in days), and a vector of covariates including age (in years) and recurrence (coded as 1 = yes; 0 = no) were recorded.The covariate recurrence is a time-dependent covariate which can be expressed as an intermediate event and modeled using the illness-death model with states "alive and disease-free", "alive with recurrence" and "dead".
By including covariates depending on the history (using a Cox proportional hazards model), we verified that the mortality transition for recurring patients is affected by the time spent in the previous state (p value < 0.001).This allowed us to conclude that the Markov assumption may be unsatisfactory for the colon cancer data set and that, consequently, Aalen-Johansen type estimators should not be used.Thus, in this section we illustrate the use of two "Markovfree" estimators (KMW and KMPW) as well as two additional estimators (IPCW and LIN) that were proposed to estimate the transition probabilities conditionally on current or past covariate measures such as age.
Below is an excerpt of the data with one row per individual.We computed the estimated values for the transition probabilities p hj (s, t) for several pairs (s, t), s < t.For illustration purposes we only report the estimated values of p hj (365, 1096) (one year and three years) for the KMW and KMPW methods with 95% bootstrap confidence intervals.

Example of application: Bladder cancer study
The methods described in Section 2.6 are illustrated using data from a bladder cancer study (Byar 1980) conducted by the Veterans Administration Cooperative Urological Research Group.In this study, patients had superficial bladder tumors that were removed by transurethral resection.Many patients had multiple recurrences (up to a maximum of 9) of tumors during the study, and new tumors were removed at each visit.For illustration purposes we re-analyze data from 85 individuals in the placebo and thiotepa treatment groups; these data are available as part of the R survival package.Here, only the first two recurrence times (in months) and the corresponding gap times, Z and Y = T − Z, are considered.Thus, we have a progressive three-state model with state "alive and disease-free", "first recurrence" and "second recurrence".From the total of 85 patients, 47 relapsed at least once and, among these, 29 experienced a new recurrence.
For large values of s and t, the transition probabilities p hj (s, t) will be difficult to estimate in a completely nonparametric way.This will be particularly true in situations where censoring percentages are high as for our data set for which we have a total amount of censoring of 66%.The location-scale method is appropriate for the bladder cancer data since this methodology is mainly relevant for estimation in the right tail of the distribution where the censoring effects are strong at those points (uncensored information is scarce).
We will calculate estimates for the transition probabilities in several points and plot these estimates.This will be done using the function transLS.
data("bladderTP", package = "TPmsm"); head(bladderTP); We computed the estimated values for the transition probabilities p hj (s, t) for several pairs (s, t), s < t.For illustration purposes we only report the estimated values of p hj (3, 8) for the LS method with 95% bootstrap confidence intervals.The success of the LS method greatly depends on the choice of an appropriate bandwidth.The selection of the optimal bandwidth is highly computationally intensive, but is crucial to the success of the location-scale method.

Conclusion
This paper discusses the implementation of some newly developed methods for the transition probabilities in the illness-death model in an R package.The TPmsm package uses seven nonparametric and semiparametric estimators.One of these estimators is the Aalen-Johansen estimator (Aalen and Johansen 1978) under the assumption of a Markovian data generating process.A modification of the Aalen-Johansen estimator (Moreira et al. 2013), based on a preliminary estimation (presmoothing) of the censoring probability for the total time, given the available information is also implemented.This method allows for a variance reduction in the presence of censoring, in particular for situations with high percentages of censored total time among the uncensored subjects in State 1.
If there is no evidence against the Markov condition then the time honored Aalen-Johansen estimator and its presmoothed version will be preferred.If the Markov property is violated, then the consistency of these estimators cannot be ensured in general.Exceptions to this are the estimator for the occupation probabilities.Alternative estimators of the transition probabilities not relying on the Markov condition were recently proposed (Meira-Machado et al. 2006;Amorim et al. 2011) and are implemented in the package.As a drawback, these alternative methods will suffer from a larger variance in estimation, particularly when the sample size is small and there is a large censoring degree.One alternative method for these scenarios was provided by Keilegom et al. (2011) for a progressive three-state model.The key of this methodology is the transfer of tail information from lightly censored areas to heavily ones.
The package also implements two methods that account for dependent censoring and allow for the inclusion of covariates.These two approaches are free from the Markov assumption.The functions implementing these methods use a kernel density and a bandwidth.We believe that the choice of the kernel density has relatively little impact on the estimation results.However, the use of different bandwidths might have a substantial effect on the performance of the estimators.To this end we implemented the use of the dpik function which is available from the R KernSmooth package (Wand 2014).It might be worthwhile to include other options and to investigate their impact on the estimation results.
A function called TPmsmOut can be used to convert an object of class 'data.frame'with the structure of the data input as described in Section 3 to the structure of the data input used in the p3state.msmpackage.Essentially, this involves a transformation of some variables and a renaming of other variables.With this function users may connect the TPmsm package with the p3state.msmpackage and perform Cox-type multi-state regression.
We plan to constantly update the TPmsm package to improve its limitations and to cope with other estimators.

Function DescriptiondgpTP
Generate data from an illness-death model (based on some known copula functions).By default returns a data set of class 'survTP'.corrTPCorrelationbetween the bivariate times for some copula distributions.survTP Set up adequate data set of class 'survTP' for implementing all the methods.transAJ Aalen-Johansen (AJ) estimates for the transition probabilities.transPAJ Presmoothed Aalen-Johansen (PAJ) estimates for the transition probabilities.transKMW Kaplan-Meier weighted (KMW) estimates for the transition probabilities.transKMPW Kaplan-Meier presmoothed weighted (KMPW) estimates for the transition probabilities.transIPCW Inverse probability of censoring weighted (IPCW) estimates for the transition probabilities.transLIN LIN-based (LIN) estimates for the transition probabilities.transLS Location-scale (LS) estimates for the transition probabilities.plotPlots for the transition probabilities.

Figure 3
Figure 3 depicts the KMW estimates of p 12 (s = 365, t) as functions of the time (for a fixed value of s = 365) together with a 95% pointwise confidence band based on simple bootstrap.The estimates shown in the main curve indicate that this probability increases until around time t = 600 and afterwards decreases.

Figure 2 :
Figure 2: Transition probability estimates using the KMW method.Colon cancer data.

Figure 3 :
Figure 3: Transition probability estimates, with bootstrap confidence bands, using the KMW method.Colon cancer data.

Figure 4 :
Figure 4: Evolution of the transition probability p 11 (365, 1096) along the covariate age with 95% bootstrap confidence bands based on the IPCW method.Colon cancer data.

Figure 5 :
Figure 5: Evolution of the transition probability p 12 (365, 1096) along the covariate age with 95% bootstrap confidence bands based on the IPCW method.Colon cancer data.

Figure 6 :
Figure 6: Evolution of the transition probabilities p hj (365, 1096) along the covariate age, based on the IPCW method.Colon cancer data.

Figure 7 :
Figure 7: Transition probability estimates using the LS method.Bladder cancer data.

Figure 8 :
Figure 8: Transition probability estimates, with bootstrap confidence bands, using the LS method.Bladder cancer data.
. There are five different transition probabilities to estimate: p 11 (s, t), p 12 (s, t), p 13 (s, t), p 22 (s, t) and p 23 (s, t).However, only three of them need to be estimated since the two other transition probabilities can be obtained from the following relations: p 11 (s, t) + p 12 (s, t) + p 13 (s, t) = 1 and p 22 (s, t) + p 23 (s, t) = 1.
and Moreira et al. (2013)we are going to obtain estimates for transition probabilities at values s = 0.5108 and t = 0.9163.The true values for the transition probabilities are: p 11 (s, t) = Each line represents the information from one individual in the study.The variable time1 denotes the sojourn time in State 1 whereas Stime is the total time of survival.event1 and event are the corresponding indicator statuses.Among the first six individuals, only individual represented by line 2 remains alive (and without having had a recurrence) at the end of the study.All the remaining individuals had a recurrence and died before the end of the study.For example, the individual represented by line 1 had a recurrence at day 968 and died at day 1521.Note that time1 < Stime means that a transition from State 1 to State 2 (i.e., recurrence) occurred.