Systematic comparison of approaches to analyze clustered competing risks data

Abstract

Background

In many clinical trials the study interest lies in the comparison of a treatment to a control group regarding a time to event endpoint like time to myocardial infarction, time to relapse, or time to a specific cause of death. Thereby, an event can occur before the primary event of interest that alters the risk for or prohibits observing the latter, i.e. a competing event. Furthermore, multi-center studies are often conducted. Hence, a cluster structure might be observed. However, commonly only the aspect of competing events or the aspect of the cluster structure is modelled within primary analysis, although both are given within the study design. Methods to adequately analyze data in such a design were recently described but were not systematically compared yet.

Methods

Within this work we provide a systematic comparison of four approaches for the analysis of competing events where a cluster structure is present based on a real life data set and a simulation study. The considered methods are the commonly applied cause-specific Cox proportional hazards model with a frailty, the Fine and Gray model for considering competing risks, and extensions of the latter model by Katsahian et al. and Zhou et al.

Results

Based on our simulation results, the model by Katsahian et al. showed the best performance in bias, square root of mean squared error, and power in nearly all scenarios. In contrast to the other three models this approach allows both unbiased effect estimation and prognosis.

Conclusion

The provided comparison and simulations help to guide applied researchers to choose an adequate method for the analysis of competing events where a cluster structure is present. Based on our simulation results the approach by Katsahian et al. can be recommended.

Background

The aim of many clinical trials is the comparison of two treatment groups regarding the occurrence of one specific event among many competing events. For example, a specific cause of death is of primary interest but another absorbing event (like another cause of death) occurs before the primary event of interest and hence the latter cannot be observed. Furthermore, clinical trials are often conducted at multiple treatment centers, which leads to a cluster structure. Commonly, the Cox proportional hazards model [1] is applied to investigate the time until the first event takes place, i.e. cause-specific analysis is applied. In this model, competing events are censored, while the treatment effect on the probability of occurrence of the event of interest is estimated [2]. Therefore, there is only one possible event during the estimation process. However, this can lead to an invalid analysis of the cumulative probabilities. Another frequently used approach is the Fine and Gray model [3] where individuals experiencing a competing event are not censored but remain in the risk set for the primary event of interest. The treatment effect is now estimated using a subdistribution hazard. Existing cluster structures are still mostly ignored in this evaluation of clinical trials with competing risks. This is considered problematic because it is assumed that existing cluster structures lead to a dependency of failure times within clusters [4, 5]. Ignoring these correlations might result in underestimation of variance of the group-specific regression parameters. Therefore, to address this problem, new approaches which analyze event times in competing risks settings that also take cluster structures into account, are needed. Usually competing risks are defined by the following properties:

• There are multiple events possible, but only one can occur at a time.

• In clinical trials, one event is always of primary interest.

• In general, researchers are interested in the first time the event occurs.

The aim of this work is therefore to compare the performance of different models used for analyzing survival data in settings with competing risks regarding their accuracy in the presence of cluster structures. In this paper, we compare two newer approaches by Katsahian et al. [4, 5] and Zhou et al. [6] which explicitly address this topic and contrast them to the commonly used Fine and Gray model [3] which addresses competing risks but not the cluster structure as well as the cause-specific Cox proportional hazards based model with a frailty term [7,8,9] which incorporates the cluster structure but ignores the presence of competing risks. However, in contrast to Katsahian et al. and Zhou et al. we are mainly interested in the effect estimation and not prognosis, i.e. we consider bias, square root of mean squared error, and empirical power for the estimated treatment effect and do not consider the performance of the models regarding prognosis. This is due to our focus on clinical trials where the main aim is to estimate a treatment effect. Hence, it is evaluated whether the models by Katsahian et al. and Zhou et al. can also be used for unbiased effect estimation and not only for prognosis. If both, effect estimation and prognosis, could be possible with one model this would be optimal.

The approach by Katsahian et al. uses an alternative estimation method, which has already been used in the context of standard survival analysis and has now been extended for competing risks settings. In addition, it uses a specific weighting technique by which individuals who have already experienced an event from a competing cause are weighted in the model. Individuals are considered to be at risk until they experience either the event of interest or censoring. The main interest of the approach by Zhou et al. is to assess marginal effects of covariates on the cumulative incidence function for the occurrence of an event of interest conditional on the covariates. The strength of this model is, that it takes into account both dependencies of failure times and censoring times.

The performance of the models is examined and compared using a Monte Carlo simulation study considering different relevant clinical settings. As mentioned above, as performance measures bias, square root of mean squared error, and empirical power for the effect estimation are considered.

Hence, we provide a neutral comparison as proposed by [10, 11]. Thereby, we highlight the properties of commonly used methods in the context of clustered competing risk data, as well as methods that have been proposed but are not used in clinical practice to analyze such data. The Cox proportional hazards model with frailty was chosen because it is the most frequently applied method in clinical trials where a time to event endpoint is of interest and a cluster structure is present due to different clinical centers involved. Although the competing risks structure is not considered in this method, it has been recommended for the primary analysis of randomized controlled trials [12]. The Fine and Gray model is of interest, since it is commonly used for the analysis of survival data where competing events are present. We wanted to investigate to which extend this method, although theoretically not fully appropriate for the setting at hand, differs from the ’correct’ methods or whether it might yield acceptable results in practice. The two ’newer’ methods were chosen because they were proposed for the analysis of clustered competing risk data and model a hazard ratio. Other methods that were described in the literature were additive models [13] or models with copulas [14]. The additive model does not result in a hazard ratio and is thus not comparable to the other methods [13]. Furthermore, additive models are not as common as multiplicative models in survival analyses. The copula approach produces an estimate which describes the degree of acceleration which was not our main interest [14].

This paper is structured as follows: We will start by the separate presentation of the statistical models. Afterwards we will present and compare the obtained simulation results of each model. In the end we discuss the methods and results and finish the article with a conclusion.

Methods

We consider a two-arm multi-center clinical study with an intervention I and a control C, where the primary endpoint (PE) is a time to event endpoint. Further, we assume that a competing event (CE) might be observed. The competing event might occur before the primary event of interest. Moreover, adminstrative censoring is assumed. The individuals $i=1,...,n$ are randomized in a 1 : 1 allocation to the two groups within each center ($k=1,...,K)$, where K is the total number of centers (i.e. clusters). We consider a two-sided test problem, i.e. test for difference with

Null-hypothesis:

$$\begin{aligned} H_0: \beta _{PE} = 0 \end{aligned}$$

(1)

and

Alternative hypothesis:

$$\begin{aligned} H_1: \beta _{PE} \ne 0. \end{aligned}$$

(2)

Thereby, $\beta _{PE}=log(HR_{PE})$ denotes the treatment effect for the primary endpoint of interest (HR = hazard ratio).

The following sections describe the methods evaluated in this work.

Cox proportional hazards model with frailty

In the Cox proportional hazards model [1] the hazard for an individual $i=1,...,n$ within a two-group comparison is modelled as follows:

$$\begin{aligned} \lambda _{i}(t)=\lambda _{0}(t)exp(\beta _{PE} X_i) \end{aligned}$$

(3)

where $\lambda _{0}(t)$ refers to a common cause-specific baseline hazard, i.e. the instantaneous baseline risk for an individual to experience an event of one specific type (primary event of interest) at time t given no prior event. $X_i$ is the treatment indicator, i.e. $X_i$ = 1 if individual i belongs to the intervention group (I), and $X_i$ = 0 for the control group (C). $\beta _{PE}$ is the corresponding coefficient.

This model can be extended to model heterogeneity in the data due to e.g. clustered data like within a multi-center trial. This extension is given as follows [8, 9]:

$$\begin{aligned} \lambda _{i}(t)=\lambda _{0}(t)exp(\beta _{PE} X_{i}+u Z_i). \end{aligned}$$

(4)

This model is also called frailty model [7,8,9], where u are independent and identically distributed from some positive scale family with expected mean value of 0 and variance $\theta$. Within this work we assume a gamma distribution with $u\sim log(\Gamma (\frac{1}{\theta },\frac{1}{\theta }))$. $Z_i$ is a vector of indicator variables, i.e. $Z_i=1$ if an individual i belongs to the cluster of interest and otherwise $Z_i=0$.

In a competing risks setting this approach is used as cause-specific analysis where the competing event time point is assumed as a censored time point. Thereby, a cause-specific treatment effect can be given, i.e. it allows one to estimate the effect of the treatment on the rate of occurrence of the primary outcome of interest in those subjects who are currently event free. Within a frailty model the treatment effect is conditional on the frailty.

For effect estimation the maximum likelihood approach can be used.

The partial log-likelihood with gamma frailty is given as [8, 9]:

$$\begin{aligned} l^{Cox}= & {} \sum \limits _{i=1}^n \delta _i\Bigg ((\beta _{PE} X_{i}+ u Z_i)\nonumber \\{} & {} - log\left. \left( \sum \limits _{j\in R_i} exp(\beta _{PE} X_{j}+u Z_j )\right) \right) \nonumber \\{} & {} +\frac{1}{\theta }\sum \limits _{k=1}^K(u_k-exp(u_k)). \end{aligned}$$

(5)

Thereby, the risk set $R_i$ is defined at the time of failure $t_i$ for the ith individual as

$$\begin{aligned} \lbrace j: (t_j\ge t_i)\rbrace . \end{aligned}$$

(6)

$\delta _i$ indicates the event time point ($\delta _i=1$) or censoring time point ($\delta _i=0$). $u_k$ is the independent and identically distributed frailty for center k.

Fine and Gray model

Fine and Gray [3] proposed a semi-parametric proportional hazards model for the subdistribution of a competing risk to assess the treatment effect on the marginal probability function. The hazard formulation for the primary event of interest is given as:

$$\begin{aligned} \lambda _{i,PE}(t)=\lambda _{0,PE}(t)exp(\beta _{PE} X_i). \end{aligned}$$

(7)

As before $X_i$ is the treatment indicator and $\beta _{PE}$ is the corresponding coefficient. The baseline subdistribution hazard for the primary event of interest is $\lambda _{0,PE}(t)$.

Within this model clustered data structure is not considered.

For effect estimation the partial log-likelihood can be used and is given as follows[3]:

$$\begin{aligned} l^{FG} =\sum \limits _{i=1}^n I(\epsilon _i=PE)\cdot \Bigg (\beta _{PE} X_i\\ \nonumber -log\left. \left( \sum \limits _{j\in R_i}exp(\beta _{PE} X_j)\right) \right) . \end{aligned}$$

(8)

Thereby, the risk set $R_i$ defined at the time of failure for the ith individual is given as

$$\begin{aligned} \lbrace j: (t_j\ge t_i)\cup (t_j\le t_i \cap \epsilon _j \ne PE)\rbrace \end{aligned}$$

(9)

with $\epsilon$ indicating the cause of failure and t the event times. Please note, that the risk sets, i.e. individuals at risk for an event, for this marginal approach and the cause-specific approach will usually not be the same.

Approach by Katsahian et al.

In 2006 Katsahian et al. [5] introduced the random effects model for the subdistribution hazard as an extension of the classical Fine and Gray model. It offers the opportunity to consider clustered event data appropriately in the model. In contrast to the standard model, frailties and random center effects are now integrated into the subdistribution hazard, resulting in the following modification of the subdistribution hazard of the individual i:

$$\begin{aligned} \lambda _{i,PE}(t)=\lambda _{0,PE}(t)exp(\beta _{PE} X_i+u Z_i). \end{aligned}$$

(10)

As before $X_i$ is the treatment indicator, $\beta _{PE}$ the corresponding coefficient, and the baseline subdistribution hazard for the primary event of interest is $\lambda _{0,PE}(t)$. u is the frailty vector, where frailties are considered independent and normally distributed with mean 0 and variance $\theta$. $Z_i$ is a vector of indicator variables, i.e. $Z_i=1$ if individual i belongs to the cluster of interest and otherwise $Z_i=0$.

For estimation the following log-likelihood function is used [4]:

$$\begin{aligned} l^{Kat}= & {} \sum \limits _{i=1}^nI(\epsilon _i=PE)\Bigg ((\beta _{PE} X_i+uZ_i)\nonumber \\{} & {} -log\left. \left( \sum \limits _{j=1}^n w_j(t_i)exp(\beta _{PE} X_j+u Z_j)\right) \right) \nonumber \\{} & {} -\frac{1}{2}\left( K\cdot log(2\pi \theta )+\frac{u^{'}u}{\theta }\right) . \end{aligned}$$

(11)

Thereby

$$\begin{aligned} \omega _i(t)=I(t_i\ge t \cup \epsilon _i \ne PE)\frac{\hat{G}(t)}{\hat{G}(t_i\wedge t)}. \end{aligned}$$

(12)

In this formula $t_i$ describes the minimum of the observed failure time and the censoring time of the individual i and $t_i\wedge t$=$min(t_i,t)$. $\omega _i(t)$ describes the weighting of the individual i at time t, according to its censoring status, where $\hat{G}(\cdot )$ describes the Kaplan-Meier estimate of the survival function of the censoring times. This formula shows that individuals who did not experience an event before time t $(t\le t_i)$ are fully considered in the model ($\omega _1(t)$=1). In case of $t>t_i$ two situations must be distinguished:

• If individual i experienced an event ($\epsilon _i=PE$) or was censored prior to time t, the weight is zero ($\omega _i=0$).

• If individual i experienced a competing event prior to time t, the individual weighting is done using $\omega _i(t)=\frac{\hat{G}(t)}{\hat{G}(t_i\wedge t)}$. In summary, individuals who have experienced a competing event ($\epsilon _i = CE$) remain at risk until some kind of censoring occurs $C_i>t$ ($C_i$ denotes the censoring time of individual i).

Approach by Zhou et al.

Zhou et al. [6] described a marginal proportional subdistribution hazards model which provides the ability to evaluate marginal effects of covariates on the cumulative incidence function. An existing correlation between individuals of the same cluster, due to unobserved factors, can be accounted for in settings of clustered competing risks. Within the marginal proportional subdistribution hazards model individuals within one cluster are considered as independent observations and the correlation between these individuals remains completely unspecified. Using the independence assumption, the Fine and Gray methodology can be used to estimate the cumulative incidence function and the effects of prognostic factors. However, the variance estimator allows not only the consideration of correlations between failure times but also the consideration of existing dependencies between the individual censoring times within a cluster. The variance estimation does not require any specification of the dependency between the individuals.

The main interest of the model is to assess the effect of the covariates on the marginal cumulative incidence function for the occurrence of event type PE conditional on the covariates. Thereby, the subdistribution hazard for PE is defined as:

$$\begin{aligned} \lambda _{i,PE}(t)=\lambda _{0,PE}(t)exp(\beta _{PE} X_{ik}). \end{aligned}$$

(13)

The baseline subdistribution hazard for the primary event of interest is $\lambda _{0,PE}(t)$. $X_{ik}$ is the indicator for an individual belonging to a specific treatment and cluster k with corresponding coefficient $\beta _{PE}$.

For effect estimation the log-likelihood is given as:

$$\begin{aligned} l^{Zhou}= & {} \sum \limits _{k=1}^K \sum \limits _{i=1}^{n_k} I(\epsilon _{ik} =PE)\Bigg (\beta _{PE} X_{ik}\nonumber \\{} & {} - log\left. \left( \sum \limits _{j\in R_{ik}} exp(\beta _{PE} X_{jk})\right) \right) . \end{aligned}$$

(14)

The number of individuals in a cluster is denoted by $n_k$. Risk set within a cluster k is defined via:

$$\begin{aligned} \lbrace j: (t_{jk}\ge t_{ik})\cup (t_{jk}\le t_{ik} \cap \epsilon _{jk} \ne PE)\rbrace . \end{aligned}$$

(15)

Simulation study

To provide a systematic comparison of the methods described in the previous sections, we conducted a simulation study with the statistic software R (Version 4.0.2) [15]. We used the packages survival [16], coxme [17], and crrSC [18] for the analysis.

R uses the Mersenne twister [19] for generating random numbers.

We considered a competing risks setting with one primary event and one competing event. Within each proportional cause-specific hazards are assumed for the two groups, i.e. a constant treatment effect over time.

To gain first insights into the performance of the four methods we considered scenarios with independent competing events.

Since we assumed a multi-center trial, different cluster counts are considered, as well as stratified randomization (i.e. stratum = center). Within each center, i.e. a cluster, the observations might be more correlated than between clusters. To model this we assumed a cluster-specific additional parameter. To be more precise a gamma-distributed frailty is added to the common baseline hazard (see formula 4).

The event times are generated as described by Bender et al. [20] for the two event processes, i.e. each individual gets an event time for the primary event of interest and the competing event. Of those the earlier time is assigned to the individual. If both times are larger than the administrative censoring time point, this censoring time point is used. Within the approach by Bender et al. the center-specific frailty can be incorporated, see e.g. [21, 22].

Table 1 Simulation Scenarios

In Table 1 the simulation scenarios are listed. In Column 2 the baseline hazard for the primary event of interest is given and in Column 3 the baseline hazard for the competing event. Columns 4 and 5 show the assumed hazard ratio for the primary event and the competing event, respectively. Columns 6 and 7 show the assumed logarithmic hazard ratio for the primary event and the competing event, respectively. In Column 8 the assumed variance for the gamma distributed frailty is given ($\theta$). In the last column the cluster count is given.

For Scenario 1 the baseline hazards for the two event types are the same. The hazard ratios differ but point into the same direction. Here, also no center-specific distribution is assumed, i.e. $\theta =0$. Nevertheless, three different cluster counts are considered for all scenarios which are depicted in scenarios a-c. Within Scenario 2 a moderate center-specific effect is assumed. In Scenario 3 the hazard ratios for the two event types change, i.e. an effect is only assumed for the primary event but not for the competing event.

Within Scenarios 4 and 5 the difference between centers increases since a higher frailty variance is assumed. For Scenarios 6-13 the baseline hazards change for either the primary event of interest or the competing event but else similar hazard ratios and frailty variances are assumed.

For all scenarios 250 individuals in total per data set were generated with about 125 in each treatment group, i.e. allocated using a binomial distribution. A follow-up of two years was assumed. For each scenario 2000 data sets were simulated and analyzed.

For performance comparison we considered the mean and standard deviation of estimated logarithmized hazard ratios for the primary event of interest as well as the corresponding absolute bias, square root of the mean squared error, and the empirical power.

Our primary interest is hence how the underlying treatment effect is estimated by the different models. Although it might also be of interest to evaluate the model performances regarding the cluster effect estimation, this cannot be done adequately. The proposed methods do not estimate the number of clusters but the cluster variance. The methods differ in their estimation approach for the cluster variance and are also different from our data simulation approach and thus it cannot be defined what a misspecification for the cluster structure would mean.

Application data set

To further illustrate the four approaches we used a data set (center) that is publicly available within the package named crrSC [23] of the statistic software R. The data set named center consists of multi-center bone marrow transplantation data and includes 400 patients from 153 transplant centers. The primary event of interest is “Acute or chronic Graft-versus-Host-Disease (GvHD)”. The competing event is “death” (without GvHD). We are interested in the influence of the source of stem cells, i.e. peripheral blood vs. bone marrow.

Results

Results of simulation study

In Tables 2 and 3 the results of the simulation study are displayed. Figure 1 further illustrates the results.

Table 2 Simulation Results

Table 3 Simulation Results

Note that the models by Fine and Gray, Katsahian et al., and Zhou et al. can only estimate the least false parameter due to our simulation set-up (which is basically based on the Cox model with frailty). This is due to the property that these models originally aim to get a prognosis taking into acount the competing event and do not focus on the estimation of the treatment effect. However, we want to evaluate whether one of the models can also produce an unbiased treatment effect estimate.

Fig. 1

Results of the simulation study: Boxplots for the estimated effects for the primary outcome. Dashed line = true underlying simulated effect; $log(HR_{PE})$ = logarithmic hazard ratio for primary endpoint

It can be seen that independent of the cluster structures, i.e. cluster count or correlation, the Cox model with frailty and the method by Katsahian et al. yield the smallest bias and square root of the mean squared error. Overall, these two models show very similar estimation results. The results of the Fine and Gray model and the model by Zhou et al. are also quite similar. Hence, bias and square root of the mean squared error coincide but are in some scenarios considerably higher compared to the results within the Cox model or the model by Katsahian et al. The bias and square root of the mean squared error within the Fine and Gray model and the model by Zhou et al. is higher in scenarios where the assumed hazard ratio for the competing event is not equal to 1.

The results seen for the bias or square root of the mean squared error are also reflected in the empirical power. The power for the Cox model and the model by Katsahian et al. coincide and exceed those of the Fine and Gray model and the model by Zhou et al. in nearly all scenarios (depending on the underlying simulated hazard ratio). The power of the model by Zhou et al. is higher than that observed for the model by Fine and Gray in all scenarios except the first where no cluster variance is assumed.

For the Cox model with frailty, 1 and 2 data sets produced an convergence error in Scenarios 10a and 11a, respectively and hence could not be included in the comparison.

Results of application

The results of the four different models for the application data set are given in Table 4. The source of stem cells was peripheral blood for 178 (44.5$\%$) patients and bone marrow for 222 (55.5$\%$) patients. For patients with peripheral blood as source of stem cells 79 (44.4$\%$) primary events of interest were observed and 36 (20.2$\%$) competing death events. For patients with bone marrow as source of stem cells 115 (51.8$\%$) primary events of interest were observed and 38 (17.1$\%$) competing death events. Median follow-up was estimated using the reverse Kaplan-Meier method and is 1639 days, the maximum observed event time for the primary event of interest is 2668 days, and maximum observed (censoring) time is 5138 days. The estimated hazard ratio for the primary event of interest is 0.85 for the Cox model with frailty but for the other approaches it is 0.82. For the approach by Fine and Gray and Zhou et al. this is in line with the results of our simulation study and hence the method by Zhou et al. seems to be an adequate extension for the Fine and Gray model when a cluster structure is present. Interestingly, the results gained by the approach by Katsahian et al. deviates more from the Cox-frailty results than seen in the simulation study. The main difference between our simulation study and the example is the number of clusters. Hence, this might be one reason but it could also be the cluster structure itself, i.e. cluster correlation and underlying distribution. To explore this, we altered the number of clusters in the example data by joining different clusters and present the results of these additional analyses in the appendix. However, the number of clusters seems not to be the reason for the difference between the Cox model results and the results by the approach by Katsahian et al. since the pattern of effect estimates remains the same. Another difference might be that the cause-specific hazards are not proportional but the subdistribution hazards.

Table 4 Application Results

Discussion

The analysis of a time to event endpoint where competing events and a cluster structure are present is a challenging task in cardiovascular or oncologic trials. Therefore, we compared four different methods that were proposed for those studies to give an overview of their properties in different clinical data situations. Here, the focus was on treatment effect estimation and not on prognosis which was already the focus of the main publications of the newly proposed methods [4, 6]. Hence, we extended the systematic comparison of the proposed methods to the effect evaluation to gain more insights in whether the proposed methods can also be used for unbiased treatment effect estimation and not only prognosis. The proposed methods differ in their performance and assumptions.

The cause-specific Cox-frailty model estimates the effect of the treatment on the rate of occurrence of the primary outcome of interest in those subjects who are currently free of any events (and conditional on the frailty). Since the simulation study does not include correlated event types, the true underlying effect for the primary endpoint can be estimated without bias, as seen in the results of the simulation study. Hence, the cause-specific Cox model is a good choice and can be recommended if the treatment effect on the rate of occurrence is of interest rather than prediction as it was already mentioned by Austin et al. and Allison [2, 24].

However, the performance of the model should be evaluated in scenarios where the endpoints are more correlated.

The Fine and Gray model allows one to estimate the effect of the treatment on the absolute risk of the primary outcome of interest over time, i.e. if one is more interested in prognosis. As described by Allison [24] the model may not be a good choice if one is interested in an unbiased treatment effect in randomized clinical trials for the primary event of interest. This is again supported by our simulation study. Moreover, it does not incorporate the cluster structure and hence we recommend to use the following model.

Katsahian et al. described a proportional subdistribution hazards model which is a frailty model. The model is also described as an extension of the Fine and Gray model to allow cluster structure, but produces unbiased treatment effects for the primary endpoint of interest for randomized clinical trials as also seen in the results of our simulation study. The model allows to estimate the cluster effect, as well as to incorporate this effect in prognostic analyses [4]. The latter is an advantage over the Cox-frailty model. Thus, this model can be recommended to estimate the treatment effect in clinical trials with competing events and a given cluster structure. However, it might be of interest to evaluate its performance in a setting where the event types are (more) correlated.

The model by Zhou et al. shows similar results in our simulation study as for the Fine and Gray model. Since this model was also described as extension of the latter, this is not surprising. We therefore conclude that if one is interested in the prognosis for the patients rather than an unbiased treatment effect on occurrence rates to use this model if a cluster structure is present in the data with competing risks.

Zhou et al. did conclude the same for their analysis. I.e. their simulation showed the potential bias and loss of power in hypothesis testing [6]. They hypothesized that this may arise from ignoring within-cluster correlations in variance estimation [6]. They also referred to the model by Katsahian et al. as an alternative. They further described that their method is particularly useful in applications with small groups of correlated observations where the correlation is mainly a confounding factor.

The application supported the results that the model by Zhou et al. is an appropriate extension of the Fine and Gray model in the case of a multi-center trial. However, it also supports that more complex simulation studies and evaluation of application might be necessary to shed light on the model by Katsahian et al.

Moreover, other approaches of interest in future work might be additive models or models using copulas [13, 14].

Conclusion

In conclusion, for clinical studies where two groups shall be compared regarding a time to event endpoint where competing events and a cluster structure are present the approach by Katsahian et al. can be recommended since it allows unbiased effect estimation and prognosis [4]. However, extended simulation studies are necessary to confirm its application in a broader range of data settings. For unbiased treatment effect estimation without the focus on prognosis, the Cox model with frailty can be used.