Efficient and flexible mediation analysis with time-varying mediators, treatments, and confounders

1 Introduction

Mediation analyses are paramount to scientific inquiry, as they provide a way of understanding the mechanisms through which effects operate [1]. For example, recent mediation analyses have helped to uncover the types of immune response that coronavirus disease 2019 vaccines trigger in order to prevent disease [2].

Multiple methods have been developed for mediation analysis in the setting of a mediator measured at a single time point, using a counterfactual or potential outcome framework. For example, the natural direct effect is defined as the effect of exposure/treatment that would have been seen in a hypothetical world where the effect operating through the mediator is disabled. Conditions for identifiability for the natural direct effect and its indirect counterpart have been given in [3] and [4].

Though the definition of the natural (in)direct effects is scientifically interesting, their identification requires assumptions about the real world, which are guaranteed to not hold in many practical situations, such as with time-varying confounding (cf. recanting witness, [5]). This restriction limits the applicability of natural direct and indirect effects in practice, as intermediate confounders are expected to be present in many applications, for example, when the effect of an intervention operates through adherence [6]. Several methods have been proposed to do away with the cross-world independence assumption and/or the assumption of no intermediate confounders, for example, partial identification methods [7–9], so-called separable effects [7,10], and effects defined in terms of stochastic interventions on the exposure [11,12]. In this article, we tackle this problem focusing on mediation effects defined in terms of contrasts between counterfactuals in hypothetical worlds in which the treatment is set to some value deterministically, whereas the mediator is drawn stochastically from its counterfactual distribution under interventions on treatment [13–15]. Efficient non-parametric estimators that leverage machine learning to alleviate model misspecification bias have been recently proposed for these parameters [16,17], but they are limited to single time-point exposures. These effects have been called interventional effects [18].

Mediation analysis methodologies for time-varying exposures and mediators are limited. Non-parametric identification formulas for interventional effects in the case of mediators and treatments measured longitudinally, assuming that the time-dependent covariates are measured either before or after the mediator, but not both, are given in refs [18,19]. The authors propose estimation methods that rely on the unlikely ability to correctly specify parametric models for the distribution of the unobservable counterfactual outcomes. Similar interventional effects where the mediator is drawn from its counterfactual distribution conditional on all the past are proposed in ref. [20], together with non-parametric efficient estimators that rely on data-adaptive regression to alleviate model misspecification bias. However, the “direct effect” defined in ref. [20] does not capture the pathway from treatment through intermediate confounder to outcome and is thus not a direct effect in the sense that we are interested in this article. Similar methods for survival analysis and treatment at a single time point have also been proposed [21]. Marginal structural models for longitudinally measured mediators under treatment at a single time point and no loss-to-follow-up are proposed in ref. [22]. More recent work includes the development of a 4-way decomposition of immediate and delayed effects using an infinite time-horizon in a reinforcement learning framework [23], and other methods considering time-varying mediators but a single exposure at baseline [21,24–27].

In this article, we develop a general longitudinal causal mediation approach that fills several gaps in the aforementioned literature. The method we propose satisfies the following: (i) the direct effect is defined in terms of the effects operating through all the pathways that do not include the mediator at any time point, (ii) the methods allow for longitudinally measured mediators and treatments, with confounders possibly measured before and after treatment, (iii) the methods allow the use of data-adaptive regression to alleviate concerns of model misspecification bias, and (iv) the methods are endowed with formulas to construct efficient estimators and computation of valid standard errors and confidence intervals, even under the use of data-adaptive regression. One limitation of prior work remains: our proposed methods can only handle categorical mediators, and the computational complexity increases with the number of categories.

The remainder of this article is organized as follows. In Section 2, we introduce the parameters of interest as well as the identification result and some simple estimators; in Section 4, we discuss efficiency theory for the interventional mediation functional, presenting estimating equations and the efficiency bound in the non-parametric model; in Section 3, we discuss simple estimators based on inverse probability weighting and sequential regression, and in Section 5, we discuss the proposed efficient estimator as well as its asymptotic properties. In Section 6, we present the results of numerical studies, and finally, in Section 7, we present the results of an illustrative study on estimating the longitudinal effect of initiating treatment for opioid-use disorder (OUD) with extended-release naltrexone (XR-NTX) vs buprenorphine-naloxone (BUP-NX) on risk of illicit opioid use during the fourth week of treatment that operates through the mediator of craving symptoms, using symptoms of depression and withdrawal as time-varying covariates. We include weekly measures for each of the first 4 weeks of treatment.

2 Notation and definition of interventional (in)direct effects

Let X 1 , … , X n denote a sample of i.i.d. observations with X = ( L 1 , A 1 , Z 1 , M 1 , L 2 , … , A τ , Z τ , M τ , L τ + 1 = Y ) ∼ P , where A t denotes a vector of intervention variables such as treatment and/or loss-to-follow-up, Z t denotes intermediate confounders, M t denotes a mediator of interest, and L t denotes the time-varying covariates. The outcome of interest is a variable Y = L τ + 1 measured at the end of the study. We let P f = ∫ f ( x ) d P ( x ) for a given function f ( x ) . We use P n to denote the empirical distribution of X 1 , … , X n and assume P is an element of the non-parametric statistical model defined as all continuous densities on X with respect to a dominating measure ν . We let E denote the expectation with respect to P , i.e., E { f ( X ) } = ∫ f ( x ) d P ( x ) . We also let ‖ f ‖ 2 denote the L 2 ( P ) norm ∫ f 2 ( x ) d P ( x ) . We use W ¯ t = ( W 1 , … , W t ) to denote the past history of a variable W , use W ̲ t = ( W t , … , W τ ) to denote the future of a variable, and use H A , t = ( L ¯ t , M ¯ t − 1 , Z ¯ t − 1 , A ¯ t − 1 ) to denote the history of all variables up until just before A t . The random variables H Z , t , H M , t , and H L , t are defined similarly as H Z , t = ( A t , H A , t ) , H M , t = ( Z t , H Z , t ) , and H L , t = ( M t − 1 , H M , t − 1 ) . For the complete history and past of a random variable, we sometimes simplify W ¯ τ and W ̲ 1 as W ¯ . By convention, variables with an index t ≤ 0 are defined as null, expectations conditioning on a null set are marginal, products of the type ∏ t = k k − 1 b t and ∏ t = 0 0 b t are equal to one, and sums of the type ∑ t = k k − 1 b t and ∑ t = 0 0 b t are equal to zero. We let g A , t ( a t ∣ h A , t ) denote the probability mass function of A t conditional on H A , t = h A , t and assume A t takes values on a finite set. The function g M , t ( m t ∣ h M , t ) is defined similarly, and we also assume that M t takes values on a finite set. The variables L t and Z t are allowed to take values on any set; i.e., they can be multivariate, continuous, etc.

We formalize the definition of the causal effects using a non-parametric structural equation model [28]. Specifically, for each time point t , we assume the existence of deterministic functions f A , t , f Z , t , f M , t , and f L , t such that A t = f A , t ( H A , t , U A , t ) , Z t = f Z , t ( H Z , t , U Z , t ) , M t = f M , t ( H M , t , U M , t ) , and L t = f L , t ( H L , t , U L , t ) . Here, U = ( U A , t , U Z , t , U M , t , U L , t , U Y : t ∈ { 1 , … , τ } ) is a vector of exogenous variables, with unrestricted joint distribution. This model can be expressed in the form of a directed acyclic graph (DAG) as in Figure 1. In this article, we are concerned with the definition and estimation of the causal effect of an intervention on A ¯ on Y , as well as its decomposition in terms of the effect through all paths involving the mediator M ¯ versus effects through all other mechanisms that do not involve any component of M ¯ . Mediation effects will be defined in terms of hypothetical interventions where the equations A t = f A , t ( H A , t , U A , t ) and M t = f M , t ( H M , t , U M , t ) are removed from the structural model, and the treatment and mediator nodes are externally assigned as follows. Let M ¯ ( a ¯ ) denote the counterfactual mediator vector observed in a hypothetical world where A ¯ = a ¯ . For a value m ¯ in the range of M ¯ , we also define a counterfactual outcome Y ( a ¯ , m ¯ ) as the outcome in a hypothetical world, where A ¯ = a ¯ and M ¯ = m ¯ . These variables are defined in terms of our assumed NPSEM in the supplementary materials.

Figure 1

DAG. For simplicity, the symbol ↦ is used to indicate arrows from all nodes in its left to all nodes in its right.

Let a ¯ ⋆ = ( a 1 ⋆ , … , a τ ⋆ ) and a ¯ ′ = ( a 1 ′ , … , a τ ′ ) denote two user-specified values in the range of A ¯ , and let G ¯ ( a ¯ ) denote a random draw from the distribution of the counterfactual variable M ¯ ( a ¯ ) conditional on V ⊆ L 1 . We define the conditional interventional effect as E [ Y ( a ¯ ′ , G ¯ ( a ¯ ′ ) ) − Y ( a ¯ ∗ , G ¯ ( a ¯ ∗ ) ) ] and decompose it as

This is the definition of interventional effect in a longitudinal setting given by ref. [18]. In what follows, we focus on the identification and estimation of the parameters E [ Y ( a ¯ ′ , G ( a ¯ ⋆ ) ) ] for fixed a ¯ ′ , a ¯ ⋆ , from which we can obtain the aforementioned effects. In a slight abuse of language, we will also refer to parameters of the type E [ Y ( a ¯ ′ , G ( a ¯ ⋆ ) ) ] as effects.

The aforementioned setup allows for the definition of causal effects for time-to-event outcomes subject to loss-to-follow-up or censoring as follows. Let A t = ( A 1 , t , A 2 , t ) , where A 1 , t denotes the exposure at time t , A 2 , t is equal to one if the unit remains uncensored at time t + 1 and zero otherwise, and Y = L τ + 1 denotes the event-free status at the end of study follow-up. Assume monotone loss-to-follow-up so that A 2 , t = 0 implies A 2 , k = 0 for all k > t , in which case all the data for k > t become degenerate. In this case, we could define the effects as above with a ¯ ′ = ( ( a ¯ 1 , 1 ′ , 1 ) , … , ( a ¯ 1 , τ ′ , 1 ) ) and a ¯ ⋆ = ( ( a ¯ 1 , 1 ⋆ , 1 ) , … , ( a ¯ 1 , τ ⋆ , 1 ) ) , contrasting regimes where treatment at time t is set to A 1 , t = a 1 , t ′ vs A 1 , t = a 1 , t ⋆ while setting censoring status A 2 , t = 1 as “not censored” for everyone. This definition of causal effect for time-to-event outcomes in terms of interventional effects bypasses some (but not all) problems that occur with natural effects due to the fact that a counterfactual longitudinal mediator may be truncated by death, which may render the counterfactual survival time undefined [29]. Furthermore, recent work has uncovered limitations of these effects as a tool to unveil mechanisms [30]. These and other issues related to the interpretation of these direct and indirect effects are discussed at length elsewhere. Interested readers are encouraged to consult the cited articles.

In what follows, we will use the notation H A , t ′ , H Z , t ′ , H M , t ′ , and H L , t ′ to refer to the intervened histories of the variables under an intervention setting A ¯ = a ¯ ′ . For example, H L , t ′ = ( L ¯ t − 1 , M ¯ t − 1 , Z ¯ t − 1 , A ¯ t − 1 = a ¯ t − 1 ′ ) . The histories H A , t ⋆ , H Z , t ⋆ , H M , t ⋆ , and H L , t ⋆ are defined analogously, e.g., H L , t ⋆ = ( L ¯ t − 1 , M ¯ t − 1 , Z ¯ t − 1 , A ¯ t − 1 = a ¯ t − 1 ⋆ ) . The following assumptions will be sufficient to prove the identification of the parameter E [ Y ( a ¯ ′ , G ( a ¯ ⋆ ) ) ] :

Assumptions 1 (i), 1(ii), and 1(iii) together with the assumed DAG in Figure 1 state that H A , t contains all the common causes of A t and Y , H M , t contains all the common causes of M t and Y , and H A , t contains all the common causes of A t and M t , respectively. These are the standard assumptions of no-unmeasured confounders for the treatment-outcome, mediator-outcome, and treatment-mediator relations. Assumptions 2(i) and 2(ii) allow the identification of E [ Y ( a ¯ ′ , m ¯ ) ∣ V = v ] for all m ¯ . Assumption 2(i) states that the value a t ′ has positive probability conditional on values ( l s , z s − 1 : s ≤ t ) that have themselves a positive conditional probability of occurring. Assumption 2(ii) states that the value m t has positive probability conditional on values ( l s , z s : s ≤ t ) that have a positive conditional probability of occurring. Assumption 2(iii) allows us to identify P ( M ¯ ( a ⋆ ) = m ¯ ∣ V = v ) and states that a t ⋆ has a positive probability conditional on values ( l s , m s − 1 , z s − 1 : s ≤ t ) that have positive conditional probability of occurring. Under these assumptions, we have the following identification result.

The aforementioned identification theorem result generalizes several identification formulas in the literature. When τ = 1 , Z = ∅ , and V = L 1 , this identification formula yields (by computing the corresponding contrasts) the identification formula for the natural direct and indirect effects as derived by ref. [4] under an additional cross-world counterfactual assumption. When τ = 1 and V = L 1 , and the confounder Z is present, this identification formula yields the identification formula for the interventional effect described by [15]. If τ > 1 , V = L 1 , and Z t = ∅ for all t , θ in Theorem 1 reduces to Formula (1) in ref. [18]. If τ > 1 and L t = ∅ for t ≤ τ , θ in Theorem 1 reduces to the identification results for Figure 5 given in page 926 of [18]. Our result is a particular case of the identification results of ref. [31].

The choice of variable V in defining the effect should be guided by the total effect to decompose. For instance, if V = L 1 then it may be the case that the total effect in (1) will be closer to the average treatment effect E [ Y ( a ¯ ′ ) − Y ( a ¯ ⋆ ) ] compared to V ⊂ L 1 . In what follows, we assume for simplicity in the presentation that V = ∅ , but the estimators presented can be easily adapted to other cases.

3 Inverse probability weighting and sequential regression

The identification formula in Theorem 1 is a complex functional of the distribution P of the observed data, and it is therefore hard to estimate. Consider, for example, a plug-in estimation strategy where the relevant components of the likelihood are estimated and then plugged in the identifying functional. This would involve finding estimates of the densities of L t and Z t conditional on the past data H L , t and H Z , t . If L t and Z t are continuous or multivariate, or if H L , t and H Z , t are high-dimensional, these densities may be hard to estimate, complicating the construction of a plug-in estimator. Furthermore, the consistency of such an estimator would depend on our ability to consistently estimate the densities involved. In general, this would require the use of flexible estimators that allow automatic selection of non-linearities, interactions, etc. However, there is no statistical theory that allows the study of the sampling distribution (both asymptotic and finite-sample) of such estimator. This results in a lack of a theoretical foundation for the construction of confidence intervals, standard errors, p -values, and other quantities that are often of interest to quantify the uncertainty of statistical estimates.

To address this problem, an estimator based on the theory of doubly robust unbiased transformations [32], targeted minimum loss-based estimation [33,34], and double machine learning [35] was proposed. Concepts in these frameworks will be used to construct estimators that will allow us to use machine learning to address model misspecification bias, while also achieving asymptotic normality, allowing us to construct formulas for computation of standard errors, confidence intervals, and p -values. The asymptotically normal estimators will be constructed in Section 5. This section focuses on describing the two main building blocks: an inverse probability weighted estimator and a sequential regression estimator.

4 Efficiency theory

The foundations of our proposed efficient estimation approach are in semi-parametric estimation theory [e.g., 39–42] and in the theory for doubly robust estimation of causal effects using sequential regression [e.g., 33,34,36,37,43–47]. Central to this theory is the study of the non-parametric efficient influence function (EIF) or canonical gradient, which characterizes the efficiency bound of the longitudinal mediation functional given in Theorem 1 and allows the development of estimators under slow convergence rates for the nuisance parameters involved [41]. Specifically, our estimators will involve finding von-Mises-type approximations for the parameters Q L , t ( m ¯ ) , Q Z , t ( m ¯ ) , and Q M , t ( m ¯ ) , which can be intuitively understood as first-order expansions with second-order error remainder terms. Because the errors in the expansion are second-order, this will mean that the resulting estimator of θ ˆ will be consistent and asymptotically normal at rate n 1 / 2 as long as the second-order error terms converge to zero at rate n 1 / 2 . This convergence rate would be satisfied, for example, if the all regression functions were used for estimation converge at rate n 1 / 4 . We will elaborate on this discussion in Section 5 when we present the asymptotic normality theorems for the proposed estimators.

Let η = { G A , t ′ , G A , t ⋆ , G M , t , Q Z , t , Q M , t , Q L , t : t = 1 , … , τ } . For t = 0 , … , τ , define

Whenever necessary, we make explicit the dependence of these functions on η using notation such as D Z , t ( η ) or D Z , t ( X ̲ t , m ̲ t ; η ) .

The terms R L , t ( η , η ˜ ) , R Z , t ( η , η ˜ ) , and R M , t ( η , η ˜ ) are the second-order error terms involving expectations of products of errors such as ( G ˜ M , s − G M , s ) ( Q ˜ L , s − Q L , s ) , ( G ˜ A , s ′ − G A , s ′ ) ( Q ˜ Z , s − Q Z , s ) , and ( G ˜ A , s ⋆ − G A , s ⋆ ) ( Q ˜ M , s − Q M , s ) , and their explicit form is given in the supplementary materials. This lemma and specifically the second-order form of these remainder terms have important implications in terms of estimation. Specifically, this lemma says that for any value η ˜ , which could represent an inconsistent estimator, regressing D Z , t + 1 ( η ˜ ) on H M , t among units with M t = m t yields a consistent estimator of Q L , t , as long as R L , t ( η , η ˜ ) is small. Because R L , t ( η , η ˜ ) is a sum of products of errors, it may be reasonable to assume that this term is small for data-adaptive regression estimators of η . Specifically, in Section 5, consistency and asymptotic normality of the estimators will require that R L , t ( η , η ˆ ) converges to zero in probability at rate n − 1 / 2 or faster. The second-order term structure of this remainder term means that this fast convergence rate can be achieved under slower convergence rates for each of the components of η ˆ . For example, it will be achievable if all the components of η ˆ converge in probability to their correct limits at rate n − 1 / 4 . This kind of rate is achievable by several data-adaptive regression methods, such as ℓ 1 regularization [48], regression trees [49], neural networks [50], or the highly adaptive lasso [51]. Analogous considerations apply to the estimation of Q Z , t and Q M , t .

An application of the delta method along with Lemma 1 yields the following EIF for the estimation of θ in the non-parametric model (see the supplementary materials for a proof).

This implies that the non-parametric efficiency bound for the estimation of θ is Var [ S ( X , η ) ] and that an efficient estimator of θ will satisfy n ( θ ˆ − θ ) = 1 n ∑ i = 1 n S ( X i , η ) + o P ( 1 ) , licensing the construction of Wald-type confidence intervals and hypothesis tests based on the central limit theorem. In the following section, we will describe an algorithm to construct such an estimator.

5 Efficient estimation using sequential doubly robust regression

Motivated by Lemma 1, we construct a sequentially doubly robust estimator as follows. Consider pseudo-outcomes D ˆ Z , t + 1 , D ˆ L , t , and D ˆ M , t + 1 . The estimation algorithm proceeds sequentially from t = τ , … , 1 performing regression of these pseudo-outcomes in datasets D t + that pool over values of m t . Specifics of this sequentially doubly robust estimator are presented in Algorithm 1.