Estimation of semi-Markov multi-state models: a comparison of the sojourn times and transition intensities approaches

2.1 Approach I: sojourn times

As defined above, an SMP, X(t), can be constructed via the sequence { ( J n , T n ) } n ≥ 0 , of states and jump times respectively. The underlying parameters of this construction involve the transition probabilities of the embedded chain, p _ij presented in (1), as well as the distributions of the sojourn times for each transition i → j, presented in (2).

The transition probabilities of the embedded chain are often organized in a matrix, P = [p _ij ] (where we set p _ii  = 0). Note that P restricted to the non-absorbing states is a stochastic matrix. Further, assuming the sojourn time distributions are continuous, they are often represented in different forms including the probability density function, the survival function or the hazard rate function. We now present the details.

The probability density function of the sojourn time is

The corresponding survival function is

(Note that S _ij (t) is a decreasing function, that is S _ij (0) = 1 and lim_t→+∞S _ij (t) = 0). The hazard function which is often thought of as the probability that a jump occurs in a specified interval (t, t + Δt) given no jump before time t, is

Here, note that by definition of conditional probability we have

It is also useful to define the probability of staying in a current state i for at least t time units. We call this the survival function of the waiting/holding time in state i, and denote it via

2.2 Approach II, intensity transition functions

The first approach required specification of the parameters using two types of objects, transition probabilities of the embedded Markov chain, and the distribution of sojourn times given a transition i → j. The second approach which we present now is more succinct in that it only requires one type of object: intensity transition functions. These functions are defined via

and are similar in nature to hazard rates. However, they should not be treated as hazard rate functions of a probability distribution. They rather indicate the instantaneous probability of making a transition from state i to state j after spending t time units in state i since the last transition. However summing up over all target states j, we do obtain a hazard rate of a probability distribution which we denote via α ˜ i ( t ) = ∑ j ≠ i α ˜ i j ( t ) . This is the hazard rate of the waiting/holding time in state i. We can use this approach to obtain an alternative expression to S _i (t) of (7).

Using Approach II, more formally, the meaning of the intensity transition functions is

where using the notation defined previously, τ_N(t)+1 is the time elapsed since the last transition, measured at time t. The conditional probability in (10) is conditional on the fact that the last transition was at time t into state i and up to time t + u there have not been any further transitions. Given this condition, it indicates the instantaneous probability of: (i) making a transition at time t + u. (ii) Making the transition into state j.

It can be shown that specification of the intensity transition functions α ˜ i j ( ⋅ ) completely specify the probability law of an SMP. See [16] for a formal description where one can also use the notation ℋ(t⁻) to define the history of the process just before time t, i.e., up to t⁻. This can formally be defined as the sigma algebra in a filtration associated with the stochastic process. Related notation often used in the literature is the (conditional on history) transition probability to be in state j at time t + u, after being in state i in time t. This is often denoted by P i j ( t , t + u ) = ℙ ( X ( t + u ) = j | X ( t ) = i , H ( t − ) ) .

We note that Approach II can be viewed as a collection of competing risk models where at the instance of arriving to a new state, a new set of cause-specific hazard rates is defined. See for example [17] where a competing risk model is presented as a multi-state model and the cause-specific hazard rates terminology is used.

2.3 Relations between the two approaches

We now wish to show how for the same SMP, one can use either the parameters of Approach I or the parameters of Approach II and interchange between them. A key relationship is the following

It is established using the conditional probabilities, defined in (1), (4), (7), and (5). This key relationship also yields an interpretation of the cumulative incidence function for the i → j transition which we denote via CIF _ij (⋅). This is a common measure used in the field of competing risks, see for example [18]; that determines the probability of transitioning into j before or at time t. By rearranging and integrating both sides of (11) we obtain the following representation of the cumulative incidence function.

We can now convert between the parameterizations of both approaches as follows.

Approach I → Approach II. Given α _ij (⋅) and p _ij , obtain α ˜ i j ( ⋅ ) as follows.

This follows directly from (11).

Approach II → Approach I. Given α ˜ i j ( ⋅ ) obtain p _ij and α _ij (⋅):

First we have,

This can also be obtained from (12) by taking t → ∞. Once p _ij values are at hand we again use (11) and isolate f _ij (t) to obtain,

from which we can obtain α _ij (t) using (6) in the standard manner.

2.4 Examples

We now present three examples illustrating the relationship between the two approaches.

Observe that α _ij (t) is independent of t and independent of the target state j.

• CTMC Approach II: Define a ℓ × ℓ generator matrix Q (non-negative entries on the off-diagonal, negative entries on the diagonal, and rows sum to 0). Denote the entries via q _ij . Then treat the process as an SMP using Approach II with

  (17) α ˜ i j ( t ) = q i j ,

Observe that α _ij (t) is independent of t.

It is well known from the theory of CTMCs the two parameterizations are equivalent.

CTMC Approach I → CTMC Approach II:

CTMC Approach II → CTMC Approach I:

We now verify that these transformations directly agree with the relationship between Approach I of the SMP and Approach II of the SMP. Using (13) and (16) we have,

Going in the other direction using (14) and remembering that since Q is a generator matrix, q i i = − ∑ k ≠ i q i j < 0 , we have,

Now using (15) we have,

This is an exponential density and hence it has a constant hazard rate α _ij (t) = λ _i .

We now show that this condition is necessary for having constant transition intensity functions. To see this, use (13) with α _ij (t) = λ _ij where for some j₁ and j₂, λ _{ij 1}  ≠ λ _{ij 2} . Now,

Since λ _{ij 1}  ≠ λ _{ij 2} we are not able to cancel out the exponent as was done in (18).

Interestingly, the converse does not hold. The previous example showed that time independent transition intensity functions always yield a CTMC and hence constant (Approach I) hazard rates.

Say now that we are using Approach I and parameterize all of the sojourn time distributions as Weibull distributions, where for transition i → j we have respective scale and shape parameters μ _ij and η _ij . Now, if we were to consider the Approach II representation, then from (13) the transition intensity functions are

where we use the notation K _ij (t) to represent all of the components of α ˜ i j ( t ) excluding t η i j − 1 . Now the form of K _ij (t) can determine if the Approach II representation is of a Weibull type or not. A Weibull type representation will follow if K _ij (t) is independent of t, otherwise it is not.

We now see that if for a fixed i and every j₁ and j₂ with j₁ ≠ j₂ we have that

then, K _ij (t) is independent of t and we obtain a Weibull type Approach II representation with,

and hence the scale parameter is modified to p _ij ^−1/ηijμ _ij and the shape parameter keeps the same form. Otherwise (if there exists j₁ and j₂ such that (21) does not hold), then α ˜ i j ( t ) cannot be of the Weibull type as in (22).

Going the other way, assume we are using Approach II with Weibull type transition intensity functions,

Then by integrating (14) we can in principle obtain p _ij . It turns out that if for a given state i, the shape parameters of all the transition intensity functions is the same (denote it via η _i ) then the integration yields an explicit form,

otherwise, there is not an explicit solution. In such a case where η _ij  = η _i for all j, we also have from Equation (15),

Further, if μ _ij  = μ _i for all j then p _ij of (23) reduces to ℓ i − 1 where ℓ _i  < ℓ is the number of possible target states from state i. Note that this form is not specific to the Weibull type intensity transition function, but will hold whenever the intensity transition functions out of state i are the same.

However, with the Weibull case we can continue further and (24) can be represented as,

which denotes a Weibull distribution with same shape parameter as the transition intensity function and a scale parameter equal to ℓ i − 1 / η i μ i .

3.1 Likelihood for Approach I

The form of ℒ^(h) based on ℋ^(h) is,

where for brevity we omit the (h) superscripts for the state and sojourn time sequences {J} and {τ}. To further understand (27) consider a recursive construction where each transition without censoring based on J_k−1 → J _k with a sojourn time of τ _k has a likelihood contribution p J k − 1 J k f J k − 1 J k ( τ k ) . Further, in case of censoring the last censored transition has likelihood contribution S J N ( h ) ( U ( h ) ) based on the survival function of the holding time (7), corresponding to the last observed state J N ( h ) .

3.2 Likelihood for Approach II

We refer to the key relationship (11) and replace p _ij  f _ij (u) in (27) with α ˜ i j ( u ) S i ( u ) to obtain,

where for brevity we omit all superscripts (h). Further, we manipulate the likelihood contribution of each subject as follows:

In the first line we use the representation in (9). For the second line, by increasing the summation from N to N + 1 we are able to remove δ by defining τ_N+1 = U, J_N+1 ≡ −1, and using the fact that α ˜ i j ( u ) ≡ 0 in cases where i is an absorbing state and 0⁰ ≡ 1.

Now, we are able to separate ℒ^(h) into the form

following similar ideas as those initially presented in [21]. To achieve such a separation, it is useful to define the indicators

We may now manipulate (29) to obtain

The form of (26), (30), and (31) allows us to treat each transition separately as if each transition intensity transition function α ˜ i j ( ⋅ ) has its own set of parameters.

3.3 Parametric forms and covariate information

In carrying out inference, we assume a parametric form for α _ij (⋅) in Approach I or α ˜ i j ( ⋅ ) in Approach II. We also allow for covariate information where the natural way to introduce covariates in a multi-state model is a Cox-like proportional hazard model which can be defined by using either the hazard function of the sojourn times or by intensity transitions, see [22]. Hence for a vector of covariates Z, we have

for Approach I. Further we have

for Approach II.

Here α i j , 0 ( θ i j ) ( ⋅ ) and α ˜ i j , 0 ( θ ˜ i j ) ( ⋅ ) are the baseline functions for the hazard rate and transition intensity functions respectively. They each follow a parametric form, determined via θ _ij and θ ˜ i j respectively. Further, β _ij and β ˜ i j are the regression parameters for transition i → j associated with the covariates of the subject. Observe that there are major differences in the interpretation of the regression coefficients β _ij of Approach I and β ˜ i j of Approach II. The former determine the hazard ratios dealing with sojourn times and bear no effect on the actual transitions of the SMP. The latter, determine the hazard ratios dealing with risks and also affect the instantaneous rate of transitioning.

When estimating the parameters for such a model, with Approach I, we also need to estimate p _ij , whereas with Approach II this is not needed as p _ij is implicitly determined via (14). Interestingly, via Approach II, we have that the covariates implicitly affect p _ij ,whereas with Approach I, we may wish to set p _ij (Z) using multinomial regression, however it isn’t common practice.

3.4 Contrasting inference for Approach I and Approach II

The key difference in the likelihood expressions between Approach I (27) and Approach II (30) and (31) is in computational tractability. As noted by [23] the optimization of the likelihood for Approach I could be challenging when the number of covariates and parameters is large due to a complex form of the sojourn time distribution. This is especially the case when the number of variables is large compared to the size of the data. It should further be noted that modeling using Approach I requires the embedded chain parameters in addition to the parameters involved in (32). This means that for a model with ℓ states, there may be up to ℓ² − ℓ more parameters to estimate when using Approach I, in comparison to Approach II, where all the parameters appear in (33).

With Approach II, the decoupling in (30) allows to optimize the likelihood for each transition i → j separately if each hazard intensity transition function (33) has its own set of parameters θ ˜ i j and β ˜ i j . In such a framework, the full likelihood of the SMP can be maximized separately by maximizing the likelihood of the different separate two-state models, see [22]. Also, the inference of the SMP can be easily handled by using survival models like the Cox model as the separate likelihood arising in survival analysis.

Hence Approach II enables fitting SMP models using survival analysis methods and software, see for example [24], [25], [26]. This also applies to the case where extensions to the proportional hazard assumption are needed. The forms in (32) and (33) are clearly based on a proportional hazard assumption with fixed covariates where there is only a multiplicative effect on the baseline hazard functions independent of time. An extension is to consider time-dependent covariates Z(t) as in [17], [27]. Further, in both proportional models the effect of variables is assumed to have a linear (or log-linear) functional form. A more flexible model is to consider non-linear effects using smooth functions of the covariates such as the generalized additive model (GAM) presented by [28]. With such extensions, using Approach II is much more straightforward than using Approach I due to the decoupling and because there is ample statistical software available for survival analysis.

One may also wish to incorporate random effects in the models using frailty models. In this case, both (32) and (33) could also include random effect to take into account the correlation between observations/subjects. This has been handled in [10] for Approach II by exploiting the decoupling of the likelihood. Such random effects may also be incorporated in Approach I, however the computational effort will be greater.

4.1 Progressive three-state model for the stanford heart transplant data

As an illustrative example, we revisit the analysis of the Stanford Heart Transplant data, freely available through the survival package. A full description of the data can be found in [35]. The dataset presents the survival of patients on the waiting list for the Stanford heart transplant program. We analyze the data similarly to [36] which use their R package p3state and propose an illness-death model similar to Figure 2. Their states represent “alive without transplant”, “alive with transplant”, and “dead”.

We consider the same dataset, where patients’ records from October 1967 (start of the program) until April 1974 includes the information of 103 patients where 69 received heart transplants, and 75 died. Three fixed covariates for each patient are available. These include the age at acceptance (age), the year of acceptance (year), and previous surgery (surgery: coded as 1 = yes; 0 = no).

We first use the data without taking into account any covariate effects to show that an exponential distribution on the baseline hazard rate of the sojourn time, α _ij (⋅) does not produce constant intensity transition functions, α ˜ i j ( ⋅ ) . This point was shown theoretically in Example 2, of Section 2.4 and we now illustrate it using data and the SemiMarkov package. Using the package we obtain both α _ij (⋅) and α ˜ i j ( ⋅ ) for (i,  j) = (1, 2) and (i,  j) = (1, 3). The results are in Figure 4 where we see that while α _ij (⋅) are constant, α ˜ i j ( ⋅ ) are clearly not constant across time. These curves in fact follow a form like (19).

Figure 4:

Baseline hazard rate of the sojourn times (Approach-I: α _ij (⋅)) and baseline hazard rate of the semi-Markov process (Approach-II: α ˜ i j ( ⋅ ) ) for transition 1 → 2 (top plot) and 1 → 3 (bottom plot).

We now continue to analyze the data fully with Approach I using the SemiMarkov package where we compare different parametric forms and covariates in the model. Table 1 presents the inference results of 6 different models investigated depending on the distribution and covariates chosen for each transition. For each model we present the β ˆ i j estimates and their estimated standard error. We also point out that the SemiMarkov package yields p-values for the associated Wald tests. We also present the estimates of the transition probabilities of the embedded chain, p ˆ i j . It is evident that these are outputs from likelihood estimation as they don’t fully agree across models (they are not proportion estimates). For the final two models, Weibull + Select and Weibull/Exp + Select are based on model selection of the significant covariates with more details in the vignette, see [12]. Observe that the final model, Weibull/Exp + Select, incorporates two different distributions.

For choosing between the different proposed models, we use the expected Kullbak-Leibler (EKL) risk as in [10]. This is done via the Akaike information criterion (AIC) which is an estimate of the EKL risk in a parametric framework using maximum likelihood method, see [37], [38]. For each of the 6 models described above, we obtain the AIC using output from SemiMarkov and present it in the first part of Table 2.

The best model according to the AIC is given by the Exponentiated Weibull model. However, due to convergence issues during the optimization (see the vignette [12], we select the second best result given by the Weibull/Exp + Select model as the most appropriate model for this dataset according to AIC using Approach I.

Moving onto inference using Approach II, we investigate and estimate different models by specifying different proportional intensity transition models as defined in Equation (33). We use flexsurv to estimate α ˜ i j ( ⋅ ) for several models which we present in the second part of Table 2. The models Exponential, Weibull, Gamma, and Generalized Gamma are all estimated with all three covariates in the model. The additional two models, Weibull + Select, and Generalized Gamma + Select, have a reduced set of covariates with a procedure consistent with that used for Approach I above. More details are in the vignette (see [12]. The AIC based performance of each of these models is presented in the second part of Table 2.

The best model according to the AIC is Generalized Gamma + select. In this model only the age covariate is used for transition 1 → 2, only the year covariate is used for transition 1 → 3, and both the surgery and age covariates are used for transition 2 → 3. For this dataset, the best model based on Approach I presents a better AIC value (1639.73) in comparison to the best model based on Approach II (1701.03).

As discussed in Section 3.3 the regression coefficients of Approach I and those of Approach II do not have the same interpretation. Indeed for the hazard rate of the sojourn time, the regression coefficients can be interpreted in terms of relative risk on the waiting time (i.e., given an i → j transition). For example the positive coefficient for ‘surgery’ for the 1 → 3 transition, as appearing in Table 1, can be interpreted as an indication that death without getting a heart transplant is less favorable for those that had previous surgery. While the hazard rate of the SMP based on transition intensity functions can be interpreted as the subject’s risk of passing from state i to state j. Interpretation of the regression coefficients using both Approach I and Approach II is further discussed in the second example below.

Since a key object of a Semi-Markov process is the transition intensity function, generally called the hazard rate of the SMP, it is also useful to visualize estimates of these functions based on both Approach I and Approach II. Further insight may be gained by plotting these against non-parametric estimators. For this, we used the Breslow estimator [39]; which yields an estimate for the baseline cumulative risk function of a cox-regression type model of each transition. Figure 5 plots the cumulative baseline (all covariates are set to 0) estimated transition intensity functions, estimated via Approach I and Approach II, against the Breslow estimator. These cumulative functions are based on the parameter estimates for the best model of Approach I (Weibull/Exp + Select) vs. the best model of Approach II (Generalized Gamma + Select). In this case, the plots in Figure 5 hint at a better fit (closer to the non-parametric estimation) of the SMP using Approach I in accordance with the results of the AIC values discussed above.

Figure 5:

Cumulative of the baseline hazard rate of the SMP (intensity transition function) estimated from best AIC model of Approach-I and Approach-II. A semi-parametric estimation is presented as benchmark.

4.2 Reversible semi-Markov model for the asthma control data

We revisit the analysis of the asthma control data which has been used to illustrate the SemiMarkov R package in [29]. See also the related papers using the same data, [40] and [41]. The data consists of the follow-up of severe asthmatic patients consulting at varied times according to their perceived needs. At each visit, several covariates were recorded and asthma was evaluated. The available data (asthma) from the SemiMarkov package presents a random selection of the 371 patients with at least two visits. In this illustration, we only use the body mass index (BMI) covariate coded: 1 if BMI ≥25 and 0 otherwise. Similarly to the analysis in [29]; we use a reversible three-state model as presented in Figure 6. This investigates the evolution of asthma based on, “optimal control” (State 1), “sub-optimal control” (State 2) and “unacceptable control” (State 3). More details about concepts of control scores can be found in [41]. Observe that in this model, none of the states are absorbing and hence the model is said to be reversible.

Figure 6:

The reversible three-state model.

We note that this data is based on scheduled visits every three months or according to patients’ perceived needs [40], and is thus potentially interval censored since the exact times of state transitions are approximated to be at visit times. This aspect of the analysis was not the focus of [29] and the related papers. While dealing with such interval censoring is important we do not focus on it here further and refer the reader to [42, 43].

For illustration purposes, we concentrate on a Weibull model as proposed in [29]. We perform a first full parametric proportional Weibull model for the sojourn times (Approach I) with the BMI covariate for each transition. This yields significant results of BMI covariate for the transitions 1 → 3 and 3 → 1. Then, we decide to run a sparse model including BMI only for these transitions to facilitate the convergence of the model which can be difficult for a model with too many parameters as reported in [29]. For this new model, BMI regression coefficients remain significant for transitions 1 → 3 and 3 → 1 with β ˆ 13 = − 0.88 and β ˆ 31 = − 0.448 , and respective p-values 0.012, and 0.023. The fact that hazard ratio of the sojourn time associated with these covariates is less than unity (estimated coefficients are negative), indicates that BMI ≥25 generally lengthens the duration of the sojourn time in state 1 when making a 1 → 3 transition and generally lengthens the duration of the sojourn time in state 3 when making a 3 → 1 transition. This can also be interpreted as a decrease of the risk of leaving “optimal control” state to “unacceptable control” as well as a decrease of the risk of leaving “unacceptable control” state to “optimal control”. However, the magnitude of the estimated coefficients cannot be used to evaluate the change in the hazard ratio on the risk (recall Equations (32) and (33) and the differences between β _ij and β ˜ i j ). We may also visualize the effect of BMI covariate on the associated risks by plotting the hazard rate of the semi-Markov model (intensity transition functions), deduced from the sojourn times. See Figure 7 where we also plot estimated transition intensity functions using Approach II estimation, which we describe now.

Figure 7:

Hazard rate of the semi-Markov process for each transition from the sojourn time and intensity transition approaches.

Approach II presented in this paper offer a direct way to model the intensity transition functions and has the advantage of a split likelihood which facilitates an efficient optimization procedure. In the vignette available in [12]; we present the results of the full Weibull proportional semi-Markov model using this Approach. For illustration, we do this in two ways: (i) Optimizing the full likelihood. (ii) Splitting the likelihood for each transition and optimizing individually. As expected, both ways yield the same likelihood estimates. This highlight the advantage of the intensity transition based approach for overcoming potential computational issues in the optimization. Moreover, in some applications, some transitions of interest could then be further explored with any sophisticated survival model of choice as the estimation for each transition is performed separately. See for example [10], [24], [26], [44].

We augment the plots in Figure 7 with the estimated transition intensity functions using Approach II. Note that similarly to Approach I, the BMI effects are in the same direction. However, the interpretation of the regression coefficients is different. For the intensity transition based estimation (Approach II), the BMI regression coefficient estimates are β ˜ ˆ 13 = − 0.5 and β ˜ ˆ 31 = − 0.67 . Here, in contrast to the Approach I estimates, the exact magnitude of the coefficients may be interpreted as the change in the hazard ratios associated with BMI ≥25.