Causal mediation analysis for time-to-event mediator and outcome in the presence of left truncation

Abstract

We propose a causal mediation approach to semi-competing risks under left truncation sampling by considering an intermediate event as a mediator and a terminal event as an outcome. We focus on the causal relationship from exposure to the terminal outcome in relation to the intermediate event. In particular, we study the direct effect, the effect of exposure on the terminal event that is not through the intermediate event, and the indirect effect—the effect of exposure on the terminal event that is mediated through the intermediate event. We propose nonparametric and semiparametric methods, both accounting for left truncation. The nonparametric estimator can be viewed as a model-free time-varying Nelson–Aalen estimator that is robust to model misspecification. The semiparametric estimator calculated with the Cox proportional hazards model enjoys flexibility in adjusting for potential confounders as covariates. The asymptotic properties for both estimators, including uniform consistency and weak convergence, were established using the martingale theorem and functional delta method. The finite sample performance of the proposed estimators was evaluated through extensive numerical studies that investigated the influences of left truncation, confounding, and sample size. The utility of the proposed methods was illustrated using a hepatitis study.

Keywords

Causal mediation analysis time-to-event outcomes time-to-event mediator left truncation

1. Introduction

Motivated by the community-based prospective cohort study REVEAL,¹ we consider the natural history of chronic hepatitis B. Scholars have established that, relative to the average person, hepatitis B patients have an increased incidence of liver cirrhosis and a higher mortality rate; however, we are particularly interested in whether hepatitis B affects mortality mediated by increasing the incidence of liver cirrhosis as a function of age. In developing an analytical methodology for this question, two statistical topics must be addressed: Semi-competing risks and left truncation. When two events are of interest, semi-competing risks are frequently encountered in biomedical research, where the intermediate event may be masked by the terminal event, but not vice versa.² In our motivating example, the fact that liver cirrhosis incidence may be masked by mortality, but having liver cirrhosis remains subject to mortality, represents a typical semi-competing risks problem. The challenge of left truncation arises from the design of most cohort studies, including the REVEAL cohort, where the participants were recruited at certain time points, and individuals who had not survived up to the time of recruitment were not included in the study. This methodology may go beyond the scope of the motivation study. We propose developing novel methods of quantifying the mediation effect of an exposure or intervention on a time-to-event outcome through a time-to-event mediator in the presence of semi-competing risks and left truncation.

Current approaches to semi-competing risks can be classified into several categories. Approaches in one category seek to characterize the joint survivor function of the times to intermediate and terminal events. For example, Day et al.³ considered the predictive hazard ratio⁴ for the dependence of two outcomes under the assumption that the ratio is constant, and the joint survivor function corresponds to the Clayton copula model,^2,5 which is also equivalent to gamma frailty models. Approaches in another category apply multistate modeling,⁶ where semi-competing risks data are modeled for transitions between different states: healthy (i.e., without any intermediate or terminal event), diseased (i.e., having an intermediate event), and dead (i.e., having a terminal event). Wang⁷ proposed a nonparametric estimator for the survival function of the sojourn time between states; by incorporating a gamma frailty in the transition hazards of the multistate model, Xu⁸ derived a new nonparametric maximum likelihood estimator for the copula parameter under the Clayton copula model. Approaches in another category focus on the survivor’s average causal effect (SACE; Frangakis and Rubin⁹ and Zhang and Rubin¹⁰). The SACE only concerns the effect on intermediate event among those who would have been survivors of terminal events, regardless of their exposure. The first two categories successfully characterized the dependence and transition probabilities of the events, but the causal interpretations of these statistical quantities remain elusive, whereas the third category focuses on a different causal effect than mediation. To this end, we propose an alternative approach to semi-competing risks using causal mediation modeling.

Causal mediation analysis provides a general framework for characterizing the mediation mechanism linking exposure and outcome. Causal mediation models originated from social science literature^11,12 and were theoretically consolidated by a potential outcome framework.¹³ Causal mediation models have been rigorously developed for counterfactual definitions, identifiability assumptions, and estimators.^14–17 Causal mediation models decompose the effect of the exposure on the outcome into an indirect effect (IE): Effect on the outcome mediated through an intermediary factor (termed a “mediator”); and a direct effect (DE), which is not through the mediator. Conceptually, formulating the semi-competing risks using a causal mediation model, where the intermediate and terminal events are the mediator and outcome, respectively, is straightforward. Although causal mediation analyses of time-to-event outcomes are available,^18–21 previously published methods require a fully observable mediator. In semi-competing risks data, the intermediate event is masked by either conventional right censoring or a terminal event. Our recent work has proposed the nonparametric estimator for mediation analyses accounts for the time to an intermediate event.²² The invited discussions by Fulcher et al.²³ and Stensrud et al.²⁴ pointed to several issues, including the ill-defined time to the intermediate event after the occurrence of the terminal event, cross-world independence, and Markov assumption. While these issues have been discussed and clarified in the rejoinder,²⁵ we propose a new estimand by adopting interventional approach^26,27 to address the former two, and an analytic approach to evaluate and weaken the Markov assumption. Moreover, we found that the proposed estimand is closely related to the separable effects of Stensrud et al.,²⁸ which is further discussed. However, the current mediation analyses for survival outcomes do not account for the bias resulting from the left truncation. Therefore, our proposed approach attempts to bridge causal mediation modeling and bivariate survival analyses under left truncation.

Left truncation can be construed as a mechanism of missing data, particularly when the time scale of interest is age.²⁹ Under the left-truncation phenomenon, researchers tend to recruit individuals with slower-than-usual progression of a disease (overall, a healthier condition than non-participants), which leads to biased sampling of the underlying population. For the motivation study on hepatitis, we were interested in the effect as a function of age because of its biological implications. From 1991 to 1992, individuals aged 30–65 years were recruited to participate in the REVEAL cohort. Consequently, those who died before the relevant age in the same recruitment period did not have any opportunity to participate in the study. Therefore, an analysis that does not account for left truncation usually leads to a substantial overestimation of the survival time.³⁰ To adjust for the bias caused by the left truncation, we adopted the conditional likelihood function approach,^31,32 where the likelihood is calculated by conditioning on the observed left truncation time. This approach enabled us to focus more on the times to liver cirrhosis and mortality without imposing fully parametric assumptions on the truncation time and was thus robust to mis-specified distributions of the truncation time.

2. Assumptions and estimands

For subject $i$ where $i = 1, \dots, m$ and $m$ is the sample size, let ${\tilde{T}}_{2 i}$ , ${\tilde{T}}_{1 i}$ ${\tilde{A}}_{i}$ , ${\tilde{Z}}_{i}$ , and ${\tilde{W}}_{i}$ be the terminal event, intermediate event, left truncation times, exposure, and confounders, respectively, in the underlying population. With these random variables, we further denote the stochastic processes of the intermediate and terminal events as ${\tilde{N}}_{1 i} (t) = I ({\tilde{T}}_{1 i} \leq t)$ and ${\tilde{N}}_{2 i} (t) = I ({\tilde{T}}_{2 i} \leq t)$ where $I (\cdot)$ is an indicator function. We made the following statistical assumptions.

(1.1)
$(\tilde{A}, \tilde{C}) ⊥ ({\tilde{T}}_{1}, {\tilde{T}}_{2}) | (\tilde{Z}, \tilde{W})$ : both time to the underlying left truncation and time to censoring are independent of the times to the underlying intermediate and terminal events given the underlying exposure and confounders.
(1.2)
$τ$ is the end time of study. If ${\tilde{T}}_{2} < {\tilde{T}}_{1}$ , then ${\tilde{T}}_{1} > τ$ ^8,33,34: if the intermediate event has not occurred prior to the terminal event, it will never occur by the end of the study.
(1.3)
For $t \in [0, τ]$ , $pr ({\tilde{T}}_{2} \leq t | Z, {\tilde{T}}_{1} \leq t^{-}, \tilde{W}) < 1, pr ({\tilde{T}}_{2} \leq t | Z, {\tilde{T}}_{1} > t^{-}, \tilde{W}) < 1$ and $G (t | Z, W) < 1$ , where $G (\cdot)$ is the cdf of the censoring time.
(1.4)
$π_{z_{a} n_{1}} (t) = pr (\tilde{Z} = z_{a}, {\tilde{N}}_{1} (t^{-}) = n_{1} | \tilde{W}) > 0$ : all possible combinations of $(\tilde{Z}, {\tilde{N}}_{1} (t^{-}))$ given $\tilde{W}$ can occur with nonzero probabilities.
(1.5)
For $s \leq t$ , $pr (d {\tilde{N}}_{2} (t) = 1 | {\tilde{N}}_{2} (t^{-}) = 0, {\tilde{N}}_{1} (t^{-}) = n_{1}, {\tilde{N}}_{1} (s^{-}) = n_{1}^{}, \tilde{Z}, \tilde{W}) = pr (d {\tilde{N}}_{2} (t) = 1 | {\tilde{N}}_{2} (t^{-}) = 0, {\tilde{N}}_{1} (t^{-}) = n_{1}, \tilde{Z}, \tilde{W})$ : the hazard of the terminal event at time $t$ only depends on the status of the intermediate event right before $t$ , not on its history.
These assumptions are used to ensure the independence mechanism for left truncation and right censoring: $(\tilde{A}, \tilde{C}) ⊥ ({\tilde{T}}_{1}, {\tilde{T}}_{2}) | (\tilde{Z}, \tilde{W})$ ^32,35,36 (Assumption (1.1)), the sequence of intermediate and terminal events (Assumption (1.2)), and stochasticity for the survival and censoring times, as well as the occurrence of joint combinations of exposure and intermediate event status (Assumptions (1.3) and (1.4)). In addition, our approach depends on a Markov assumption (Assumption (1.5)) that is assumed in the multistate and frailty approaches.^8,37 The Markov assumption states that given the same status of the intermediate event at time $t^{-}$ and the fact that subjects with the same $\tilde{Z}$ survive until $t^{-}$ , the risk of having a terminal event at $t$ is independent of when the intermediate event occurs.

Note that the Markov assumption is evaluated in the data application, and the evaluation was conducted using the Cox model, which includes the duration since the occurrence of the intermediate event as a covariate. Such diagnostic analysis supports the plausibility of the Markov assumption (details in the Web Appendix). In addition to assumption diagnostics, inclusion of the duration-related covariates may be used as an analytic approach to satisfy or weaken the Markov assumption, which can easily be implemented by our semiparametric approach (details are in the “Semiparametric estimator” section). In addition, the Markov assumption and Assumption (2.1) in our study are related to the assumption required for causal interpretations of the Kaplan–Meier estimates of the time-varying covariate by Sjölander³⁸ (details are provided in the Discussion section).

We introduce potential outcome notations for intermediate and terminal event processes. We denote $n_{2} (t^{-})$ as the setting of ${\tilde{N}}_{2} (t^{-})$ , ${\tilde{N}}_{1} (t; z)$ is the counterfactual process of ${\tilde{N}}_{1} (\cdot)$ at $t$ had $\tilde{Z}$ been set to $z$ . ${\tilde{N}}_{1} (t; z, n_{2} (t^{-}) = 0)$ is the counterfactual process of ${\tilde{N}}_{1} (\cdot)$ at $t$ , had $\tilde{Z}$ and ${\tilde{N}}_{2} (t^{-})$ been set to $z$ and 0. ${\tilde{N}}_{2} (t; z)$ is the counterfactual process of ${\tilde{N}}_{2} (\cdot)$ at $t$ had $\tilde{Z}$ been set to $z$ , and ${\tilde{N}}_{2} {t; z, n_{1} (t^{-}), n_{2} (t^{-}) = 0}$ is the counterfactual process for ${\tilde{N}}_{2} (\cdot)$ at $t$ had $\tilde{Z}$ , ${\tilde{N}}_{1} (t; z^{}, n_{2} (t^{-}) = 0)$ and ${\tilde{N}}_{2} (t^{-})$ been set to $z$ (not necessarily equal to $z^{}$ ), $n_{1} (t^{-})$ and 0, respectively. To identify the proposed counterfactual hazard defined in the “Direct and indirect effects” section, under the interventional approach,^26,27 we make the following causal assumptions. (2.1)
${\tilde{N}}_{2} {t; z, n_{1} (t^{-}), n_{2} (t^{-})} ⊥ {\tilde{Z}, {\tilde{N}}_{1} (t^{-}), {\tilde{N}}_{2} (t^{-})} | \tilde{W}$ : there is no unmeasured confounding for the association of ${\tilde{N}}_{2} (t)$ and ${\tilde{Z}, {\tilde{N}}_{1} (t^{-}), {\tilde{N}}_{2} (t^{-})}$ given $\tilde{W}$ .
(2.2)
${\tilde{N}}_{1} {t; z, n_{2} (t)} ⊥ {\tilde{Z}, {\tilde{N}}_{2} (t)} | \tilde{W}$ : there is no unmeasured confounding for the association of ${\tilde{N}}_{1} (t)$ and ${\tilde{Z}, {\tilde{N}}_{2} (t)}$ given $\tilde{W}$ .
(2.3)
${\tilde{N}}_{1}^{} {t; z, n_{2} (t) = 0}$ is a random draw from the distribution of ${\tilde{N}}_{1} {t; z, n_{2} (t) = 0}$ .
We stress that the above assumptions do not include the cross-world independence assumption, and are much weaker than those used by Huang.²² Assumptions (2.1) and (2.2) can be illustrated by using a single world intervention graph (SWIG; Richardson and Robins;³⁹ Figure 1). For $t_{1} < t_{2}$ , we denote the time-specific error terms: $ϵ (t_{1})$ for ${\tilde{N}}_{1} (t_{1})$ , $ϵ^{} (t_{1})$ for ${\tilde{N}}_{1}^{} (t_{1})$ , $e (t_{1})$ , $e (t_{2})$ , $ε (t_{1})$ and $ε (t_{2})$ for ${\tilde{N}}_{2} (\cdot)$ , and the following confounders ${\tilde{W}}_{d}$ , ${\tilde{W}}_{b}$ , ${\tilde{W}}_{c}$ , and ${\tilde{W}}_{a}$ need to be adjusted to ensure that
$\begin{aligned} {\tilde{N}}_{1} {t; z, n_{2} (t)} ⊥ Z | {\tilde{W}}_{d}, \\ {\tilde{N}}_{1} {t; z, n_{2} (t)} ⊥ {\tilde{N}}_{2} (t) | {\tilde{W}}_{b}, \\ {\tilde{N}}_{2} {t; z_{a}, {\tilde{N}}_{1}^{} {t^{-}; z_{b}, n_{2} (t^{-})}, n_{2} (t^{-})} ⊥ \tilde{Z} | {\tilde{W}}_{c}, \end{aligned}$
(1)

$\begin{aligned} {\tilde{N}}_{2} {t; z_{a}, {\tilde{N}}_{1}^{} {t^{-}; z_{b}, n_{2} (t^{-})}, n_{2} (t^{-})} ⊥ {\tilde{N}}_{2} (t^{-}) | {\tilde{W}}_{a} . \end{aligned}$
(2)
Therefore, outcome-outcome ( ${\tilde{N}}_{2} (t_{1})$ - ${\tilde{N}}_{2} (t_{2})$ ) confounder ( ${\tilde{W}}_{a}$ ), mediator-outcome confounder ( ${\tilde{W}}_{b}$ ), exposure-outcome confounder ( ${\tilde{W}}_{c}$ ) and exposure-mediator confounder ( ${\tilde{W}}_{d}$ ) need to be adjusted to ensure Assumptions (2.1) and (2.2). We propose sensitivity analyses to evaluate Assumptions (2.1) and (2.2) in the Web Appendix. Moreover, as mentioned in equations (1) and (2), we use $z_{a}$ and $z_{b}$ to describe the potential intervention values for ${\tilde{N}}_{2} (\cdot)$ and ${\tilde{N}}_{1} (\cdot)$ , respectively, received from exposure, which can be set at distinct levels.

Figure 1.
Single world intervention graph (SWIG) that describe the intervention of $Z = z_{b}$ and ${\tilde{N}}_{2}^{} (t) = n_{2} (t) = 0$ in ${\tilde{N}}_{1}^{} (t)$ as well as intervention of $Z = z_{a}$ , ${\tilde{N}}_{2} (t^{-}) = n_{2} (t^{-}) = 0$ and ${\tilde{N}}_{1}^{*} (t^{-}; z_{b}, n_{2} (t^{-}) = 0)$ in ${\tilde{N}}_{2} (t)$ .
2.1. Direct and indirect effects

Our key estimand is denoted by

\frac{d}{d t} Λ (t; z_{a}, z_{b}) \overset{d e f}{=} pr [{\tilde{N}}_{2} {t; z_{a}, {\tilde{N}}_{1}^{*} (t^{-}; z_{b}, n_{2} (t^{-}) = 0), n_{2} (t^{-}) = 0} = 1],

where

{\tilde{N}}_{1}^{*} (t^{-}; z_{b}, n_{2} (t^{-}) = 0)

is a random draw from an identical probabilistic distribution of

{\tilde{N}}_{1} (t^{-}; z_{b}, n_{2} (t^{-}) = 0)

, using an idea analogous to the interventional approach.^26,27 We express the counterfactual hazard without conditioning the covariates

\tilde{W}

to simplify the notations, but it is trivial to stratify by

\tilde{W}

throughout the development. Under Assumptions (2.1)–(2.3) prove that (detailed in the Appendix)

\frac{d}{d t} Λ (t; z_{a}, z_{b}) = \sum_{n_{1} \in {0, 1}} \frac{d}{d t} Λ_{n_{1}} (t | z_{a}) w_{n_{1}} (t | z_{b}),

(3)

where

\frac{d}{d t} Λ_{n_{1}} (t | z_{a})

pr {d {\tilde{N}}_{2} (t) = 1 | \tilde{Z} = z_{a}, {\tilde{N}}_{1} (t^{-}) = n_{1}, {\tilde{N}}_{2} (t^{-}) = 0}

and

w_{n_{1}} (t | z_{b}) = pr {{\tilde{N}}_{1} (t^{-}) = n_{1} | \tilde{Z} = z_{b}, {\tilde{N}}_{2} (t^{-}) = 0}

. The identification formula is identical to that by Huang,²² but is obtained under weaker assumptions without cross-world independence. Using the counterfactual hazard, we can easily obtain the cumulative counterfactual hazard:

Λ (t; z_{a}, z_{b}) = \sum_{n_{1} \in {0, 1}} \int_{0}^{t} w_{n_{1}} (s | z_{b}) d Λ_{n_{1}} (s | z_{a}) = \sum_{n_{1} \in {0, 1}} Λ_{n_{1}} (t; z_{a}, z_{b}),

and the counterfactual survival probabilities

S (t; z_{a}, z_{b}) = \exp {- Λ (t; z_{a}, z_{b})},

where

Λ_{n_{1}} (t; z_{a}, z_{b}) = \int_{0}^{t} w_{n_{1}} (s | z_{b}) d Λ_{n_{1}} (s | z_{a})

In the context of the definitions by Pearl¹⁵ and Robins,¹⁶ DE is defined as the effect of $\tilde{Z}$ on the terminal event not mediated by the intermediate event, as $ρ_{D E} (t) = S (t; z_{a} = 1, z_{b} = 0) / S (t; z_{a} = 0, z_{b} = 0) = \exp {- Δ_{D E} (t)}$ , where

\begin{aligned} Δ_{D E} (t) & = Λ (t; z_{a} = 1, z_{b} = 0) - Λ (t; z_{a} = 0, z_{b} = 0) \\ = Λ (t; 1, 0) - Λ (t; 0, 0), \end{aligned}

(4)

IE is defined as the effect of

\tilde{Z}

on the terminal event mediated by affecting the intermediate event, as

ρ_{I E} (t) = S (t; z_{a} = 1, z_{b} = 1) / S (t; z_{a} = 1, z_{b} = 0) = \exp {- Δ_{I E} (t)}

, where

\begin{aligned} Δ_{I E} (t) & = Λ (t; z_{a} = 1, z_{b} = 1) - Λ (t; z_{a} = 1, z_{b} = 0) \\ = Λ (t; 1, 1) - Λ (t; 1, 0) . \end{aligned}

(5)

Additionally, the total effect (TE) can be obtained using

Δ_{T E} (t) = Λ (t; 1, 1) - Λ (t; 0, 0) = Δ_{D E} (t) + Δ_{I E} (t)

. Note that the effect measure is not limited to the cumulative hazard difference. Because our estimation focuses on the cumulative hazard, other effect measures such as the cumulative hazard ratio,

Λ (t; 1, 0) / Λ (t; 0, 0)

for

D E

Λ (t; 1, 1) / Λ (t; 1, 0)

for

I E

and

Λ (t; 1, 1) / Λ (t; 0, 0)

for

T E

; survival probability ratios,

\exp {- Δ_{D E} (t)}

\exp {- Δ_{I E} (t)}

and

\exp {- Δ_{T E} (t)}

, can be easily obtained.

3. Correspondence of separable effect

Separable effect has recently been proposed to identify causal effects in the setting of competing risks.²⁸ In the separable effect approach, the treatment $\tilde{Z}$ can be decomposed into ${\tilde{Z}}_{1}$ and ${\tilde{Z}}_{2}$ , where ${\tilde{Z}}_{1}$ affects ${\tilde{N}}_{1} (\cdot)$ without any direct effect on ${\tilde{N}}_{2} (\cdot)$ and ${\tilde{Z}}_{2}$ affects ${\tilde{N}}_{2} (\cdot)$ without any direct effect on ${\tilde{N}}_{1} (\cdot)$ .

We denote the potential values for separable exposures ${\tilde{Z}}_{1}$ and ${\tilde{Z}}_{2}$ as $z_{1}$ and $z_{2}$ , respectively, to avoid confusion with $z_{a}$ and $z_{b}$ in the proposed method. Although both $z_{1}$ and $z_{b}$ affect the mediator ${\tilde{N}}_{1} (\cdot)$ , and $z_{2}$ and $z_{a}$ affect the outcome ${\tilde{N}}_{2} (\cdot)$ , their influence is on different effect measures. In our proposed method, $z_{b}$ affects the prevalence of the mediator, and $z_{a}$ affects the hazard of the outcome conditional on the mediator, as shown in (3). In the separable effect, $z_{1}$ and $z_{2}$ affect cause-specific hazard of the mediator and outcome, respectively, as will be shown in the following dismissible component conditions.

To establish a connection with Stensrud et al.,²⁸ we derived the counterfactual hazard function under competing risks, using the assumptions specified by Stensrud et al.²⁸ The counterfactual counting processes of ${\tilde{T}}_{1}$ and ${\tilde{T}}_{2}$ in the separable model are ${\tilde{N}}_{1} (\cdot; {\tilde{Z}}_{1} = z_{1}, {\tilde{Z}}_{2} = z_{2})$ and ${\tilde{N}}_{2} (\cdot; {\tilde{Z}}_{1} = z_{1}, {\tilde{Z}}_{2} = z_{2})$ , respectively. The key to identification of separable effect is the two dismissible component conditions, which can be expressed as:

\begin{aligned} pr {{\tilde{N}}_{2} (t; {\tilde{Z}}_{1} = 1, {\tilde{Z}}_{2} = z_{2}) = 1 | {\tilde{N}}_{1} (t^{-}; {\tilde{Z}}_{1} = 1, {\tilde{Z}}_{2} = z_{2}) = {\tilde{N}}_{2} (t^{-}; {\tilde{Z}}_{1} = 1, {\tilde{Z}}_{2} = z_{2}) = 0} \\ = pr {{\tilde{N}}_{2} (t; {\tilde{Z}}_{1} = 0, {\tilde{Z}}_{2} = z_{2}) = 1 | {\tilde{N}}_{1} (t^{-}; {\tilde{Z}}_{1} = 0, {\tilde{Z}}_{2} = z_{2}) = {\tilde{N}}_{2} (t^{-}; {\tilde{Z}}_{1} = 0, {\tilde{Z}}_{2} = z_{2}) = 0}, \\ pr {{\tilde{N}}_{1} (t; {\tilde{Z}}_{1} = z_{1}, {\tilde{Z}}_{2} = 1) = 1 | {\tilde{N}}_{1} (t^{-}; {\tilde{Z}}_{1} = z_{1}, {\tilde{Z}}_{2} = 1) = {\tilde{N}}_{2} (t^{-}; {\tilde{Z}}_{1} = z_{1}, {\tilde{Z}}_{2} = 1) = 0} \\ = pr {{\tilde{N}}_{1} (t; {\tilde{Z}}_{1} = z_{1}, {\tilde{Z}}_{2} = 0) = 1 | {\tilde{N}}_{1} (t^{-}; {\tilde{Z}}_{1} = z_{1}, {\tilde{Z}}_{2} = 0) = {\tilde{N}}_{2} (t^{-}; {\tilde{Z}}_{1} = z_{1}, {\tilde{Z}}_{2} = 0) = 0} . \end{aligned}

Thus, with the above notations and assumptions, we define the following counterfactual hazard function of the separable effect method, denoted by

Λ^{s e p} (t; z_{1}, z_{2})

, which can be estimated using the following second equality under the above dismissible component conditions (details in the Web Appendix):

\begin{aligned} Λ^{s e p} (t; z_{1}, z_{2}) \overset{d e f}{=} & \int_{0}^{t} pr {{\tilde{N}}_{2} (s; {\tilde{Z}}_{1} = 1, {\tilde{Z}}_{2} = z_{2}) = 1 | {\tilde{N}}_{2} (s^{-}; {\tilde{Z}}_{1} = 1, {\tilde{Z}}_{2} = z_{2}) = 0} d s \\ = & \int_{0}^{t} w_{0}^{s e p} (s | z_{1}, z_{2}) d Λ_{0} (s | z_{2}), \end{aligned}

(6)

where

w_{0}^{s e p} (s | z_{1}, z_{2}) = \frac{e^{- \int_{0}^{s} \sum_{i \in {1, 2}} h_{i} (u | Z = z_{i}) d u}}{\int_{s^{-}}^{\infty} h_{2} (u | Z = z_{2}) e^{- \int_{0}^{u} \sum_{i \in {1, 2}} h_{i} (v | Z = z_{i}) d v} d u},

and

h_{i} (u | Z = z) = pr {{\tilde{N}}_{i} (u) = 1 | {\tilde{N}}_{1} (u^{-}) = {\tilde{N}}_{2} (u^{-}) = 0, \tilde{Z} = z}

. On the other hand, our proposed estimand

Λ (t; z_{a}, z_{b})

under competing risks leads to the following identification formula under Assumptions (2.1)–(2.3)

\int_{0}^{t} w_{0} (s | z_{b}) d Λ_{0} (s | z_{a}),

(7)

where

w_{0} (s | z_{b})

can be expressed as

\frac{e^{- \int_{0}^{s} \sum_{i \in {1, 2}} h_{i} (u | Z = z_{b}) d u}}{\int_{s^{-}}^{\infty} h_{2} (u | Z = z_{b}) e^{- \int_{0}^{u} \sum_{i \in {1, 2}} h_{i} (v | Z = z_{b}) d v} d u} .

Therefore, although the causal mechanisms differ between the proposed approach and the separable effects approach, their identifications are similar, as indicated by equations (6) and (7). Furthermore, treating

z_{a} = z_{2}

and

z_{b} = z_{1}

, the expression for equation (6) is similar to (7). The two equations are identical in setting

z_{1} = z_{2} = z_{a} = z_{b}

. Moreover, both (6) and (7) can be viewed as weighted integrals for

d Λ_{0} (s | z_{a})

with different weight functions

w_{0}^{s e p} (\cdot | z_{1}, z_{2})

and

w_{0} (\cdot | z_{b})

, which describes the probability of

min {{\tilde{T}}_{1}, {\tilde{T}}_{2}}

being greater than

s^{-}

, conditional on

{\tilde{T}}_{2}

being greater than

s^{-}

under different intervention mechanisms for the exposure.

4. Proposed estimators

Here, we construct the estimation for the proposed method based on the likelihood approach (details are in Section S.3 of the Web Appendix). Moreover, from the data with left truncation, we observe only individuals with information on $(T_{2}, T_{1}, A, Z)$ , where $(T_{2}, T_{1}, A, Z) = ({\tilde{T}}_{2}, {\tilde{T}}_{1}, \tilde{A}, \tilde{Z})$ for those satisfying $A \leq T_{2}$ . The corresponding processes under left truncation and right censoring are given by $N_{1} (t) I (T_{1} \leq t, T_{1} \leq X)$ and $I (X \leq t, δ = 1)$ , respectively, where $X = min {T_{2}, C}$ and $δ = I (T_{2} \leq C)$ .

4.1. Nonparametric estimator

We consider a nonparametric approach where potential confounding factors are adjusted using stratification. According to the conditional likelihood function mentioned in S.3.1 of the Web Appendix (equation (A8)), the $z$ -stratified conditional log-likelihood, denoted as $\log L_{c} {Λ_{n_{1}} (\cdot | z)}$ , of the entire data set under right censoring and left truncation is equal to:

\begin{aligned} \sum_{i = 1}^{m} { & \int_{0}^{\infty} \log d Λ_{n_{1}} (s | z) I {N_{1 i} (s^{-}) = n_{1}, Z_{i} = z} d N_{2 i} (s) \\ - \int_{A_{i} \lor T_{n_{1} i}}^{X_{i}} I {N_{1 i} (s^{-}) = n_{1}, Z_{i} = z} d Λ_{n_{1}} (s | z)} . \end{aligned}

Let

d {\hat{Λ}}_{n_{1}} (s | z) = {\hat{Λ}}_{n_{1}} (s | z) - {\hat{Λ}}_{n_{1}} (s^{-} | z)

be the estimated jump of

Λ_{n_{1}} (\cdot | z)

s

from this likelihood function. Then, we obtain

d {\hat{Λ}}_{n_{1}} (X_{j} | z)

by solving the root of

\frac{\partial \log L_{c} {Λ_{n_{1}} (\cdot | z)}}{\partial d Λ_{n_{1}} (X_{j} | z)} = 0

. Thus, we obtain the nonparametric estimator

d {\hat{Λ}}_{n_{1}} (X_{j} | z) = \frac{\sum_{i = 1}^{m} d N_{2 n_{1} i}^{N P} (X_{j} | z)}{\sum_{i = 1}^{m} Y_{n_{1} i}^{N P} (X_{j} | z)} I {\sum_{i = 1}^{m} Y_{n_{1} i}^{N P} (X_{j} | z) > 0},

in which

d N_{2 n_{1} i}^{N P} (s | z) = d N_{2 i} (s) I {N_{1 i} (s^{-}) = n_{1}, Z_{i} = z}, Y_{n_{1} i}^{N P} (s | z) = I {(A_{i} \lor T_{n_{1} i}) \leq s \leq X_{i}} I (N_{1 i} (s^{-}) = n_{1}, Z_{i} = z), T_{n_{1} i} = I (n_{1} = 1) T_{1 i} and δ_{n_{1} i} = I (T_{2 i} \leq C_{n_{1} i}) are n_{1} -specific censoring indicators with C_{n_{1} i} =

C_{i} I (n_{1} = 1) + min (T_{1 i}, C_{i}) I (n_{1} = 0)

For the estimand $w_{n_{1}} (\cdot | z)$ , according to equation (8), we formalize $L_{w}$ as

\begin{aligned} \prod_{i = 1}^{m} \prod_{j = 1}^{m} & [{1 - w_{1} (X_{j} | z)}^{I {A_{i} \leq X_{j} \leq (T_{1 i} \land X_{i}), Z_{i} = z}} \\ \times w_{1} (X_{j} | z)^{I {T_{1 i} \leq X_{j} \leq X_{i}, Z_{i} = z}} I (A_{i} \leq X_{j})], \end{aligned}

(8)

Then, for each

X_{j}

, we can estimate

w_{1} (X_{j} | z)

by finding the root

{\hat{w}}_{1} (X_{j} | z_{b})

of equation

\partial \log L_{w} / \partial w_{1} (X_{j} | z) = 0

. It follows that our nonparametric estimators of

w_{n_{1}} (X_{j} | z)

and

n_{1} \in {0, 1}

are

{\hat{w}}_{1} (X_{j} | z) = \frac{\sum_{i = 1}^{m} Y_{1 i}^{N P} (X_{j} | z)}{\sum_{i = 1}^{m} {Y_{1 i}^{N P} (X_{j} | z) + Y_{0 i}^{N P} (X_{j} | z)}},

and

{\hat{w}}_{0} (X_{j} | z) = 1 - {\hat{w}}_{1} (X_{j} | z)

, where

T_{1 i} \land X_{i} = min {T_{1 i}, X_{i}}

4.2. Semiparametric estimator

We further develop model-based semiparametric estimators for $Λ_{n_{1}} (\cdot | z)$ and $w_{n_{1}} (\cdot | z)$ . Given the status of the intermediate event $N_{1} (t^{-}) = n_{1}$ , we propose a Cox proportional hazards model for $Λ_{n_{1}} (t | z)$ :

Λ_{n_{1}} (t | z) = Λ_{n_{1}} (t) \exp (β_{Z n_{1}} z),

(9)

where

Λ_{n_{1}} (t)

is the baseline cumulative hazard function and

β_{Z n_{1}}

is the regression coefficient for

Z

, both conditional on status

N_{1} (t^{-}) = n_{1}

. Then, according to equation (A8) of the Web Appendix and model (9), the conditional log-likelihood can be formulated as

\begin{aligned} \log L_{c} {β_{Z n_{1}}, Λ_{n_{1}} (\cdot)} = \sum_{i = 1}^{m} & [\int_{0}^{\infty} [β_{Z n_{1}} Z_{i} + \log d Λ_{n_{1}} (s)] I {N_{1 i} (s^{-}) = n_{1}} d N_{2 i} (s) \\ - \int_{A_{i} \lor T_{n_{1} i}}^{X_{i}} \exp (β_{Z n_{1}} Z_{i}) I {N_{1 i} (s^{-}) = n_{1}} d Λ_{n_{1}} (s)] \end{aligned}

(10)

Therefore,

d {\hat{Λ}}_{n_{1}} (X_{j})

can be obtained by solving

\frac{\partial \log L_{c} {β_{Z n_{1}}, Λ_{n_{1}} (\cdot)}}{\partial d Λ_{n_{1}} (X_{j})} = 0

, and thus

d {\hat{Λ}}_{n_{1}} (X_{j}) = \frac{\sum_{i = 1}^{m} d N_{2 n_{1} i}^{S P} (X_{j})}{\sum_{i = 1}^{m} Y_{n_{1} i}^{S P} (X_{j}) \exp ({\hat{β}}_{Z n_{1}} Z_{i})} I {\sum_{i = 1}^{m} Y_{n_{1} i}^{S P} (X_{j}) > 0},

where

d N_{2 n_{1} i}^{S P} (s) = d N_{2 i} (s) I {N_{1 i} (s^{-}) = n_{1}}

and

Y_{n_{1} i}^{S P} (s) = I {(A_{i} \lor T_{n_{1} i}) \leq s \leq X_{i}} I {N_{1 i} (s^{-}) = n_{1}}

and

{\hat{β}}_{Z n 1}

is explained in the following. We note that the estimator for the cumulative baseline hazard

{\hat{Λ}}_{n_{1}} (t) = \sum_{j : X_{j} \leq t} d {\hat{Λ}}_{n_{1}} (X_{j})

is the Breslow-type estimator.

{\hat{β}}_{Z n_{1}}

can be obtained by solving the following partial likelihood score:

\sum_{i = 1}^{m} \int_{0}^{\infty} {Z_{i} - \frac{\sum_{l = 1}^{m} z_{l} Y_{n_{1} l}^{S P} (s) \exp (β_{Z n_{1}} z_{l})}{\sum_{l = 1}^{m} Y_{n_{1} l}^{S P} (s) \exp (β_{Z n_{1}} z_{l})}} d N_{2 n_{1} i}^{S P} (s) = 0,

that depends on the status of

N_{1} (t^{-})

For $w_{n_{1}} (\cdot | z)$ , we relate $w_{1} (t | z)$ at a given time $t$ to $z$ using a logistic regression model, where $t = X_{j}$ and $j = 1, \dots, m$ :

logit w_{1} (X_{j} | Z_{i}) = α_{0} (X_{j}) + α (X_{j}) Z_{i} .

(11)

Then, we estimate

{α_{0} (X_{j}), α (X_{j})}

by using the maximum likelihood estimators

{{\hat{α}}_{0} (X_{j}), \hat{α} (X_{j})}

which maximizes

\prod_{i = 1}^{m} \prod_{j = 1}^{m} \frac{\exp [{α_{0} (X_{j}) + α (X_{j}) Z_{i}} I (A_{i} \lor T_{1 i} \leq X_{j} \leq X_{i})]}{{[1 + \exp {α_{0} (X_{j}) + α (X_{j}) Z_{i}}]}^{I (A_{i} \leq X_{j} \leq X_{i})}} .

The model-based semiparametric approach can efficiently adjust for confounding by simply including potential confounders as covariates in the models (9) and (11). We also note that to estimate DE and IE from the counterfactual hazard (3),

w_{n_{1}} (t)

only needs to be estimated at time points with

δ_{j} = 1

4.3. Asymptotic results

Our nonparametric and semiparametric estimators are based on the counting process of the terminal event, which depends on the status of the intermediate event $N_{2 n_{1}}^{N P} (\cdot | z)$ and $N_{2 n_{1}}^{S P} (\cdot)$ . To study the asymptotic properties of the proposed estimators, we first establish the martingality of newly defined processes. Let ${N_{2 n_{1} i} (\cdot), Y_{n_{1}} (\cdot)}$ and $A_{n_{1}} (t)$ denote the generalized notation of the counting processes and compensators, respectively, in both methods: ${N_{2 n_{1}}^{N P} (\cdot | z), Y_{2 n_{1}}^{N P} (\cdot | z)}$ and $\int_{0}^{t} Y_{2 n_{1}}^{N P} (s | z) d Λ_{n_{1}} (s | z)$ for the nonparametric method and ${N_{2 n_{1}}^{S P} (\cdot), Y_{2 n_{1}}^{S P} (\cdot)}$ and $\int_{0}^{t} Y_{2 n_{1}}^{S P} (s) \exp (β_{Z_{n_{1}}} Z) d Λ_{n_{1}} (s)$ for the semiparametric method.

Theorem 1
Under assumptions (1.1), (1.2), and (1.5), $N_{2 n_{1}} (t) - A_{n_{1}} (t)$ is a zero-mean martingale under filtration:
$F_{t} = σ {N_{2 n_{1} i} (s), Y_{n_{1} i} (s^{+}), N_{1 i} (s), Z_{i}; i = 1, \dots, m, n_{1} \in {0, 1}, 0 \leq s \leq t} .$

We simplify the notation by letting $η$ denote the associated parameters in the nonparametric method, $η = {w_{0} (\cdot | z), w_{1} (\cdot | z), Λ_{0} (\cdot | z), Λ_{1} (\cdot | z)}$ or those in the semiparametric method ${α_{0} (\cdot), α_{1} (\cdot), β_{Z_{0}}, Λ_{0} (\cdot), β_{Z_{1}}, Λ_{1} (\cdot)}$ , and $\hat{η}$ and $η^{}$ denote the estimator and the theoretical counterpart, respectively. From the properties of the aforementioned martingale and the result of Andersen and Gill,⁴⁰ one can show that estimators for the counterfactual hazard ${\hat{Λ}}_{n_{1}} (t; z_{a}, z_{b})$ are asymptotically unbiased and consistent for a given $t$ . From propositions II.5.2 and II.5.3 in Andersen et al.,⁴¹ we further prove that unbiasedness and consistency hold uniformly.
Theorem 2
For nonparametric and semiparametric methods, $m \to \infty$ and cardinality of set ${i | T_{1 i} \leq X_{i}}$ goes to infinity,
$sup_{0 \leq t \leq τ} | {\hat{Λ}}_{n_{1}} (t; z_{a}, z_{b}) - Λ_{n_{1}} (t; z_{a}, z_{b}) |$
(12)
converges to zero in probability, and
$E {sup_{0 \leq t \leq τ} | {\hat{Λ}}_{n_{1}} (t; z_{a}, z_{b}) - Λ_{n_{1}} (t; z_{a}, z_{b}) |}$
(13)
converges to zero, where $τ$ is the end time of the study.

Applying the triangle inequality, we can prove that the proposed estimators for DE and IE also enjoy uniform asymptotic unbiasedness: Both
$E {sup_{0 \leq t \leq τ} | Δ_{D E} (t | \hat{η}) - Δ_{D E} (t | η^{}) |} and E {sup_{0 \leq t \leq τ} | Δ_{I E} (t | \hat{η}) - Δ_{I E} (t | η^{}) |}$
converge to 0, and uniform consistency
$sup_{0 \leq t \leq τ} | Δ_{D E} (t | \hat{η}) - Δ_{D E} (t | η^{}) | and sup_{0 \leq t \leq τ} | Δ_{I E} (t | \hat{η}) - Δ_{I E} (t | η^{}) |$
converge to 0 in probability. To study the weak convergence of the proposed estimators, we first establish that of $\hat{η}$ .
Proposition 1
As $m \to \infty$ , $\sqrt{m} {\hat{η} - η^{}}$ converge weakly to a zero-mean Gaussian process denoted by $W$ .

The aforementioned result can be verified by considering the joint likelihood function $L (η)$ as a special case in Zeng and Lin⁴² Because the estimand $Λ (t; z_{a}, z_{b})$ is a function of $η$ , the weak convergence results of $Δ_{D E} (\cdot | \hat{η})$ and $Δ_{I E} (\cdot | \hat{η})$ can be derived using the functional delta method.⁴³
Theorem 3
As $m \to \infty$ and cardinality of set ${i | T_{1 i} \leq X_{i}}$ goes to infinity, $\sqrt{m} {Δ_{D E} (\cdot | \hat{η}) - Δ_{D E} (\cdot | η^{})}$ and $\sqrt{m} {Δ_{I E} (\cdot | \hat{η}) - Δ_{I E} (\cdot | η^{})}$ converge weakly to zero-mean Gaussian processes, denoted by $Δ_{D E η^{}}^{'} (W)$ and $Δ_{I E η^{}}^{'} (W)$ , with variances at time $t$
$\begin{aligned} σ_{D E}^{2} (t) = {\frac{\partial Δ_{D E} (t | η)}{\partial η^{⊤}} cov (\hat{η}) \frac{\partial Δ_{D E} (t | η)}{\partial η} |}_{η = η^{}}; \end{aligned}$

$\begin{aligned} σ_{I E}^{2} (t) = {\frac{\partial Δ_{I E} (t | η)}{\partial η^{⊤}} cov (\hat{η}) \frac{\partial Δ_{I E} (t | η)}{\partial η} |}_{η = η^{}}, \end{aligned}$
respectively, where $Δ_{D E η^{}}^{'}$ and $Δ_{I E η^{}}^{'}$ are the derivatives of $Δ_{D E}$ and $Δ_{I E}$ at $η^{}$ , respectively.

To estimate the pointwise variance in practice, we plug in $\partial Δ_{D E} (t | η) / \partial η |_{η = \hat{η}}$ and $\partial Δ_{I E} (t | η) / \partial η |_{η = \hat{η}}$ , where the forms of $\partial Δ_{D E} (t | η)^{⊤} / \partial η$ and $\partial Δ_{I E} (t | η)^{⊤} / \partial η$ and the limiting expressions for $cov (\hat{η})$ are provided in the Web Appendix.

Theorem 2 also indicates that $Δ_{T E} (\cdot | \hat{η})$ is both asymptotically unbiased and uniformly consistent. Moreover, $\sqrt{m} {Δ_{T E} (\cdot | \hat{η}) - Δ_{T E} (\cdot | η^{})}$ converges weakly to zero mean Gaussian process with variance
$σ_{T E}^{2} (t) = {[\frac{\partial Δ_{D E} (t | η)}{\partial η^{⊤}} + \frac{\partial Δ_{I E} (t | η)}{\partial η^{⊤}}] cov (\hat{η}) [\frac{\partial Δ_{D E} (t | η)}{\partial η} + \frac{\partial Δ_{I E} (t | η)}{\partial η}] |}_{η = η^{}}$
at time $t$ .

By treating $Λ_{0} (X_{i}^{}), Λ_{1} (X_{i}^{}), α_{0} (X_{i}^{})$ and $α_{1} (X_{i}^{})$ for ${i | δ_{i} = 1}$ as parameters. $cov (\hat{η})$ can be obtained by the block diagonal matrix $M = diag (M_{A}, M_{B_{0}}, M_{B_{1}})$ , where $M_{A}$ denotes the covariance matrix of ${{\hat{α}}_{0} (X_{i}^{}), {\hat{α}}_{1} (X_{i}^{})}^{'} s$ . $M_{B_{0}}$ is the covariance matrix of ${{\hat{β}}_{Z_{0}}, {\hat{Λ}}_{0} (X_{i}^{})}^{'} s$ and $M_{B_{1}}$ is the covariance matrix of ${{\hat{β}}_{Z_{1}}, {\hat{Λ}}_{1} (X_{i}^{*})}^{'} s$ (the details of the expression can be found in the Web Appendix).
5. Numerical studies

5.1. Simulation

We conducted simulation experiments to evaluate the finite-sample performance of the proposed methods. Each experiment was conducted based on a total of 1,000 replicates under the sample size $m$ = 500 or 1,000. In particular, the simulation data were generated as follows: the data of $\tilde{Z} = I ({\tilde{Z}}_{C} + {\tilde{E}}_{C} > 0) + 1$ , where ${\tilde{Z}}_{C}$ , the confounder, follow a Bernoulli distribution with $pr ({\tilde{Z}}_{C} = 0) = pr ({\tilde{Z}}_{C} = 1) = 0.5$ and ${\tilde{E}}_{C}$ follow a standard normal distribution; the data of $\tilde{A}$ follow a uniform distribution over the time interval $[0, 100]$ ; censoring time $C$ follows a Weibull distribution with the scale of 2 and shape of 5; the data of ${\tilde{T}}_{1}$ follow a Weibull distribution with the scale $\exp (α_{z} \tilde{Z} + α_{c} {\tilde{Z}}_{C})$ and shape 1, where $α_{z} = 0.5$ and 0 represent the alternative hypothesis $H_{a}$ (with IE) and null hypothesis $H_{0}$ (without IE), respectively, and $α_{c} = 1$ or 0 for the setting with or without the confounding. ${\tilde{T}}_{2} = {\tilde{T}}_{1} + 0.5 \times {\tilde{E}}_{W}$ , where ${\tilde{E}}_{W}$ follows a Weibull distribution with the scale $\exp (β_{z} \tilde{Z} + β_{c} {\tilde{Z}}_{C})$ and shape 1 where $β_{z} = 0.25$ and 0 under $H_{a}$ and $H_{0}$ , respectively, and $β_{c} = 1$ and 0 under settings with and without confounding, respectively. Moreover, we introduced additional simulations to evaluate the performance of the proposed semiparametric approach when the proportional hazards assumption is violated (detailed results illustrated in Figure S3 of the Web Appendix). Specifically, in these simulations, we replaced the Weibull distribution with a log-normal distribution. The new parameters were set with a variance of 1 and means of $- α_{z} \tilde{Z}$ and $- β_{z} \tilde{Z}$ for the scale parameters that were originally used in the Weibull scenarios.

To demonstrate the bias resulting only from the left truncation, but not confounding ( $α_{c} = β_{c} = 0$ ), Figure 2 illustrates that the averages of the estimates from both nonparametric and semiparametric methods, without accounting for left truncation, deviate from the true values if nonzero effects are present. However, no bias exists under the null hypothesis. As depicted in Figure 3 ( $m = 500$ ) and S1 ( $m = 1000$ in the Web Appendix), our proposed methods properly correct the bias caused by left truncation, and the averages of our method’s DE and IE estimates align accurately with the true values. The performance of the proposed variance estimators was evaluated by comparing the pointwise averages of the estimated asymptotic standard errors (NP_asy and SP_asy) at each time point with those of the empirical standard error calculated by averaging across 1,000 replicates (NP_emp and SP_emp). The averages of the proposed variance estimators match well with their empirical counterparts, and the variances in the semiparametric approach are smaller than those in the nonparametric approach (Figures 3 and S1). Moreover, the confidence intervals (CIs) calculated with the proposed variance estimators in nonparametric and semiparametric approaches achieve proper (approximately 95%) coverage rates for both DE and IE (Table 1 [ $m = 500$ ] and ST1 [ $m = 1, 000$ , in the Web Appendix]). On the other hand, as illustrated in Figure S3, the estimates from the proposed semiparametric approach align with the true values under the $H_{0}$ condition. However, under $H_{a}$ , they may introduce potential bias when the proportional hazards assumption fails. To demonstrate the performance of the proposed semiparametric method in adjusting for confounders, we set $α_{c} = β_{c} = 1$ . Figure 4 illustrates that without accounting for the confounder ${\tilde{Z}}_{C}$ , the averages of the estimates deviate from the true values, and the bias can be corrected by adjusting for ${\tilde{Z}}_{C}$ .

Table 1.

CP of proposed methods under sample size $m = 500$ .

	Percentile of the follow-up time.
	H $_{0}$					H $_{a}$
CP ( $%$ )	20 $%$	40 $%$	50 $%$	60 $%$	80 $%$	20 $%$	40 $%$	50 $%$	60 $%$	80 $%$
Direct effect
NP	96.7	96.7	96.2	95.4	96.3	97.4	96.5	96.8	97.1	97.0
SP	96.8	96.3	96.3	96.2	96.0	96.5	97.2	96.9	96.9	97.2
Indirect effect
NP	97.5	97.0	96.1	95.3	93.2	94.7	95.4	95.6	96.2	95.4
SP	97.0	96.8	96.5	95.5	93.7	97.8	97.3	96.9	96.8	95.6

CP: coverage probability; NP: nonparametric method; SP: semiparametric method.

5.2. Data application

In the REVEAL motivation study, data were collected from a community-based prospective cohort of 23,820 individuals recruited during 1991 and 1992. At enrollment, the status of hepatitis B was determined by hepatitis B surface antigen (HBsAg) from the collected serum samples. Studies have argued that most HBsAg-seropositive Taiwanese residents have hepatitis B perinatally before 3 years of age.^44,45 Based on the epidemiological characteristics of Taiwan, we assume that hepatitis B in the REVEAL study occurred at birth.^46,47 The REVEAL study was followed to the end of 2008; during the follow-up, liver cirrhosis and death were ascertained through data linkage with the National Cancer Registry and Death Certification Profile, respectively. In other words, it is highly likely that both the start time of left truncation and the exposure-receiving time are at the same time, that is, birth. It is more biologically meaningful to characterize hepatitis B-induced mortality in relation to the incidence of liver cirrhosis as a function of age. However, using age as the time scale is subject to bias from left truncation, where individuals entering the study at an older age may be healthier than those at a younger age. We used the proposed methods to estimate the DE and IE by age and adjusted for bias from the left truncation. To evaluate the performance of the variance estimators, we calculated the CIs using asymptotic variances and bootstrapping. Two analyses were conducted in this study. First, for a better comparison between the proposed nonparametric and semiparametric methods, we minimized potential confounding by focusing on a more homogeneous subpopulation consisting of 4,895 male participants aged less than 50 years at study entry and without any history of alcohol consumption (analysis I). Second, a total of 12,789 participants aged less than 50 years with confounders were included: Sex (male vs. female) history of alcohol consumption (yes vs. no) and alanine transaminase levels at study entry ( $< 15$ unit/liter vs. $15 - 44$ unit/liter and $\geq 45$ unit/liter) were analyzed using the semiparametric method (analysis II).

As displayed in Figures 5 and 6, both nonparametric and semiparametric methods in analysis I and the semiparametric method in analysis II (setting sex as male and other covariates as their sample means) provide similar results in point estimates. For the CIs, we found that the $95 %$ pointwise intervals calculated with asymptotic variances were similar to those obtained by bootstrapping with 1,000 replicates (constructed by $2.5 %$ and $97.5 %$ quantiles of the estimator for each time point) and that the intervals in the semiparametric method were slightly narrower than those in the nonparametric method. Both methods reveal that hepatitis B increased mortality through liver cirrhosis incidence (particularly after the age of 50 years) and that the effect not through liver cirrhosis was nonsignificant across different ages. The results suggest that most hepatitis B-related mortality occurs due to liver cirrhosis. Such analysis may also have implications for health policies. For example, to reduce preventable deaths among hepatitis B carriers, frequent liver cirrhosis screening after the age of 45 years may be implemented for this high-risk population. Assumption (1.5) is evaluated by comparing the survival curves for those with $N_{1} (t^{-}) = n_{1}$ and $N_{1} (s^{-}) = n_{1}^{*}$ and those with $N_{1} (t^{-}) = n_{1}$ and $s \leq t$ , which supports the plausibility of the assumption (details in the Web Appendix); to evaluate Assumptions (2.1)–(2.3), we conduct sensitivity analyses for the proposed nonparametric and semiparametric methods in the Web Appendix, which supports the validity of our finding. We also applied copula and multistate models to analyze the REVEAL data under analysis I. The Clayton copula model was adopted herein using the estimators of Jiang et al.⁴⁸ Kendall’s tau of 0.895 suggests a high correlation between the event times of liver cirrhosis and mortality, 95 $%$ bootstrap CI $= (0.873, 0.934)$ for those without hepatitis B, and 0.893 $(0.881, 0.920)$ for patients with hepatitis B. The joint survival functions (Figure S4 in the Web Appendix) demonstrate that hepatitis B leads to a relatively lower joint survival probability. For analysis using the multistate model, we analyzed Path A from a healthy state to incidence of liver cirrhosis Path B from a healthy state to mortality, Path C from liver cirrhosis incidence to mortality via the R packages “mstate”.⁴⁹ Multistate modeling suggests that hepatitis B significantly affects the liver cirrhosis incidence from the healthy state (Path A) with a log-hazard ratio of 2.38 (95 $%$ CI $= (2.06, 2.69)$ ); the effects of hepatitis B on Path B (log-hazard ratio $=$ 0.12 (95 $%$ CI $= (- 0.08, 0.33)$ ) and Path C (log-hazard ratio=0.28 (95 $%$ CI $= (- 0.21, 0.76)$ ) are not significant.

6. Discussion

Here, we develop novel approaches to analyze semi-competing risks data in the presence of left truncation by causal mediation modeling. Although existing methods focus more on joint survivor functions (copula models) or transition probabilities between states (multistate models), we aim to characterize DE and IE, which may be more meaningful and scientifically interpretable. For statistical inference, the proposed methods are based on conditional likelihood; we calculate the likelihood of intermediate and terminal events conditioning on the truncation time.

To establish the estimand, the random-draw approach is adopted. The main idea of the random draw approach can be realistically expressed by considering the intervene $\tilde{Z} = z_{b}$ and ${\tilde{N}}_{2} (t) = n_{2} (t)$ in an independent random copy (such as an experiment controlled under similar conditions) of ${\tilde{N}}_{1}^{*} (t)$ to obtain ${\tilde{N}}_{1}^{*} (t; z_{b}, n_{2} (t))$ . Thus, ${\tilde{N}}_{1}^{*} (t; z_{b}, n_{2} (t))$ would have an identical distribution to ${\tilde{N}}_{1} {t; z_{b}, n_{2} (t)}$ but ${\tilde{N}}_{1}^{*} {t; z_{b}, n_{2} (t)} ⊥ {\tilde{N}}_{1} {t; z_{b}, n_{2} (t)}$ . Thus, as illustrated in Figure 1, we borrow the mechanism of ${\tilde{N}}_{1}^{*} {t; z_{b}, n_{2} (t)}$ instead of the original ${\tilde{N}}_{1} {t; z_{b}, n_{2} (t)}$ in ${\tilde{N}}_{2} {t; z_{a}, {\tilde{N}}_{1} {t; z_{b}, n_{2} (t)}, n_{2} (t)}$ , and then the cross-world and $(B_{5})$ assumptions in Huang²² can be escaped.

Two major differences exist between the proposed nonparametric and semiparametric methods. First, adjustment for confounders is more efficient in the semiparametric method. The semiparametric method is model based and allows covariate adjustment in regression models to account for potential confounding factors. For the nonparametric method, confounders can be adjusted by stratification, which, however only applies to categorical confounders. Although stratum-specific effects may be more robust, they lose efficiency and provide no pooled effect. Second, the semiparametric method depends on the proportional hazard assumption and model assumptions, such as linearity for the included covariates; the nonparametric method is free of these assumptions, but its efficiency is traded off for robustness.

Moreover, if time-varying confounders are present, additional research is needed to fully address these complexities when confounders affected by the exposure are as shown in Figure 7(a). However, when the confounder mechanism follows the structure depicted in Figure 7(b), our semiparametric approach, which utilizes the Cox proportional hazards model, can effectively adjust for time-varying exposure-outcome confounders ( $Z ⟵ \tilde{W} (t) ⟶ {\tilde{N}}_{2} (t)$ ) by treating them as covariates. Additionally, this approach employs logistic regression to estimate model parameters at each event time, allowing for the adjustment of time-varying exposure-mediator confounders ( $Z ⟵ \tilde{W} (t) ⟶ {\tilde{N}}_{1} (t)$ ) by incorporating them as covariates.

In the presence of left truncation, the estimation of DE and IE is unbiased if the effects are null, that is, the exposure is independent of the terminal event (Figure 2). This is because the exposure and survival times are independent in the underlying population, and the independence would still be preserved under biased sampling. However, if the effects are not null (i.e., the exposure is associated with the survival time), either through DE or IE, then biased sampling due to left truncation tends to recruit individuals with longer-than-usual survival. Moreover, the distributions of the exposure and mediator (intermediate event) in these survivors may differ from those in the underlying population.

The occurrence of an event, such as birth or diagnosis of a disease, is usually the start time for survival analysis in the presence of left truncation. The event of birth or disease occurrence can be accurately determined. However, causal inference research concerns the effect of exposure or treatment, which makes the initiation of exposure or treatment a natural time of origin. The start time of a treatment is well defined in a clinical trial setting, but it may not be a trivial task to define when exposure to a pathogen or environmental risk factor occurs in observational studies. In our data application, we exploit the epidemiologic characteristics that most hepatitis B carriers are exposed to the virus at birth.^44–47 For observational studies with a less well-defined timing of exposure, careful consideration needs to be made when specifying the time of origin.

The Markov assumption and Assumption (2.1) in our study are related to the assumption required for causal interpretations of the Kaplan–Meier estimates of the time-varying covariate by Sjölander.³⁸ As illustrated in Figure 1, the presence of time-specific variability for ${\tilde{N}}_{2} (\cdot)$ : $e (t_{1})$ and $ε (t_{1})$ for ${\tilde{N}}_{2} (t_{1})$ and $e (t_{2})$ and $ε (t_{2})$ for ${\tilde{N}}_{2} (t_{2})$ , where $e (t_{1})$ and $e (t_{2})$ are independent errors, but $ε (t_{1})$ and $ε (t_{2})$ share a common cause ${\tilde{W}}_{a}$ . Under our Assumption (2.1), ${\tilde{W}}_{a}$ must be adjusted to ensure ${\tilde{N}}_{2} (t^{-}) ⊥ {\tilde{N}}_{2} {t; z, n_{1} (t^{-}), n_{2} (t^{-}) = 0} | {\tilde{W}}_{a}$ for valid causal interpretation. The adjustment of ${\tilde{W}}_{a}$ may also satisfy Assumption (1) for Sjölander³⁸ for causal interpretation. Assumption (2.1) may be weaker than Assumption (1) for Sjölander,³⁸ because we still allow time-specific variability $e (t_{1})$ and $e (t_{2})$ rather than adjusting for all the variability of ${\tilde{T}}_{2}$ .

Supplemental Material

sj-pdf-1-smm-10.1177_09622802241313291 - Supplemental material for Causal mediation analysis for time-to-event mediator and outcome in the presence of left truncation

Supplemental material, sj-pdf-1-smm-10.1177_09622802241313291 for Causal mediation analysis for time-to-event mediator and outcome in the presence of left truncation by Jih-Chang Yu and Yen-Tsung Huang in Statistical Methods in Medical Research

Footnotes

Declaration of conflicting interests

The authors 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: This work was supported by the National Science and Technology Council, Taiwan (grant 112-2118-M-305-005-MY2).

ORCID iDs

Jih-Chang Yu

Yen-Tsung Huang

Supplemental material

Supplemental material for this article is available online.

Appendix A. Prove the equation (1)

References

Chen

Yang

, et al. Risk of hepatocellular carcinoma across a biological gradient of serum hepatitis b virus dna level. J Am Med Assoc 2006; 295: 65–73.

Fine

Jiang

Chappell

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

Day

Bryant

Lefkopoulou

. Adaptation of bivariate frailty models for prediction, with application to biological markers as prognostic indicators. Biometrika 1997; 84: 45–56.

Oakes

. Bivariate survival models induced by frailties. J Am Stat Assoc 1989; 84: 487–493.

Clayton

. A model for association in bivariate life tables and its application to epidemiological studies of familiar tendency in chronic disease epidemiology. Biometrika 1978; 65: 141–151.

Kalbfleisch

Prentice

. The statistical analysis of failure time data. Inc., New Jersy: John Wiley and Sons, 2002.

Wang

. Nonparametric estimation of the sojourn time distributions for a multipath model. J R Stat Soc Ser B 2003; 65: 921–935.

Kalbfleisch

Tai

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

Frangakis

Rubin

. The defining role of principal stratification and effects for comparing treatments adjusted for post-treatment variables: from treatment noncompliance to surrogate endpoints. Biometrics 2002; 58: 191–199.

10.

Zhang

Rubin

. Estimation of causal effects via principal stratification when some outcomes are truncated by “death”. J Educ Behav Stat 2003; 28: 353–368.

11.

Baron

Kenny

. The moderator-mediator variable distinction in social psychological research: conceptual, strategic, and statistical consideration. J Pers Soc Psychol 1986; 51: 1173–1182.

12.

MacKinnon

. Introduction to statistical mediation analysis. NY: Taylor & Francis, 2008.

13.

Rubin

. Bayesian inference for causal effects. Ann Stat 1978; 6: 34–58.

14.

Imai

Keele

Yamamoto

. Identification, inference and sensitivity analysis for causal mediation effects. Stat Sci 2010; 25: 51–71.

15.

Pearl

. Direct and indirect effects. In: Proceedings of the seventeenth conference on uncertainty and artificial intelligence, 2001, pp.411–420. Morgan Kaufmann, San Francisco.

16.

Robins

. Semantics of causal DAG models and the identification of direct and indirect effects. New York, NY: Oxford University Press, 2003.

17.

VanderWeele

Vansteelandt

. Conceptual issues concerning mediation, intervention and composition. Stat Inference 2009; 2: 457–468.

18.

Huang

Cai

. Mediation analysis for survival data using semiparametric probit models. Biometrics 2016; 72: 563–574.

19.

Lange

Hansen

. Direct and indirect effects in a survival context. Epidemiology 2011; 22: 575–581.

20.

Tchetgen Tchetgen

. On causal mediation analysis with a survival outcome. Int J Biostat 2011; 7: 33.

21.

VanderWeele

. Causal mediation analysis with survival data. Epidemiology 2011; 22: 582–585.

22.

Huang

. Causal mediation of semicompeting risks. Biometrics 2021a; 77: 1143–1154.

23.

Fulcher

Shpitser

Didelez

, et al. Discussion on “causal mediation of semicompeting risks” by yen-tsung huang. Biometrics 2021; 77: 1165–1169.

24.

Stensrud

Young

Martinussen

. Discussion on “causal mediation of semicompeting risks” by yen-tsung huang. Biometrics 2021; 77: 1160–1164.

25.

Huang

. Rejoinder to “causal mediation of semicompeting risks”. Biometrics 2021b; 77: 1170–1174.

26.

Lin

VanderWeele

. Interventional approach for path-specific effects. J Causal Inference 2016; 5: 1–10.

27.

Vansteelandt

Daniel

. Interventional effects for mediation analysis with multiple mediators. Epidemiology 2017; 28: 258–265.

28.

Stensrud

Young

Didelez

, et al. Separable effects for causal inference in the presence of competing events. J Am Stat Assoc 2020; 117: 1–9.

29.

Tsai

Jewell

Wang

. A note on the product-limit estimator under right censoring and left truncation. Biometrika 1987; 74: 883–886.

30.

Huang

Ning

Qin

. Semiparametric likelihood inference for left-truncated and right-censored data. Biostatistics 2015; 16: 785–798.

31.

Lynden-Bell

. A method of allowing for known observational selection in small samples applied to 3cr quasars. Mon Not R Astron Soc 1971; 155: 95–118.

32.

Wang

. Nonparametric estimation from cross-sectional survival data. J Am Stat Assoc 1991; 86: 130–143.

33.

Haneuse

Lee

. Semi-competing risks data analysis: accounting for death as a competing risk when the outcome of interest is nonterminal. Circ: Cardiovasc Qual Outcomes 2016; 9: 322–331.

34.

Lee

Gilsanz

Haneuse

. Fitting a shared frailty illness-death model to left-truncated semi-competing risks data to examine the impact of education level on incident dementia. BMC Med Res Methodol 2021; 21: 1–13.

35.

Cheng

Huang

. Combined estimating equation approaches for semiparametric transformation models with length-biased survival data. Biometrics 2014; 70: 608–618.

36.

Peng

Fine

. Nonparametric estimation with left-truncated semicompeting risks data. Biometrika 2006; 93: 367–383.

37.

Lee

Haneuse

Schrag

, et al. Bayesian semi-parametric analysis of semi-competing risks data: investigating hospital readmission after a pancreatic cancer diagnosis. J R Stat Soc Ser C 2015; 64: 253–273.

38.

Sjölander

. A cautionary note on extended kaplan–meier curves for time-varying covariates. Epidemiology 2020; 31: 517–522.

39.

Richardson

Robins

. Single world intervention graphs (swigs): a unification of theh counterfactual and graphical approaches to causality, 2013.

40.

Andersen

Gill

. Cox’s regression model for counting processes: a large sample study. The Ann Stat 1982; 1100–1120.

41.

Andersen

Borgan

Gill

, et al. Statistical models based on counting processes. New York : Springer Science & Business Media, 2012.

42.

Zeng

Lin

. Maximum likelihood estimation in semiparametric regression models with censored data. J R Stat Soc: Ser B (Stat Methodol) 2007; 69: 507–564.

43.

Van der Vaart

. Asymptotic statistics. Cambridge university press, Vol. 3, 2000.

44.

Chen

Wang

. Epidemiology of hepatitis b virus infection in the asia-pacific region. J Gastroenterol Hepatol 2000; 15: E3–E6.

45.

Hsu

Chang

Chen

, et al. Baseline seroepidemiology of hepatitis b virus infection in children in taipei, 1984: a study just before mass hepatitis b vaccination program in taiwan. J Med Virol 1986; 18: 301–307.

46.

Stevens

Beasley

Tsui

, et al. Vertical transmission of hepatitis b antigen in taiwan. N Engl J Med 1975; 292: 771–774.

47.

Stevens

Neurath

Beasley

, et al. Hbeag and anti-hbe detection by radioimmunoassay: correlation with vertical transmission of hepatitis b virus in taiwan. J Med Virol 1979; 3: 237–241.

48.

Jiang

Fine

Chappell

. Semiparametric analysis of survival data with left truncation and dependent right censoring. Biometrics 2005; 61: 567–575.

49.

de Wreede

Fiocco

Putter

, et al. MSTATE: an r package for the analysis of competing risks and multi-state models. J Stat Softw 2011; 38: 1–30.

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

8.37 MB