Joint AFT random-effect modeling approach for clustered competing-risks data

Competing risks data arise when occurrence of an event hinders observation of other types of events, and they are encountered in various research areas including biomedical research. These data have been usually analyzed using the hazard-based models, not survival times themselves. In this paper, we propose a joint accelerated failure time (AFT) modeling approach to model clustered competing risks data. Times to competing events are assumed to be log-linear with normal errors and correlated through a scaled random effect that follows a zero-mean normal distribution. Inference on the model parameters is based on the h-likelihood. Performance of the proposed method is evaluated through extensive simulation studies. The simulation results show that the estimated regression parameters are robust against the violation of the assumed parametric distributions. The proposed method is illustrated with three real competing risks data sets.


Introduction
Competing risks data arise when occurrence of an event hinders observation of other types of events, and they are encountered in various research areas such as medicine, biology, econometrics, and engineering [1].Usually the event of interest is the primary event (e.g.time to recurrence), and other types of events (e.g.time to death) are competing events.Two broad classes of models for analyzing the competing risks data have been used based on the Cox's proportional hazards model: one is to model the cause-specific hazard of each event type separately [2] and the other is to model the subhazard (i.e. the hazard function of a subdistribution) for a particular event of interest [3].
The two classes of models have been extended to clustered competing risks data, which are correlated within a cluster, via frailty.Huang and Wolfe [4] proposed a cause-specific frailty model that is an extension of the ordinary frailty model using expectation-maximization (EM) algorithm.Katsahian et al. [5] extended the Fine-Gray model [3] to a subdistribution hazard frailty model with a random center effect.Gorfine and Hsu [6] studied a flexible cause-specific hazard frailty model, which is an extension of the bivariate frailty model [7].Christian et al. [8] and Ha et al. [9] have proposed the hierarchical likelihood (h-likelihood) methods for the inferences of the cause-specific hazard and subdistribution hazard frailty models, respectively.
For the competing risks analyses, however, most attentions have been drawn to the hazard-based models rather than to the AFT model in survival analysis literature.It is well known that the estimated regression parameters in the AFT model are robust against violation of the model assumption, and it directly describes the relationship between the logarithm of failure time T and covariates which is easy to interpret [10,11].Thus, in this paper we propose a new AFT joint modeling approach based on the random effects for analyzing jointly clustered competing risks data.The model inference is based on the h-likelihood method [1,10,12].Unlike the classical likelihood for fixed parameters only, the h-likelihood is constructed for both fixed parameters and unobserved random effects at the same time.The h-likelihood avoids the integration with respect to the random effects themselves, whereas the marginal likelihood approach often involves intractable integration to eliminate the random effects [1,12].
The rest of the paper is organized as follows.In Section 2, we propose a joint AFT competing-risks model with random effects and outline the construction of the hlikelihood.Section 3 presents the estimation procedure using the h-likelihood method.Simulation studies for evaluating the performance of the proposed method are presented in Section 4. In Section 5, we illustrate the proposed method using two clustered competingrisk data sets.In Section 6, we present an extension of the proposed method, with a motivated example data set, and the likelihood ratio test for association between the event times.The discussion is presented in Section 7. Finally, the technical details are provided in the Supplementary Materials.

Joint AFT competing-risks models
When a cluster has two outcomes such as main and competing event times, there may be correlations between the two outcomes because of an unobserved shared cluster effect [1,5,6].As ignoring such associations, biases may occur in marginal models.However, the joint model describes the association between the two outcomes via an unobserved random effect that makes the joint model more informative [13,14].
In clinical studies, we consider a clustered data set with a total sample size n = q i=1 n i , where q is the number of clusters (centers), and n i is the number of observations (subjects) in the ith (i = 1, . . ., q) cluster.For the jth (j = 1, . . ., n i ) subject in the ith cluster, let T ijk be time to event from the kth (k = 1, . . ., K) cause.Let T ij be time to the first event which is denoted by T ij = min(T ij1 , T ij2 , . . ., T ijK ), and let ξ ij ∈ {1, . . ., K} be the corresponding cause of event (type of event).Let C ij denote the independent censoring time.Denote by V i a shared unobserved random effect of the ith cluster.We assume that given V i = v i , C ij is conditionally independent and non-informative of (T ij , ξ ij ) for j = 1, . . ., n i [1,8].
In this paper, for simplicity, we first propose a cause-specific joint AFT model with a shared random effect based on two types of events (K = 2).
The proposed model (1) gives an explicit expression for correlation between log T ij1 and log T ij2 : which motivates the interpretation of the proposed model (1).Here, the association parameter γ allows the magnitude of the association to be different between two event times (i.e.log T ij1 and log T ij2 ).For example, if γ > 0 (γ < 0), there is a positive (negative) correlation between two event times, and there is no correlation between two event times if γ = 0.The proof of (2) is given in the Supplementary Material A. Besides, we plot the relationship between the association parameter γ and the correlation coefficient ρ between log T ij1 and log T ij2 for the proposed model which is shown in Figure 1.As shown in Figure 1, the correlation (ρ) has the same sign as the association parameter γ , and the magnitude of the correlation increases as |γ | increases.Several authors have studied the joint competing risks frailty model with a shared random effect [1,4,[15][16][17] which corresponds to the proposed AFT model (1).In particular, [1,15] have mentioned that when γ = 0, the hazard rate in Type 2 event does not depend on random effects v i (i = 1, . . ., q) and is noninformative for the hazard rate in Type 1 event, so that the two hazard rates are not associated for the cause-specific frailty model with shared frailties.

H-likelihood construction
The observed event times and event indicator are, respectively, defined by where ) is time to the first event, C ij is the censoring time, δ ijk is the censoring indicator of the kth (k = 1, 2) type of event, and I(•) is the indicator function.Based on the two assumptions mentioned at the beginning of Section 2, following Ha et al. [1] and Christian et al. [8], the h-likelihood for the proposed model ( 1) is defined by: where where 1ijk (k = 1, 2) are the conditional log-likelihoods for (Y ij , δ ijk ) (k = 1, 2) given V i = v i , respectively, 2i is the log-likelihood for V i , and (•) is the cumulative distribution function of N(0, 1).Note here that where Then, the h-likelihood (4) can be expressed as

Estimation procedure
Below we present the h-likelihood estimation procedure [1,8] for the proposed model ( 1). (

1) Estimation of the fixed and random effects
Given dispersion parameters θ = (φ 1 , φ 2 , α, γ ) T , the score equations of fixed and random effects τ = (β T 1 , β T 2 , v T ) T are given by where s = 0, 1, . . ., p − 1, and i = 1, . . ., q.Here, is the hazard function of N(0, 1), where ϕ(•) and (•) are the density and cumulative distribution functions of N(0, 1), respectively.Because of censoring, where and μ ij1 and μ ij2 are given in (6).Therefore, according to Ha et al. [10] and Buckley and James [18], we consider pseudo-response variables Y * ijk (k = 1, 2) for the joint AFT model with censored data.The definitions are as follows: Then, under the two assumptions (i.e.conditional independence and non-informativeness) in Section 2, we have that [10].Let y * ijk be the observed value of Y * ijk .Furthermore, it can be easily shown that when log Substituting y * ijk into the score equations ( 8) reduces them, respectively, to Let , where z ij = (z ij1 , . . ., z ijq ) T is the q × 1 group indicator vector whose rth (r = 1, . . ., q) element is ∂μ ij1 /∂v r , and v = (v 1 , . . ., v q ) T is the q × 1 vector of random effects v i 's.Then, the corresponding second derivatives are, respectively, given by where s, l = 0, 1, . . ., p − 1, and r, t = 1, . . ., q, and Thus, the (2p + q) × (2p + q) negative Hessian matrix H = −∂ 2 h/∂τ ∂τ T can be expressed as where ) are diagonal weight matrices [10] with the elements w ijk /φ k , and Q = I q /α with the q × q identity matrix I q .Here, X and Z are the model matrices of fixed effects β k and random effects v, respectively.
Let β = ( β1 T , β2 T ) T .Following Ha et al. [1] and Christian et al. [8], we can easily show that the maximum h-likelihood (MHL) estimates can be obtained by solving the following AFT iterative weighted least squares (IWLS) equations: where Here, μ k and y * k are n × 1 vectors of μ ijk and y * ijk (k = 1, 2), respectively.Note that there is an inverse of W k (i.e.W −1 k , k = 1, 2) in w * k which can lead to convergence problems due to the computation of complex weight elements w ijk with ζ(m ijk ).Thus, for further convergence, in this paper we use an alternative procedure without calculating the inverse of W k (k = 1, 2) [1].That is, we use the following w * * , instead of computing directly Ww * in (15): (
Secondly, for the estimation of the association parameter γ , we use a grid search method based on p τ (h).Given ( φ1 , φ2 , α), let p * τ (h) be a profile likelihood of γ , defined by Then, we find an optimal value γ that maximizes p * τ (h), which is given by The proposed fitting algorithm is summarized as follows: Step 0: Obtain the initial estimates of (β 1 , β 2 ) under the model (1) without random effects v.The initial estimates of association parameter γ and v are all used as zero values, and the initial estimates of variance components (φ 1 , φ 2 , α) are used as (0.5, 0.5, 0.1).
Step 3: Repeat Step 1 and Step 2 until the maximum absolute difference of the previous and current estimates of (τ , θ) is less than 10 −8 .
After convergence has been achieved, we calculate the standard error of τ using the inverse of the negative Hessian matrix (14) [1,8,9].Note that the initial estimates in Seep 0 are based on models without both random effects and association (i.e.v = 0 and γ = 0) [1,19].We have found from simulation studies that there is no impact on the estimated results even if they are chosen differently.

Simulation study
As the proposed cause-specific joint AFT model ( 1) is parametric, we conduct two kinds of numerical studies to evaluate the performance of proposed estimation procedures for the proposed AFT model (1) based on 500 replications of simulated data.
The first one is when the proposed model is correctly specified.That is, the simulated data are generated from normal distributions for both random effect and error terms.The second one is when the proposed model is incorrectly specified, which is useful for investigating the robustness of the proposed method.That is, the simulated data are generated from non-normal distributions for the random effects, whereas the random error terms follow normal or non-normal distributions.The proposed model ( 1) is then fitted to the two simulation cases using the estimation procedure given in Section 3.
From 500 replications of simulated data, for fixed effects β k (k = 1, 2) of the both types of events we calculate the mean, the standard deviation (SD), the mean of the estimated standard errors (SE) and the mean squared error (MSE) of βk (k = 1, 2).In addition, we calculate the empirical coverage probability (CP) for a nominal 95% confidence interval for β k (k = 1, 2).For the dispersion parameters θ, the mean, SD, and MSE for θ are also obtained.
Below are the details on simulation designs and results for those two kinds of studies, respectively.

Simulation design (1) Simulation design I
For the cause-specific joint model that includes two types of events, we use a technique that is similar to [8,20] to generate survival time T. The conditional cause-specific survival functions for Type 1 and Type 2 events given v i are, respectively, given by: where The corresponding hazard functions are, respectively, given by Similar to [20], the conditional survival function of all-cause survival time T is as follows: where The procedure of simulation design is as follows: (i) Sample sizes are n = q i=1 n i = 200, 400, 800, where (q, n i ) = (20, 10), (20,20), (80, 10), respectively.(ii) Consider a single binary covariate x which is generated from a Bernoulli distribution with a success probability of 0.5.(iii) The true values of parameters are based on the resulting estimates from the multicenter bladder cancer data in Table 6.Here, the true fixed effects (regression parameters) are similar to the estimates of the main covariate (CHEMO).(iv) Set the true fixed effects of two types of events as β 1 = (β 10 , β 11 ) T = (6.1,0.9) T and β 2 = (β 20 , β 21 ) T = (8.3,−0.4) T , respectively.(v) Set the true dispersion parameters as (φ 1 , φ 2 ) = (3.1,1.5), α = 0.1, 0.5, and γ = 0.2, 1, −1, respectively.(vi) Generate v from the normal distribution N(0, α) with α = 0.1 and 0.5.(vii) For the generation of survival time T, we use the inverse transformation method [21], i.e. let S T (t ij |v i ) = U, where U is a random variable that follows the standard uniform distribution.For the calculation of survival time T, we use a numerical rootfinding algorithm [22], i.e. find the roots of the equation S T (t ij |v i ) − U = 0 to obtain survival time T. To obtain the proportion of Type 1 and Type 2 events, a binomial experiment is conducted with probability λ 1 (T)/(λ 1 (T) + λ 2 (T)) on cause 1. (vii) Generate censoring time C from a uniform distribution U(0, c), where c is used to adjust the right censoring rate to 20% (lower) and 40% (higher). (

2) Simulation design II
To investigate the robustness of the proposed model, we set a sample size n = q i=1 n i = 400 with (q, n i ) = (20,20).Consider one covariate which is generated from a Bernoulli distribution with probability 0.5.The true values of fixed effects and dispersion parameters are the same as before.
Here, we consider two scenarios for evaluating the robustness of the proposed model.The first one specifies the distribution of random effects v i (i = 1, . . ., q) as the misspecified distribution (i.e.non-normal distributions).For this purpose, we use the following three distributions for v i , i.e. two points (2P), extreme value (EV), or mixture normal (MN) distributions: (i) For the 2P distribution, generate v i with equal probability 0.5 at the points -1 and 1. (ii) For the EV distribution, generate v i through a transformation such that the mean and variance are 0 and 1, respectively.(iii) For the MN distribution of v i , the unimodal (MN-1) and bimodal (MN-2) distributions are, respectively, considered.The unimodal distribution is 0.5N(− √ 0.5, 0.5) + 0.5N( √ 0.5, 0.5), and the bimodal distribution is 0.
The true association parameter is set to γ = 0.5.The generations of survival time T and censoring time C are similar to those for the design I.The right censoring rates are the same as before as 20% and 40%.
The second scenario specifies the distributions of the error terms ij1 and ij2 as the EV, given that random effects v i are normal or non-normal distributions.The details are as follows: (i) Set the true association parameter γ = 0.5.(ii) Generate v i from N(0, 1) or non-normal (the same as the first scenario) distributions.(iii) Generate ij1 and ij2 from the EV distributions EV(0, 1/1.5) and EV(0, 1), respectively, such that given v i , Type 1 or Type 2 event times follow Weibull distributions.(iv) Survival time T is generated similarly as in design I.
Specifically, when ijk ∼ EV(0, σ k ) with location parameter 0 and scale parameters Here, λ ijk = exp(−μ ijk ψ k ) is scale parameter and ψ k = 1/σ k is shape parameter, (k = 1, 2).Furthermore, the conditional survival function of survival time T can be expressed as: (v) Generation of censoring time C and the censoring rates are the same as before.

(1) Correctly specified case
We first present the simulation results when the proposed model ( 1) is correctly specified under the simulation design I. Tables 1 and 2 show the simulation results of fixed effects (β 10 , β 11 , β 20 , β 21 ) for the both types of events when the variances of random effects are α = 0.1 and α = 0.5, respectively.The simulation results for the corresponding dispersion parameters (φ 1 , φ 2 , α, γ ) are summarized in Table 3 .When the censoring rate is 20%, the proportions of Type 1 and Type 2 events are about, 64% and 16%, respectively, and when the censoring rate is 40%, they are about 50% and 10%.
First, from Table 1 under α = 0.1 we find that the estimated fixed effects of Type 1 and Type 2 events overall look reasonable as sample size (q and/or n i ) increases.In particular, the estimated fixed effects are also close to the true values for the three different association cases (γ = 0.2, 1, −1).Note that the SEs in Tables 1 and 2 are the average of 500 estimated standard errors of βkj (k = 1, 2; j = 0, 1), and the SDs are the empirical estimates of the true {var( βkj )} 1/2 ((k = 1, 2; j = 0, 1).We also see that there is a good agreement between SEs and SDs, except for the case of the estimated intercept ( β20 ) in the Type 2 event.A possible reason is that the observed proportions of Type 2 events are relatively small (i.e.16% for 20% censoring and 10% for 40% censoring) as compared to those of Type 1 events.Accordingly, the CPs of the 95% CI are also in good agreement with the nominal value of 0.95, except for lower CPs of the intercept β 20 due to the underestimated SEs.As expected, we see that the biases and variations (SEs, SDs and MSEs) tend to decrease as the sample size increases, but that they tend to increase as the censoring rate increases.The results in Table 2 when α = 0.5 are overall similar to those in Table 1, except for slightly larger variations.
Table 3 presents the simulation results of dispersion parameters (φ 1 , φ 2 , α, γ ) T when the variances of random effects are α = 0.1 and α = 0.5, respectively.Here, φ1 and φ2 are the estimates of the variances of random error terms for two types of events, and γ is the estimate of the association parameter.We observe that the estimates ( φ1 , φ2 , α, γ ) of all dispersion parameters are close to their true values.As expected, the SDs and MSEs of the estimated dispersion parameters decrease as the sample size increases, and they increase when the censoring rate becomes larger.In addition, we calculate the estimate ρ of correlation ρ in (2).We find from Table 3 that the estimates of ρ are overall improved with the sample size.Note that α is the variance of the random effect v which is nonnegative and Table 1.Simulation results on the proposed estimation method over 500 replications for the regression parameters in the joint AFT competing risks model under the variance of random effect α = 0.1 and the variances of random errors (φ 1 , φ 2 ) = (3.1,1.5); true regression parameters are β 1 = (β 10 , β 11 ) = (6.1,0.9),    bounded on zero, leading that the estimates of γ seem to be unstable when α decreases.Moreover, Table 3 shows that the estimates of γ is less biased as α increases.This is a possible reason why the variations (i.e.SDs and MSEs) of γ decrease as α increases.In summary, the results from Tables 1-3 indicate that the proposed h-likelihood method overall works well when the assumed models are correctly specified.
(2) Incorrectly specified cases Now, we investigate the performance of the proposed method under simulation design II when the distribution of the random effect v is misspecified, and the random error terms k (k = 1, 2) are both misspecified or not.Table 4 shows the simulation results under the model mis-specification.We observe that when the distributions of the random errors are correctly specified as N(0, φ k ) (k = 1, 2), but that of the random effect is mis-specified (i.e.non-normal), the estimates of all fixed effects tend to be robust.The CPs of β are also in good agreement with 95% nominal level, except for the intercept β 20 as shown in Table 1.As expected, the biases and variations (SEs, SDs and MSEs) increase as the censoring rate increases.When the random errors follow the extreme value distribution EV(0, σ k ) (k = 1, 2) (i.e.non-normal), even if the random-effect distribution is normal, the estimates of the slopes (β 11 , β 21 ) for both types are overall robust in terms of biases and CPs, except for the estimates of the intercepts (β 10 , β 20 ).
We find that the estimates of the dispersion parameters (α, γ ) related to random effect are somewhat robust when censoring is 20%, even if the assumed distributions of random errors and random effects, or only random effects are misspecified.However, they are slightly sensitive when the censoring is 40%.We observe that the estimates of the dispersion parameters (φ 1 , φ 2 ) related to random errors seem to be very sensitive to the true values of two scale parameters (σ 1 , σ 2 ) in the extreme value distribution of the error terms.
Accordingly, the results presented in Table 4 imply that the proposed h-likelihood method overall performs well when the assumed model is incorrectly specified, especially under the low censoring rate.

(3) Cases of ignoring association
In addition, we investigate the behaviors of the proposed method when the association between the two event (Types 1 and 2) times in the proposed model ( 1) is ignored (i.e.γ = 0) under the simulation design I.The simulation schemes are based on the design I above under n = 400 with (q, n i ) = (20, 20) and 20% censoring rate.Here, we consider the four combination according to sizes of (α, γ ), with (α, γ ) = (0.1, 0.2), (1, 0.2), (0.1, 1), (1, 1), which correspond to Cases 1-4 in Table 5.We also compute the relative biases (Rbias) of all available estimated parameters, defined by Rbias = {(O ς − ς)/ς} × 100, where ς = (β, θ) with θ = (φ 1 , φ 2 , α, γ ) T .The simulation results are presented in Table 5.Following [4], we compare the full model (proposed model) with the reduced model (i.e. proposed model with γ = 0).In terms of the Rbias, regardless of the size of α, the full model performs well overall, but the reduced model ignoring the association induces bias, especially when γ and/or α are large (i.e.Cases 2-4) under Type 2. These results confirm the simulation results of the joint frailty model by [15].Moreover, the SDs and MSEs of the estimated parameters in the full model are overall slightly smaller than those in the reduced model.Notes: The simulation studies are conducted with 500 replications under the sample size n = 400 with (q, n i ) = (20,20), and the association parameter γ = 0.5 and variance of random effects α = 1; The estimates of (φ 1 , φ 2 ) under the EV distribution of the random error indicate those of scale parameters (σ 1 , σ 2 ) of the EV distribution.The remaining true parameters are the same as Tables 1-3; : distribution of random errors, V: distribution of random effects, CR: censoring rate, N: normal distribution, 2P: two Points distribution, EV: extreme value distribution, MN-1: unimodal mixed normal distribution, MN-2: bimodal mixed normal distribution.
In particular, the reduced model gives relatively lower CPs for β 20 and β 21 when both γ and α are large (i.e.Case 4) under Type 2 event.
In conclusion, the simulation results from Tables 1-5 indicate that the proposed method overall performs well when the assumed model is correctly or incorrectly specified.

Illustration
For illustration of the proposed method, we consider two practical examples.

Multi-center bladder cancer data
The multi-center bladder cancer data set is from a multi-center clinical trial conducted by the European Organization for Research and Treatment of Cancer (EORTC) [23].The data set consists of 392 bladder cancer patients from 21 centers, and the number of patients per center varied from 3 to 78, with the mean of 18.9 and the median of 14.The event time of interest (Type 1 event) is time to the first bladder recurrence, and the competing event (Type 2 event) is a death prior to recurrence.Of the 392 patients, 200 (51.02%)patients had recurrence of bladder cancer (Type 1 event), 81 (20.66%) patients died prior to the recurrence (Type 2 event).In addition, 111 (28.32%) patients who were still alive and without recurrence were censored at the date of the last available follow-up.Notes: The simulation studies are conducted with 500 replications under sample size n = 400 with (q, n i ) = (20, 20) and 20% censoring rate.Rbias= 100( ς − ς)/ς is the relative bias, where ς is the true value of some parameter.
As in Section 4.2, we fit the full model (1) (M1; proposed model) and the reduced model (γ = 0 in (1), M2) with the multi-center bladder cancer data using the estimation procedure presented in Section 3. The estimates and the corresponding standard errors (SEs) of the fixed effects from both types of events in M1 and M2 are summarized in Table 6.We can see that the effect of the main covariate (CHEMO) is positively significant at a 5% significant level in both models in terms of t-test (t − value = Estimate/SE).That is, in M1 the estimate of the effect of CHEMO is 1.011 (SE = 0.250) and in M2 it is 1.008 (SE = 0.250), which means that the CHEMO significantly prolongs time to the first bladder recurrence.In terms of the estimates and SEs of the covariates, we can see that the results are overall very similar in both M1 and M2 models, especially for the Type 1 event, so we cannot determine which model is the best choice.It may be due to the fact that the estimates of the association parameter γ and the variance α of the random effect are somewhat small under full model (i.e. ( γ , α) = (0.170, 0.112) in M1), as shown in the first case of simulation results of Table 5.
For model selection, we considered the conditional Akaike information criterion (cAIC [1,[24][25][26]) which is defined as: Here, h 1 = ijk 1ijk is the conditional log-likelihood given in (4), and df c = trace(H −1 H * ) is an effective degree of freedom with H = −∂ 2 h/∂τ 2 and H * = −∂ 2 h 1 /∂τ 2 .Note that a smaller cAIC indicates a better model.The cAIC values of M1 and M2 are 1412.154and 1412.326,respectively.As mentioned above, the difference is very small, so that we may select M2 in terms of the parsimony principle of the model.
In addition, we have found out that the results of the proposed AFT model (M1) in Table 6 are similar to those of the cause-specific frailty model with shared log-normal frailties [17], especially for the significance of the estimated regression parameters: see Table S1 in Supplementary Material C.

Multi-center BMT data
Bone marrow transplantation (BMT) is one of the most commonly used treatments for leukaemia patients.A multi-center study of BMT for leukaemia patients was conducted at four hospitals in Australia and the United States, from March 1, 1984 to June 30, 1989 [27].As the graft-versus-host disease (GVHD) may develop during this period, the recovery after transplantation is quite a complex process [28].
One hundred and thirty seven (137) patients from four hospitals with either acute myelocytic leukaemia (AML) or acute lymphoblastic leukaemia (ALL) were enrolled for the study.The number of patients per hospital varied from 17 to 76, with the mean of 34.25 and the median of 22.Total 82 patients died or relapsed.Here, the event of interest (main event or Type 1 event) is death or relapse, and the competing event (Type 2 event) is GVHD.The details of this study can be found in Copelan et al. [29].The numbers of Type 1 and Type 2 events and no event (censoring) are 46 (33.6%), 71 (51.8%), and 20 (14.6%), respectively.Here, 46 was calculated as 82-36 = 46, where 36 cases had GVHD before death or relapse.The survival data within the same hospital may be correlated due to unobservable shared commonality among patients, thus the cluster effect in the competing risk analysis cannot be ignored [27].
The potential risk factors (covariates) we consider are as follows [27]: -AML risk group (AML.Low= 1, AML.High = 2, otherwise = 0) -Donor age in years (D.age) -Donor sex (D.sex; Male = 1, Female = 0) -Donor cytomegalovirus immune status (D.CMV; CMV Positive = 1, CMV Negative = 0) -French-American-British classification based on standard morphological criteria for AML patients (FAB; FAB grade 4 or 5 and AML = 1, otherwise = 0) -The methotrexate used as a Graft-Versus-Host-Prophylactic (MTX; Yes = 1, No = 0) Similar to the multi-center bladder cancer data of Section 5.1, for the multi-center BMT data we also fitted the full model M1 (proposed model) and the reduced model M2 (γ = 0).The estimation results are given in Table 7.We can see that the effect of AML.Low is significant in both types of events (Types 1 and 2) under the both models.Furthermore, D.age is significant in the competing event (Type 2) for both models.The estimate of the association parameter γ is 0.384 in our proposed model, and the estimates of the random effects' variance α are 1.063 and 1.039 in M1 and M2, respectively.
For model selection, the cAIC values for M1 and M2 are 645.289and 646.880, respectively, so that M1 is chosen.This indicates that there may be correlations among the outcomes within hospitals, leading that our proposed model is more reasonable.
Furthermore, the results of the proposed AFT model (M1) in Table 7 are also similar to those of the cause-specific frailty model with shared log-normal frailties [17], particularly for the significance of the estimated regression parameters: see Table S2 in Supplementary Material C.

Extension to three types of events
In this session we present an extension of our procedure with two types of events to three types of events.We consider a multi-center breast cancer data set (B-14) from a randomized double-blinded multi-center clinical trial, conducted by the National Surgical Adjuvant Breast and Bowel Project (NSABP) [30].The aim of this study was to investigate the effect of tamoxifen against placebo for patients who had negative axillary lymph nodes and estrogen receptor positive breast cancer after surgery.As in Christian et al. [8], we consider a high risk subset of patients from the B-14 study, with tumor size greater than 2.5 centimeters.Seven hundred and thirty seven (731) women (placebo: 371, tamoxifen: 360) with followup were eligible for the study.The median age for women in either placebo or treatment group was 55 years.A series of multiple types of treatment failure were recorded; local, regional, or distant recurrence of the original cancer as well as a new second primary cancer or death.
We consider the three types of events as follows: -Type 1: A local or regional recurrence.
-Type 2: A new second primary cancer in the contralateral breast.
-Type 3: A distance recurrence, other new second primary cancer or death.The numbers of Types 1-3 events and no events (censoring) until the last follow-up were 113 (13.45%), 64 (7.62%), 388 (46.19%) and 275 (32.74%), respectively.The covariates considered are age and treatment (1: tamoxifen, 0: placebo); and the treatment is a main covariate in this trial.For modeling of the three types of events, the joint AFT model ( 1) is extended to a joint model with a shared random effect v i (i = 1, . . ., q), composed of three interlinked AFT submodels, where i = 1, . . ., q, j = 1, . . ., n i .Here, T ij1 is the event time of interest, and T ij2 and T ij3 are different competing event times.The definitions of parameters are similar to those of Section 2. Note here that γ 1 [γ 2 ] are association parameters that are used to measure the correlation between Type 1 and Type 2 [Type 1 and Type 3], respectively.Let T ij = min(T ij1 , T ij2 , T ij3 ) be time to the first event.Then, the corresponding observed random variables, Y ij and δ ij , are the same as (3), with k = 1, 2, 3.Under the two assumptions mentioned in Section 2, the h-likelihood for model (27) with three types of events is defined by where (29) Similar to (15), the corresponding IWLS equations for τ where are diagonal matrices, with the elements are w ijk /φ k , Here Q = I q /α with the q × q identity matrix I q .w * where H is the (Kp + q) × (Kp + q) negative Hessian matrix with K = 3.We use the multi-center breast cancer data set to fit the extended model (27).The estimates and SEs of fixed effects and estimates of dispersion parameters are summarized in Table 8.We find that in full model (M1; extended model) the estimate (SE) of the tamoxifen effect in Type 1 event are 0.700 (0.181), implying that the tamoxifen effect for Type 1 event is significant in terms of t-test (t − value = Estimate/SE) at a 5% significant level.In other words, tamoxifen significantly prolongs time to local or regional recurrence.Similarly, the age effect in Type 1 is also significant.The estimates of the association parameters γ 1 [γ 2 ] are 1.084 [1.590], respectively, which implies a somewhat strong positive association between Type 1 and Type 2 [Type 2 and Type 3]; for example, we can see that patients who experienced a local or regional recurrence (Type 1) will be at a greater risk for developing a new second primary cancer (Type 2) in the contralateral breast.The estimate of the random effect's variance α is 0.662, indicating a heterogeneity among the patients.Furthermore, in terms of the cAIC, our extended (full) model (cAIC = 2889.613)performs better than the reduced model with γ 1 = γ 2 = 0 (i.e.cAIC = 3367.299).
In addition, similarly to two real data sets in Section 5, the results of the proposed AFT model (M1) in Table 8 are again similar to those of the cause-specific frailty model with shared log-normal frailties [17], except for the significance of age effect in Type 3: see Table S3 in Supplementary Material C.

The likelihood ratio test (LRT) for association
We are interested in the inference of correlation ρ in Equation ( 2), but ρ is not a model parameter in the joint model (1).As shown in Figure 1, the behavior of ρ in Equation ( 2) is similar to that of association parameter γ [4].Thus, following [1,31], we present a profile LRT for γ based on Equation (19), i.e. p * τ (h).That is, let L(γ ) = p * τ (h).Then, under the null hypothesis H 0 : γ = 0, the LRT statistic (LRT s ) defined by follows asymptotically chi-square with degree of freedom 1.The LRT statistic's critical value is 3.84 at level 5%.Note here that L( γ ) and L(γ = 0) are calculated under full and reduced models, respectively.Below we illustrate the LRT method (32) with a small numerical study and the three real data used in previous sections.Note: CR: censoring rate, q: number of clusters (centers), n i : cluster size. (
From the observed power on Table 9, the LRT statistic under γ = 0 is somewhat conservative for a small sample size n = 400, but the observed power at γ = 0 (i.e. the observable rate of rejecting H 0 at the 5% level under γ = 0) becomes overall closer to the nominal level 5% with a larger sample size n = 1, 600.As expected, the power of all tests under γ > 0 increases as γ , α or n increases, but, as expected, it decreases when the censoring rate increases from 20% to 40%.

(2) Real data sets
We also calculate the value of LRT statistic based on the three real data sets, i.e. bladder cancer, BMT, and breast cancer data sets in Tables 6-8.
For two real data sets in Tables 6 and 7, L( γ ) and L(γ = 0) are calculated under M1 (Full) and M2 (Reduced), respectively.In Table 6, γ = 0.170 and the corresponding correlation ρ = 0.009 which is very small; the value of LRT s is, as expected, also small with 0.114.In Table 7, γ = 0.384 and the corresponding ρ = 0.099 which is slightly larger; the value of LRT s is, as expected, slightly large with 1.082.Thus, H 0 : γ = 0 is not rejected at level 5% for the two data sets (bladder cancer and BMT data).
For the breast cancer data set in Table 8, let η = (γ 1 , γ 2 ).Then, under null hypothesis H 0 : η = 0 (i.e.γ 1 = γ 2 = 0), the LRT statistic (LRT s ) defined by follows asymptotically chi-square with degree of freedom 2. The LRT statistic's critical value is 5.99 at level 5%.In Table 8, ( γ1 , γ2 ) = (1.084,1.590) and the corresponding ( ρ1 , ρ2 ) = (0.229, 0.427) which are relatively larger.As expected, the value of LRT s in (33) is very large with 45.194, which gives the rejection of H 0 : η = 0 at level 5%.In summary, the significance results by the LRT above are identical to the model selection results by cAIC, except for the BMT data in Table 7.These LRT results may indicate that the estimates of ρ for the three data sets support the model selection performed and quantify very intuitively the degree of the (linear) dependence between the (log) competing times.However, the relationship between both results would require a further theoretical research.

Discussion
We have successfully applied the h-likelihood method to the joint AFT model for clustered competing risks data.Simulation results demonstrate that our proposed estimation procedure performs well regardless of model misspecification.In addition, we compare the proposed AFT model with the cause-specific frailty model with log-normal frailties.Both models provide similar results in terms of the significance of covariates' effects based on three real data sets.However, the AFT model directly describes the relationship between survival time and covariates which is easy to interpret.
The main advantage of the proposed modeling approach based on the h-likelihood method is that it avoids intractable or high dimensional integration over random effects in computing the marginal likelihood [1,8,19].Therefore, our method provides an efficient procedure for fitting the cause-specific joint AFT model for clustered competing risks data.
It would be a further work to find directly the SE for the dispersion parameters θ = (φ 1 , φ 2 , α, γ ) T in the current procedure because it requires extensive computations for the second derivatives of p τ (h) with respect to (φ 1 , φ 2 , α), and p * τ (h) with respect to γ , respectively.Extending the cause-specific joint AFT model with shared random effects to that with correlated random effects would be also an interesting further work.

Disclosure statement
No potential conflict of interest was reported by the author(s).

Figure 1 .
Figure 1.The relationship of γ and the correlation coefficient between log T ij1 and log T ij2 for the proposed model.

1 T , w * 2 T , w * 3 T
2,3) are n × 1 vectors, and w * = (w * ) T .For the estimation of dispersion parameters θ = (φ 1 , φ 2 , φ 3 , α, γ 1 , γ 2 ) T , we use the procedures similar to Section 3. The corresponding adjusted profile h-likelihood p τ (h) with τ The research of Lin Hao was supported by the Science and Technology Project of Weifang University of Science and Technology [grant number KJRC2023012].The research of Il Do Ha was supported by the National Research Foundation of Korea (NRF) grant funded by the Korea government (MSIT) [grant number NRF-2020R1F1A1A01056987].The research of Youngjo Lee was supported by the National Research Foundation of Korea (NRF) grant funded by the Korea government (MSIT) [grant number 2019R1A2C1002408].Jong-Hyeon Jeong's research was supported in part by the National Institute of Health (NIH) [grant numbers 5-U10-CA69974-09 and 5-U10-CA69651-11].
Note: CR: censoring rate, q: number of clusters(centers), n i : cluster size.
Note: CR: censoring rate, q: number of clusters (centers), n i : cluster size.

Table 4 .
(Continued)Simulation results of the estimated dispersion parameters when the distributions of random effect V and random errors k (k = 1, 2) or only the random effect V are misspecified; true dispersion parameters are α = 1, γ = 0.5, φ 1 = 3.1 and φ 2 = 1.5.

Table 5 .
Comparison of simulation results between the proposed (full) and reduced (γ = 0) models under four combinations of α and γ .

Table 6 .
Estimated parameters and standard errors (SEs) from fitting the joint AFT model (full) and reduced model for the multi-center bladder cancer data.

Table 7 .
Estimated parameters and standard errors (SEs) from fitting the joint AFT model (full) and reduced model for the multi-center BMT data.

Table 8 .
Estimated parameters and standard errors (SEs) from fitting the joint AFT model (full) and reduced model for the multi-center breast cancer data.

Table 9 .
Simulation results for the power of the likelihood ratio test over 500 replications under the proposed model with censoring rates of 20% and 40%.