Joint analysis of longitudinal data and competing terminal events in the presence of dependent observation times with application to chronic kidney disease

ABSTRACT In many prospective clinical and biomedical studies, longitudinal biomarkers are repeatedly measured as health indicators to evaluate disease progression when patients are followed up over a period of time. Patient visiting times can be referred to as informative observation times if they are assumed to carry information in addition to that of the longitudinal biomarker measures alone. Irregular visiting times may reflect compliance with physician instruction, disease progression and symptom severity. When the follow-up time may be stopped by competing terminal events, it is possible that patient observation times may correlate with the competing terminal events themselves, thus making the observation times difficult to assess. To explicitly account for the impact of competing terminal events and dependent observation times on the longitudinal data analysis in the context of such complex data, we propose a joint model using latent random effects to describe the association among them. A likelihood-based approach is derived for statistical inference. Extensive simulation studies reveal that the proposed approach performs well for practical situations, and an analysis of patients with chronic kidney disease in a cohort study is presented to illustrate the proposed method.


Introduction
Longitudinal data arise commonly in many clinical and biomedical studies. The repeated measurements in studies are usually biomarkers which can be used as powerful health indicators, e.g. prostate-specific antigen for patients with prostate cancer [41], estimated glomerular filtration rate (eGFR) values for chronic kidney disease (CKD) patients [32,38] and CD4 counts for AIDS patients [40]. Such follow-up subjects typically have additional associated features and complications, e.g. failure resulting from multiple causes during the study period and right censoring. There is an abundance of literature on the analysis of such complex data structure based on joint modeling. This work is motivated by a clinical study of patients with CKD treated at Division of Nephrology of Taichung Veterans General Hospital in Taiwan. Patients with CKD recruited from December 2001 to December 2011 were included in the study, and a follow-up was conducted until July 2012. The eGFR were repeatedly examined for each patient by the modification of diet in renal disease formula at each doctor's visit to assess the progression of the disease. The National Kidney Foundation of the United States considers the eGFR to be the best overall index of kidney function for diagnosing kidney damage and classifying disease severity. In this article, we focus on patients diagnosed with stage 5 CKD (eGFR smaller than 15 ml/min/1.73 m 2 ) at their first visit. The effective sample size was 1096.
The main interest here is to study the decline rate of eGFR. For each patient, there are three primary outcome variables, including eGFR values, observation times and competing terminal events. Details of the data structure are described as follows. As the disease progresses, renal function may deteriorate to a failure situation in which the patient requires renal replacement therapy (RRT), either dialysis or renal transplantation. However, the deterioration of kidney function may also accompany a higher risk of mortality. Patients were followed until either the initiation of RRT or death, i.e. each patient had two competing terminal events, RRT and death. Each patient with stage 5 CKD was required to visit a nephrologist once a month. However, the proportions of gap times in the intervals (0, 1.5], (1.5, 4.5] and (4.5, ∞) were 16.05%, 71.65% and 12.31%, respectively. Most patients' visiting times were irregular, which not only reflects the patients' compliance with physician instruction, but might also depend on the progression of the disease and its symptoms. Thus, visiting times may carry additional information as compared with the longitudinal eGFR measures alone, and are referred to as informative observation times. Moreover, patient observation times may be correlated to the aforementioned two competing terminal events through their underlying health conditions. For such complex data, the typical methodology of analyzing longitudinal measurements, competing risks and recurrent visit times separately is inadequate. For example, Liang and Zeger's [18] generalized estimating equations (GEEs) are a typical method to infer longitudinal data. However, this method could yield biased inferences when the longitudinal measurements are correlated with the observation process [21,33,35]. Consequently, it would be better to handle these data sets simultaneously to explore the relationship of these outcomes, and determine the factors related to the deterioration of kidney function as well as the risk factors for terminal events.
A class of joint models with random effects is proposed for longitudinal data with a univariate survival endpoint [26,36,40,41]. Joint modeling and analysis of longitudinal data with possibly informative observation times has been developed in the literature [2,12,19,21]. Also, extensive statistical literature has been devoted to modeling recurrent events in the presence of a terminal event [11,39,42]. However, joint models to handle longitudinal data in the presence of informative observation times and a dependent terminal event are sparse in the literature. Recently, Liu et al. [24], Liu and Huang [23] and Kim et al. [14] dealt with such complex data by using joint models with random effects and proposed likelihood approaches for inference. Under the framework of joint modeling without specifying the distributions of random effects, Sun et al. [34] developed an estimating equation approach for parameter estimation. However, these methods did not consider the presence of competing terminal events.
Competing risks arise in many fields, including biostatistics, economics as well as engineering and the social sciences. However, as the identification of all cause latent failure times is a well-known problem, the cause-specific hazard functions and cumulative incidence functions [8,9,13,37] are often of primary consideration. When covariates are available, the joint distribution is identifiable under some regularity conditions [1,10]. Consequently, the frailty model and copula model are popular and attractive [3,5,27,30,31] since they assume the existence of latent correlation structures. This paper motivated by the CKD study proposes a joint model for longitudinal measurements of eGFR values and competing terminal event times in the presence of dependent recurrent visit times. We utilize random effects to characterize the association between these processes. That is, the correlation between three types of data, longitudinal measurements data, visiting times and competing risks data, are attributed to unobserved frailties, which may reflect a patient's compliance, attitudes of hygiene and underlying health condition. To this end, we adopt the frailty model for competing terminal events instead of the copula model. We specify a bivariate Weibull model based on correlated frailties for each competing terminal event. The longitudinal outcome process is fitted with a linear mixed-effects model and observation times are described by an intensity function with random effects. Due to the existence of unobserved random effects, a likelihood approach with expectation-maximization (EM) algorithm is often adopted for estimation. However, it would be very computationally time-consuming, especially under the complicated joint model considered here. In this article, we use Gaussian quadrature technique to calculate the likelihood function and obtain maximum likelihood estimates (MLEs), which can be implemented by many statistical software packages.
The paper is organized as follows. In Section 2, the data are described and a new joint frailty model is formulated. A likelihood-based method is developed in Section 3 and extensive simulation studies to assess the performance of the proposed approach are conducted in Section 4. The method applied to the motivating study is presented in Section 5. This article concludes with some discussions and remarks in Section 6.

Notation and model specification
For the sake of simplicity and clarity in developing the model with proper assumptions, we utilize the CKD data from the motivating study to serve as an instructive example throughout this section. The observation time for measuring eGFR is calculated from the patient's first visit to a nephrologist. The eGFR value for each patient was repeatedly measured at each visit and followed up until he reached the initiation of RRT or death. Consequently, besides the patients' background records, the data consist of longitudinal eGFR values, observation times with one of the terminal events. Let t ij , j = 1, . . . , m i , denote the observation times for the individual i starting t i0 = 0, and Y ij denote the longitudinal measurement at time t ij . Here and in what follows let the subscript i be the index for a subject, i = 1, 2, . . . , n. Let Z i be a p-dimensional external covariate. In our motivating example, the longitudinal measurement Y ij is the natural logarithm transformation of eGFR at time t ij . Since the observation times t ij , j = 1, . . . , m i , may be correlated with the longitudinal measurements due to the progression and symptoms of the disease, we have to take them into account via modeling. Let T 1i and T 2i denote the corresponding elapsed times from recruitment until the initiation of RRT and the time of death of a patient, respectively. Since the observation process is stopped by one of two terminal events and might be correlated with both T 1i and T 2i we adopt a random effects model to explain the correlation.
In this article, we treat longitudinal measurements, observation times and competing terminal event times as outcomes of the study, and propose a joint model for inference. In the joint model, two random effects U 1i and U 2i , which are assumed to independently follow N(0, σ 2 1 ) and N(0, σ 2 2 ), respectively, are employed to characterize the correlation of multiple outcomes in the present complex data structure. Given (U 1i , U 2i , Z i ), longitudinal measurements, observation times and competing terminal event times are assumed to be independent. A linear mixed-effect model is specified for longitudinal measurements as follows: where the error term ij is distributed as i.i.d. N(0, σ 2 ) and independent of (U 1i , U 2i ), and γ 0 and γ 1 are fixed-effects parameters. In this model, the correlation between repeated measurements (Y i1 , . . . , Y im i ) is measured by U 1i . Similarly, conditional on U 1i , U 2i and Z i , the intensity function of the observation process at time t is denoted by λ R (t|Z, in which λ R0 (t) is an unknown baseline intensity function and α is the covariate effect. Two competing times T 1i and T 2i , which stop the follow-up, are specified to follow frailty Weibull distributions with the following cumulative hazard functions where exp{η 1 U 1i + η 2 U 2i } and exp{η * 1 U 1i + η * 2 U 2i } are correlated frailties. For k = 1,2, the scale and shape parameters of T ki are exp{Z i β (k) } and κ −1 ki , respectively. To allow the shape parameters to be freely adjusted by covariates, we perform regression on κ ki by This idea of constructing frailty competing risks models follows from Abbring and van den Berg [1] and Heckman and Honoré [10]. Following Lambert et al. [16], Models (3) and (4) can be represented in the following accelerated failure time (AFT) frailty model: where * 1i and * 2i are independent of (U 1i , U 2i ) and follow an extreme value distribution with Pr( * ki > u) = exp{− exp(u)}, k = 1, 2, independently. Moreover, the AFT model presentation gives more direct interpretations on T 1i and T 2i and describes the effects of covariates directly on the length of event times. If ρ (1) (3) and (4) reduce to mixed proportional hazards model [1]. The proposed joint model (1), (2), (3) and (4) is identifiable [1,10,24].
With regard to our CKD data, Model (1) describes the log(eGFR) trajectory for a patient. Model (2) specifies the intensity function of a patient's observation times. Formulae (3)-(5) present the survival models for the competing times T 1i and T 2i , time to receiving RRT and time to death/out of reach, respectively. The random effect U 1i accounts for the withinindividual dependency for the repeatedly measured log(eGFR). The random effect U 2i together with φU 1i specify the within-individual dependence between observation times. Moreover, the sign and magnitude of correlation between log(eGFR), observation times and competing times are explicitly measured by φ, η 1 , η 2 , η * 1 and η * 2 . Let τ be an administrative censoring due to the end of the study and C i be an independent censoring time, e.g. inter-hospital transfer for conveniently continuing CKD treatment in our motivating example. Then the observed follow-up time In this definition, variable i is a censoring indicator, where i = 1 means that one of the terminal events is observed, and δ i indicates the terminal case, where δ i = 1 (or δ i = 2) indicates that the event 1 (or 2) occurs first. For a random sample of n independent subjects, the observed data are In summary, the parameters needed to be estimated are (γ 0 , γ 1 , σ 2 , α, β (1) , ρ (1) , all of which are denoted by θ.

Statistical inference
Under the assumption that all outcomes of interest are independent given the latent random effects and covariates, the complete likelihood function of the observed data and random effects {(u 1i , u 2i ), i = 1, . . . , n} is (2) 1/κ 2i , f (u 1i ; σ 2 1 ) and f (u 2i ; σ 2 2 ) are the probability densities of U 1i and U 2i , respectively. Since U 1i and U 2i are unobservable, it is common to use the EM algorithm [6] to obtain MLEs by treating the random effects as missing data. It can be seen that the EM algorithm is not suitable here, since with this complicated joint model it is necessary to calculate many conditional expectations of random effects given observed data, which usually have no closed forms. Although some Monte Carlo methods (e.g. Metropolis-Hastings algorithm) can be used to approximate such terms in the E-step, the computational cost of numerical methods usually restricts their applicability.
In this article, we calculate the full likelihood function of θ given the observed data to obtain MLEs. As a consequence, the full likelihood function is Here, we calculate Equation (8) directly via Gauss-Hermite quadrature and maximize it to obtain MLEs. The Gaussian quadrature technique is a useful numerical integration tool and advocated in the literature [4,[22][23][24][25]38]. SAS Proc NLMIXED or free software, aML ( [20]; freely available from www.applied-ml.com) can be used to implement Gauss-Hermite quadrature for calculating MLEs. Nonparametric term λ R0 (t) is of infinite-dimensional function. Although we can convert the problem to a parametric one by treating λ R0 (t) as having discrete mass at each observation time, the number of mass points is in the order of n. Thus computing the inversion of a high-dimensional square matrix will result in instability of MLEs and the standard asymptotic theory for MLEs may not apply. To avoid these problems and to obtain stable standard error estimates, we adopt a piecewise constant intensity function to approximate the unknown function λ R0 (t). This piecewise constant intensity is weakly parametric and approximates various shapes of unknown functions well. An abundant body of literature has proposed using piecewise constant functions for approximating unknown functions when analyzing various types of data [7,17,[23][24][25]29]. We adopt a piecewise constant function with J intervals for approximating λ R0 (t). To construct such a step function, divide the space of observation times into J intervals of approximately equal probability: Here, we follow Lawless and Zhan [17], Feng et al. [7] and Liu et al. [24] to set J as 8-10. All the parameters needed to be estimated are (γ 0 , γ 1 , σ 2 , α, β (1) , 1 , β (2) , ρ (2) . We use the hat-accent (∧) on a parameter to indicate the MLE of the parameter. For example,γ 0 is the MLE of γ 0 . The simulation code written in SAS is available upon request.

Simulation study
In this section, we conduct extensive simulations to evaluate the performance of the proposed method. The covariate Z i is generated from Bernoulli distribution with success probability .5. The error term ij , j = 1, . . . , m i , in Equation (1) follows a normal distribution with mean 0 and variance σ 2 = 1 independently, and the baseline intensity function λ R0 (t) is specified as 0.08t. For any subject i, two random effects U 1i and U 2i independently follow normal distributions with mean 0 and variances set to be σ 2 1 = 1 and σ 2 2 = 0.5, respectively. The censoring time C i follows a uniform distribution in (10,15), and the administrative censoring time is chosen as τ = 15.
Two different values for φ, −1 and 1, are considered, and are denoted as Case I and Case II, respectively. Other parameters are listed in Table 1.
Case I indicates that the correlation between longitudinal response outcomes and the intensity of observation process is negative. That is, smaller the Y ij values correlate with higher intensities of repeated measures. Lower measured values are also associated with higher risks of the two terminal events. Moreover, a higher intensity of observation process accompanies higher hazards for two competing terminal events. The setting for (η 1 , η 2 ) and (η * 1 , η * 2 ) suggests a positive association between T 1i and T 2i . However, for Case II, the lower Y ij value is associated with a lower intensity of recurrent event process, but higher hazards for T 1i and T 2i . A lower intensity of recurrent events accompanies higher hazards for the two competing terminal events. But, the two competing risk event times are also positively correlated, as in Case I. From the simulation results, for case I, each subject on the average has 7 repeated measures, but for case II the average number is about 15. About 40% of the total number of subjects received RRT, about 30% died, and the remaining subjects were censored.
We adopt the non-adaptive Gaussian quadrature with 13 quadrature points to calculate the integration (8) and choose a piecewise constant baseline intensity function λ * R0 (t) with J = 10 intervals for approximating λ R0 (t). These 10 intervals are chosen by every 10th quantile of {t ij , j = 1, . . . , m i , i = 1, . . . , n}, and thus each interval (s j−1 , s j ] contains roughly equal numbers of t ij 's. To facilitate good convergence rates and speed up the computing, suitable initial values are necessary. We temporarily assume that all U 1i and U 2i are zeros. By treating (1) as a fixed effect linear model, the ordinary least-squares estimate can be set as an initial value for (γ 0 , γ 1 ). The estimates based on homogeneous Poisson process assumption for recurrent event data {(t i1 , . . . , t im i ,T i , Z i ), i = 1, . . . , n} can be set as initial values of α and λ * j , j = 1, . . . , 10. By temporarily assuming independence between two terminal event times, and setting the initial values of ρ (1) 1 and ρ (2) 1 as zero, the initial values for ρ (1) 0 , β (1) , ρ (2) 0 and β (2) can be obtained by conducting a survival analysis based on the Weibull distribution assumption. All these initial values can be quickly achieved by using commands in many existing statistical software packages, such as the procedures reg and lifereg in SAS for linear regression analysis and survival analysis, respectively. However, for the initial values for φ, η k and η * k , k = 1,2, we recommend using various starting values to check whether the estimates are trapped into local maxima. Totally, 500 replicate data sets are generated for each case with sample sizes n = 200 and n = 400. The simulation results are summarized in Tables 2 and 3. In both tables, the parameter estimates (EST), corresponding sample standard deviations (SE), averages of estimated standard errors (ASE) as well as the empirical coverage probabilities of 95% Wald's confidence intervals (CP) are reported. It can be seen that the overall performance of the proposed estimation method appears to be satisfactory. All parameter estimates are nearly unbiased and the standard error estimates agree well with the simulation standard errors of the parameter estimates. Moreover, as the sample size increases, the biases and variances decrease, and the coverage percentages get closer to the nominal level. For investigating the biases induced by ignoring the other competing terminal event, all the generated data sets are also analyzed by an approach similar to that described above. Thus, the competing terminal event 2 is treated as an independent censoring event and model (3) is still fitted to T 1i . It can be seen that there exists biases onφ andσ 2 1 in the approach ignoring the other competing terminal event and that the CP values are significantly lower than the nominal value. In some cases,β (1) 0 andρ (1) 0 also have poor CP performance.

An illustrative example
In this section, we analyze the CKD data set described in Section 1 by the proposed method.
The study is a historically prospective clinical investigation of CKD patients assessing the effect of risk factors on the deterioration of kidney function, need for RRT and mortality. One of the main objectives of the present paper is to determine the rate of eGFR decline among patients who were diagnosed with stage 5 CKD at their first visit, as they have a high risk for receiving RRT. In total, there were 1096 patients included. Since nephrology nurses contacted all the patients to ensure that they attended their scheduled follow-ups and provided complete information, we treat the patients who went out of reach as informative dropouts. In this study, we treat out of reach and death as the same cause of failure. Among these patients, 407 (37.14%) were observed to receive RRT and 136 (12.41%) experienced another terminal event, mortality/out of reach, where 38 died and 98 went out of reach. A total of 22 patients transferred to other hospitals for continuing CKD treatment, and it is reasonable to treat hospital transfer as an independent source of censoring. For all the 1096 patients, the average follow-up time was 17 months with a standard deviation of 14.00. Among the 407 patients who received RRT, the average number of visits per subject was 3.99 with a standard error of 3.43 (the median and the range were 3 and 24, respectively). However, among the 136 patients who either died or went out of reach, the average number of visits per subject was 2.52 with a standard error of 2.27 (the median and the range were 2 and 11, respectively). The Welch two sample t-test resulted in a p -value of 3.1 × 10 −8 , indicating that the cause of terminating may be significantly related to the pattern of visits in the observation process. For all the 543 uncensored patients, Figure 1 shows the eGFR measures against observation times from their enrolment until terminating follow-up. It can be seen that the pattern for patients finally receiving RRT is different from that for patients who died or went out of reach in respect to their follow-up times and observation times. Consequently, we have to consider the pattern of visits and competing terminal events simultaneously.
The main risk factors of interest are the binary indicators for diabetes (comorbidity rate 41%), hypertension (comorbidity rate 76%), hyperlipidemia (comorbidity rate 11%), male gender (50%) and advanced age at entry (age ≥ 65, 53%), all of which were collected at the patient's first visit. All of these five baseline covariates are included in models (1), (2), (3) and (4). As mentioned in the beginning, logarithmic transformation of eGFR is taken as the longitudinal outcome to reduce extreme values of eGFR. Since the emphasis is placed in the rate of eGFR decline, the observation time in months is also included as one of the covariates in model (1). The interaction of gender and time is also considered in model (1) to evaluate the effect of gender on the rate of decline, which is one of medical interest. As in the simulation studies, the likelihood value is calculated by using non-adaptive Gaussian quadrature with 13 quadrature points and adopting piecewise constant function with 10 intervals for approaching λ R0 (t).
The results of the analysis based on our model are summarized in Tables 4 and 5. The linear mixed-effect model of longitudinal data provided sufficient evidence to show that the log(eGFR) value decreases with time (slope = −0.011, p-value< 0.001) and the rate of decline in log(eGFR) for females is different from that for males (p-value = 0.007). These findings imply that females experienced a 12.06% decrease in eGFR per year, while males experienced a 9.69% annual decrease. After adjusting for other factors, patients with hyperlipidemia were found to have significantly larger eGFR values than patients without hyperlipidemia, and female patients were found to have significantly lower eGFR values than male patients. Table 4 also indicates that patients diagnosed with stage 5 CKD at younger ages (age < 65) had lower eGFR values than patients who were diagnosed at older ages (age > 65). This means that the patients who were diagnosed with stage 5 CKD at a younger age (age < 65) suffered from more severe damage to their kidney function than patients who reached stage 5 at an older age (age > 65). However, the effects of concomitant diabetes and hypertension on eGFR values were nonsignificant. None of the covariates had a significant effect on either the intensity function of the observation process or hazard function of death/out of reach. Nevertheless, for T 1i , time to initiation of RRT, the scale component results show that the patients concomitant with diabetes, hypertension or diagnosed with stage 5 CKD at younger ages had significantly smaller scales in T 1i . However, the corresponding shape component results indicate that patients with hypertension or diagnosed with stage 5 CKD at younger ages had larger shape parameter estimates for T 1i . According to Equations (6) and (7), patients with diabetes have a significant shift in the distribution of their log(T 1i ) to the left when compared with patients without diabetes. Similarly, the distribution of log(T 1i ) for patients with hypertension also has a shift to the left when compared with patients without hypertension. The same conclusion can be drawn for patients diagnosed with stage 5 CKD at younger ages than those diagnosed at older ages. However, the log(T 1i ) of patients with hypertension has a smaller dispersion than that of those without hypertension. Patients diagnosed with stage 5 CKD at younger ages also had smaller dispersions in log(T 1i ) than those diagnosed at older ages. Consequently, from the perspective of the linear model for log(T 1i ) in Equation (6), diabetes, hypertension and having been diagnosed with stage 5 CKD at younger ages resulted in significantly larger shifts to the left, and hypertension and having been diagnosed with stage 5 CKD at younger ages correspond significantly with smaller variations in log(T 1i ). The risk factor diabetes may tend to engender smaller variation since its p-value is 0.069. Since diabetes is a serious chronic disease, and one of principal concern of Taiwanese physicians, we discuss this disease in more detail in the following. There were 446 patients with diabetes in the dataset used in this study, among whom 39%, 14% and 47% were receiving RRT, dead and censored, respectively. Among the 650 patients who did not have diabetes, the respective proportions of patients receiving RRT, dead and censored were 36%, 11% and 53%. The joint model analysis results reveal that although diabetes was not associated with any significant change in eGFR value, patients with diabetes did receive RRT earlier than those without diabetes. Regarding mortality, the estimated effects of diabetes implicitly reveal that the diabetic patients tended to have smaller log(T 2i ) with smaller variation, although the effects are nonsignificant. This finding can be demonstrated by depicting the eGFR measurement trajectories of all the patients with or without diabetes which appears as supplementary material on the Journal of Applied Statistics website. It reveals that although there is no significant difference in the pattern of eGFR measurement The association analysis did not reveal a significant correlation between the log(eGFR) values and observation process (φ = 0.059, p-value = .527). However, the smaller log(eGFR) values tended to increase the risk of receiving RRT (η 1 = −2.579, p-value < .001 ). The nonsignificant estimateη * 1 means that the log(eGFR) values might be uncorrelated with the risk of death/out of reach. The significantly positive estimateη 2 implies a higher intensity for the observation process corresponding to a higher risk for receiving RRT. But the significantly negative estimateη * 2 means that a higher intensity of visits corresponds with a smaller risk of mortality/out of reach. Finally, the two competing terminal event times T 1i and T 2i were demonstrated to have a significant negative correlation. One possible reason for this correlation is that a patient with a higher hazard for receiving RRT tends to visit nephrologists more often. As a consequence, the patient receives more care, which may decrease the hazard of mortality/out of reach. The proposed joint model gives investigators a comprehensive understanding of the relationship between disease development and other risk factors.
In order to check our model adequacy, we performed a residual analysis for the proposed model. By calculating empirical Bayes' estimates for random effects, denoted by (Û 1i ,Û 2i ), i = 1, . . . , n, we obtained residuals for the longitudinal model (1). Plotting the residuals against each covariate can be a model-checking method. Similarly, under the empirical Bayes' estimates for random effects, the martingale residuals for the observation process (9) are used to assess the goodness-of-fit of model (2). A graphical method is proposed to assess the appropriateness of the fitted model for competing risk data. By utilizing simple algebraic manipulation of models (6) and (7) and the independence assumption for ( * 1i , * 2i , U 1i , U 2i ), * 1i can be demonstrated to be subject to usual independent random censoring by which is a function of ( * 2i , U 1i , U 2i ). As a consequence, we have two right-censored data sets, {(r 1i , i I (δ i = 1)), i = 1, . . . , n} and {(r 2i , i I(δ i = 2)), i = 1, . . . , n}, where κ 1i +η 1iÛ1i +η 2iÛ2i , κ 2i +η * 1iÛ 1i +η * 2iÛ 2i , The Kaplan-Meier survival curve of {(r 1i , i I (δ i = 1)), i = 1, . . . , n} is a reasonable graphical display for model checking. If the Kaplan-Meier survival curve is close to the assumed survival distribution of * 1i , the fitted model (3) should be adequate. The same logic is applied to model (4) based on {(r 2i , i I(δ i = 2)), i = 1, . . . , n}. Some graphics of the residuals are shown in Figure 2. For assessing the adequacy of the linear mixed-effects model, Figure 2 (a) shows the residual plot for longitudinal measurements against observation times, and Figure 2 (b) presents the corresponding lowess smooth curve (solid black) and horizontal line at 0 (dashed red). Both lines indicate that our longitudinal model fits the data well. The boxplots of martingale residuals (9) against diabetes are shown in Figure 2 (c), which indicates that the fitted model (2) is adequate. Boxplots against other covariates (not reported here) display similar results. Figure 2 (d) examines the goodness-of-fit for model (3), which shows the Kaplan-Meier survival curve of {(r 1i , i I(δ i = 1)), i = 1, . . . , n} by a solid curve and the assumed survival function of an extreme value random variable by a dashed one. The closeness of both curves suggests that the proposed model is adequate. However, the discrepancy between the Kaplan-Meier survival curve of {(r 2i , i I(δi = 2)), i = 1, . . . , n} and the assumed extreme value survival curve appearing as supplementary material in the Journal of Applied Statistics website is noted, which means that the assumed model (4) may not fit the second-caused (mortality/out of reach) competing risk data well. The discrepancy may result from pooling mortality and out of reach as a cause. To remedy this, we may consider the CKD data with three competing terminal events, receiving RRT, mortality and out of reach. However, the number of expired patients is too small to provide enough information for estimation. It would be worthwhile to continuously collect data for providing more information for the inference of a mortality event.
In order to understand the effectiveness of the proposed model, we analyze the data with two reduced models as well. First, we analyze longitudinal eGFR data only fitted by model (1), not considering dependent observation times and competing terminal events. Since E(Y ij ) = γ 0 + Z i γ 1 and Cov(Y ij , Y ij ) = Var(U 1i ), j = j , we analyze the data with a generalized estimating equation (GEE) approach with an exchangeable correlation structure [18]. The results shown in Table 4 are similar to those obtained through our proposed model, except for the interaction of gender and time which was found to be nonsignificant. This means that the difference in the rates of decline between males and females cannot be detected if the information of dependence of observation times and competing terminal events is ignored. We also fit longitudinal eGFR data by a reduced joint model without considering death/out of reach, i.e. treated as an independent censoring. The results shown in Tables 4 and 5 reveal that the covariate effect estimates and the corresponding p-value in longitudinal model (1) are similar to that of our joint model. The major difference is that this reduced model could not detect the significant effect of hypertension on the time of initiation of RRT, neither in the scale nor shape components.

Concluding remarks
The longitudinal data with two competing terminal events in the presence of dependent observation times are studied in this article. The observation times may depend on a patient's health status. Different terminal events are also associated with the observation patterns and longitudinal measurements. This type of longitudinal data set arises commonly in many clinical studies, as well as in engineering and the social sciences. For example, when we are interested in the maintenance costs of a vehicle, the owner may not maintain his car regularly, which may reflect not only the owner's attitude but also the underlying condition of the car. Consequently, we cannot ignore the potential correlation between the maintenance process and the longitudinal cost process. Furthermore, each car eventually experiences one of several competing terminal events, e.g. sold or scraped, which may also relate to the car's underlying condition. Our proposed model and method can be applied for this type of data. This joint model distinguishes the mutual relationship and stresses the importance of awareness of different dependent terminal events, as emphasized by Koller et al. [15]. To avoid computational complexity induced by the EM algorithm, the Gauss-Hermite quadrature technique is adopted for calculating the full likelihood to obtain MLEs, which can be simply conducted with existing statistical software. The simulation studies demonstrate that the proposed estimation approach performs well and that the method is suitable for the longitudinal biomarker measurements of CKD patients.
Although only two competing causes of failure are considered in the joint model, the approach derived in this article can be easily extended to that with an arbitrary, but known finite number of competing terminal events. The proposed model can also easily incorporate time-dependent external covariates. Let Z i (t) be an external covariate process. The model for longitudinal data is Y ij = γ 0 + Z(t ij ) γ 1 + U 1i + ij , j = 1, . . . , m i . For Equation (2), we consider the cumulative intensity function of the observation process as follows: where R0 (t) is an unknown cumulative baseline intensity function. Similarly, Equations (3) and (4) can be extended to e Z i (t) β (2) 1/κ 2i (t) , in which κ 1i (t) = exp{ρ (1) 0 + Z i (t) ρ (1) 1 } and κ 2i (t) = exp{ρ (2) 0 + Z i (t) ρ (2) 1 }. Then our approach is still valid.
In the present article, we adopt normal random effects for characterizing the association between three types of outcomes. In fact, there are other types of distributions to be considered in the joint models, for example, the log-gamma distribution. For identifiability, the mean of each random effect is required to be one. However, one commonly used random effect in the literature is the normal distribution. Based on our approach, practitioners and clinicians can assume other distributions for U 1i and U 2i and apply the Gaussian quadrature technique with the probability integral transform [28] to obtain the likelihood function (8) and the corresponding MLEs. The developed model-checking method in Section 5 can also be adopted to assess the model adequacy. Although specifying semiparametric structures for terminal time models (3) and (4) is feasible, it would lead to three unknown infinite-dimensional parameters in the joint model, which would certainly induce instability in the estimation. To this end, a Weibull model is proposed instead and flexibly allows investigators to adjust for covariates on the scale and shape components. But how to extend our model for more flexible classes and develop a feasible estimating approach will be a topic for future research.
In many situations, each patient has more than one biomarker to be measured simultaneously, which together reflect a patient's disease progression and should be treated as internal variables. This leads to multivariate longitudinal data. Consequently, there are numerous statistical challenges for such complex data sets containing multivariate longitudinal measurements, dependent observation times and competing terminal events. In this case, constructing a joint model with random effects is a feasible way to handle multiple outcomes simultaneously. It is worthy of future research.