Two-Step Estimation Procedure for Parametric Copula-Based Regression Models for Semi-Competing Risks Data

TL;DR

This paper introduces a new two-step method for analyzing complex data involving non-terminal and terminal events using copula-based regression.

Contribution

A novel two-step algorithm combined with BOBYQA is proposed for robust estimation in semi-competing risks regression models.

Findings

The proposed method effectively reduces the influence of initial values in complex copula structures.

Simulation results confirm the method's performance improvements.

Application to the Amsterdam Cohort Study showed practical improvements.

Abstract

Non-terminal and terminal events in semi-competing risks data are typically associated and may be influenced by covariates. We employed regression modeling for semi-competing risks data under a copula-based framework to evaluate the effects of covariates on the two events and the association between them. Due to the complexity of the copula structure, we propose a new method that integrates a novel two-step algorithm with the Bound Optimization by Quadratic Approximation (BOBYQA) method. This approach effectively mitigates the influence of initial values and demonstrates greater robustness. The simulations validate the performance of the proposed method. We further applied our proposed method to the Amsterdam Cohort Study (ACS) real data, where some improvements could be found.

Full text

1. Introduction

In follow-up biomedical studies, objects may experience two types of events: a non-terminal event (e.g., disease relapse) or a terminal event (death). The occurrence of a non-terminal event can preclude the occurrence of a terminal event, but not vice versa. This is referred to as semi-competing risks (Fine et al. [1]).

There are two conventional statistical methods for analyzing semi-competing risks data: the random effects model and the copula-based model. As in the study by Xu et al. [2], the random effects model usually adopts a shared frailty non-negative random variable to generate the joint distribution. In contrast, Shih et al. [3] showed that the copula-based model constructs the joint distribution by connecting two marginal distributions, and the copula parameters can explicitly estimate the dependence strength. For the copula-based model, the marginal distribution of times to events conditioned on covariates can be specified independently, and the covariate effects on marginal events can be accessed directly. Moreover, Ramadan et al. [4] highlighted that the copula-based model can capture asymmetric and nonlinear connections between variables and simplify the computation of certain probability measures. Using the copula function, Fayomi et al. [5] proposed the BFGMLG family to model bivariate continuous data with heavy-tailed and skewed distributions; they applied this to environmental, medical, and computer science data. Haj Ahmad et al. [6] proposed a bivariate modified extended exponential (BMExE) distribution and adopted this model to study the relationship between processor and memory reliability in data science. Ramadan et al. [4] introduced the BAPB-XII distribution and applied it to engineering data. In this article, we consider the copula model for semi-competing risks data.

Fine et al. [1] proposed estimators of the association parameter from the concordance-estimating function and used a novel plug-in estimator for the unspecified marginal distribution of a non-terminal event. Jiang et al. [7] proposed another copula-based estimator with better properties based on the principle of pseudo-self-consistency. Peng and Fine [8] considered a class of time-varying coefficient semiparametric regression models for marginals and copulas, and they proposed an estimator based on the moment method. Hsieh and Huang [9] analyzed semi-competing risks data based on the conditional likelihood approach. They adopted semiparametric transformation models as the marginal regressions, which include the popular proportional hazards and proportional odds models as special cases.

For the copula-based model, Shih and Louis [3] proposed a two-stage procedure for the estimation of parameters in the marginal distribution and copula function. In the first stage, the parameters for the marginal distributions are estimated independently for each event, disregarding the dependencies between events; in the second stage, the parameter associated with the copula is estimated via maximizing the log-likelihood using the plug-in estimation of the marginals. However, this approach may not be suitable for the copula-based models for semi-competing risks data, as a terminal event may dependently censor the non-terminal event. Chen [10] considered a nonparametric maximum likelihood estimation (MLE) approach, which involves simultaneously estimating the parameters that are associated with the marginals and the copula via maximization of the log-likelihood function. Arachchige et al. [11] provided a two-stage pseudo-maximum likelihood estimation approach (PMLE). More specifically, in the first stage, they obtained a consistent estimate of the marginal of the terminal event; in the second stage, they estimated the marginal of the non-terminal event and the copula parameters simultaneously by maximizing the pseudo-likelihood function with the plug-in estimation of the terminal event. Sun et al. [12] proposed a two-step method: in the first step, they obtained the parameter estimation separately as the initial value while, in the second step, they obtained all parameters of the marginal distribution and copula function simultaneously through a joint distribution. They called this procedure a two-step procedure. This method can mitigate the influence of initial values, and is adopted in our article.

In addition, due to the complex structure of a copula-based objective function, the derivative information is unavailable, unreliable, or computationally expensive to obtain. In such cases, derivative-free optimization (DFO) methods serve as a viable alternative. This article adopts the Bound Optimization by Quadratic Approximation (BOBYQA) algorithm proposed by Powell [13]. As summarized by Ragonneau [14], Powell pioneered his first DFO method using conjugate directions in 1964. From the 1990s onward, he developed five model-based DFO algorithms, each tailored to distinct problem classes: COBYLA for nonlinearly constrained problems, UOBYQA and NEWUOA for unconstrained problems, BOBYQA for bound-constrained problems, and LINCOA for linearly constrained problems. However, a key limitation of BOBYQA is its limited scalability to high-dimensional problems.

In this article, we consider the parametric copula-based regression model for semi-competing risks data. Fully parametric approaches are appealing because they are generally identifiable, can produce smooth survival functions, and give a direct interpretation [3,15]. Within the realm of parametric regression copula-based models, traditional survival distributions—particularly the exponential, Weibull, and Gompertz distributions—remain significant, both theoretically and practically [16]. The primary challenge in parametric copula-based modeling lies in the optimization process, which is complicated by the intricate likelihood structure of semi-competing risks data. This complexity makes the estimation particularly sensitive to the choice of initial values. We implement a two-step estimation approach to address this issue and reduce the dependence on initial value selection. First, we employed the method developed by Shih and Louis [3] to obtain preliminary parameter estimates. These estimates are then used as initial values for the second step maximum likelihood estimation based on the joint likelihood of the two marginal events, thereby improving the stability and reliability of the final parameter estimates. We employ the BOBYQA optimization algorithm and demonstrate the competitive performance of our proposed method in comparison to other existing methods through simulations and real data on HIV/AIDS from the Amsterdam Cohort Study (ACS).

The remainder of this article is organized as follows. In Section 2, we provide a brief description of the motivating data. In Section 3, we introduce the fully parametric copula-based regression models that utilize various copula functions and survival functions. We also explain the main procedure of the state-of-the-art BOBYQA optimization method. Subsequently, we applied our method to the real data on HIV/AIDS from the Amsterdam Cohort Study (ACS) (Section 4) and compared our method with other optimization methods in Section 5. The discussion is given by concluding the article in Section 6.

2. Data

The Amsterdam Cohort Study (ACS) on HIV/AIDS among homosexual men is a progressive cohort study that commenced in 1984 to investigate HIV/AIDS. This dataset consists of the time (in years) from HIV infection to AIDS, virus phenotype switching from non-syncytium-inducing to syncytium-inducing (SI switch), and mortality data for 329 men who have sex with men (MSM). Data from the period until combination anti-retroviral therapy are available (1996). The censoring rates for the SI switch event and mortality are approximately 65.7% and 45.9%, respectively. After removing the five partially observed cases, the final dataset comprises 324 samples, and a summary of the variables is listed in Table 1.

This dataset is included in the mstate library and has been analyzed by Geskus et al. [17,18] and Sorrell et al. [19]. However, they only consider the association between events. In this article, we are also interested in the effect of CCR5 and age on the time to SI switch and death events. Due to the sample size limitation, the parametric model may be more suitable, and we specify the marginal distribution for SI switch event and death event.

3. Methods

3.1. The Copula-Based Model Specifications

In this section, we introduce the model for semi-competing risks events based on the conditional copula function proposed by Patton [20]. For the semi-competing event, we denote $[eqn]$ as the time to a non-terminal event and $[eqn]$ as the time to a terminal event. Moreover, the times $[eqn]$ are subject to right censoring $[eqn]$ . Then, the time to first event is $[eqn]$ , with the censoring indicator $[eqn]$ if $[eqn]$ ; otherwise, it is $[eqn]$ . The time to the second event is $[eqn]$ , with the censoring indicator $[eqn]$ if $[eqn]$ ; otherwise, it is $[eqn]$ .

We assume the parameters of the copula and marginal survival distribution depend on the covariate vector $[eqn]$ , where p is the dimension of covariates. The marginal survival function is denoted as $[eqn]$ for the non-terminal event and $[eqn]$ for the terminal event. For the joint survival distribution, we employ a bivariate copula function $[eqn]$ with a univariate parameter $[eqn]$ :

$$[eqn]$$

Subsequently, the density function for the copula is as follows:

$$[eqn]$$

and the density corresponding to the survival function can be written as

$$[eqn]$$

where $[eqn]$ , $[eqn]$ and $[eqn]$ are the marginal densities of $[eqn]$ and $[eqn]$ , respectively.

In this article, we mainly focus on four commonly used univariate-parameter copulas: Clayton, Frank, Gumbel, and normal. To incorporate covariate $[eqn]$ effects while maintaining parameter constraints, we employ distinct link functions corresponding to each copula’s parameter space. The four copula families and their respective parameter link functions are formally specified as follows:

• (i)Clayton copula:

$$[eqn]$$

• (ii)Frank copula:

$$[eqn]$$

• (iii)Gumbel copula:

$$[eqn]$$

• (iv)Normal copula:

$$[eqn]$$

where $[eqn]$ is the inverse of the cumulative distribution function of the standard univariate Gaussian distribution.

The above three copula functions (i–iii) belong to the Archimedean copula family, and the normal copula (iv) belongs to the elliptical copula family. For more details, we refer the reader to Nelsen [21].

3.2. The Parametric Model for Event to Time

For the marginal distribution of non-terminal and terminal events, we consider three popular parametric distributions: exponential, Weibull, and Gompertz. The covariates’ effects on marginal distribution are specified as an exponential linear regression model.

For the exponential marginal distribution, the survival function is given as $[eqn]$ for non-terminal and terminal events, respectively. The relationship between parameter $[eqn]$ and covariates $[eqn]$ is characterized by the following function:

$$[eqn]$$

where $[eqn]$ are regression coefficients for non-terminal and terminal events, respectively, for $[eqn]$ . We apply the exponential link function to guarantee that $[eqn]$ is positive.

Concerning the Weibull distribution, the survival function can be written as $[eqn]$ , where $[eqn]$ and $[eqn]$ are scale and shape parameters, respectively, for $[eqn]$ . We assume the shape parameters $[eqn]$ and $[eqn]$ are constant, and the scale parameters $[eqn]$ and $[eqn]$ are determined by $[eqn]$ as

$$[eqn]$$

for $[eqn]$ , where $[eqn]$ are regression coefficients.

As for the Gompertz marginal distribution, the survival function can be written as $[eqn]$ for $[eqn]$ , where $[eqn]$ is the rate parameter, and $[eqn]$ is the shape parameter. As is the case for the Weibull distribution, let the shape parameter $[eqn]$ be constant and the rate parameter $[eqn]$ depend on $[eqn]$ in the following way:

$$[eqn]$$

for $[eqn]$ , where $[eqn]$ are regression coefficients.

3.3. Joint Likelihood Function

For the semi-competing bivariate right-censored data, the likelihood based on the copula function $[eqn]$ is expressed as

$$[eqn]$$

where n is the sample size, and $[eqn]$ is the density function for the copula defined in Equation (2). We assume that the censoring times are independent of the survival times and are not informative.

Here,

$$[eqn]$$

and

$$[eqn]$$

are called h-functions in the literature [22], and the censoring indicator is $[eqn]$ if $[eqn]$ ; otherwise, it is $[eqn]$ . The censoring indicator is $[eqn]$ if $[eqn]$ ; otherwise, it is $[eqn]$ . The exact form of the likelihood for the four types of copula functions previously mentioned is given by Wei et al. [23].

3.4. Estimation

From Equation (11), we can see that the structure of the log-likelihood function is relatively complex; therefore, it is often a computationally challenging task to estimate the unknown parameters. We define $[eqn]$ as a vector of the parameters, where $[eqn]$ ( $[eqn]$ ) is the constant parameter of marginal distribution, $[eqn]$ ( $[eqn]$ ) are regression coefficients, and b denotes the regression coefficients in the linear predictor of the link function associated with the copula association parameter $[eqn]$ mentioned in Section 3.3. It is known that the terminal event can censor the non-terminal event, so the parameter estimation of the non-terminal event may be biased when using the two-stage method because it treats the censored time as non-informational. Thus, we leverage the popular two-step estimation procedures by Sun et al. [12,24]. The method is shown in Algorithm 1. Algorithm 1 Two-step estimation procedure.Step 1: Estimate the parameters of marginal distributions and copula function separately.

• (1)First obtain the initial values of $[eqn]$ through grid search method. Then estimate these parameters based on

$$[eqn]$$

where $[eqn]$ denotes the likelihood for the marginal survival distribution, $[eqn]$ .

• (2)Obtain the initial values of $[eqn]$ through grid search method. Then estimate it by

$$[eqn]$$

where $[eqn]$ and $[eqn]$ are the estimated values by above procedure (1), $[eqn]$ , and $[eqn]$ is the joint likelihood based on the copula function. Step 2: Simultaneously maximize the joint likelihood

$$[eqn]$$

with the initial values of $[eqn]$ , $[eqn]$ , $[eqn]$ and $[eqn]$ attained from Step 1.

In Algorithm 1, the grid search method establishes initial values for parameters $[eqn]$ and $[eqn]$ through the following systematic procedure. First, define search ranges for each parameter: $[eqn]$ and $[eqn]$ . Second, discretize each interval into n equally spaced points: $[eqn]$ and $[eqn]$ . Third, generate all possible parameter combinations from the Cartesian product of $[eqn]$ and $[eqn]$ , resulting in $[eqn]$ distinct pairs $[eqn]$ , where p denotes the dimension of $[eqn]$ . Finally, evaluate the marginal survival distribution likelihood $[eqn]$ , $[eqn]$ and $[eqn]$ and identify the optimal parameters $[eqn]$ . Apply an analogous procedure to determine the initial values for b. The parameter estimates derived from Step 1, formally expressed as $[eqn]$ serve as the initialized parameter vector for the joint optimization process in Step 2.

The accuracy and stability of parameter estimation are critically dependent on appropriate initialization. Direct likelihood maximization in single-step procedures frequently encounters convergence failures when using arbitrary initial values, particularly in complex models. In contrast, the two-stage procedure can be an effective initial method. Given the complexity of the copula itself, the estimation of parameters for the copula and marginal distributions, as outlined in the maximization Equations (12)–(14), is a challenging task. Therefore, we adopt the BOBYQA optimization algorithm instead of the novel limited memory Broyden–Fletcher–Goldfarb–Shanno bounded (L-BFGS-B) method.

The BOBYQA algorithm is a model-based derivative-free trust-region method proposed by Michael J.D. Powell [13]. It is an iterative algorithm used to obtain a minimum of a function $[eqn]$ , where p is the dimension of variable $[eqn]$ , subject to some box bounds on the variables. It mainly consists of building or updating the interpolation set, constructing or updating the interpolation-based quadratic approximation function, and solving the trust-region subproblem.

At iteration k, it employs the interpolation method to construct or update the quadratic model $[eqn]$ ( $[eqn]$ ) approximation to the objective function $[eqn]$ ( $[eqn]$ ) without using its derivatives. For a selected collection of points $[eqn]$ , $[eqn]$ , where m is a fixed integer from the interval $[eqn]$ , and a local approximation $[eqn]$ of objective $[eqn]$ , is constructed as follows:

$$[eqn]$$

and $[eqn]$ satisfies the interpolation conditions

$$[eqn]$$

As the number m of interpolation points is less than $[eqn]$ (the necessary number for a pure interpolation method) and $[eqn]$ is recommended by Powell [13,25], the method uses the remaining degree of freedom by solving

$$[eqn]$$

where $[eqn]$ denotes the square of the Frobenius norm, i.e., $[eqn]$ [26]. Further, based on the established quadratic approximation function $[eqn]$ , it solves its trust-region subproblem. Equally important, the set of m interpolation points is also updated during the iterations. The BOBYQA algorithm can be applied to solve functions with large p parameters, as mentioned by Ragonneau [27]. When compared with the L-BFGS-B method, the strength of the BOBYQA algorithm is that it can approximate functions whose Hessian matrix is a singular matrix, as also manifested in Section 4. For more information, see [13]. The algorithm is included in the R “optimx” library [28], a replacement and extension of the “optim” function in the “stats” library [29].

As mentioned above, the two-stage procedure refers to an algorithm that optimizes parameters separately, and the two-step procedure denotes an algorithm that optimizes parameters jointly. To evaluate the performance of our proposed TStep-BOBYQA, we compared it with three other methods: TStage-BFGS, TStep-BFGS, and TStage-BOBYQA. TStage-BFGS refers to the optimization of parameters through a two-stage procedure that combines the L-BFGS-B optimization algorithm with initial values obtained via grid search. TStep-BFGS refers to the optimization of parameters using a two-step procedure that integrates the L-BFGS-B optimization algorithm with initial values derived from TStage-BFGS. TStage-BOBYQA refers to the optimization of parameters through a two-stage procedure that combines the BOBYQA optimization algorithm with initial values obtained through grid search. In contrast, TStep-BOBYQA utilizes a two-step procedure coupled with the BOBYQA optimization algorithm, employing initial values sourced from TStage-BOBYQA. These methods are summarized in Table 2.

4. Application

In this section, we apply the methods described in Section 3 to the ACS data presented in Section 2. This dataset includes the time to SI switch, time to death from AIDS, the continuous covariate age, and the categorical covariate CCR5. As mentioned in Section 3.2, we utilized three various parametric distributions (i.e., exponential, Weibull, and Gompertz distributions) to model the marginal distribution of survival time. We employed the four different copula functions: Clayton, Frank, Gumbel, and normal, to assess the association between the two survival endpoints. Finally, we selected the best-fitting model based on the Akaike information criterion (AIC). We also assume that the hazard rate and the association between the two events all depend on covariates. We estimate the hazard ratios of the non-terminal and terminal events as well as the regression coefficients of the covariates for the copula function parameters. To eliminate the influence of units, we standardized the times to events and age using a min-max normalization method.

The estimated hazard ratios for each covariate and the regression coefficients of the covariates for the association parameter are reported in Table 3, Table 4, Table 5 and Table 6. The hazard ratios and regression coefficients for the association parameter are similar between the L-BFGS-B methods (TStage-BFGS; TStep-BFGS) and the BOBYQA methods (TStage-BOBYQA; TStep-BOBYQA), with some exceptions. For the Clayton–Gompertz model, parameter estimation can not be obtained using TStep-BFGS (see Table 3). For the Frank–Gompertz model, the variance estimation of regression coefficients of association can not be obtained using TStage-BOBYQA (see Table 4). Additionally, for the Weibull distribution, although the estimated parameters are comparable, the confidence intervals for parameters estimated by the L-BFGS-B methods (TStage-BFGS; TStep-BFGS) are wider than those obtained through the BOBYQA methods (TStage-BOBYQA; TStep-BOBYQA), as illustrated in Table 3, Table 4, Table 5 and Table 6. The methods mentioned above are briefly outlined in Table 2.

From Table 3, Table 4, Table 5 and Table 6, we can see that the two-step methods (TStep-BFGS and TStep-BOBYQA) provide lower AIC values than the two-stage methods (TStage-BFGS and TStage-BOBYQA), with the exception of the copulas combined with the Weibull distribution models. The parameters estimated by the two-stage procedure are identical to those estimated by the two-step procedure, resulting in equivalent AIC values.

4.1. Hazard Ratios and Associations

SI following HIV infection. Based on the estimation of four different copula functions combined with three different marginal distribution models, no association was found between age and SI switch. Similarly, no association was observed between CCR5 and SI switch, except in two models: the Frank–exponential model (HR:0.705, 95% CI: 0.425 to 0.986, obtained through TStep-BFGS and TStep-BOBYQA; see Table 4) and the normal–Gompertz model (HR: 0.698, 95% CI: 0.402 to 0.995, using TStep-BFGS and TStep-BOBYQA; see Table 6). Thus, the “WM” groups are shown to have a lower risk of SI switch.

Death following HIV infection. According to the results of all models, age is not found to be associated with death, as the value of 1 is within their HR confidence intervals. We can see that the CCR5 HR is less than the value of 1, and their HR confidence intervals do not include a value of 1; this means that it is associated with death, and the “WM” groups have a lower risk of death compared to the “WW” groups. These conclusions are consistent for the four different copula models.

Association between SI and death. From the results of all models, the association between SI switch and death is stronger for the patients with the CCR5 “WM” genotype than for those with the CCR5 “WW” genotype, except for cases including the Frank–Weibull model, Frank–exponential model, and normal–exponential model. Age does not appear to influence the association for both SI switch and death events, as the value of 0 is within the CI of the age parameter coefficient estimation.

4.2. Results for the Preferred Model

For the data from the Amsterdam Cohort Study on HIV infection and AIDS, we employed the AIC to select the optimal model. Our comparative analysis reveals that the two-step estimation method consistently demonstrates superior performance, yielding lower AIC values compared to the two-stage optimization approach in most situations. Among the copula models, the normal copula exhibits the most favorable AIC outcomes, as per Table 3, Table 4, Table 5 and Table 6. Moreover, the Gompertz distribution presents the lowest AIC values within each copula model. These systematic comparisons conclusively identify the normal–Gompertz copula model as the optimal choice for the ACS data.

As demonstrated by the normal–Gompertz survival model in Table 6, individuals with the CCR5 “WM” genotype exhibit significantly reduced risks of both seroconversion illness (SI switch) and mortality compared to the “WW” genotype group. Specifically, the hazard ratios are 0.698 (95% CI: 0.402–0.995) for SI switch and 0.406 (95% CI: 0.235–0.578) for death. This result aligns with the genetic mechanism proposed by Eugen-Olsen et al. [30].

Age demonstrated no statistically significant association with clinical outcomes in the normal–Gompertz survival model. The hazard ratio (HR) for seroconversion illness (SI switch) was 2.065 (95% CI: 0.164–3.967), and the mortality HR reached 3.434 (95% CI: 0.654–6.215), indicating insufficient evidence to guarantee their relationship.

5. Simulation Study

In this section, we investigate the performance of our proposed methods through simulations based on different copula functions. On the one hand, we focus on evaluating the estimation accuracy and efficiency of the proposed two-step process and comparing its performances with the two-stage method. On the other hand, we compare the parameter estimation precision achieved by the optimizer BOBYQA (TStage-BOBYQA; TStep-BOBYQA) with that of the popular optimizer L-BFGS-B (TStage-BFGS; TStep-BFGS). The comprehensive simulation results are presented in Table 7, Table 8, Table 9 and Table 10.

5.1. Design

The data were generated from various copula models with survival distributions. The main procedure follows the methodology outlined by Wei et al. [23]. We generated 3000 individuals, mimicking the standardized real data analysis from the ACS dataset. The covariate vector is Z = $[eqn]$ , with $[eqn]$ from Uniform(0, 1) and $[eqn]$ from Bernoulli(0.3). The parameters of two survival distributions were generated using the exponential regression model, as illustrated in Section 3.2. The times to non-terminal and terminal events were simulated from the exponential, Weibull, and Gompertz survival distributions. The censoring time $[eqn]$ was generated independently from a uniform distribution. The true regression coefficients are specified as in Table 7, Table 8, Table 9 and Table 10. We generated 500 datasets with 3000 individuals.

5.2. Performance Measures

To estimate the regression coefficients for the marginal and association parameters, we maximized the likelihood described in Section 3.2 by using the L-BFGS-B optimization algorithm (implemented in the optim function in R) and the BOBYQA algorithm (implemented in the optimx function in R). The performance of our proposed methods was measured using mean squared error (MSE) and coverage probability (CP).

5.3. Results

The results are presented in Table 7, Table 8, Table 9 and Table 10. As expected, in the two-stage optimization procedure, the estimation of parameters for the non-terminal regression model is inferior to that of the terminal regression model in terms of MSE and CP, as the non-terminal time to event is censored by the terminal time to event. Furthermore, we can see that the estimation of the parameters for the copula function is suboptimal due to inaccuracies in estimating the parameters for the non-terminal time-to-event regression model. This highlights the limitations of the two-stage optimization method when applied to semi-competing risks data.

When comparing the two-stage and two-step optimization procedures, the coverage probability achieved by the two-step method (TStep-BFGS and TStep-BOBYQA) is closer to the nominal level than that of the two-stage method (TStage-BFGS and TStage-BOBYQA) in most scenarios. Regarding the MSE criteria, the two-step method (TStep-BFGS and TStep-BOBYQA) demonstrates a slight reduction for the exponential distribution with all copula functions when compared to the two-stage method (TStage-BFGS and TStage-BOBYQA) (Table 7). Therefore, the two-step method outperformed the two-stage method.

When comparing the results of the optimization methods, TStage-BFGS yields results similar to those of TStage-BOBYQA with some exceptions. Specifically, for the Clayton–Gompertz and Frank–Gompertz models, TStage-BOBYQA demonstrates a reduction in MSE compared to TStage-BFGS (see Table 7 and Table 8). TStage-BFGS achieves 0.0% coverage probabilities for all coefficients, except for the parameter $[eqn]$ in the Frank–Weibull model (see Table 8). Furthermore, TStep-BFGS produces results comparable to those of TStep-BOBYQA. It is critical to note that the bounds of parameter values must be specified properly for TStep-BFGS; otherwise, parameter estimation cannot be obtained, rendering it impractical. The BOBYQA optimization method demonstrates greater reliability than the L-BFGS-B optimization method. Consequently, our proposed TStep-BOBYQA optimization method outperforms the other optimization methods in these simulations.

6. Discussion

In this article, we present copula-based regression models for semi-competing events with right censoring. We mainly consider four well-known univariate copula functions (i.e., Clayton, Frank, Gumbel, and normal) with three parametric distributions (i.e., exponential, Weibull, and Gompertz) as the marginal distributions. To evaluate the effects of covariates on the time to events, we incorporated these covariates into an interpretable regression model. Due to the complexity of copula-based likelihood functions and the censoring data, estimating the parameters of copula functions and marginal distributions is a great challenge. The initial values of parameters always have a heavy impact on the accuracy of parameter estimation. We employed a two-step procedure method to alleviate the influence and adopted the BOBYQA optimization solver.

In our analysis of real data on HIV/AIDS from the ACS, we employed copula-based models to construct the association between two events. We obtained consistent results, with the exception of the variable CCR5 and its effect on the SI switch event. Specifically, CCR5 shows no relation to the SI switch event in most models, except for two cases: the Frank–exponential model (HR:0.705, 95% CI: 0.425 to 0.986, using TStep-BFGS and TStep-BOBYQA; see Table 4) and the normal–Gompertz model (HR: 0.698, 95% CI: 0.402 to 0.995, using TStep-BFGS and TStep-BOBYQA; see Table 6), which show that CCR5 “WM” is associated with a lower risk of SI switch. Additionally, we conducted several simulations to evaluate the performance of different models based on MSE and CP criteria. In general, in view of MSE criteria, the two-step method performs slightly better than the two-stage method when utilizing both the L-BFGS-B and BOBYQA optimizers. From the aspect of CP, the two-step method consistently outperforms the two-stage method. For the BOBYQA optimizer, it yields nearly the same accurate results as the L-BFGS-B optimizer in certain simulations, including the Gumbel–exponential and Gumbel–Gompertz models in Table 9, as well as the normal–exponential model in Table 10. However, in terms of CP, the BOBYQA optimizer outperforms the L-BFGS-B optimizer in some other simulations. Specifically, the TStage-BFGS method yields a CP further from 95% for the variable coefficients of the terminal regression models in the Clayton–Weibull, Clayton–Gompertz (Table 7), and Frank–Weibull (Table 8) models. Furthermore, in some simulations, careful specification of the bounds is essential for the L-BFGS-B optimizer; otherwise, parameter estimation may not be achievable.

As it is known, one appealing feature of the copula model is that it allows for separate modeling and estimation of the association parameter and marginal distributions, such that the marginal distributions do not depend on the dependency structure. Consequently, the dependency and margins can be estimated separately. However, different distributional choices within the family may induce significantly different dependence structures, as studied by Shih and Louis [3] and Wang and Wells [31]. The methods for selecting appropriate copula models are critical to avoid incorrect association modeling, which can lead to misleading or incorrect conclusions. In addition to the AIC and BIC criteria, Zhu et al. [32] applied the diagnostic method proposed by Chen and Bandeen-Roche [33] to choose the appropriate copula model based on the observed bivariate failure times $[eqn]$ . Further work could see the use of asymmetric and heavy-tailed distributions for parametric copulas [34]. Therefore, research is warranted to develop a method for copula model selection.

In this article, we mainly consider the BOBYQA derivative-free optimization algorithm. More relevant algorithms surveyed by Rios and Sahinidis [35] need to be investigated, and they may yield better results. In practical applications, the selection of a copula model and the choice of optimizer must be considered for semi-competing risks data. At present, there is no established criterion for selecting an appropriate optimization algorithm for copula models across different semi-competing risks datasets. Optimization methodologies for semi-competing risks models remain a significant and understudied research area, warranting further systematic investigation.

Moreover, there may be issues related to the misspecification of the parametric marginal survival function. To enhance flexibility, semi-parametric or nonparametric models can be considered. Sun and Ding [36] proposed a class of copula-based semiparametric transformation models for bivariate data subject to general interval censoring. For the marginal model, they used the semi-parametric transformation model and approximate the infinite-dimensional nuisance parameters using sieves with Bernstein polynomials. The major challenges in analyzing semi-competing risks data are parameter identifiability, the stability of parameter estimation, and efficiency. Fine [37] used an iteration strategy between estimating parameters of a terminal event and estimating the parameters of copula and non-terminal events. Finally, in the context of real data, there are only two factors: age and CCR5 for AIDS data. Some other unmeasured covariates may influence the time-to-event, and the frailty term can be added for future work, as mentioned by Wei et al. [23].