The competing risks multistate model
Figure 1 depicts the multistate model of a competing risks process (X
_{
t
} )_{
t≥0 }with initial state 0 and two competing event states 1 and 2. X
_{
t
} denotes the state that an individual is in at time t. The restriction to two competing events is for ease of presentation only. Initially, every individual is in state 0 at time origin, i.e., X
_{0} = 0. The individual stays in this state, i.e., X
_{
t
} = 0 until occurrence of any first event. Usually, there is one event of interest, modelled by transitions into state 1, and all other first event types may be subsumed into the competing event state 2.
The competing risks process moves out of the initial state 0 at the event time T. At time T , the process enters one of the competing event states 1 or 2. Hence, the state X
_{
T
} that an individual is in at time T denotes the type of the first event, often called cause of failure. X
_{
T
} equals either 1 or 2.
Key quantities in competing risks are the causespecific hazards (CSHs) α
_{0j
}(t), j = 1, 2,
The dt notation in (1) is a short, but more intuitive form of the more common α
_{0j
}(t) = lim_{Δt↘0 }P(T ∈ [t, t + Δt), X
_{
T
} = j  T ≥ t)/Δt. The interpretation of (1) is that the CSH at time t times the length dt of a very (infinitesimally) small time interval equals the conditional probability of making an 0 → j transition within that very small time interval [t, t + dt). The CSHs can be envisaged as forces of transition, which work along the arrows in Figure 1. Analogous to cumulative distribution functions, basic statistical inference addresses cumulative quantities. We therefore also write A
_{0j
}(t) for the cumulative CSHs, A
_{0j
}(t) = .
The CSHs are key, because, as seen below, the pertinent probabilities are deterministic functions of the CSHs, and because both nonparametric estimation and Coxtype regression modelling build on them. The fundamental nonparametric estimator is the NelsonAalen estimator of the cumulative CSHs
j = 1, 2, where the summation is over all observed event times s ≤ t. Variance estimation and, properly standardised, asymptotic normality of is discussed in detail in Chapter IV of [7]. The most important regression approach are Cox models for the CSHs. If, as is usual, we assume that covariates have different, i.e., causespecific effects on the CSHs, the partial likelihood factorises such that a Cox model for α
_{01}(t) may be fitted by additionally censoring type 2 events, and vice versa for a Cox model for α
_{02}(t), see, e.g., Chapter VII of [7]. The R [15] packages mvna[16] and survival[17] provide for convenient computation of the NelsonAalen estimator and the Cox model, respectively.
The α
_{0j
}'s sum up to the allcause hazard α
_{0·}(t)dt = P (T ∈ dt  T ≥ t) with cumulative allcause hazard A
_{0·}(t). The survival function of T is therefore a function of both α
_{0j
}'s, . As usual, P (T > t) is estimated using the KaplanMeier estimator, which is a deterministic function of the causespecific NelsonAalen estimators.
The socalled cumulative incidence functions (CIFs) are the expected proportion of individuals experiencing a certain competing event over the course of time,
j = 1, 2. The interpretation of the right hand side of (2) is that the CIF at time t equals the integral over all 'probabilities' to make an 0 → j transition precisely at time u. For this, an individual has to stay in state 0 until just before time u, which is reflected by the term P (T > u). Conditional on being in state 0 just before time u, an 0 → j transition at time u is reflected by the term α
_{0j
}(u) du.
One consequence is that the CIFs are involved functions of all CSHs via P (T > u). The CIFs P (T ≤ t, X
_{
T
} = 1) and P (T ≤ t, X
_{
T
} = 2) add up to the allcause distribution function 1  P (T > t). The AalenJohansen estimators of the CIFs can be obtained from (2) by substituting P (T > u) with the KaplanMeier estimator and α
_{0j
}(u) du with the increment of the causespecific NelsonAalen estimator. A detailed discussion of the AalenJohansen estimator is in Chapter IV of [7]; variance estimation is assessed in [18]. In R, survival is typically used for estimating P (T > t). The AalenJohansen estimator may be computed using etm[19], and prediction of the CIFs based on Cox models for the CSHs is implemented in mstate[20].
It is crucial to any competing risks analysis that both causespecific hazards completely determine the stochastic behaviour of the competing risks process [7, Chapter II.6]. This is mirrored above by the CIFs and the allcause distribution function being deterministic functions of all CSHs. However, these deterministic relationships are involved. Therefore, we will pursue an algorithmic interpretation of the CSHs below.
Algorithmic interpretation of the causespecific hazards
Thinking of the CSHs as momentary forces of transition which move along the arrows in Figure 1 suggests that competing risks data are generated over the course of time as follows:

1.
The event time T is generated with distribution function 1  P (T > t), i.e., with hazard α
_{01}(t) + α
_{02}(t) = α
_{0·}(t).

2.
At time T , event type j occurs with probability α
_{0j
}(T )/α
_{0·}(T ), j = 1, 2.
Using this algorithm for simulation studies in competing risks has been discussed in [13]. It is important to note, however, that the algorithm goes beyond the computational question of how to implement simulations. Rather, the algorithm reflects the probabilistic question of how to build a probability measure based on the CSHs. This aspect is discussed in detail by [12] in the more general context of multistate models which are realized as a nested series of competing risks experiments.
The algorithmic perspective of this paper then implies that the task of statistical inference is to detect the ingredients of the above algorithm. We illustrate this approach in the data example below. To this end, we note that the analysis of a combined endpoint is restricted to step 1 of the above algorithm. Here, the effect of a treatment, say, on both CSHs determines whether the occurrence of an event (of any type) is delayed or accelerated. In step 2, the type of an event again depends on the treatment effect on both CSHs. We illustrate below that interpretation is straightforward if the treatment effects on the CSHs work in opposite directions or if one CSHs remains unaffected. However, interpretation will become more challenging if there are unidirectional effects on both CSHs. We will find that, in general, it is also mandatory to consider causespecific baseline hazards in the interpretation.
The 4D study
The background of the 4D study was that statins are known to be protective with respect to cardiovascular events for persons with type 2 diabetes mellitus without kidney disease, but that a potential benefit of statins in patients receiving hemodialysis had until then not been assessed. Patients undergoing hemodialysis are at high risk for cardiovascular events. The 4D study was a prospective randomised controlled trial evaluating the effect of lipid lowering with atorvastatin in 1255 diabetic patients receiving hemodialysis. Patients with type 2 diabetes mellitus, age 1880 years, and on hemodialysis for less than 2 years were enrolled between March 1998 and October 2002. Patients were randomly assigned to doubleblinded treatment with either atorvastatin (619 patients) or placebo (636 patients) and were followed until death, loss to followup, or end of the study in March 2004.
The 4D study was planned [21] and analysed [14] for an event of interest in the presence of competing risks. The event of interest was defined as a composite of death from cardiac causes, stroke and nonfatal myocardial infarction, whichever occurred first. The other competing event was death from other causes. Wanner et al. reported a CSH ratio of 0.92 (95%confidence interval [0.77, 1.10]) for the event of interest. There was essentially no effect on the competing CSH.
The simulation study below will use the data of the placebo group. In this group, 243 (38.2% of 636) events of interest and 129 (20.3%) competing events were observed. There were 264 (41.5%) censored patients. The NelsonAalen estimators for the event of interest and for the competing event are displayed in Figure 2.
'Empirical simulation'
Simulations were based on the empirical probabilities defined by . To avoid the issue of model misspecification, which is outside the scope of the present investigations, 'ideal' Cox models were used for generating data in the treatment group as explained below. Event times and event types were generated following the twostep algorithm described earlier.
To be specific, event times in the placebo group were drawn from the distribution defined by the KaplanMeier estimator derived from . Here, a practical complication arose in that the KaplanMeier estimator only spent about 76% of the probability mass. I.e., the KaplanMeier curve did not drop down to 0 because of rightcensoring, which is a common phenomenon in clinical trials. This was handled by putting point mass beyond the largest observed time; corresponding realisations were always censored. In other words, on average about 100%  76% = 24% of the simulated event times equalled some time that was larger than the largest observed time in the original data set. The corresponding observations were always censored with censoring times generated as explained below. Event types were generated by the binomial experiment of the twostep algorithm, substituting the CSHs by the increments of the cumulative CSHs. The algorithm was analogously applied, when transformed cumulative CSHs were used for the placebo group.
Data in the treatment group were generated based on 'ideal' Cox models. That is, CSH ratios exp(β
_{1}) for the event of interest and exp(β
_{2}) for the competing event were specified. If the original empirical baseline hazards were used, data were drawn based on the probabilities defined by . If transformations of were used for generating placebo data, the CSH ratios acted on the transformed causespecific baseline hazards. This reflects the situation that the true underlying baseline CSHs are α
_{01}(t) and α
_{02}(t) (or transformations thereof) in the control group, while the CSHs of the treatment group are exp(β
_{1})α
_{01}(t) and exp(β
_{2})α
_{02}(t) for the case of untransformed baseline hazards. Random censoring times were generated for all individuals based on the KaplanMeier estimator of the censoring survival function in the placebo group. Note that this estimator did spend 100% of the probability mass.
Overview of simulation scenarios
The scenarios investigated in the next Section differed both with respect to the choice of β
_{1} and β
_{2} and in terms of the cumulative baseline hazards, i.e., the cumulative CSHs in the placebo group. The choice of the β's broadly falls into three categories.
One category is characterised by β
_{2} = 0, i.e., there is no effect on the competing CSH. A prime example are implantable cardioverter defibrillators [11], which display a beneficial effect on the CSH of sudden cardiac death but no effect on the CSH for death from other causes. The assumption of β
_{2} = 0 straightforwardly implies the direction of the treatment effect on the CIF: If β
_{1} < 0, the CIF of interest in the treatment group is always less than the one in placebo group. The relationship is reversed, if β
_{1} > 0. This is intuitively understood thinking of the CSHs as momentary forces of transition, and it is reflected in step 2 of the simulation algorithm. E.g., if β
_{1} < 0 and β
_{2} = 0, the binomial event type 1 probabilities are reduced for the treatment group.
A second category is characterised by opposite treatment effects on the CSHs. This category straightforwardly implies the direction of the treatment effect on the CIF, too: The constellation β
_{1} < 0 and β
_{2} > 0 implies a smaller CIF of interest in the treatment group but also a larger competing CIF. These relations are reversed for β
_{1} > 0 and β
_{2} < 0. This is again intuitively implied by thinking of the CSHs as momentary forces of transition, and it is also reflected in step 2 of the simulation algorithm. E.g., if β
_{1} < 0 and β
_{2} > 0, the binomial event type 1 probabilities are reduced for the treatment group.
Finally, unidirectional treatment effects on the CSHs constitute the third category. Interestingly, the interpretation of unidirectional effects is straightforward when both competing events are fatal in the sense that a treatment with β
_{1} < 0 and β
_{2} < 0, say, is beneficial. But unidirectional effects also present the most challenging scenario in terms of understanding the resulting course of the CIFs. The interpretation for step 1 of the simulation algorithm is straightforward. If, e.g., β
_{1} < 0 and β
_{2} < 0, events of any type will happen later. The interpretational challenge becomes apparent in the second step of the algorithm: If, e.g., β
_{2} < β
_{1} < 0, the relative magnitude of the CSHs changes such that the binomial event of interest probabilities are increased. The interpretational difficulty is the increase of this probability, although β
_{1} is negative. This constellation may result in an (eventually) increased CIF of interest.
All three categories are encountered in practice. Examples from clinical trials in hospital epidemiology are given in [22].
Difference to the latent failure time model
Some readers may be more familiar with competing risks as arising from riskspecific latent times, say, T
^{(1)} and T
^{(2)}. The connection to our multistate framework is T = min(T
^{(1)}, T
^{(2)}) with event type X
_{
T
} = 1, if T
^{(1)} < T
^{(2)}, and X
_{
T
} = 2 otherwise. The latent failure time model imposes an additional structure, which has been heavily criticised mainly for three reasons. The dependence structure of T
^{(1)} and T
^{(2)} is, in general, not identifiable [23]. Since the latent times are unobservable, there is something hypothetical about them, which questions their plausibility [24]. Perhaps most importantly, it has been disputed whether the latent failure time point of view constitutes a fruitful approach to answer questions of the original subject matter [25].
Despite of this critique, latent times are the predominant approach for simulating competing risks data [13]. Assuming, for tractability, T
^{(1)} and T
^{(2)} to be independent with hazards equal to the causespecific hazards α
_{01}(t) and α
_{02}(t), respectively, is computationally correct in that simulation based on this model yields the right data structure.
However, nothing is gained from assuming the additional latent structure, either. As a consequence, we will emphasise simulation and interpretation along the lines of the Algorithmic interpretation of the causespecific hazards outlined earlier. This avoids the concerns on identifiability, plausibility and usefulness.
E.g., in the 4D study, the typical interpretation of the latent times would be that T
^{(1)} is the time until death from cardiac causes, stroke or nonfatal myocardial infarction, while T
^{(2)} is the time until death from other causes. Such an interpretation has given rise to debating whether, say, a patient may conceptually still die from other causes after having died because of a cardiac event. In contrast to this, our approach only assumes that a patient is conceptually at risk of experiencing any of these events, provided that none of these events has happened so far.