Semiparametric copula method for semi-competing risks data subject to interval censoring and left truncation: Application to disability in elderly

Abstract

We aim to evaluate the marginal effects of covariates on time-to-disability in the elderly under the semi-competing risks framework, as death dependently censors disability, not vice versa. It becomes particularly challenging when time-to-disability is subject to interval censoring due to intermittent assessments. A left truncation issue arises when the age time scale is applied. We develop a flexible two-parameter copula-based semiparametric transformation model for semi-competing risks data subject to interval censoring and left truncation. The two-parameter copula quantifies both upper and lower tail dependence between two margins. The semiparametric transformation models incorporate proportional hazards and proportional odds models in both margins. We propose a two-step sieve maximum likelihood estimation procedure and study the sieve estimators’ asymptotic properties. Simulations show that the proposed method corrects biases in the marginal method. We demonstrate the proposed method in a large-scale Chinese Longitudinal Healthy Longevity Study and provide new insights into preventing disability in the elderly. The proposed method could be applied to the general semi-competing risks data with intermittently assessed disease status.

Keywords

Copula cumulative incidence function disability in elderly interval censoring left truncation semi-competing risks semiparametric transformation model

1 Introduction

Disability in the elderly is a long-term impairment largely attributable to disease conditions such as dementia, cardiovascular diseases, respiratory diseases, cancer, arthritis, and diabetes.^1–3According to World Health Organization, more than $46 %$ of elderly aged 60 years and over have disabilities, and more than $250$ million older people experience moderate to severe disability. Multiple large-scale longitudinal studies (i.e. the Health and Retirement Study,⁴ English Longitudinal Study of Ageing,⁵ and Chinese Longitudinal Healthy Longevity Study⁶) are investigating key topics relevant to the ageing process including disability. An important challenge for studying disability in the elderly is the high mortality. The time to disability is heavily censored by death. Because the time-to-death is possibly correlated with the time-to-disability, conventional survival methods, such as the Cox regression model for estimating covariate effects, are not applicable. Valid estimation and inference need to account for the dependence between disability and death.

Motivated by the Chinese Longitudinal Healthy Longevity Study (CLHLS), with 12,969 elderly being recruited in 2002, this work aims to quantify the dependence strength between disability and death, estimate covariate effects on disability, as well as to characterize the cumulative incidence function of disability. For example, although the elderly are at risk of developing disability and progressing to death,⁷ the explicit dependence strength between the time to disability and the time to death is unknown. Also, despite the wealth of evidence that sex plays an important role in the risk of disability,⁸ its precise effect on time-to-disability among the Chinese elderly is unavailable in the literature. Moreover, the cumulative incidence function is desired to depict the absolute risk profile of disability, which is more appropriate than the survival function for decision making.⁹

The disability data belong to the category of semi-competing risks data,¹⁰ where one individual can experience both the non-terminal event of interest (i.e. disease) and the terminal event (i.e. death). The non-terminal event is dependently censored by the terminal event, but not vice versa. Semi-competing risks data arise in a broad range of applications including cancer,¹¹ cardiovascular diseases,⁹ dementia,¹² and diabetes.¹³ The semi-competing risks data are usually handled with two approaches: the illness-death model^14–16 and the copula model.^17–19 The illness-death model focuses on the transition intensities between states (i.e. health, disease, and death). A frailty effect can be introduced into the illness-death model to account for the dependence among transitions.^20,12 The copula model constructs the joint distribution of disease and death by connecting marginal distributions via the copula function. The dependence strength can be explicitly estimated by the copula parameter. Because the marginal distribution in the copula model only involves time and covariates, the covariate effect on the disease can be interpreted on the marginal level. We choose the copula model in our motivating example because of its straightforward interpretation of dependence strength and covariate effects.

Modeling the semi-competing risks data of disability in CLHLS must address two important data attributes. First, the time scale by age is more preferred than by study entry in ageing studies like CLHLS.²¹ However, the age scale induces the left truncation issue since CLHLS only recruits individuals aged 60 years or above. Second, since the recruited individuals were intermittently assessed every 3–4 years, the time-to-disability is interval-censored. Ignoring these data attributes may lead to biased estimation results.²² Few works on disability have fully considered semi-competing (i.e. dependent-censoring), left-truncation, and interval-censoring. In this work, we aim to develop a copula model for semi-competing risks data subject to left truncation and interval censoring and apply it to the CLHLS data.

Despite the wealthy literature of copula models for semi-competing risks data subject to right censoring, there are few copula methods under interval censoring. Several works use copula to model bivariate interval-censored data.^23–25 However, these earlier works could not deal with interval-censored semi-competing risks data. Ma et al.²⁶ develop a copula model for dependent current status data where the non-terminal event is only examined once at either the terminal event time or the study stopping time. The method does not apply to the general interval censoring situation, where each individual can be potentially examined multiple times before the terminal event or the study stopping time.

We propose a broad class of copula-based semiparametric regression models. Specifically, we use a two-parameter copula function that can flexibly account for dependence between the non-terminal and terminal events in both upper and lower tails. The dependence strength is explicitly expressed by Kendall’s $τ$ . We use separate semiparametric transformation models for the non-terminal and terminal events, which incorporate the proportional hazards (PH) and proportional odds (PO) models. We approximate the infinite-dimensional parameters using sieves with Bernstein polynomials and develop a computationally efficient two-step sieve maximum likelihood estimation procedure. We establish the asymptotic properties for the sieve estimators. In addition, we estimate the cumulative incidence function for the non-terminal event based on the sieve estimators. We demonstrate the advantage of the proposed method through intensive simulation studies and the CLHLS data.

The remaining of this article is organized as follows. Section 2 introduces the model and the joint likelihood functions. Section 3 presents the sieve maximum likelihood estimation procedure, the estimation of the cumulative incidence function, and the asymptotic properties. Section 4 illustrates extensive simulation studies for the estimation performance of the proposed method. We analyze the CLHLS data and present the findings in Section 5. Finally, we discuss and conclude in Section 6. Additional results, the regularity conditions, and proofs are provided in the Supplemental Appendix.

2 Notation and likelihood

2.1 Copula model for semi-competing risks data

Assume that there are $n$ independent subjects in the study. For subject $i$ , denote $T_{i 1}$ as the non-terminal event time, $T_{i 2}$ the terminal event time, and $Z_{i j}$ a $p_{j}$ -dimensional covariate vector, $j = 1, 2$ . We define the marginal survival function for $T_{i j}$ as $S_{j} (t_{i j} | Z_{i j}) = Pr (T_{i j} > t_{i j} | Z_{i j})$ and the joint survival function for subject $i$ as $S (t_{i 1}, t_{i 2} | Z_{i 1}, Z_{i 2}) = Pr (T_{i 1} > t_{i 1}, T_{i 2} > t_{i 2} | Z_{i 1}, Z_{i 2})$ .

According to the Sklar’s theorem,²⁷ the joint survival function of $T_{i 1}$ and $T_{i 2}$ can be expressed by a unique copula function $C_{η}$ connecting two continuous marginal survival functions:

S (t_{i 1}, t_{i 2} | Z_{i 1}, Z_{i 2}) = C_{η} {S_{1} (t_{i 1} | Z_{i 1}), S_{2} (t_{i 2} | Z_{i 2})},

(1)

where the copula parameter

η

measures the dependence strength between the two marginal distributions and is explicitly connected to Kendall’s

τ

. As shown in equation (1), the marginal functions are modeled separately from the dependence structure, and therefore the covariate effects are naturally interpreted at the marginal level. One popular copula group for semi-competing risks data is the Archimedean family. Two commonly used Archimedean copulas are the Clayton and Gumbel copulas, which account for lower or upper tail dependence between two margins using one single copula parameter.²⁸ We consider a flexible two-parameter Archimedean copula model (“Copula2”)²⁹:

C_{α, κ} (u, v) = [1 + {(u^{- 1 / κ} - 1)^{1 / α} + (v^{- 1 / κ} - 1)^{1 / α}}^{α}]^{- κ}, α \in (0, 1], κ \in (0, \infty),

(2)

where

u

and

v

are uniformly distributed margins,

α

and

κ

are the two dependence parameters. One distinct feature of Copula2 is that it can simultaneously account for both upper and lower tail dependence of two margins. The dependence strength can be explicitly expressed via Kendall’s

τ = 1 - 2 α κ / (2 κ + 1)

. The lower and upper tail dependence can be measured by

τ_{L} = 2^{- α κ}

and

τ_{U} = 2 - 2^{α}

, respectively.³⁰ In particular, as seen in equation (2), Copula2 incorporates Clayton (

α = 1

) and Gumbel (

κ \to \infty

) copulas as special cases and thus offers a flexible framework for modeling semi-competing risks data.

2.2 Observed data and joint likelihood without left truncation

In many longitudinal studies like CLHLS, the non-terminal event is interval-censored due to intermittent assessments. Denote ${C_{i, k} : k = 0, 1, \dots, K_{i}, 0 = C_{i, 0} < C_{i, 1} < \dots < C_{i, K_{i}} < C_{A_{i}}}$ as assessment times, where $C_{A_{i}}$ is the administrative censoring time (e.g. end of study). We assume ${C_{i, k}}$ and $C_{A_{i}}$ are non-informative for $T_{i j}$ conditional on $Z_{i j}$ , which is reasonable in CLHLS as the assessment schedule does not depend on disability and mortality. The terminal event can be observed precisely in CLHLS and is subject to right censoring by $C_{A_{i}}$ . Thus, we observe $Y_{i} = min (T_{i 2}, C_{A_{i}})$ and $δ_{i} = I (T_{i 2} \leq C_{A_{i}})$ . The non-terminal event is observed as $(L_{i}, R_{i}]$ , which is the smallest interval containing $T_{i 1}$ . The construction of $(L_{i}, R_{i}]$ is complicated in CLHLS as it depends on $C_{i, k}$ , $Y_{i}$ and whether we know the status of the non-terminal event if the terminal event is observed. We will discuss the details below. Overall, we observe $D_{i} = {L_{i}, R_{i}, Y_{i}, δ_{i}, Z_{i j}}$ for subject $i$ .

Next, we will describe four censoring scenarios in CLHLS as illustrated in Figure 1 and construct the likelihood function under each scenario.

Scenario 1: $T_{i 1}$ and $T_{i 2}$ both right-censored

Figure 1.

Graphical representation of four censoring scenarios for disability ( $T_{1}$ ) and death ( $T_{2}$ ) in Chinese Longitudinal Healthy Longevity Study (CLHLS).

Scenario 1 contains two different sub-scenarios: both Scenarios 1a and 1b represent the situation where neither disability nor death occurs by the study end time $C_{A_{i}}$ , and thus both events are right censored by $C_{A_{i}}$ . The difference between Scenarios 1a and 1b lies in whether disability occurs before death (1a) or not (1b). Although the CLHLS data cannot distinguish the two sub-scenarios, we illustrate both in Figure 1 for completeness. The observed data for both Scenarios 1a and 1b are $(L_{i}, R_{i}] = (C_{A_{i}}, \infty], Y_{i} = C_{A_{i}}$ and $δ_{i} = 0$ . Therefore, Scenarios 1a and 1b share the same likelihood expressed as

Pr (T_{i 1} > L_{i}, T_{i 2} > Y_{i} | Z_{i 1}, Z_{i 2}) = C_{α, κ} {S_{1} (L_{i} | Z_{i 1}), S_{2} (Y_{i} | Z_{i 2})} .

Scenario 2: $T_{i 1}$ interval-censored and $T_{i 2}$ right-censored

Disability is observed between two adjacent assessment times $C_{i, k}$ and $C_{i, k + 1}$ , while death is right censored by the study end time $C_{A_{i}}$ . The observed data are $(L_{i}, R_{i}] = (C_{i, k}, C_{i, k + 1}], Y_{i} = C_{A_{i}}$ and $δ_{i} = 0$ . The contribution to the likelihood is

Pr (L_{i} < T_{i 1} \leq R_{i}, T_{i 2} > Y_{i} | Z_{i 1}, Z_{i 2}) = C_{α, κ} {S_{1} (L_{i} | Z_{i 1}), S_{2} (Y_{i} | Z_{i 2})} - C_{α, κ} {S_{1} (R_{i} | Z_{i 1}), S_{2} (Y_{i} | Z_{i 2})} .

Scenario 3: $T_{i 1}$ right-censored and $T_{i 2}$ exactly observed

Scenario 3 involves three different sub-scenarios. Scenarios 3a and 3b represent situations where disability is not observed at the assessment time $C_{i, k}$ , but death occurs before the next assessment $C_{i, k + 1}$ . Therefore, disability is right-censored by $C_{i, k}$ , and death time is observed at $T_{i 2}$ . The difference betweenScenarios 3a and 3b lies in whether disability occurs before death (3a) or not (3b). While the CLHLS data cannot distinguish the two sub-scenarios, we illustrate both in Figure 1 for completeness. The observed data for both Scenarios 3a and 3b are $(L_{i}, R_{i}] = (C_{i, k}, \infty], Y_{i} = T_{i 2}$ and $δ_{i} = 1$ . Scenario 3c is similar to 3b, except that we know disability does not occur by the time of death, leading to the right censoring of disability by the death time $T_{i 2}$ . For example, for some elderly who passed away, their disability statuses before death can be determined through carefully designed interviews with their descendants. If disability did not happen by the time of death based on the interview results, we observed scenario 3c with $(L_{i}, R_{i}] = (T_{i 2}, \infty], Y_{i} = T_{i 2}$ and $δ_{i} = 1$ . In summary, the general likelihood for all Scenarios 3a, 3b, and 3c can all be expressed as

Pr (T_{i 1} > L_{i}, T_{i 2} = Y_{i} | Z_{i 1}, Z_{i 2}) = - {\dot{S}}_{2} (Y_{i} | Z_{i 2}) \times C_{α, κ}^{(0, 1)} {S_{1} (L_{i} | Z_{i 1}), S_{2} (Y_{i} | Z_{i 2})},

where

{\dot{S}}_{2}

is the first-order derivative of

S_{2}

and

C_{α, κ}^{(0, 1)} (u, v) = \partial C_{α, κ} (u, v) / \partial v

Scenario 4: $T_{i 1}$ interval-censored and $T_{i 2}$ exactly observed

Scenario 4 includes two different sub-scenarios. Scenario 4a represents the situation where disability is observed between two adjacent assessment times $C_{i, k}$ and $C_{i, k + 1}$ , and death is observed after $C_{i, k + 1}$ . The observed data for Scenario 4a are $(L_{i}, R_{i}] = (C_{i, k}, C_{i, k + 1}], Y_{i} = T_{i 2}$ and $δ_{i} = 1$ . Scenario 4b represents the situation where disability is not observed at time $C_{i, k}$ and death happens before the next visit $C_{i, k + 1}$ , but we know disability occurs between $C_{i, k}$ and death through interviewing the descendants of the elderly. The observed data for Scenario 4b are $(L_{i}, R_{i}] = (C_{i, k}, T_{i 2}], Y_{i} = T_{i 2}$ and $δ_{i} = 1$ . In summary, the general likelihood for both Scenarios 4a and 4b can be expressed as

\begin{aligned} Pr (L_{i} < T_{i 1} \leq R_{i}, T_{i 2} = Y_{i} | Z_{i 1}, Z_{i 2}) = & - {\dot{S}}_{2} (Y_{i} | Z_{i 2}) \times [C_{α, κ}^{(0, 1)} {S_{1} (L_{i} | Z_{i 1}), S_{2} (Y_{i} | Z_{i 2})} \\ - C_{α, κ}^{(0, 1)} {S_{1} (R_{i} | Z_{i 1}), S_{2} (Y_{i} | Z_{i 2})}] . \end{aligned}

Overall, the joint likelihood of combining all scenarios can be summarized as

\begin{aligned} L_{n} (S_{1}, S_{2}, α, κ | D) & = \prod_{i = 1}^{n} {- {\dot{S}}_{2} (Y_{i} | Z_{i 2})}^{δ_{i}} \times [C_{α, κ}^{(0, δ_{i})} {S_{1} (L_{i} | Z_{i 1}), S_{2} (Y_{i} | Z_{i 2})} \\ - C_{α, κ}^{(0, δ_{i})} {S_{1} (R_{i} | Z_{i 1}), S_{2} (Y_{i} | Z_{i 2})}], \end{aligned}

(3)

where

C_{α, κ}^{(0, δ_{i})} (u, v) = \partial C_{α, κ}^{(0 + δ_{i})} (u, v) / \partial v^{δ_{i}}

. The fact that

T_{1}

is a non-terminal event and

T_{2}

is a terminal event is indeed reflected in the observed data and the likelihood function in formula (3). When

T_{i 1}

is right censored with

R_{i} = \infty

, the term

C_{α, κ}^{(0, δ_{i})}

involving

R_{i}

becomes

0

. In particular, the likelihood in (3) incorporates the special case of dependent current status data with only one assessment time (i.e. either

C_{A_{i}}

T_{i 2}

) for the non-terminal event, where we observe either

L_{i} = 0

R_{i} = \infty

for all subjects. Next, we will extend the likelihood function (3) to handle the left truncation issue in CLHLS.

2.3 Observed data and joint likelihood with left truncation

Suppose a subject $i$ without disability is recruited at time $A_{i}$ , and then intermittently assessed by ${C_{i, k} : k = 1, \dots, K_{i}, A_{i} < C_{i, 1} < \dots < C_{i, K_{i}}}$ till the terminal event time $T_{i 2}$ or administrative stopping time $C_{A_{i}}$ . Similar to Section 2.2., we assume $A_{i}$ , ${C_{i, k}}$ and $C_{A_{i}}$ are independent of $T_{i j}$ conditional on $Z_{i j}$ . The censoring scenarios are the same as the scenarios depicted in Figure 1. The observed data become $D_{i} = {A_{i}, L_{i}, R_{i}, Y_{i}, δ_{i}, Z_{i j}}$ . The likelihood adjusting for left truncation is expressed as

L_{n} (S_{1}, S_{2}, α, κ | D) = \prod_{i = 1}^{n} {- {\dot{S}}_{2} (Y_{i} | Z_{i 2})}^{δ_{i}} \times \frac{C_{α, κ}^{(0, δ_{i})} {S_{1} (L_{i} | Z_{i 1}), S_{2} (Y_{i} | Z_{i 2})} - C_{α, κ}^{(0, δ_{i})} {S_{1} (R_{i} | Z_{i 1}), S_{2} (Y_{i} | Z_{i 2})}}{C_{α, κ} {S_{1} (A_{i} | Z_{i 1}), S_{2} (A_{i} | Z_{i 2})}} .

(4)

The likelihood in (4) incorporates the function in (3) when

A_{i} = 0

(i.e. no left truncation). Moreover, Figure 1 covers all possible censoring scenarios when the disease is intermittently assessed and dependently censored by death. Therefore, the likelihood function in (4) is applicable to the general interval-censored and possibly left-truncated semi-competing risks data in a broad range of applications, such as cancer, dementia, and diabetes.

2.4 Semiparametric transformation model for marginal distributions

We consider the semiparametric transformation model for the marginal survival function of $T_{i j}$ :

S_{j} (t_{i j} | Z_{i j}) = \exp [- G_{j} {e^{Z_{i j}^{T} β_{j}} Λ_{j} (t_{i j})}],

(5)

where

G_{j} (\cdot)

is a pre-specified increasing function,

β_{j}

is an unknown regression coefficient vector, and

Λ_{j} (\cdot)

is an unknown non-decreasing function. The model in (5) incorporates a broad class of survival models. In particular, model (5) becomes a PH or PO model when

G_{j} (x) = x

\log (1 + x)

. We choose the Box-Cox transformation

G_{j} (x) = {(1 + x)^{r_{j}} - 1} / r_{j}

, which includes PH (

r_{j} = 1

) and PO (

r_{j} = 0

). For the derivative of

S_{j} (t_{i j} | Z_{i j})

, we have

{\dot{S}}_{j} (t_{i j} | Z_{i j}) = - \exp [- G_{j} {e^{Z_{i j}^{T} β_{j}} Λ_{j} (t_{i j})}] \times {\dot{G}}_{j} {e^{Z_{i j}^{T} β_{j}} Λ_{j} (t_{i j})} \times e^{Z_{i j}^{T} β_{j}} \times λ_{j} (t_{i j})

, where

{\dot{G}}_{j}

is the derivative of

G_{j}

and

\int_{0}^{t_{i j}} λ_{j} (s) d s = Λ_{j} (t_{i j})

3 Estimation and inference

3.1 Sieve likelihood with Bernstein polynomials

In the likelihood function (4), we are interested in estimating the unknown parameter $θ \in Θ$ :

Θ = {θ = (β_{1}^{T}, β_{2}^{T}, α, κ, Λ_{1}, Λ_{2})^{T} \in B \otimes M^{1} \otimes M^{2}} .

Here

B = {(β = (β_{1}^{T}, β_{2}^{T})^{T}, α, κ) \in R^{p} \times R^{(0, 1]} \times R^{+}, ‖ β ‖ + ‖ α ‖ + ‖ κ ‖ \leq M}

with

p

being the dimension of

β

and

M

being a positive constant. Denote

M^{j}

as the collection of all bounded, continuous, nondecreasing and nonnegative functions over

[c, u]

, where

0 \leq c < u < \infty

. In practice,

[c, u]

can be chosen as the range of all observed finite times (i.e.

L_{i}

R_{i}

Y_{i}, i = 1, \dots, n

In the log-likelihood function $l_{n} (θ; D) = \log L_{n} (θ; D) = \sum_{i = 1}^{n} \log L (θ; D_{i}) = \sum_{i = 1}^{n} l (θ; D_{i})$ , there are finite-dimensional parameters of interest $(β, α, κ)$ and two infinite-dimensional nuisance parameters $(Λ_{1}, Λ_{2})$ , which need to be estimated simultaneously. Unlike the right-censored data, tools such as partial likelihood and martingale cannot be applied to the interval-censored data due to the absence of exact event times. Instead, we employ the sieve approach and form a sieve likelihood. The sieve approach is one of the most well-developed methods for estimating infinite-dimensional parameters.^31–33 Specifically, we use Bernstein polynomials to build a sieve space $Θ_{n} = {θ_{n} = (β^{T}, α, κ, Λ_{1 n}, Λ_{2 n})^{T} \in B \otimes M_{n}^{1} \otimes M_{n}^{2}}$ . Here, $M_{n}^{j}$ is the space defined by Bernstein polynomials:

M_{n}^{j} = {Λ_{j n} (t) = \sum_{k = 0}^{m_{j n}} ϕ_{j k} B_{j, k} (t, m_{j n}, c, u) : \sum_{k = 0}^{m_{j n}} | ϕ_{j k} | \leq M_{j n}; 0 \leq ϕ_{j 0} \leq \cdot \cdot \cdot \leq ϕ_{j m_{j n}}; j = 1, 2},

where

B_{j, k} (t, m_{j n}, c, u)

represents the Bernstein basis polynomial defined as:

B_{j, k} (t, m_{j n}, c, u) = (_{k}^{m_{j n}}) {(\frac{t - c}{u - c})}^{k} {(1 - \frac{t - c}{u - c})}^{m_{j n} - k}; k = 0, \dots, m_{j n},

with degree

m_{j n} = o (n^{ν_{j}})

for some

ν_{j} \in (0, 1)

In practice, we choose $m_{j n}$ with the smallest AIC (i.e. Akaike information criterion) value. With pre-specified $m_{j n}$ , we solve $ϕ_{j k}$ together with other parameters $(β, α, κ)$ . One advantage of Bernstein polynomials is that they can naturally achieve the non-negativity and monotonicity properties of $Λ_{j} (t)$ through re-parameterization.³² Another advantage of Bernstein polynomials is that they do not require the specification of interior knots, making them flexible for use. With the sieve space defined above, $Λ_{j} (t)$ will be approximated by $Λ_{j n} (t) \in M_{n}^{j}$ and $λ_{j} (t)$ can be estimated by $λ_{j n} (t) = \sum_{k = 0}^{m_{j n}} ϕ_{j k} {\dot{B}}_{j, k} (t, m_{j n}, c, u)$ with ${\dot{B}}_{j, k} (t, m_{j n}, c, u)$ being the first order derivative of $B_{j, k} (t, m_{j n}, c, u)$ .

In the next section, we describe an estimation procedure to maximize $l_{n} (θ; D)$ over the sieve space $Θ_{n}$ for obtaining the sieve maximum likelihood estimators ${\hat{θ}}_{n} = ({\hat{β}}_{n}^{T}, {\hat{α}}_{n}, {\hat{κ}}_{n}, {\hat{Λ}}_{1 n}, {\hat{Λ}}_{2 n})^{T}$ .

3.2 Estimation procedure for sieve maximum likelihood estimators ${\hat{θ}}_{n}$

We develop a two-step sieve maximum likelihood estimation procedure that ensures computational efficiency in semi-competing risks data. The procedure generally applies to any choice of Archimedean copulas and marginal models. In principle, one may obtain the sieve estimators by directly maximizing the log-likelihood function in one step. However, due to the complex structure of the likelihood function, directly maximizing the log-likelihood in a one-step procedure using randomly selected initial values will cause a high convergence failure rate. Therefore, by following previous works on copula model estimation,^34,24 we propose a two-step estimation procedure by introducing a separate first step. The goal of the first step is to obtain good initial values for the optimization of the complex likelihood in the second step, which will greatly guarantee the convergence performance of the sieve estimation. All the parameters are updated through the full likelihood in the second step. Therefore, the coverage probabilities and variance estimates are all correct, as illustrated in simulations in Section 4.

Specifically, we first estimate $(β_{2}, Λ_{2 n})$ marginally in step 1(a). Their consistent estimators are then plugged into the joint likelihood to form a pseudo-likelihood. In step 1(b), the rest parameters $(β_{1}, α, κ, Λ_{1 n})$ are estimated through maximizing the pseudo-likelihood function. Finally, using initial values from steps 1(a) and 1(b), we simultaneously update all the unknown parameters under the joint log-likelihood function in step 2. The reason for estimating $(β_{1}, Λ_{1 n})$ in step 1(b), not step 1(a), is that the non-terminal event is dependently censored by the terminal event, leading to biased estimation in the marginal model in step 1(a). The estimation procedure is described below:

Obtain initial estimates of $θ_{n}$ :

$({\hat{β}}_{2 n}^{(1)}, {\hat{Λ}}_{2 n}^{(1)}) = {\arg \max}_{(β_{2}, Λ_{2 n})} l_{2 n} (β_{2}, Λ_{2 n})$ , where $l_{2 n}$ denotes the marginal sieve log-likelihood accounting for left truncation and right censoring in the terminal event $T_{i 2}$ ;

$({\hat{β}}_{1 n}^{(1)}, {\hat{α}}_{n}^{(1)}, {\hat{κ}}_{n}^{(1)}, {\hat{Λ}}_{1 n}^{(1)}) = {\arg \max}_{(β_{1}, α, κ, Λ_{1 n})} l_{n} (β_{1}, {\hat{β}}_{2 n}^{(1)}, α, κ, Λ_{1 n}, {\hat{Λ}}_{2 n}^{(1)})$ , where ${\hat{β}}_{2 n}^{(1)}$ and ${\hat{Λ}}_{2 n}^{(1)}$ are the initial estimates from (a), and $l_{n}$ is the joint sieve log-likelihood.

Simultaneously maximize the joint sieve log-likelihood to get final estimates:

${\hat{θ}}_{n} = ({\hat{β}}_{n}, {\hat{α}}_{n}, {\hat{κ}}_{n}, {\hat{Λ}}_{1 n}, {\hat{Λ}}_{2 n}) = {\arg \max}_{(β, α, κ, Λ_{1 n}, Λ_{2 n})} l_{n} (β, α, κ, Λ_{1 n}, Λ_{2 n})$ , with initial values $({\hat{β}}_{n}^{(1)} = ({\hat{β}}_{1 n}^{(1)}, {\hat{β}}_{2 n}^{(1)}), {\hat{α}}_{n}^{(1)}, {\hat{κ}}_{n}^{(1)}, {\hat{Λ}}_{1 n}^{(1)}, {\hat{Λ}}_{2 n}^{(1)})$ obtained from step 1.

The two-step estimation procedure guarantees almost

100 %

convergence in our simulations. We obtain the estimates for the parameters of interest (

β, α, κ

) from step 2. The initial values from step 1 are already consistent estimators.^35,36 Step 2 yields correct variance-covariance estimates for all the parameters using the joint likelihood. Standard optimization algorithms, such as Newton-Raphson or Broyden–Fletcher–Goldfarb–Shanno (BFGS), can be employed to obtain the estimators and observed information matrix. For the variance estimation, we invert the observed information matrix of all parameters (including the nuisance parameters).³⁷ An alternative variance estimation approach is bootstrap. However, bootstrap is quite time-consuming, especially with a large sample size (e.g. n = 12682 in CLHLS). Therefore, we estimate variances through the observed Fisher information, which produces correct coverage probabilities as shown in simulations in Section 4. Due to the complex structure of the joint sieve log-likelihood, we use Richardson’s extrapolation³⁸ to approximate the score function and observed Fisher information numerically, which works well in the simulated and real-data analysis.

3.3 Estimation of cumulative incidence function

In the semi-competing risks data analysis, instead of using the marginal survival function, it is more reasonable to use the cumulative incidence function to characterize the risk profile of the non-terminal event of interest.^9,13 Given that a subject is alive and free of disability by age $a$ , the cumulative incidence function for disability is expressed as

\begin{aligned} Pr (T_{1} \leq t, T_{1} \leq T_{2} | T_{1} > a, T_{2} > a, Z_{1}, Z_{2}) & = \frac{Pr (a < T_{1} \leq t, T_{1} \leq T_{2} | Z_{1}, Z_{2})}{Pr (T_{1} > a, T_{2} > a | Z_{1}, Z_{2})} \\ = \frac{Pr (a < T_{1} \leq T_{2} < t | Z_{1}, Z_{2}) + Pr (a < T_{1} \leq t \leq T_{2} | Z_{1}, Z_{2})}{Pr (T_{1} > a, T_{2} > a | Z_{1}, Z_{2})} \\ = \frac{\int_{a}^{t} [C_{α, κ}^{(1, 0)} {S_{1} (t_{1} | Z_{1}), S_{2} (t | Z_{2})} - C_{α, κ}^{(1, 0)} {S_{1} (t_{1} | Z_{1}), S_{2} (t_{1} | Z_{2})}] \times {\dot{S}}_{1} (t_{1} | Z_{1}) d t_{1}}{C_{α, κ} {S_{1} (a | Z_{1}), S_{2} (a | Z_{2})}} \\ + \frac{C_{α, κ} {S_{1} (a | Z_{1}), S_{2} (t | Z_{2})} - C_{α, κ} {S_{1} (t | Z_{1}), S_{2} (t | Z_{2})}}{C_{α, κ} {S_{1} (a | Z_{1}), S_{2} (a | Z_{2})}} . \end{aligned}

(6)

As shown in formula (6), the cumulative incidence function measures the absolute cumulative risks of the non-terminal event

T_{1}

over time. It can be used to evaluate the cumulative impact of a covariate on the population as a whole. The cumulative incidence function is analogous to the distribution function with the range from 0 to 1, thus making it easy to interpret in practice. The cumulative incidence function provides natural interpretations of results from a clinical or epidemiological point of view. Therefore, it is a popular effect measure in (semi-)competing risks data analysis.^{39,40,9,41,16} We estimate the cumulative incidence function using the sieve maximum likelihood estimators

{\hat{θ}}_{n}

. The function is evaluated on

t \in [c, u]

. For the integral in the formula, we approximate it numerically using Gauss–Legendre quadratures.

3.4 Asymptotic properties of sieve estimators

We study the asymptotic properties of the sieve maximum likelihood estimators ${\hat{θ}}_{n}$ . Let $θ_{0} = (β_{0}^{T}, α_{0}, κ_{0}, Λ_{10}, Λ_{20})^{T}$ denote the true value of $θ \in Θ$ . Denote $P$ as the true probability measure and $P_{n}$ as the empirical measure. Let $| v |$ be the Euclidean norm for a vector $v$ , $‖ f ‖_{\infty} = s u p_{t} | f (t) |$ be the supremum norm for a function $f (t)$ , and $‖ f ‖_{2} = (\int | f |^{2} d P)^{1 / 2}$ be the $L_{2} (P)$ norm for a function $f$ under the probability measure $P$ . Define the distance between $θ^{(1)} = (β^{(1) T}, α^{(1)}, κ^{(1)}, Λ_{1}^{(1)}, Λ_{2}^{(1)})^{T} \in Θ$ and $θ^{(2)} = (β^{(2) T}, α^{(2)}, κ^{(2)}, Λ_{1}^{(2)}, Λ_{2}^{(2)})^{T} \in Θ$ as

d (θ^{(1)}, θ^{(2)}) = (| β^{(1)} - β^{(2)} |^{2} + | α^{(1)} - α^{(2)} |^{2} + | κ^{(1)} - κ^{(2)} |^{2} + ‖ Λ_{1}^{(1)} - Λ_{1}^{(2)} ‖_{2}^{2} + ‖ Λ_{2}^{(1)} - Λ_{2}^{(2)} ‖_{2}^{2})^{1 / 2} .

The following theorems present the asymptotic results. The proofs are supplied in the Supplemental Appendix.

Theorem 3.1
(Convergence rate). Assume that Conditions 1-2 and 4-5 in the Supplemental Appendix hold. Let $m_{n} = o (n^{ν})$ , where $ν \in (0, 1)$ and $q$ be the smoothness parameters of $Λ_{j}$ as defined in Condition 4. Then, we have
$d ({\hat{θ}}_{n}, θ_{0}) = O_{p} (n^{- \min {q ν / 2, (1 - ν) / 2}}) .$

Theorem 3.2
(Asymptotic normality and efficiency). Suppose Conditions 1–5 in the Supplemental Appendix hold. Define ${\hat{b}}_{n} = ({\hat{β}}_{n}^{T}, {\hat{α}}_{n}, {\hat{κ}}_{n})^{T}$ and $b_{0} = (β_{0}^{T}, α_{0}, κ_{0})^{T}$ . If $1 / (2 + q) < ν < 1 / 2$ , then
$n^{1 / 2} ({\hat{b}}_{n} - b_{0}) = I^{- 1} (b_{0}) n^{1 / 2} P_{n} l^{} (b_{0}, Λ_{10}, Λ_{20}; D) + o_{p} (1) \to_{d} N {0, I^{- 1} (b_{0})},$
where $I (b_{0}) = P l^{} (b_{0}, Λ_{10}, Λ_{20}; D)^{\otimes 2}$ and $l^{*} (b_{0}, Λ_{10}, Λ_{20}; D)$ is the efficient score function with respect to the log-likelihood function under equation (4) as defined in the Supplemental Appendix. Therefore, ${\hat{b}}_{n}$ is asymptotically normal and efficient.

Theorem 3.1 states that the sieve estimator ${\hat{θ}}_{n}$ has a polynomial convergence rate. Although this rate is lower than $n^{- 1 / 2}$ , Theorem 3.2 states that the estimators of the finite-dimensional parameters ( $β, α, κ$ ) are asymptotically normal and semiparametrically efficient.
4 Simulation studies

We conduct simulations to evaluate the performance of the proposed method, including parameter estimation and cumulative incidence function estimation.

4.1 Data generation

The true event times are generated from various copula models (i.e. Clayton, Gumbel, AMH, Frank, and Gaussian) with Weibull (PH model, scale $λ = 0.1$ and shape $k = 2$ ) or Loglogistic (PO model, scale $λ = 1$ and shape $k = 2$ ) margins. For the dependence strength between margins, we set Kendall’s $τ = 0.2$ or $0.6$ , corresponding to weak or strong dependence. The covariate vector $Z$ includes three covariates: one from a Normal distribution $N (6, 2)$ and two from Bernoulli distributions with $p = 0.5$ . Two margins share the same $Z$ . The regression coefficients are set as $β_{1}^{T} = (β_{11}, β_{12}, β_{13})^{T} = (0.1, 0.1, 0), β_{2}^{T} = (β_{21}, β_{22}, β_{23})^{T} = (0.1, 0.1, 0)$ . The parameters (i.e. shape, scale, and regression coefficients) are selected to mimic the real-data analysis. The truncation times $A$ are obtained from uniform $(0, a)$ . We generate the intermittent examination times by following the procedure in Sun and Ding,²⁴ assuming the length between two adjacent times follows an exponential distribution. The administrative censoring times $C_{A}$ are generated from uniform distribution. Overall, the right-censoring rates for $T_{1}$ and $T_{2}$ are $60 %$ and $25 %$ , respectively. The left-truncation rates are $20 %$ for both events. Under all scenarios, we use $n = 500, m_{n 1} = m_{n 2} = 3$ , and perform $1, 000$ replications.

4.2 Parameter estimation under correct copula specification

We evaluate the estimation performance of the proposed method. The true model is Clayton copula with Weibull (PH) or Loglogistic (PO) margins under Kendall’s $τ = 0.6$ . We fit three models: the true parametric model, the proposed model, and a marginal semiparametric model for the non-terminal event accounting for interval censoring and left truncation. Specifically, the marginal semiparametric model is a univariate model for the non-terminal event (i.e. $T_{1}$ ), which does not account for the dependent censoring by the terminal event. The proposed method considers the joint distribution of the non-terminal and terminal events. The true parametric model uses the true copula and parametric margins, which serves as the gold standard model. For the dependence estimation, we report Kendall’s $τ$ estimates for both Copula2 and true models to make direct comparisons between the two models. Table 1 summarizes the estimation results of regression coefficients and Kendall’s $τ$ . Similar to the estimators from the true model, the proposed sieve maximum likelihood estimators are all virtually unbiased, and all empirical coverage probabilities are close to the nominal level. The marginal model does not account for dependent censoring and thus produces biased estimators and wrong coverage probabilities. Moreover, the marginal method could not estimate the dependence strength between the two events. We also consider the situation of weaker dependence by setting Kendall’s $τ = 0.2$ . As shown in Table S1 in the Supplemental Appendix, the proposed method performs well, whereas the marginal method produced biased estimation results. We also conduct simulations with a smaller sample size ( $n = 300$ ). As shown in Table S2 in the Supplemental Appendix, we observe similar findings as simulation results with $n = 500$ . We observe two convergence failures out of 1000 replications (failure rate $0.2 %$ ), probably due to optimization’s numerical instability with a small sample of $n = 300$ .

Table 1.
Estimation results for interval-censored and left-truncated semi-competing risks data under three models: the true Clayton copula with parametric margins (Weibull for proportional hazards and Loglogistic for proportional odds; denoted as Clayton-PM), the proposed Copula2 with semiparametric margins (Copula2-S), and the marginal semiparametric model (Marginal-S) for the non-terminal event. Kendall’s $τ$ is set as 0.6.

Clayton-PM Copula2-S Marginal-S

Param Bias SE SEE (CP) Bias SE SEE (CP) Bias SE SEE (CP)

Proportional hazards

$β_{11}$ 0.0009 0.0275 0.0272 (0.947) $-$ 0.0014 0.0278 0.0270 (0.946) $-$ 0.0199 0.0412 0.0382 (0.890)

$β_{12}$ $-$ 0.0021 0.1070 0.1060 (0.949) $-$ 0.0003 0.1070 0.1060 (0.943) $-$ 0.0177 0.1590 0.1540 (0.928)

$β_{13}$ 0.0019 0.0782 0.0768 (0.939) 0.0011 0.0753 0.0769 (0.955) $-$ 0.0047 0.1140 0.1110 (0.941)

$β_{21}$ 0.0009 0.0237 0.0236 (0.949) $-$ 0.0001 0.0237 0.0234 (0.950) – – –

$β_{22}$ $-$ 0.0006 0.0945 0.0919 (0.943) 0.0001 0.0930 0.0921 (0.950) – – –

$β_{23}$ 0.0014 0.0679 0.0664 (0.939) 0.0018 0.0649 0.0665 (0.958) – – –

$τ$ 0.0054 0.0258 0.0275 (0.955) 0.0091 0.0280 0.0278 (0.935) – – –

Proportional odds

$β_{11}$ $-$ 0.0028 0.0504 0.0507 (0.954) $-$ 0.0028 0.0502 0.0499 (0.948) $-$ 0.0243 0.0638 0.0547 (0.877)

$β_{12}$ $-$ 0.0087 0.2010 0.2000 (0.949) $-$ 0.0063 0.2020 0.1990 (0.944) $-$ 0.0438 0.2520 0.2130 (0.872)

$β_{13}$ 0.0023 0.1460 0.1450 (0.946) 0.0001 0.1440 0.1440 (0.949) 0.0063 0.1820 0.1570 (0.906)

$β_{21}$ $-$ 0.0008 0.0449 0.0452 (0.962) $-$ 0.0016 0.0450 0.0448 (0.955) – – –

$β_{22}$ $-$ 0.0113 0.1810 0.1790 (0.949) $-$ 0.0077 0.1800 0.1780 (0.945) – – –

$β_{23}$ 0.0031 0.1310 0.1290 (0.949) $-$ 0.0027 0.1310 0.1290 (0.949) – – –

$τ$ 0.0028 0.0249 0.0264 (0.963) 0.0123 0.0255 0.0260 (0.945) – – –

	Clayton-PM	Copula2-S	Marginal-S
Proportional hazards
$β_{11}$	0.0009	0.0275	0.0272 (0.947)	$-$ 0.0014	0.0278	0.0270 (0.946)	$-$ 0.0199	0.0412	0.0382 (0.890)
$β_{12}$	$-$ 0.0021	0.1070	0.1060 (0.949)	$-$ 0.0003	0.1070	0.1060 (0.943)	$-$ 0.0177	0.1590	0.1540 (0.928)
$β_{13}$	0.0019	0.0782	0.0768 (0.939)	0.0011	0.0753	0.0769 (0.955)	$-$ 0.0047	0.1140	0.1110 (0.941)
$β_{21}$	0.0009	0.0237	0.0236 (0.949)	$-$ 0.0001	0.0237	0.0234 (0.950)	–	–	–
$β_{22}$	$-$ 0.0006	0.0945	0.0919 (0.943)	0.0001	0.0930	0.0921 (0.950)	–	–	–
$β_{23}$	0.0014	0.0679	0.0664 (0.939)	0.0018	0.0649	0.0665 (0.958)	–	–	–
$τ$	0.0054	0.0258	0.0275 (0.955)	0.0091	0.0280	0.0278 (0.935)	–	–	–
Proportional odds
$β_{11}$	$-$ 0.0028	0.0504	0.0507 (0.954)	$-$ 0.0028	0.0502	0.0499 (0.948)	$-$ 0.0243	0.0638	0.0547 (0.877)
$β_{12}$	$-$ 0.0087	0.2010	0.2000 (0.949)	$-$ 0.0063	0.2020	0.1990 (0.944)	$-$ 0.0438	0.2520	0.2130 (0.872)
$β_{13}$	0.0023	0.1460	0.1450 (0.946)	0.0001	0.1440	0.1440 (0.949)	0.0063	0.1820	0.1570 (0.906)
$β_{21}$	$-$ 0.0008	0.0449	0.0452 (0.962)	$-$ 0.0016	0.0450	0.0448 (0.955)	–	–	–
$β_{22}$	$-$ 0.0113	0.1810	0.1790 (0.949)	$-$ 0.0077	0.1800	0.1780 (0.945)	–	–	–
$β_{23}$	0.0031	0.1310	0.1290 (0.949)	$-$ 0.0027	0.1310	0.1290 (0.949)	–	–	–
$τ$	0.0028	0.0249	0.0264 (0.963)	0.0123	0.0255	0.0260 (0.945)	–	–	–

SE: standard deviation of the point estimate; SEE: mean of the standard error estimates; CP: $95 %$ coverage probability.

Moreover, we evaluate the three models’ performance in two special cases: (1) no left truncation; (2) dependent current status data where the non-terminal event is only censored by the study stopping time or the terminal event time. The results are summarized in Tables S3 and S4 in the Supplemental Appendix. Overall, the performance of our proposed and the true models are close, whereas the marginal model leads to biased estimation and incorrect coverage probabilities.

4.3 Parameter estimation under copula misspecification

We evaluate the estimation performance of the proposed model on data generated from copula models that do not belong to the Copula2 family, such as AMH copula with $τ = 0.2$ ( $τ$ is always $< 1 / 3$ for AMH), Frank copula with $τ = 0.6$ , and Gaussian copula with $τ = 0.6$ . In particular, AMH and Frank copulas belong to the Archimedean family, whereas the Gaussian copula does not. The margins follow Loglogistic distributions (PO model). In Table 2, the regression coefficient estimates using Copula2 are unbiased with coverage probabilities close to 95%. The estimation biases for Kendall’s $τ$ are minimal with good coverage probabilities except in the scenario where data are generated from Frank (coverage probability = 89%). Overall, the two-parameter copula model demonstrates good robustness to misspecification in copula functions.

Table 2.
Estimation results for interval-censored and left-truncated semi-competing risks data fitted with the proposed Copula2-S model (misspecified copula) where the true data are generated from AMH, Frank, and Gaussian copulas.

AMH Frank Gaussian

Param Bias SE SEE (CP) Bias SE SEE (CP) Bias SE SEE (CP)

$β_{11}$ 0.0053 0.0363 0.0351 (0.937) 0.0084 0.0309 0.0300 (0.928) 0.0022 0.0306 0.0293 (0.952)

$β_{12}$ 0.0071 0.1360 0.1370 (0.945) 0.0057 0.1200 0.1160 (0.932) 0.0008 0.1180 0.1160 (0.947)

$β_{13}$ 0.0022 0.1030 0.0989 (0.939) −0.0005 0.0832 0.0831 (0.951) 0.0013 0.0848 0.0837 (0.944)

$β_{21}$ 0.0008 0.0269 0.0264 (0.946) 0.0044 0.0270 0.0266 (0.933) 0.0006 0.0261 0.0256 (0.948)

$β_{22}$ 0.0002 0.1030 0.1030 (0.958) 0.0040 0.1060 0.1030 (0.936) −0.0016 0.1020 0.1010 (0.948)

$β_{23}$ 0.0021 0.0750 0.0742 (0.952) 0.0005 0.0752 0.0738 (0.940) −0.0006 0.0762 0.0729 (0.942)

$τ$ $-$ 0.0023 0.0531 0.0550 (0.940) −0.0231 0.0458 0.0346 (0.894) −0.0140 0.0374 0.0349 (0.943)

	AMH	Frank	Gaussian
$β_{11}$	0.0053	0.0363	0.0351 (0.937)	0.0084	0.0309	0.0300 (0.928)	0.0022	0.0306	0.0293 (0.952)
$β_{12}$	0.0071	0.1360	0.1370 (0.945)	0.0057	0.1200	0.1160 (0.932)	0.0008	0.1180	0.1160 (0.947)
$β_{13}$	0.0022	0.1030	0.0989 (0.939)	−0.0005	0.0832	0.0831 (0.951)	0.0013	0.0848	0.0837 (0.944)
$β_{21}$	0.0008	0.0269	0.0264 (0.946)	0.0044	0.0270	0.0266 (0.933)	0.0006	0.0261	0.0256 (0.948)
$β_{22}$	0.0002	0.1030	0.1030 (0.958)	0.0040	0.1060	0.1030 (0.936)	−0.0016	0.1020	0.1010 (0.948)
$β_{23}$	0.0021	0.0750	0.0742 (0.952)	0.0005	0.0752	0.0738 (0.940)	−0.0006	0.0762	0.0729 (0.942)
$τ$	$-$ 0.0023	0.0531	0.0550 (0.940)	−0.0231	0.0458	0.0346 (0.894)	−0.0140	0.0374	0.0349 (0.943)

SE: standard deviation of the point estimate; SEE: mean of the standard error estimates; CP: $95 %$ coverage probability.

4.4 Cumulative incidence function estimation

We evaluate the estimation accuracy of the cumulative incidence function by the proposed method. Data are from the Clayton copula with Weibull margins. We fit the true parametric model, the proposed method, and the marginal model. For the true and proposed methods, we estimate the cumulative incidence functions $P r (T_{1} \leq t, T_{1} \leq T_{2} | T_{1} > a, T_{2} > a, Z)$ on a sequence of pre-specified time points $t$ beyond age $a = 2$ and average the values across 1000 replications. The covariate vector $Z$ is set as $(6, 1, 1)$ , corresponding to a representative individual in the simulated population. For the marginal model, we use $P r (T_{1} \leq t | T_{1} > a, Z)$ as the cumulative incidence function by following Austin et al.⁹ Figure 2 illustrates that the proposed method produces almost identical cumulative incidence profiles as the true model under both Kendall’s $τ = 0.6$ and $0.2$ . The marginal model results in upward biases even when the dependence strength is weak at Kendall’s $τ = 0.2$ . Our observations are consistent with earlier works comparing copula and marginal models for cumulative incidence estimation under right censoring.⁹ In addition, we quantify the estimation error between the estimated and true cumulative incidence by the mean square errors (MSE) averaged over all time points and replications, which are $0.0010$ (SD $= 0.0016$ ) under $τ = 0.6$ and $0.0011$ (SD $= 0.0017$ ) under $τ = 0.2$ for the proposed method, indicating that the cumulative incidences are well estimated. We also consider the situation where data are generated from the Clayton copula with Loglogistic margins. As shown in Figure S1 in the Supplemental Appendix, the proposed method performs as well as the true parametric model, while the marginal method causes obvious overestimation.

Figure 2.

Simulation results for cumulative incidence function estimation when Kendall’s $τ = 0.6$ (left) and $0.2$ (right) under three models: the proposed Copula2 with semiparametric margins (Copula2-S, bottom dashed curve), the true Clayton copula with Weibull margins (Clayton-WB, solid curve), and the marginal semiparametric model (Marginal-S, top dashed curve). For the copula models, the cumulative incidence is $P r (T_{1} \leq t, T_{1} \leq T_{2} | T_{1} > a, T_{2} > a, Z)$ ; for the marginal model, the incidence is $P r (T_{1} \leq t | T_{1} > a, Z)$ .

5 Application to the CLHLS data

One major objective of CLHLS is to investigate the development of disability related to the activities of daily life in the elderly. A total of 12,969 participants free of disability over 60 years of age were enrolled in 2002 and underwent five examinations in 2005, 2008, 2011, 2014, and 2018. Six tasks essential to independent living were measured to evaluate disability status at each examination, including eating, toileting, bathing, dressing, indoor activity, and continence. The outcome of interest is time-to-disability, where disability is defined as the inability to perform any task independently.⁴² The time-to-disability is subject to dependent censoring by death, leading to semi-competing risks data. Due to the recruitment and intermittent assessment design, the time-to-disability is left-truncated and interval-censored. To the best of our knowledge, most existing CLHLS studies do not properly account for these complexities simultaneously.

We calculate time-to-disability from age 60. The truncation time is the enrollment age minus 60 years. We consider several covariates available at the baseline, including sex, education, and marital status. We exclude 180 subjects with missing disability status and 107 subjects with missing covariate values. The final analysis cohort contains 12,682 subjects. The mean follow-up time is 3.8 years with SD $= 4.7$ . Overall, 5863 (46.2%) subjects developed disability and 8128 (64.1%) subjects died. The basic characteristics of this final cohort are summarized in Table 3.

Table 3.
Characteristics of subjects in the cohort of CLHLS.

CLHLS (n = 12682)

Age

Mean (SD) 84.2 (11.4)

Median (range) 84.0 (61.0–120.0)

Sex (n, %)

Male (0) 5813 (45.8%)

Female (1) 6869 (54.2%)

Education (n, %)

$< 5$ years (0) 10193 (80.4%)

$\geq 5$ years (1) 2489 (19.6%)

Marital status (n, %)

Unmarried (0) 8382 (66.1%)

Married (1) 4300 (33.9%)

Censoring (n, %)

Scenario 1 4009 (31.6%)

Scenario 2 545 (4.3%)

Scenario 3 2810 (22.2%)

Scenario 4 5318 (41.9%)

	CLHLS (n = 12682)
Age
Mean (SD)	84.2 (11.4)
Median (range)	84.0 (61.0–120.0)
Sex (n, %)
Male (0)	5813 (45.8%)
Female (1)	6869 (54.2%)
Education (n, %)
$< 5$ years (0)	10193 (80.4%)
$\geq 5$ years (1)	2489 (19.6%)
Marital status (n, %)
Unmarried (0)	8382 (66.1%)
Married (1)	4300 (33.9%)
Censoring (n, %)
Scenario 1	4009 (31.6%)
Scenario 2	545 (4.3%)
Scenario 3	2810 (22.2%)
Scenario 4	5318 (41.9%)

CLHLS: Chinese Longitudinal Healthy Longevity Study.

We analyze the CLHLS data with the proposed method. We check various combinations of transformation functions (PH and PO models) and Bernstein polynomial degrees (3 to 6) and choose the model that produces the smallest AIC, which corresponds to the PH models with Bernstein polynomial degrees 3 for both margins. The AIC results are summarized in Table S6 in the Supplemental Appendix. The proposed method estimates the dependence parameters as $\hat{α} = 0.19$ and $\hat{κ} = 3.60$ , corresponding to $\hat{τ} = 0.83$ , indicating strong dependence between disability and death. This result suggests that the dependence pattern is neither Clayton (i.e. $α = 1$ ) nor Gumbel (i.e. $κ \to \infty$ ). We estimate the lower tail dependence as ${\hat{τ}}_{L} = 2^{- \hat{α} \hat{κ}} = 0.62$ and the upper tail dependence as ${\hat{τ}}_{U} = 2 - 2^{\hat{α}} = 0.86$ , suggesting that Copula2 captures both lower and upper tail dependence. We further illustrate the estimated Copula2 density function (i.e. $C_{\hat{α}, \hat{κ}}^{(1, 1)} (S_{1}, S_{2}) = \partial^{2} C_{\hat{α}, \hat{κ}} (S_{1}, S_{2}) / \partial S_{1} \partial S_{2}$ ) in Figure 3. Each point represents the estimated density value at the corresponding $S_{1}$ and $S_{2}$ values. The color reflects the estimated density, with the green color representing low density and the white color representing high density. Figure 3 demonstrates that Copula2 captures strong dependence in both upper and lower tails. Such strong dependence highlights the importance of accounting for the dependence between disability and death in modeling disability among the elderly.

Figure 3.

The Contour plot of the estimated Copula2 density function (i.e. $C_{\hat{α}, \hat{κ}}^{(1, 1)} (S_{1}, S_{2})$ ) in the joint modeling of disability and death. The x-axis and y-axis represent $1 - S_{1} (t)$ and $1 - S_{2} (t)$ , respectively. The values of $S_{1} (t)$ and $S_{2} (t)$ are sampled from the uniform distribution.

For the covariate effects on disability as shown in Table 4, by using the proposed method, sex (hazard ratio (HR) 0.814, confidence interval (CI) 0.779–0.850), education (HR 0.914, CI 0.862–0.969), and marital status (HR 0.815, CI 0.777–0.854) are significantly associated with reduced risks of developing disability. Specifically, females have a lower risk than males; people with higher education have a lower risk than people without; married people have a lower risk than unmarried people. These findings provide important insights into the prevention and management of disability. We also report results from a marginal semiparametric model on disability in Table 4. Although the estimated HRs in the marginal model are all smaller than 1, indicating the protective effects of these covariates, the coefficient values are dramatically different from that in the copula model. Specifically, the HRs are $0.947, 0.977$ , and $0.882$ for sex, education, and marital status in the marginal model, larger than the HRs ( $0.814, 0.914$ , and $0.815$ ) from the copula model. Moreover, the standard error estimates in the marginal model are larger than that of the copula model. Therefore, the marginal model reports non-significant associations ( $p$ -values $> 0.05$ ) between sex/education and disability.

Table 4.

Covariate effects on disability and death in CLHLS using the proposed Copula2-S model and the Marginal-S model.

		Copula2-S		Marginal-S for disability		Marginal-S for death
	Covariates	HR ( $95 %$ CI)	$p$ -value	HR ( $95 %$ CI)	$p$ -value	HR ( $95 %$ CI)	$p$ -value
Disability	Sex	$0.814$ $(0.779, 0.850)$	$4.2 \times 10^{- 20}$	$0.947$ $(0.894, 1.003)$	$6.5 \times 10^{- 2}$	–	–
	Education	$0.914$ $(0.862, 0.969)$	$2.6 \times 10^{- 3}$	$0.977$ $(0.906, 1.054)$	$5.5 \times 10^{- 1}$	–	–
	Marital	$0.815$ $(0.777, 0.854)$	$6.4 \times 10^{- 17}$	$0.882$ $(0.830, 0.937)$	$4.9 \times 10^{- 5}$	–	–
Death	Sex	$0.782$ $(0.749, 0.816)$	$4.5 \times 10^{- 29}$	–	–	$0.772$ $(0.735, 0.809)$	$3.0 \times 10^{- 26}$
	Education	$0.890$ $(0.839, 0.942)$	$7.1 \times 10^{- 5}$	–	–	$0.900$ $(0.855, 0.949)$	$1.3 \times 10^{- 3}$
	Marital	$0.781$ $(0.745, 0.819)$	$2.8 \times 10^{- 24}$	–	–	$0.861$ $(0.817, 0.906)$	$1.5 \times 10^{- 8}$

CLHLS: Chinese Longitudinal Healthy Longevity Study; HR: hazard ratio; CI: confidence interval.

For the modeling of death, the findings are consistent between the copula and marginal methods, as shown in Table 4. Specifically, both methods report significant protective effects for sex, education, and marital status. The $p$ -values from the copula method are generally smaller than the marginal model, probably due to the large sample size of CLHLS.

We estimate the cumulative incidence functions of disability for females and males under both Copula2-S and Marginal-S models. The rest covariates are set as median values. For Copula2-S, the cumulative incidence function is $P r (T_{1} \leq t, T_{1} \leq T_{2} | T_{1} > 60, T_{2} > 60, Z_{1}, Z_{2})$ . For Marginal-S, the cumulative incidence function is $P r (T_{1} \leq t | T_{1} > 60, Z_{1})$ . As shown in Figure 4, the proposed method finds that females below age 87 have lower cumulative incidence than males while females above 87 years have a higher incidence. The crossover at age 87 years is because males have higher mortality risks than females. About $60 %$ of males and $65 %$ of females develop disability by age 100 years. These findings are consistent with existing literature studying the relationship between sex and disability.^43,44 However, the marginal method does not reveal the crossover and generally overestimates the cumulative incidences, particularly after age 80 years. Overall, the estimated cumulative incidence functions by the proposed method provide valuable information to characterize the progression profiles for disability in the elderly.

Figure 4.

Estimation of cumulative incidence functions for disability in females and males by the proposed Copula2-S and the Marginal-S methods.

6 Discussion and conclusion

Semi-competing risks data frequently arise in clinical and epidemiological studies. With a high mortality risk for a disease, there could be incorrect estimation and inference results using a marginal method that ignores the dependence between disease and death. The situation becomes further complicated when the time of disease is subject to interval censoring and left truncation, which is commonly seen in many studies with the disease being intermittently assessed. We propose a flexible two-parameter copula-based semiparametric transformation model for semi-competing risks data subject to interval censoring and left truncation. The proposed method handles all possible censoring scenarios as depicted in Figure 1. By quantifying the dependence between disease and death, the proposed method provides correct estimations and highlights cumulative incidence profiles of disease progression associated with different baseline covariates. We develop a two-step sieve maximum likelihood estimation procedure and prove the sieve estimators’ asymptotic properties. All the methods have been built into an R package {CopulaCenR}, which is available on CRAN at https://cran.r-project.org/package=CopulaCenR.

Several copula selection methods are developed for semi-competing risks data under right censoring.^45,46 However, there is a lack of goodness-of-fit tests under interval censoring. In CLHLS, we use AIC to guide the model selection for simplicity. A formal goodness-of-fit test is highly desirable.

We use the semi-parametric transformation model for the marginal distributions instead of parametric models (e.g. Weibull). When the parametric assumption is satisfied, the copula model with correct parametric margins produces correct estimation and inference results. However, when the parametric assumption is unsatisfied, the copula model with misspecified parametric margins will lead to biased estimation and inference results, even when the copula structure is correctly specified.^47,34 We choose the semiparametric transformation model because of two advantages. First, the semiparametric model does not require a parametric assumption, which renders great flexibility in practice. In particular, its estimators of the finite-dimensional parameters (i.e. copula parameter(s) and regression coefficients) are proved to be as efficient as the estimators under the copula model with parametric margins. Second, the transformation model flexibly incorporates different survival assumptions (i.e. PH and PO). In contrast, the parametric model is usually restricted to one specific assumption (e.g. Weibull is PH, Loglogistic is PO). Because of its advantages, the semiparametric transformation model has been widely used in the literature on survival analysis.^32,48,24

It is well recognized that the PH/PO assumptions may not hold in practice. The semiparametric transformation model incorporates various model assumptions through the transformation function, which is usually selected by AIC or BIC (i.e. Bayesian information criterion). For example, Li et al.⁴⁹ find that the transformation function $g (x) = \log (1 + 0.6 x) / 0.6$ gives smaller BIC than transformation functions corresponding to PH/PO in the modeling of Alzheimer’s disease progression. One drawback of using a non-PH/PO transformation function is the lack of interpretation of covariate effects. Because one major goal of this paper is to estimate and interpret covariate effects on disability progression, we choose the marginal models among PH and PO transformation functions. We find that the copula model with PH transformation functions in both margins has the smallest AIC among all candidate models in CLHLS (as shown in Table S6 of Supplemental Appendix).

There are other possible approaches along the direction of the illness-death model, which can account for dependence between transitions by introducing a frailty effect. As discussed in several existing works,^50,51,24 the marginal survival functions are modeled differently between copula and frailty models. Specifically, the marginal function under the copula model only involves the time and covariate effects, whereas the marginal function under the frailty model includes time, covariate effects, and the frailty parameter. As a result, the model parameters (i.e. covariate effects) are usually not directly comparable between the two methods. In this article, the objectives of our real study lead us to pursue the copula model instead of the frailty-based illness-death model because the copula model offers a straightforward interpretation of covariate effects and dependence strength. In addition, the copula model accounts for the potential differential mortality between diseased and non-diseased groups through covariate effects and the dependence between $S_{1}$ and $S_{2}$ induced by the copula. In contrast, the illness-death model accounts for the potential mortality via modeling the transitions (i.e. diseased to death, non-diseased to death) and covariate effects. More details are discussed in the Supplemental Appendix.

We apply the proposed method to CLHLS, a large-scale longitudinal study of disability in the elderly, and successfully identify significant covariates associated with disability development, some of which are not detected in the marginal model. We find that the dependence between the time to disability and the time to death is extremely strong with Kendall’s $\hat{τ} = 0.83$ . In particular, the dependence lies in both upper and lower tails (i.e. ${\hat{τ}}_{L} = 0.62, {\hat{τ}}_{U} = 0.86$ ). It would be interesting to study the effects of covariates (i.e. age, sex) on the dependence in each tail in future work. In addition, to model disability development on the age scale, we set the baseline time at age 60 years, by which many covariates in CLHLS were not collected. One future direction is incorporating time-dependent covariates (i.e. BMI, exercise, disease burden) into the proposed method.

Our work is the first research on disability that adopts a solid statistical model that appropriately handles interval-censored and left-truncated semi-competing risks data. Our findings provide new insights and tools for the prevention and management of disability in the elderly. The proposed method can be applied to the general interval-censored and possibly left-truncated semi-competing risks data in clinical and epidemiological studies with intermittently assessed disease status.

Supplemental Material

sj-pdf-1-smm-10.1177_09622802221133552 - Supplemental material for Semiparametric copula method for semi-competing risks data subject to interval censoring and left truncation: Application to disability in elderly

Supplemental material, sj-pdf-1-smm-10.1177_09622802221133552 for Semiparametric copula method for semi-competing risks data subject to interval censoring and left truncation: Application to disability in elderly by Tao Sun, Yunlong Li, Zhengyan Xiao, Ying Ding and Xiaojun Wang in Statistical Methods in Medical Research

Footnotes

Acknowledgements

The authors thank the editor(s) and referees for their insightful comments and suggestions that have led to a significant improvement of this paper. The authors gratefully thank Peking University Open Access Research Database and National Archive of Computerized Data on Aging sponsored by National Institute of Aging (NIA/NIH) for providing the CLHLS data.

Declaration of conflicting interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The authors disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: The authors acknowledge the support (to XJW and TS) from the National Natural Science Foundation of China (grant nos. 72101261, 72141306), the Ministry of Education of China (grant no. 20JZD023), National Bureau of Statistics of China (grant no. 2021LZ18), and the Fund for Building World-class Universities (Disciplines) of Renmin University of China.

ORCID iDs

Tao Sun

Xiaojun Wang

Supplemental material

Supplemental material for this article is available online.

References

Sousa

Ferri

Acosta

et al. Contribution of chronic diseases to disability in elderly people in countries with low and middle incomes: a 10/66 dementia research group population-based survey. Lancet 2009; 374: 1821–1830.

Griffith

Raina

et al. Population attributable risk for functional disability associated with chronic conditions in Canadian older adults. Age Ageing 2010; 39: 738–745.

Fong

. Disability incidence and functional decline among older adults with major chronic diseases. BMC Geriatr 2019; 19: 1–9.

Juster

Suzman

. An overview of the health and retirement study. J Hum Res 1995; 30: 7–56.

Steptoe

Breeze

Banks

et al. Cohort profile: the English longitudinal study of ageing. Int J Epidemiol 2013; 42: 1640–1648.

Zeng

. Reliability of age reporting among the chinese oldest-old in the CLHLS datasets. In: Healthy longevity in China. Springer, 2008. pp. 61–78.

Majer

Nusselder

Mackenbach

et al. Mortality risk associated with disability: a population-based record linkage study. Am J Public Health 2011; 101: 9–15.

Jacob

Marron

Boudreau

et al. Age, race, and gender factors in incident disability. J Gerontol: Ser A 2018; 73: 194–197.

Austin

Lee

Fine

. Introduction to the analysis of survival data in the presence of competing risks. Circulation 2016; 133: 601–609.

10.

Fine

Jiang

Chappell

. On semi-competing risks data. Biometrika 2001; 88: 907–919.

11.

Emura

Sofeu

Rondeau

. Conditional copula models for correlated survival endpoints: individual patient data meta-analysis of randomized controlled trials. Stat Methods Med Res 2021; 30: 2634–2650.

12.

Lee

Rondeau

Haneuse

. Accelerated failure time models for semi-competing risks data in the presence of complex censoring. Biometrics 2017; 73: 1401–1412.

13.

Gao

Zeng

Lin

. Semiparametric regression analysis of interval-censored data with informative dropout. Biometrics 2018; 74: 1213–1222.

14.

Joly

Commenges

Helmer

et al. A penalized likelihood approach for an illness-death model with interval-censored data: application to age-specific incidence of dementia. Biostatistics 2002; 3: 433–443.

15.

Rouanet

Joly

Dartigues

et al. Joint latent class model for longitudinal data and interval-censored semi-competing events: application to dementia. Biometrics 2016; 72: 1123–1135.

16.

Zhou

Zhang

McLain

et al. Semiparametric regression of the illness-death model with interval censored disease incidence time: an application to the ACLS data. Stat Methods Med Res 2020; 29: 3707–3720.

17.

Peng

Fine

. Regression modeling of semicompeting risks data. Biometrics 2007; 63: 96–108.

18.

Chen

. Maximum likelihood analysis of semicompeting risks data with semiparametric regression models. Lifetime Data Anal 2012; 18: 36–57.

19.

Wang

. Evaluating association between two event times with observations subject to informative censoring. J Am Stat Assoc 2021; 0: 1–13.

20.

Kalbfleisch

Tai

. Statistical analysis of illness-death processes and semicompeting risks data. Biometrics 2010; 66: 716–725.

21.

Lamarca

Alonso

Gomez

et al. Left-truncated data with age as time scale: an alternative for survival analysis in the elderly population. J Gerontol Ser A: Biol Sci Med Sci 1998; 53: 337–343.

22.

Wang

Sun

. A pairwise pseudo-likelihood approach for left-truncated and interval-censored data under the Cox model. Biometrics 2021; 77: 1303–1314.

23.

Bogaerts

Lesaffre

. A new, fast algorithm to find the regions of possible support for bivariate interval-censored data. J Comput Graph Stat 2012; 13: 330–340.

24.

Sun

Ding

. Copula-based semiparametric regression method for bivariate data under general interval censoring. Biostatistics 2021; 22: 315–330.

25.

Petti

Eletti

Marra

et al. Copula link-based additive models for bivariate time-to-event outcomes with general censoring scheme. Comput Stat Data Anal 2022; 175: 107550.

26.

Sun

. Sieve maximum likelihood regression analysis of dependent current status data. Biometrika 2015; 102: 731–738.

27.

Sklar

. Fonctions de répartition à n dimensions et leurs marges. Publications de L’Institut de Statistique de L’Université de Paris 1959; 8: 229–231.

28.

Nelsen

. An introduction to copulas. New York: Springer-Verlag, 2006.

29.

Joe

. Multivariate models and dependence concepts. London: Chapman & Hall, 1997.

30.

Sun

Ding

. copulacenr: copula based regression models for bivariate censored data in R. R Journal 2020; 12: 266.

31.

Huang

Rossini

. Sieve estimation for the proportional-odds failure-time regression model with interval censoring. J Am Stat Assoc 1997; 92: 960–967.

32.

Zhou

Sun

. A sieve semiparametric maximum likelihood approach for regression analysis of bivariate interval-censored failure time data. J Am Stat Assoc 2017; 112: 664–672.

33.

Lee

Wong

Lam

et al. Analysis of clustered interval-censored data using a class of semiparametric partly linear frailty transformation models. Biometrics 2022; 78: 165–178.

34.

Sun

Liu

Cook

et al. Copula-based score test for bivariate time-to-event data, with application to a genetic study of AMD progression. Lifetime Data Anal 2019; 25: 546–568.

35.

Sun

Wang

Sun

. Estimation of the association for bivariate interval-censored failure time data. Scand J Stat 2006; 33: 637–649.

36.

Cui

Zhao

Sun

. A new copula model-based method for regression analysis of dependent current status data. Stat Interface 2018; 11: 463–471.

37.

Ding

Nan

. A sieve M-theorem for bundled parameters in semiparametric models, with application to the efficient estimation in a linear model for censored data. Ann Stat 2011; 39: 2795–3443.

38.

Lindfield

Penny

JET

. Microcomputers in numerical analysis. New York: Halsted Press, 1989.

39.

Andersen

Pohar Perme

. Inference for outcome probabilities in multi-state models. Lifetime Data Anal 2008; 14: 405–431.

40.

Varadhan

Xue

Bandeen-Roche

. Semicompeting risks in aging research: methods, issues and needs. Lifetime Data Anal 2014; 20: 538–562.

41.

Meira-Machado

Sestelo

. Estimation in the progressive illness-death model: a nonexhaustive review. Biom J 2019; 61: 245–263.

42.

Yuan

Mao

et al. Association of body mass index with disability in activities of daily living among chinese adults 80 years of age or older. JAMA Network Open 2018; 1: e181915.

43.

Fong

Sherris

Yap

. Forecasting disability: application of a frailty model. Scand Actuar J 2017; 2017: 125–147.

44.

Fuino

Wagner

. Duration of long-term care: socio-economic factors, type of care interactions and evolution. Insur: Math Econ 2020; 90: 151–168.

45.

Hsieh

Wang

Adam Ding

. Regression analysis based on semicompeting risks data. J R Stat Soc: Ser B (Statistical Methodology) 2008; 70: 3–20.

46.

Hsieh

Tsai

et al. Dependence model selection for semi-competing risks data. Open J Stat 2020; 10: 228.

47.

Kim

Silvapulle

. Comparison of semiparametric and parametric methods for estimating copulas. Comput Stat Data Anal 2007; 51: 2836–2850.

48.

Zeng

Mao

Lin

. Maximum likelihood estimation for semiparametric transformation models with interval-censored data. Biometrika 2016; 103: 253–271.

49.

Sun

. Penalized estimation of semiparametric transformation models with interval-censored data and application to Alzheimer’s disease. Stat Methods Med Res 2020; 29: 2151–2166.

50.

Goethals

Janssen

Duchateau

. Frailty models and copulas: similarities and differences. J Appl Stat 2008; 35: 1071–1079.

51.

Wienke

. Frailty models in survival analysis. New York: Chapman and Hall/CRC, 2010.

Supplementary Material

Please find the following supplemental material available below.

For Open Access articles published under a Creative Commons License, all supplemental material carries the same license as the article it is associated with.

For non-Open Access articles published, all supplemental material carries a non-exclusive license, and permission requests for re-use of supplemental material or any part of supplemental material shall be sent directly to the copyright owner as specified in the copyright notice associated with the article.

0.00 MB

0.69 MB

	Clayton-PM			Copula2-S			Marginal-S
Param	Bias	SE	SEE (CP)	Bias	SE	SEE (CP)	Bias	SE	SEE (CP)
Proportional hazards
$β_{11}$	0.0009	0.0275	0.0272 (0.947)	$-$ 0.0014	0.0278	0.0270 (0.946)	$-$ 0.0199	0.0412	0.0382 (0.890)
$β_{12}$	$-$ 0.0021	0.1070	0.1060 (0.949)	$-$ 0.0003	0.1070	0.1060 (0.943)	$-$ 0.0177	0.1590	0.1540 (0.928)
$β_{13}$	0.0019	0.0782	0.0768 (0.939)	0.0011	0.0753	0.0769 (0.955)	$-$ 0.0047	0.1140	0.1110 (0.941)
$β_{21}$	0.0009	0.0237	0.0236 (0.949)	$-$ 0.0001	0.0237	0.0234 (0.950)	–	–	–
$β_{22}$	$-$ 0.0006	0.0945	0.0919 (0.943)	0.0001	0.0930	0.0921 (0.950)	–	–	–
$β_{23}$	0.0014	0.0679	0.0664 (0.939)	0.0018	0.0649	0.0665 (0.958)	–	–	–
$τ$	0.0054	0.0258	0.0275 (0.955)	0.0091	0.0280	0.0278 (0.935)	–	–	–
Proportional odds
$β_{11}$	$-$ 0.0028	0.0504	0.0507 (0.954)	$-$ 0.0028	0.0502	0.0499 (0.948)	$-$ 0.0243	0.0638	0.0547 (0.877)
$β_{12}$	$-$ 0.0087	0.2010	0.2000 (0.949)	$-$ 0.0063	0.2020	0.1990 (0.944)	$-$ 0.0438	0.2520	0.2130 (0.872)
$β_{13}$	0.0023	0.1460	0.1450 (0.946)	0.0001	0.1440	0.1440 (0.949)	0.0063	0.1820	0.1570 (0.906)
$β_{21}$	$-$ 0.0008	0.0449	0.0452 (0.962)	$-$ 0.0016	0.0450	0.0448 (0.955)	–	–	–
$β_{22}$	$-$ 0.0113	0.1810	0.1790 (0.949)	$-$ 0.0077	0.1800	0.1780 (0.945)	–	–	–
$β_{23}$	0.0031	0.1310	0.1290 (0.949)	$-$ 0.0027	0.1310	0.1290 (0.949)	–	–	–
$τ$	0.0028	0.0249	0.0264 (0.963)	0.0123	0.0255	0.0260 (0.945)	–	–	–

Semiparametric copula method for semi-competing risks data subject to interval censoring and left truncation: Application to disability in elderly

Abstract

Keywords

1 Introduction

2 Notation and likelihood

2.1 Copula model for semi-competing risks data

3.1 Sieve likelihood with Bernstein polynomials

3.2 Estimation procedure for sieve maximum likelihood estimators θ ^ n

3.3 Estimation of cumulative incidence function

4.1 Data generation

4.2 Parameter estimation under correct copula specification

Supplemental Material

sj-pdf-1-smm-10.1177_09622802221133552 - Supplemental material for Semiparametric copula method for semi-competing risks data subject to interval censoring and left truncation: Application to disability in elderly

Footnotes

Acknowledgements

Declaration of conflicting interests

Funding

ORCID iDs

Supplemental material

References

Supplementary Material

3.2 Estimation procedure for sieve maximum likelihood estimators ${\hat{θ}}_{n}$