Simple analysis of non-Markov models: A case study on heart transplant data

Abstract

In transplantation studies, often several response measurements are collected for patients while they are on the waiting list. In this settings it is often of primary interest to assess whether the available history of a patient can be used for predicting patient survival as well as further performance on the list. In this work, we use a multi-state model approach to analyze the performance of patients described by their urgency status that changes in time while waiting for a new organ. We use the pseudo-value approach introduced by Andersen et al. (2003) and apply it for the Aalen-Johansen estimator of the state occupation probabilities since the transition probabilities were found to depend on the history. This approach allows us to study the impact of baseline information on the occupation probabilities treating the dependence on the history as a nuisance. It was found that the previous state, the current state and time from the moment of entering the waiting list had an impact on the future performance of the patient. Depending on those, patients were more likely to come back to the particular status in which they were before, die or get a transplant. To address the problem of those competing events, a multinomial approach was used for the next state given the previous state observed.

Keywords

measurement error correction Monte Carlo sampling nonlinear mixed models survival data

1 Introduction

1.1 Description of the data

In transplantation studies, the urgency status is measured over time for patients waiting for an organ transplant reflecting their changing disease state. In these settings, it is often of primary interest to investigate whether the available history on the urgency status of a patient can be used for predicting survival as well as the patient’s future status on the waiting list. We use a multi-state model approach to address this problem. In particular, we consider data from the Eurotransplant heart transplantation waiting list, where patients are classified as U (Urgent), HU (High Urgent), T (Transplantable) and NT (Not Transplantable) and due to the nature of the last option, these categories can be only partially ordered. The absorbing states were death (D), getting a transplant (TT) or removal from other reasons than the previous two (R). A graph depicting all possible transitions between the states is presented in Figure 1. The evaluation time was different for each patient and depended on the previous classification meaning that more severe patients were evaluated more frequently.

Figure 1

Source: Authors’ own.

1.2 Standard approaches using multi-state models

This type of data is typically analyzed using multi-state models. A standard assumption in this framework is that the process is Markov and homogeneous in time. That implies a constant transition (intensity) matrix that does not depend on time or past history of the process (but may depend on time-fixed baseline covariates), and allows the likelihood methods to be used. The results for the heart data suggested that the homogeneity assumption does not seem to be satisfied. A standard piecewise constant intensity model relaxing this assumption was still not sufficient, and therefore a non-parametric Aalen-Johansen (A-J) estimator was used to estimate the transition probabilities. Applied to our data, the multinomial approach model revealed a dependence on the history, in particular on the state visited before the current one. In such a non-Markov setting it is known that A-J estimators for the transition probabilities are not consistent and may produce biased results. Alternatively, Datta and Satten (2001) showed that for non-Markovian processes, the state occupation probabilities may still be consistently estimated, provided that data are not subject to stage-dependent censoring. Therefore, we employ the pseudo-value approach of Andersen et al. (2003) to directly model the effect of the baseline covariates on the occupation probabilities, treating the dependence on the history as a nuisance. It uses the pseudo-values from a jackknife statistic constructed for a non-parametric estimator of the transition probabilities. This approach has been applied for cumulative incidence functions in a competing risk model (Klein and Andersen, 2005; Klein, 2006), to the restricted mean survival time (Andersen, Hansen and Klein, 2004) in simple multistate models (Andersen and Klein, 2006), and further studied by Graw et al. (2009) providing more theoretical justification. However, the consistency and asymptotic normality was proven only in the special case of competing risk models. We extend the univariate pseudo-value approach for the multivariate model on occupation probabilities as suggested by Andersen et al. (2003) and present the simulation results for more general non-Markov models.

1.3 Interval-censoring problem

An additional complexity in the practical use of multi-state model designs is the possible non-ignorability of the observation process and the interval-censoring problem. Regarding the former issue, in our study more severely ill patients are monitored more close and the next sampling time is chosen on the basis of the current disease state. In particular, given a current state re-evaluation by Eurotransplant audit group is mandatory every fixed number of days and the length of this period depends on the particular current state. For example, for status HU it is 7 days and for state U, 28 days. Therefore, we have the so-called doctors care sampling scheme. Grüger et al. (1991) showed that this scheme is not informative when likelihood methods are used since the likelihood given by this examination scheme is proportional to the likelihood for the fixed examination scheme. For the non-parametric approach, the A-J estimator is also valid even if the process is not Markov since the censoring process is not stage-dependent and does not depend on any part of the history prior the current state, nor on the covariates.

The second issue of interval-censoring comes from the fact that the exact transition times are often only observed for the final but not for the intermediate states. In our analysis, it was taken into account only in the initial analysis assuming a Markov model, with parametric or eventually piecewise-constant intensities. That was possible due to the fact that for the multi-state model with continuous time the inference problem can be decoupled into several survival problems, and assuming the ignorability of the observation process, the likelihood for the whole observation of the trajectory can be written (Commenges, 2002). For the pseudo-value approach, as originally proposed in Andersen et al. (2003), we assumed the exact transition times ignoring the interval-censoring and allowing only for possible right censoring. As an alternative, in the interval-censoring situation, usually the kernel-based methods are applied. Datta and Sundaram (2006) suggested smoothed product limit estimators for the occupation probabilities for the current status data. Moreover, for handling the interval-censoring problem in simple non-homogenous Markov models, a penalized likelihood approach was also proposed. The possible consequences of ignoring the interval-censoring in a simple three-state ‘illness–death’ model have been discussed by Joly et al. (2002).

The rest of the article is organized as follows. In Section 2, we describe the methodology for non-Markov models that was applied for the analysis of the real heart transplant data. The obtained results can be found in Section 3. In Section 4, we present results from a small simulation study. Finally, in Section 5 we discuss the applied methodology and provide some final conclusions. General framework for multi-state models is provided in Supplementary material. All the codes and an example of the simulated data set are available at http://www.statmod.org/smij.

2 Non-Markov models

In multi-state models, one usually assumes that the process is Markov and homogenous in time. Then the transition probabilities can be derived from the transition intensities which can be modelled using baseline covariates by means of standard hazard-based models. These may include both multiplicative hazard models and additive hazard models. Provided that there are no loops, i.e., that there is no way to come back to a given state, the explicit expressions for transition intensities are available and yield plug-in methods for the corresponding probabilities, both for Markov and semi-Markov homogenous models. Even though the plug-in methods seem to be simple and straightforward, they often lead to complicated relations between the covariates and transition probabilities. Moreover, they can be used only when the transition intensities are constant or piecewise-constant in time. This implies that in regression models for transition intensities only time-fixed covariates are allowed. An attractive alternative involves the direct modelling on transition probabilities using non-parametric A-J estimator.

Relaxing the Markov assumption makes the standard parametric approaches inapplicable and restricts the use of non-parametric A-J estimator only for the occupation probabilities. In general, literature on non-Markov multi-state models is scarce and involves mainly non-parametric methods for a direct estimation of the transition probabilities. For a non-Markov illness–death process without recovery Meira-Machado et al. (2006) derived nonparametric estimators of the transition probabilities, which they compared with the A-J estimator under different scenarios in a simulation study. The results presented systematic bias of the A-J estimator when the Markov assumption is violated. Nevertheless, the state occupation probabilities can still be estimated using the A-J estimator since the product-integral estimator for those probabilities is consistent without requiring the Markov assumption. The variance of the occupation probability estimator in a non-Markov setting is not readily available and thus the Bootstrap method is required to estimate it (Glidden, 2002). The above apply for the case of non-Markov process with random censoring. The stage-dependent censoring can be handled via inverse probability of censoring weighting, as proposed by Datta and Satten (2002). In general, as noted by Commenges (2002), the likelihood or Bayesian approaches seem to be safer than marginal models because of the selection and censoring problems. However, in non-Markov situation, as for our data, the likelihood-based approaches do not apply and one must rely on marginal inference that usually ignores the interval-censoring problem. Some proposed approaches (Cook et al., 2004; Sutradhar and Cook, 2008) allow us to handle conditionally the Markov progressive model under interval-censoring using the random effects approach.

2.1 Pseudo-values approach for the non-parametric Aalen-Johansen estimator

In order to assess the effect of covariates on the occupation probabilities, we will apply the pseudo-values approach proposed by Andersen et al. (2003) to the A-J estimator of the occupation probabilities. This approach was originally proposed to model the cumulative incidence function in a competing risk model. To introduce the pseudo-value approach, let K = {1, 2, …, N} denote the finite state space for the considered multi-state process X with the time interval $[0,]$ , $<$ . First, we calculate the so-called pseudo-values for the function of interest F(X), which in our case is the A-J estimator of the occupation probability. More specifically, the pseudo-values are defined as:

{\overset{}{}}_{i} (t) M \overset{}{F} (X)(t) (M 1) \overset{}{F} (X)^{i} (t),

(2.1)

where $\overset{}{F} (X)(t)$ and $\overset{}{F} (X)^{i} (t)$ are the estimators of F(X(ρ)) at time ρ calculated for all individuals and for the subset without individual i, respectively, i = 1, …, M. The idea of pseudo-values is based on the leave-one-out estimator $\overset{}{F}_{i 1}^{M} {\overset{}{F}}^{i}$ (Shao and Tu, 1995). If $\overset{}{F}$ is unbiased, then $E (\overset{}{F}) F$ . In particular, it can be shown that ${\overset{}{F}}_{Jack} \overset{}{F} b_{Jack}$ , where $b_{Jack} (M 1)(\overset{}{F} \overset{}{F})$ is an unbiased estimate of F up to the second order. That holds, in particular, for the mean of the pseudo-values ${\overset{}{}}_{i}$ .

In the pseudo-value approach of Andersen et al. (2003), ${\overset{}{}}_{i} (t)$ are calculated for all individuals i at some arbitrary chosen time points ρ₁, …, t_k. Therefore, we obtain k pseudo-values for each subject. Note that due to the fact that ${\overset{}{}}_{i} (t_{1}),, {\overset{}{}}_{i} (t_{k})$ are evaluated for the same subject, they will be correlated. Hence, in the next step, and in order to measure the effects of covariates on F(X), the generalized estimating equations (GEE) approach is utilized with ${\overset{}{}}_{i}$ as a response. Let Z_i denote a vector of covariates of interest. We model the expected value of the pseudo-values as:

E (_{i}) g^{1} (^{T} Z_{i}),

(2.2)

with g(.) denoting a monotonic link function. Estimates of β are based on the unbiased estimating equations:

\underset{i}{} \{\frac{}{} g^{1} (^{T} Z_{i})}^{T} V^{1} \{{\overset{}{}}_{i} g (^{T} Z_{i})} \underset{i}{} U_{i} () U () 0,

(2.3)

where V^–1 is a working covariance matrix. The covariance matrix of $\overset{}{}$ is obtained using the sandwich estimator (Aalen and Johansen, 1978).

Following Graw et al. (2009) here we require that:

E (_{i} | Z_{i}) g^{1} (^{T} Z_{i}) O_{P} (1).

(2.4)

This condition of conditional unbiasedness of the pseudo-values, given the covariates formulated by Graw et al. (2009), relaxes unbiasedness of the pseudo-values, as formulated in Andersen et al. (2003). This relaxation is required in order to be able to use an appropriate theorem of the GEE approach (Liang and Zeger, 1986) to prove the large sample properties of the GEE solution for b. Condition (2.4) holds trivially when the observations are uncensored, and it was shown by Graw et al. (2009) to also hold for the pseudo-values in the right censored situation. The central argument was a second-order von Mises expansion of the A-J estimate which leads to an appropriate representation of the jackknife pseudo-values (Knott and Frangos, 1983). The interval-censoring was not considered.

We apply the pseudo-values approach to the A-J estimator of the occupation probabilities. We define first the A-J estimator ${\overset{}{F}}_{h} (t)$ for the probability p_h(ρ) of occupying state h at time ρ as:

{\overset{}{F}}_{h} (t)_{k 1}^{N} p_{k} (0) {\overset{}{p}}_{kh} (0, t), 2 mmh K; t,

(2.5)

where p_k(0) is the initial distribution of the process at time 0 and ${\overset{}{p}}_{kh} (0, t)$ are the estimates of the transition probabilities between states k and h from time 0 to ρ, that are obtained from the A-J product-integral estimator for the transition probability matrix P(s, ρ) = (P_kh(s, ρ)) for s = 0 as:

\overset{}{P} (0, t) \underset{(0, t]}{} (I d \overset{}{A} (u)).

(2.6)

In equation (2.6), $\overset{}{A}$ denotes the standard Nelson–Aalen estimator for the cumulative transition intensity matrix (Aalen and Johansen, 1978). Note, as mentioned in the Introduction, that in this non-parametric approach we ignore interval-censoring allowing only for potential right censoring when calculating the A-J estimator.

${\overset{}{F}}_{h} (t)$ is calculated for all individuals and we denote by ${\overset{}{F}}_{h}^{i} (t)$ the A-J estimate for the subset without individual i, i = 1, …, M. Then for subject i at time ρ a pseudo-value is calculated:

{\overset{}{}}_{i} (t) M {\overset{}{F}}_{h} (t) (M 1) {\overset{}{F}}_{h}^{i} (t).

(2.7)

Then, to assess the effects of covariates on the occupation probabilities, we fit a regression model on pseudo-values using the GEE approach according to (2.2) and (2.3), under the assumption (2.4).

Note that in the complete case, i.e., when there is no censoring, F_h(ρ) could be estimated unbiasedly as the proportion of subjects in state h at time ρ. Then i_i(ρ) is simply F_h(ρ). With censoring, we replace possibly unobserved F_h(ρ) by the pseudo-value ${\overset{}{}}_{i} (t)$ . In the regression model for i, in which we model E(i_i), the ${\overset{}{}}_{i}$ is used as a response. The validity of the assumption (2.4) was analyzed in Graw et al. (2009) only in the situation of competing risks. However, following their suggestion, this condition can also be proved in a more general model by constructing the pseudo-values using the smoothed mappings of the Nelson-Aalen and the Kaplan–Meier functionals.

We considered the univariate pseudo-value approach, in which we model each occupation probability separately as well as the multivariate pseudo-values approach where all occupation probabilities are modelled together, as suggested by Andersen and Klein (2006). Therefore, in the multivariate approach we model the vector ${\overset{}{}}_{i} (t) ({\overset{}{}}_{1 i} (t),, {\overset{}{}}_{N_{i}} (t))$ of the pseudo-values for the occupation probabilities for all N states.

Apart from the formal results of Graw et al. (2009), the pseudo-value approach was mainly evaluated using the simulation methods. The choice of numbers of time points and their location as well as the choice of a link function were shown to have a moderate impact on the results (Klein and Andersen, 2005).

2.2 Multinomial models approach

The regression models based on pseudo-values, presented in the previous section, allow us to study the effect of baseline covariates on the occupation probabilities in a non-Markov model, when the dependence on the history is a nuisance. However, it does not allow us to assess the impact of the history on those probabilities. To overcome this limitation we propose here an alternative approach, under which we can study how the observed history affects the occupation probabilities. This is based on viewing a homogenous Markov model as series of multinomial models for each observation conditionally on the previous observation. These models may be fitted using standard software for multinomial logistic regression. However, when the time is not discrete and the Markovian assumption is violated, we may still fit the multinomial model on the transition level, taking as a response the next state observed and adjusting for the time of transition as well as for baseline covariates.

In particular, for each individual we consider the triples (s_, s, t) such that the person visited the particular state p at time s before the current state c, i.e., X(s_) = p, X(s) = c, X(t) = h, h, c, p ∈ K, s_, s, t ∈ Γ; s_ < s < t. For these triples we considered the following model:

\Pr (X (t) h | X (s) c; X (s_{}) p) \frac{\exp (_{h})}{1 \underset{k r}{} \exp (_{k})}, h r;

_{h}_{h 0}_{h 1} t_{h 2} Z,

(2.8)

where X(t) is the next state, t-time of transition, Z-baseline covariates. For the reference category r we have:

\Pr (X (t) r | X (s) c; X (s_{}) p) \frac{1}{1 \underset{k r}{} \exp (_{k})} .

(2.9)

In the general setting, current and previous states could be included as covariates in a model (2.8–2.9). In our application, conditioning on c and p in (2.8) and (2.9) is realized by considering the corresponding subset of all transitions for the reasons explained in the next section. We work now on the transition level, instead of the individual level. Since we pool all individuals together, we need to account for the correlation introduced by using the same patient more than once. For that reason standard errors (SEs) for the multinomial model estimates were calculated by means of leave-one-out jackknife estimators (Shao and Tu, 1995). In particular, let r denote the estimated probability of interest from the multinomial model ((2.8–2.9) and $_{i}$ denote the same probability estimated excluding individual i. Then we calculate the jackknife SE as simple standard deviation of the ${\overset{}{}}_{1}$ estimates:

{SE}^{Jack} (\overset{}{}) \sqrt{\frac{N}{N 1}_{i 1}^{N} ({\overset{}{}}_{i} {\overset{}{\overset{}{}}}_{i})^{2}}, {\overset{}{\overset{}{}}}_{i} \frac{1}{N}_{i 1}^{N} {\overset{}{}}_{i}

(2.10)

3 Analysis of the heart transplant data

We apply the described methodology to the transplantation study introduced in Section 1. The data were taken from an international data base of the Eurotransplant heart recipient waiting list. We consider 2921 recipients who entered the list in the period from 1 January 2006 to 31 December 2008. Recipients’ observation was censored at 31 March 2010. At that time 528 patients died (18%), 1566 (54%) received a transplant, 239 (8%) were removed from the list because of other reasons and 588 (20%) were still on the waiting list. At the moment of entering the list some of the baseline information was also collected, namely: country of origin (seven countries), blood group and cardiovascular disease (categorized into dilated cardiomyopathy (DCM), coronary artery disease (CAD) and others).

Following a standard approach in the multi-state model framework, we started our analysis assuming a Markov process and modelled each intensity using the parametric approach. We considered only patients with more than one transition. In addition, due to numerical reasons and in order to be able to fit the model, we also imposed the constraints for the intensities with small number of transitions setting them to zero. To evaluate the fit of the Markov model, we graphically compared the observed and expected prevalence for each state for different models. We did not expect to observe substantial differences when assuming the exact transition times and without that assumption. Plausible reasons for the observed discrepancies include the dependence of transition rates on omitted covariates or on time (non-homogenous Markov process). The other possible reason for the discrepancies is the dependence of transition intensities on the time spent in the current state (semi-Markov) or the history of the process (non-Markov process). We estimated transition probabilities P_hr(0, t) using the A-J estimator, where h is not an absorbing state and h ! r. P_hr(s, t) were relatively close to P_hr(0, t – s) indicating that these probabilities seem to depend on the time interval t – s. This suggests that the process is homogenous. Nonetheless, as we observed, the time-homogenous Markov model had a poor fit. Therefore, the main issue for the data at hand seems to be that the Markov assumption is not satisfied.

Further examination of the history of individual profiles suggested common patterns such as (T, NT, T, NT, …) or (T, HU, T, HU, T, HU, …). Since the A-J estimator cannot be used to investigate the Markov assumption, we decided to apply the multinomial model approach presented in Section 2.2 adjusted for blood group, cardiovascular disease and informed consent (IC) law (binary). In our setting, the transition to the same state is not allowed and additionally for some transitions in the real data we observed very low frequency. Therefore, conditioning on c and p in (2.8) and (2.9) is realized by considering the corresponding subset of all transitions. Death was chosen as the baseline category for most of the models, unless there were none or very few transitions to death from a given current state c, for a given previous state p. Separate models were fitted for given states c and p and jackknife SE were calculated according to formula (2.10). In practice, in most of the cases due to the small subgroup sample sizes the adjustment for baseline covariates was not possible.

The results of this analysis are presented in Tables 1 and 2. None of the effect of baseline covariates, except the effect of blood group for a model with the previous state T and current HU in Table 1, was found to be significant. This is probably attributed to the lack of power due to the small size of the subset sample. Therefore, we present only the model adjusted for time and not for baseline characteristics (except the model in Table 1). Further results can be found in Supplementary material. Table 1 presents the results from the multinomial model for the previous state Transplantable and the current state HU. Because there were only three transitions to state R from HU, they were excluded from the analysis. Among the patients with the previous state Transplantable and the current state U there was only one case of transition to the state R and it was excluded from the analysis. As can be seen from Table 2, a multinomial model grabs the peak in the probability of being in state HU in the early times. This results in the highest estimate for log(OR) = 1.52, which is almost significant (p = 0.058). This is due to the fact that most of the transitions to state HU from other states took place in early times of being on the waiting list. For the current state NT there were only two transitions to the state HU and one to the state U and they were excluded from the analysis. There were no transitions to TT. As can be concluded, from NT patients were most likely to go back to T. The probability of death was greater than the probability of being removed. We analyzed in the same manner other histories. The results for those models are provided in Supplementary material.

Table 1
Estimates of log(OR) for the effect of current state and time on the probability for the next transition from the current state HU based on the multinomial model with the previous state Transplantable, adjusted for time and baseline covariates. Baseline category is probability of death. Jackknife SEs are calculated. The estimates for the time effect are not listed.

Previous = T, Current = HU

Intercept Blood B vs A Blood AB vs A

Pr(NT)/Pr(D) 0.42(1.13) –0.54(0.58) –0.90(1.30)

Pr(HU)/Pr(D) — — —

Pr(U)/Pr(D) –13.76(6.81) –1.15(0.62) –1.44(1.26)

Pr(T)/Pr(D) –0.23(1.01) –1.28(0.54) –1.13(1.01)

Pr(R)/Pr(D) — — —

Pr(TT)/Pr(D) 2.58(0.68) –0.5(0.48) 0.01(0.76)

	Previous = T, Current = HU
Pr(NT)/Pr(D)	0.42(1.13)	–0.54(0.58)	–0.90(1.30)
Pr(HU)/Pr(D)	—	—	—
Pr(U)/Pr(D)	–13.76(6.81)	–1.15(0.62)	–1.44(1.26)
Pr(T)/Pr(D)	–0.23(1.01)	–1.28(0.54)	–1.13(1.01)
Pr(R)/Pr(D)	—	—	—
Pr(TT)/Pr(D)	2.58(0.68)	–0.5(0.48)	0.01(0.76)

Source: Authors’ own.

Table 2

Estimates of log(OR) for the effect of the current state and time on the probability for the next transition from the current state U and NT based on the multinomial model with the previous state Transplantable, adjusted for time. Baseline category is probability of death. Jackknife SEs are calculated. The estimates for the time effect are not listed.

	Previous = T
	Current = U	Current = NT
	Intercept	Intercept
Pr(NT)/Pr(D)	–1.25(0.93)	—
Pr(HU)/Pr(D)	1.52(0.80)	—
Pr(U)/Pr(D)	—	—
Pr(T)/Pr(D)	0.52(0.81)	1.52(0.22)
Pr(R)/Pr(D)	—	–0.91(0.25)
Pr(TT)/Pr(D)	–0.23(0.82)	—

Source: Authors’ own.

As mentioned in the previous section, the A-J estimates for the transition probabilities are consistent only if the Markov assumption is satisfied. However, for our data, and as shown in the previous analysis using the multinomial approach, the Markov assumption seems to be violated. Therefore, in the final analysis we use the A-J estimates for the state occupation probability, which are consistent regardless of whether the Markov assumption is satisfied. To measure the effect of covariates in these occupation probabilities we employed the pseudo-value approach introduced in Section 2.1. To apply this approach we need to first choose the time points, at which the estimators are to be calculated. In the related framework of the competing risks models, it has been shown that the choice of the number of grid points can be proven crucial (Klein and Andersen, 2005). To be flexible, here we have chosen 10 time points based on the quantiles of the observed time distribution. To take into account the within-individual correlation (for each patient we calculate 10 pseudo-values), we use a GEE model with an unstructured correlation matrix and identity link. This model contained as covariates time, age, blood group, type of the disease (CAD, DCM and other (OD)) and IC indicator. The blood group A and CAD disease were set as the reference levels. The GEE model with time and baseline covariates was fitted for each occupation probability separately in univariate pseudo-values model (2.7) and the multivariate version for all occupation probabilities modelled together. In multivariate GEE model, the independence correlation matrix was used because of convergence problems for any other option. The results for the univariate model for state Transplanted are presented in Table 3 and illustrated in Figure 2. We have analyzed in the same univariate manner the occupation probabilities for the remaining states.

Table 3

Estimates of the effect of baseline covariates on occupation probability for state Transplanted. Results from the univariate regression model on pseudo-values.

	Estimate	SE	p-value
Intercept	–1.090	0.051	< 0.001
Time (Days)	–0.723	0.001	< 0.001
Age	–0.018	0.001	< 0.001
Blood B vs A	0.034	0.027	0.208
Blood AB vs A	0.154	0.036	< 0.001
Blood 0 vs A	–0.526	0.019	0.006
DCM vs CAD	–0.096	0.042	0.021
Other disease vs CAD	–0.046	0.044	0.293
IConsent (Y vs No)	–0.182	0.019	< 0.001

Source: Authors’ own.

Using the univariate pseudo-values approach, the probabilities of being in the given state can be compared between the different groups defined in terms of baseline characteristics such as blood group. However, we cannot compare the two complementary probabilities of being in state h and r. For example, we can compare if patients with AB blood group have higher probability of getting a transplant than patients with A blood group. Based on the univariate pseudo-values approach we can also compare the probability of death between the patients with A and AB blood group. However, we cannot assess what is more probable: death or getting a transplant, having the specified blood group. Therefore, we also fitted the multivariate pseudo-values model. The results regarding the impact of baseline covariates were similar as for the univariate approach and can be found in Supplementary material. However, based on the multivariate model, clinician could conclude for example that a patient with blood group AB is more likely to get a transplant than die as comparing to patient with blood group A. Similar questions could be also addressed using the multinomial approach.

Figure 2

Source: Authors’ own.

All computations have been performed in R. In particular, we used function msm() from msm package for fitting the Markov model. Piecewise-constant intensity models were also fitted using this function. For the A-J estimator calculation the etm() function from the etm package was used. Transitions to the same state are not allowed in the etm() function and therefore the times for those transitions were combined. For the multinomial models on the transition level, the multinom() function in the nnet package was used with the additional code written for the jackknife standard errors calculation. To apply pseudo-value approach for the occupation probabilities pseudo-values calculation needed to be implemented. Then the geese() function from the geepack package was used to fit the GEE models on the pseudo-values. All the codes used in our analyses are available at http://www.statmod.org/smij.

4 Simulation study

4.1 Design

To investigate the performance of the pseudo-values approach for non-Markov models we performed a limited simulation study of four different scenarios with 500 samples per scenario. In all scenarios, four transient and three final states corresponding to the analyzed real data set were considered. For each data set, and for N = 1500 subjects, we simulated the waiting times T_hr for a transition from state h to r from an exponential distribution with rate λ_hr. If the simulated waiting time T_hr was smaller than the simulated waiting time T_hs, then the individual moved from h to r. Each subject was initially in state T and could have at maximum 10 transitions. Individuals being in any of the non-final states after 10th transitions were considered as censored. In Scenario I, we simulated data from a homogenous Markov model with no group effect. Individual waiting times were drawn from an exponential distribution with equal rate parameters λ = exp(–1) for all possible transitions so that the average waiting time for any transition was around exp(1) = 2.7. In Scenario II, we simulated data from the Markov model with a group effect on the intensity of transition to state TT. To simulate the group effect in each data set the individual was randomly chosen to be in one of the two groups with probability 0.5. Then for all individuals in group 1 the transition times were drawn from an exponential distribution with equal rate parameters λ = exp(–1) for all possible transitions, whereas in group 2 the waiting times for the transition to TT were twice shorter. In Scenario III, the non-Markov model was simulated in which the group effect reflected the dependence on the history. The transition times were drawn from an exponential distribution with equal rate parameters λ = exp(–1) for all possible transitions except that in group 2 the waiting times for a transition to state r were twice shorter if the individual visited state r just prior the current state. Finally, in Scenario IV, we simulated both group effects together, the group effect on the transition intensity to state TT and the group effect related to the history. Therefore, for group 1, the waiting times to state TT were simulated to be twice shorter than in the other group and additionally for the subjects from the same, group 1, the waiting times to the previously visited state were twice shorter than to the other states.

The regression models on univariate and multivariate pseudo-values were fitted for each data set as well as multinomial models. To estimate a true group effect on occupation probabilities within the pseudo-values approach framework we fitted GEE models on the true occupation probabilities. Since we estimated twice longer waiting times as a group effect on the history or on the transition to state TT, we use log(2) as the true group effect in a multinomial model on a given intensity. For a non-Markov model from Scenario IV we compare the results with and without conditioning on the previous state for the estimated effect on the transition to state TT. As a summary for each scenario, bias and root mean squared error (RMSE) were reported. In the regression models on pseudo-values we used k = 10 time points in a GEE model placed in the quantiles of the transition times distribution and unstructured variance–covariance matrix. Additionally, for the Scenarios II and IV we have examined the impact of different number of time points in GEE model when increasing the sample size. In particular, we have fitted the Markov model from Scenario II and non-Markov model from Scenario IV with k = 5 and k = 10 time points for the simulated samples of size N = 500 and N = 1000.

4.2 Results

The pseudo-values approach did not reveal any difference between the groups when the effect was simulated independently on the response as in Scenario I or when the effect was only dependent on the history in Scenario III (results not shown). Figure 3 illustrates the median of the A-J estimates for the occupation probability for state Transplanted, Death and Removed for all simulated data sets for the Scenarios II and IV. As can be observed, the group effect of different transition intensities q_iTT on the probability of occupying state TT is reflected on occupation probabilities for all final states when the process is Markov (Scenario II) or not Markov (Scenario IV). For the state TT we could observe slightly stronger group effect for the Markov process with somewhat larger variability. This is related to lack of the history effect preventing from the transition to state TT in the Markov model. The pointwise confidence bounds are based on the simulations. In contrast, as noted by Glidden (2002), the pointwise confidence bounds using the recursive formula derived by Andersenb et al. (1993) for the occupation probabilities do not depend on the covariance structure $\sqrt{n} {(p_{h} (t) {\overset{}{F}}_{h} (t)}$ and depend only on $var {p_{h} (t) {\overset{}{F}}_{h} (t)}$ , which is the same for Markov and non-Markov processes.

We did not observe any substantial effect on the occupation probabilities of any non-absorbing state, neither for the Markov nor for the non-Markov scenario.

Table 4
Bias and RMSE of the group effect on occupation probability of state Transplanted, Removal and Death for a Markov Model (Scenario II) and non-Markov model (Scenario IV). Results from univariate and multivariate pseudo-values approaches for μ = 500 and μ = 1000, k = 5 and 10 time points.

Markov model
Non-Markov model

k = 100
k = 5
k = 10
k = 5

n = 500 Bias RMSE Bias RMSE Bias RMSE Bias RMSE

Transplanted

Univariate –0.012 0.037 –0.015 0.055 –0.015 0.021 –0.033 0.065

Multivariate –0.025 0.017 –0.019 0.023 0.013 0.008 –0.021 0.034

Death

Univariate –0.003 0.006 0.007 0.015 –0.004 0.005 0.014 0.066

Multivariate –0.001 0.003 0.008 0.005 –0.008 0.002 0.009 0.028

Removal

Univariate 0.003 0.018 0.006 0.047 –0.009 0.019 0.017 0.032

Multivariate
–0.016
0.011
0.045
0.014
–0.018
0.017
0.027
0.020

k = 10
k = 5
k = 10
k = 5

n = 500
Bias
RMSE
Bias
RMSE
Bias
RMSE
Bias
RMSE

Transplanted

Univariate −0.004 0.009 −0.011 0.035 0.014 0.005 −0.013 0.010

Multivariate −0.013 0.006 −0.009 0.020 0.022 0.004 −0.003 0.008

Death

Univariate −0.002 0.003 0.006 0.008 −0.012 0.006 0.009 0.015

Multivariate −0.001 0.002 0.007 0.010 −0.001 0.001 0.003 0.010

Removal

Univariate −0.014 0.012 0.004 0.024 −0.008 0.010 0.012 0.021

Multivariate −0.030 0.008 0.020 0.016 −0.016 0.007 0.019 0.013

	Markov model	Non-Markov model
Transplanted
Univariate	–0.012	0.037	–0.015	0.055	–0.015	0.021	–0.033	0.065
Multivariate	–0.025	0.017	–0.019	0.023	0.013	0.008	–0.021	0.034
Death
Univariate	–0.003	0.006	0.007	0.015	–0.004	0.005	0.014	0.066
Multivariate	–0.001	0.003	0.008	0.005	–0.008	0.002	0.009	0.028
Removal
Univariate	0.003	0.018	0.006	0.047	–0.009	0.019	0.017	0.032
Multivariate	–0.016	0.011	0.045	0.014	–0.018	0.017	0.027	0.020
	k = 10	k = 5	k = 10	k = 5
n = 500	Bias	RMSE	Bias	RMSE	Bias	RMSE	Bias	RMSE
Transplanted
Univariate	−0.004	0.009	−0.011	0.035	0.014	0.005	−0.013	0.010
Multivariate	−0.013	0.006	−0.009	0.020	0.022	0.004	−0.003	0.008
Death
Univariate	−0.002	0.003	0.006	0.008	−0.012	0.006	0.009	0.015
Multivariate	−0.001	0.002	0.007	0.010	−0.001	0.001	0.003	0.010
Removal
Univariate	−0.014	0.012	0.004	0.024	−0.008	0.010	0.012	0.021
Multivariate	−0.030	0.008	0.020	0.016	−0.016	0.007	0.019	0.013

Source: Authors’ own.

When comparing univariate and multivariate pseudo-values approaches, the latter seems to lead to smaller RMSE. When investigating the impact of the number of time points in the GEE model for the different sample sizes for Scenarios II and IV we observe that increasing the number of time points k and sample size μ leads to smaller RMSE (Table 4). We observed rather smaller RMSE in the non-Markov model, especially for the multivariate pseudo-value approach, and more time points.

Table 5 illustrates the group effect from the multinomial model conditional on the previous state Transplantable adjusted for time and current state. For the data simulated from two non-Markov models we observe smaller bias and RMSE for the effect on the history P(t)/P(D) for the Scenario III when there is no effect on the probability of the transition to TT, namely on P(TT)/P(D). For P(TT)/P(D) in the non-Markov model from Scenario IV obtained with no conditioning on the previous state we observed comparable bias and RMSE as in the Markov model. Conditioning on the previous state when estimating P(TT)/P(D) results in larger bias and RMSE which is something to be expected since we take only subset of all transitions. Even when conditioning on the previous state the RMSE and bias for P(TT)/P(D) are smaller than for the effect on the history P(T)/P(D).

Figure 3

Source: Authors’ own.

In an analogous manner, a group effect on the transition intensity to non-final state was simulated. For this scenario, the pseudo-value approach detected only the group effect on the occupational probability of absorbing states in the non-Markov model. The multinomial model was able to capture the simulated group effect on transition to T both in the Markov and non-Markov models (results not shown).

Table 5

Bias and RMSE of the group effect on the probability of transition to state Transplanted (TT) and the effect on history in multinomial models for the Markov model (Scenario II) and non-Markov model (Scenario IV), adjusted for time and current state. Baseline category is probability of death.

	n = 500		n = 1000
	Bias	RMSE	Bias	RMSE
Markov model
P(TT)/P(D)	−0.054	−0.209	0.042	0.163
Previous = T
P(TT)/P(D)	−0.068	−0.255	0.045	0.223
Non-Markov model
P(TT)/P(D)	−0.055	0.174	−0.018	0.157
Previous = T
P(TT)/P(D)	−0.072	0.249	−0.036	0.208
P(T)/P(D)	−0.292	0.366	−0.201	0.352
Non-Markov model
Previous v = T
P(T)/P(D)	−0.194	0.351	−0.172	0.277

Source: Authors’ own.

We have simulated the dependence on the history through the last state since that was suspected for the analyzed real data. However, the true dependence of the process could be more complicated. The simulation study demonstrated that the pseudo-value approach can capture the group effect of the transition intensity to one final state on the occupational probability of all final states of such a non-Markov model. Additionally, the multinomial model was able to capture also the effect related to the history.

5 Discussion

Motivated by study on patients on the waiting list for heart transplantation, in this article we presented two approaches for analyzing multi-state data for which the Markov assumption is not defensible. In particular, we use the pseudo-value approach to estimate the effects of baseline covariates on the occupation probabilities when the history dependence is a nuisance. When the interest is in a dependence history, the additional information can be derived from the multinomial model approach. Both approaches are simple in principle and can be easily implemented for any general non-Markov multi-state model. We observed less efficiency in the multinomial approach in estimating the history effects than for the effects on the transition to the final state. We should consider rather series of multinomial models with the fixed time point in each model and for each of these models we should consider the subset of individuals that have made a transition in that fixed time point. Since in our case the transition times are very dense and change rather continuously, it is very hard to find the reasonable large group of individuals with the same transition time. The pseudo-value approach requires independent censoring. One can however relax that assumption to conditional independence of censoring and event times given covariates and then apply an inverse probability of censoring weighting techniques (Scheike and Zhang, 2007). The method was recently extended for the case of clustered time-to-event data in the competing risk model by Logan et al. (2011). So far, the pseudo-values approach was mainly evaluated using the simulation results (Andersen et al., 2003; Klein and Andersen, 2005; Graw et al., 2009). Based on the simulation studies the choice of the number of time points, their location and also the choice of the correlation matrix in a GEE model was shown to have an impact on the results. The impact of the number of time points was also demonstrated in our simulation study. We have observed an increased efficiency for the multivariate pseudo-values approach comparing to the univariate model. The efficiency also increased for larger sample size and more time points in the GEE model. Adjusting for the optimal number of time points could be more problematic for the occupation probabilities of non-absorbing states since they usually exhibit a non-monotonic behaviour. These aspects need to be further studied along with the assessment of the goodness-of-fit and the choice of link function in GEE models.

Supplementary materials

Supplementary materials are available online at http://www.statmod.org/smij.

Footnotes

Acknowledgements

We thank Jacqueline M Smits from Eurotransplant International Foundation in Leiden for providing the data set analyzed in the article and for the medical consultation.

References

Aalen

Johansen

(1978) An empirical transition matrix for non-homogenous Markov chains based on censored observations. Scandinavian Journal of Statistics, 5, 141–50.

Andersen

Borgan

Gill

Keiding

(1993) Statistical models based on counting processes. New York: Springer-Verlag.

Andersen

Hansen

Klein

(2004) Regression analysis of restricted mean survival time based on pseudo-observations. Life-time Data Analysis, 10, 335–50.

Andersen

and Klein (2006) Regression analysis for multistate models based on pseudo-value approach with application to bone marrow tranplantation studies. Scandinavian Journal of Statistics, 34, 3–16.

Andersen

Klein

Rosthoj

(2003) Generalized linear models for correlated pseudo-observations with applications to multi-state models. Biometrika, 90, 15–27.

Commenges

(2002) Likelihood for interval-censored observations from multistate models. SORT, 27, 1–12.

Cook

Grace

Lee

K-A

Gladman

(2004) A conditional Markov model for clustered progressive multistate processes under incomplete observation. Biometrics, 60, 436–43.

Datta

Satten

(2001) Validity of the Aalen-Johansen estimators of stage occupation probabilities and Nelso-Aalen estimators of integrated transition hazards for non-Markov models. Statistics and Probability Letters, 55, 403–11.

Datta

Satten

(2002) Estimation of integrated transition hazards and stage occupation probabilities for non-Markov systems under dependent censoring. Biometrics, 58, 792–802.

10.

Datta

Sundaram

(2006) Nonparametric estimation of stage occupation probabilities in a multistage model with current status data. Biometrics, 62, 829–37.

11.

Glidden

(2002) Robust inference fpr event probabilities with non-Markov event data. Biometrics, 58, 361–68.

12.

Graw

Gerds

Schumacher

(2009) On pseudo-values for regression analysis in competing risks models. Lifetime Data Analysis, 15, 241–55.

13.

Grüger

Kay

Schumacher

(1991) The validity of inferences based on incomplete observations in disease state models. Biometrics, 47, 595–605.

14.

Joly

Commenges

Helmer

Letenneur

(2002) A penalized likelihood approach for an illness death model with interval-censored data: Application to age-specific incidence of dementia. Biostatistics, 3, 433–43.

15.

Klein

(2006) Modeling competing risks in cancer studies. Statistics in Medicine, 25, 1015–34.

16.

Klein

Andersen

(2005) Regression modeling of competing risks data based on pseudo-values of the cumulative incidence function. Biometrics, 61, 223–29.

17.

Knott

Frangos

(1983) Variance estimation for the jackknife using von Mises expansions. Biometrika, 70, 501–4.

18.

Liang

Zeger

(1986) Longitudinal data analysis using generalized linear models. Biometrika, 78, 13–22.

19.

Logan

Zhang

Klein

(2011) Marginal models for clustered time-to-event data with competing risks using pseudo-values. Biometrics, 67, 1–7.

20.

Meira-Machado

Una-Alvarez

Cadarsao-Suarez

(2006) Nonparametric estimation of transition probabilities in a non-Markov illness-death model. Lifetime Data Analysis, 12, 325–44.

21.

Scheike

Zhang

(2007) Direct modelling of regression effects for transition probabilities in multistate models. Scandinavian Journal of Statistics, 12, 325–44.

22.

Shao

(1995) The jackknife and bootstrap. New York: Springer-Verlag.

23.

Sutradhar

Cook

(2008) Analysis of interval-censored data from clustered multi-state processes: Application to joint damage in psoriatic arthritis. Journal of the Royal Statistical Society Series C, 57, 553–66.

	Markov model				Non-Markov model
	k = 100		k = 5		k = 10		k = 5
	n = 500	Bias	RMSE	Bias	RMSE	Bias	RMSE	Bias	RMSE
Transplanted
Univariate	–0.012	0.037	–0.015	0.055	–0.015	0.021	–0.033	0.065
Multivariate	–0.025	0.017	–0.019	0.023	0.013	0.008	–0.021	0.034
Death
Univariate	–0.003	0.006	0.007	0.015	–0.004	0.005	0.014	0.066
Multivariate	–0.001	0.003	0.008	0.005	–0.008	0.002	0.009	0.028
Removal
Univariate	0.003	0.018	0.006	0.047	–0.009	0.019	0.017	0.032
Multivariate	–0.016	0.011	0.045	0.014	–0.018	0.017	0.027	0.020
	k = 10		k = 5		k = 10		k = 5
	n = 500	Bias	RMSE	Bias	RMSE	Bias	RMSE	Bias	RMSE
Transplanted
Univariate	−0.004	0.009	−0.011	0.035	0.014	0.005	−0.013	0.010
Multivariate	−0.013	0.006	−0.009	0.020	0.022	0.004	−0.003	0.008
Death
Univariate	−0.002	0.003	0.006	0.008	−0.012	0.006	0.009	0.015
Multivariate	−0.001	0.002	0.007	0.010	−0.001	0.001	0.003	0.010
Removal
Univariate	−0.014	0.012	0.004	0.024	−0.008	0.010	0.012	0.021
Multivariate	−0.030	0.008	0.020	0.016	−0.016	0.007	0.019	0.013