Targeted maximum likelihood based estimation for longitudinal mediation analysis

2.1 Structural causal model and random interventions

Consider the following structural causal model [13,34], where the randomness of the observed data O is captured by the so-called exogenous nodes U X ’s, and the f X mappings are deterministic for each of the endogenous variables:

For some static treatment A ¯ K = a ¯ K , we define counterfactuals by static interventions on the structural causal model.

For a pair of static treatment and control interventions, a ¯ and a ¯ ′ , denote by Γ t a ¯ ′ the conditional density of the control intervention counterfactual Z t ( a ¯ ′ ) given parent nodes under a ¯ ′ , that is, Γ t a ¯ ′ ( z t ∣ r ¯ t , z ¯ t − 1 , l ¯ t − 1 ) = p ( Z t ( a ¯ ′ ) = z t ∣ R ¯ t ( a ¯ ′ ) = r ¯ t , Z ¯ t − 1 ( a ¯ ′ ) = z ¯ t − 1 , L ¯ t − 1 ( a ¯ ′ ) = l ¯ t − 1 ) . The vector Γ ¯ a ¯ ′ = ( Γ t a ¯ ′ : t = 1 , … , K ) thus represents a sequence of counterfactual conditional densities for the mediator process under the control intervention. We define random intervention counterfactuals under ( a ¯ , Γ ¯ a ¯ ′ ) by forcing mediator variables to follow the distributions of the control intervention counterfactuals, Z t ( a ¯ , Γ ¯ a ¯ ′ ) ∼ Γ t a ¯ ′ , and for nodes X t ∈ R ¯ ∪ L ¯ inserting A ¯ t = a ¯ t so that X t ( a ¯ , Γ ¯ a ¯ ′ ) = f X t ( Pa ( X t ( a ¯ , Γ ¯ a ¯ ′ ) ∣ a ¯ ) , U X ) . That is,

In this framework, a causal target parameter that describes the mediation of treatment effects on an outcome Y is denoted by E [ Y ( a ¯ , Γ ¯ a ¯ ′ ) ] . Then the decomposition of the total effect into a NIE and a NDE is given by

To explain the cross-world assumptions, we also define counterfactuals X ( a ¯ , z ¯ ) with static interventions A ¯ = a ¯ and Z ¯ = z ¯ . Note that if after a random draw under ( a ¯ , Γ ¯ a ¯ ′ ) we have that Z ¯ ( a ¯ , Γ ¯ a ¯ ′ ) = z ¯ , then the structural causal model implies X ( a ¯ , Γ ¯ a ¯ ′ ) = X ( a ¯ , z ¯ ) . However, defining NIE and NDE without random intervention usually requires additional assumptions such as conditional independence between L ( a ¯ , z ¯ ) and Z ( a ¯ ′ ) across the two counterfactual worlds. Such “cross-world” assumptions are not desirable because they are not verifiable, neither empirically nor conceptually [35]. Cross-world assumptions are also incompatible with time-varying covariates R ¯ . On the other hand, the random intervention framework does not require “cross-world” assumptions. We focus on static treatment rules A ¯ = a ¯ or a ¯ ′ , but the methodology can be generalized to dynamic treatment regimes, where the treatment and control interventions on A t are decided by deterministic functions of the available history Pa ( A t ) .

2.2 Right censored survival outcomes

Suppose T is the time point where an event of interest happens. Let Y t = I T > t ∈ L t be the monotone process of staying event free in the study at the t th time point. The counting process Y t starts with 1 at t = 0 and only jumps to 0 if an event happens at the time point t . If an event happens at the time point t , then for all X ∈ Ch ( Y t ) \ Y ¯ , conditional on Y t = 0 , we set X to be a degenerated discrete variable such that X = ∅ with conditional probability 1. The outcome of interest can be the survival beyond the study length K , in which case the target parameters take the form of

In real applications, one typically only observes T ˜ = min { T , C } , where C is the censoring time. To incorporate censored data, we create bivariate treatment nodes A t = ( A t C , A t E ) that consist of the monotone process of remaining uncensored, A t C = I C > t , and the treatment assignment, A t E . The process A t C starts with value 1 and jumps to 0 at the time point where censoring occurs. Conditioning on C t = 0 , for all X ∈ Ch ( C t ) , we set X to be a discrete variable that equals a degenerated value ∅ with conditional probability 1; i.e., no information is available for the observed data likelihood after censoring occurs.

2.3 Competing risks

Consider a competing events framework [36–38], where at the time T where an individual reaches one of several absorbing states the event type Δ ∈ { 1 , 2 , … , J } is observed. For example, Δ = 1 may indicate the onset of a cancer, and Δ = 2 death due to other causes. One can define a multi-dimensional counting process in discrete time, e.g., when J = 2 by { Y t = ( N t ( 1 ) , N t ( 2 ) ) : t = 1 , … , K } , such that N t ( 1 ) = I T ≤ t , Δ = 1 , N t ( 2 ) = I T ≤ t , Δ = 2 . Note that 1 − N t ( 1 ) − N t ( 2 ) = I T > t now indicates that the individual is alive and event-free. Suppose Y t ∈ L t . For Y t = ( 0 , 0 ) , we set X = ∅ for any X ∈ Ch ( Y t ) \ Y ¯ with probability 1. If N t ( j ) = 1 for the j th type of event, we additionally set N t ′ ( j ) = 1 and N t ′ ( j ′ ) = ∅ for all j ′ ≠ j , t ′ ≥ t . The target parameter in a competing risk framework is typically multidimensional and can be the risks of all events under ( a ¯ , Γ ¯ a ¯ ′ ) across all time points, E [ Y ¯ ( a ¯ , Γ ¯ a ¯ ′ ) ] = ( E [ N t ( j ) ( a ¯ , Γ ¯ a ¯ ′ ) ] : t = 1 , … , K ; j = 1 , … , J ) .

3.1 NDE and NIE

We define NDE and NIEs for a multi-dimensional outcome of interest Y ¯ = ψ ¯ ( L ¯ K ) with random interventions as defined in Section 2:

Examples for the choice of Y ¯ include: Y ∈ L K for an univariate outcome variable, see Section 2.1, Y K = I T > K ∈ L K for censored survival endpoints, see Section 2.2, and Y ¯ = ( N t ( 1 ) , N t ( 2 ) : t = 1 , … , K ) for competing risk outcomes, see Section 2.3.

NDE is the direct effect of a treatment while forcing mediators to have the same distribution as their control group counterfactuals. NIE is the indirect treatment effect achieved by not changing treatment values but by changing mediator distributions. The structural causal model of Section 2.1 implies that NDE and NIE decompose the total effect, which can be defined with or without random intervention:

3.2 Identification

To identify the causal mediation targets as statistical parameters of the observed data distribution, we adopt the following assumptions from [8]. For any random intervention ( a ¯ , Γ ¯ a ¯ ′ ) of interest:

• (A1) Sequential exchangeability: R ¯ t K ( a ¯ ′ ) , Z ¯ t K ( a ¯ ′ ) , L ¯ t K ( a ¯ ′ ) , R ¯ t K ( a ¯ , z ¯ ) , L ¯ t K ( a ¯ , z ¯ ) ⊥ A t ∣ Pa ( A t ) .

• (A2) Mediator randomization: R ¯ t + 1 K ( a ¯ , z ¯ ) , L ¯ t K ( a ¯ , z ¯ ) ⊥ Z t ∣ Pa ( Z t ) .

• (A3) (Strong) Positivity: p 0 ( a t * ∣ Pa ( a t * ∣ a ¯ * ) ) > 0 if p 0 ( Pa ( a t * ∣ a ¯ * ) ) > 0 for a ¯ * = a ¯ or a ¯ ′ ; also, p 0 ( r t ∣ Pa ( r t ∣ a ¯ t ′ ) ) > 0 if p 0 ( r t ∣ Pa ( r t ∣ a ¯ t ) ) > 0 , p 0 ( l t ∣ Pa ( l t ∣ a ¯ t ′ ) ) > 0 if p 0 ( l t ∣ Pa ( l t ∣ a ¯ t ) ) > 0 , and p 0 ( z t ∣ Pa ( z t ∣ a ¯ t ) ) > 0 if p 0 ( z t ∣ Pa ( z t ∣ a ¯ t ′ ) ) > 0 .

3.3 Comparison with sequential regression

We have defined and identified our target parameters as functions of the observed data likelihood p , and we will focus on such full likelihood representation in the following sections. In this subsection, we compare with conditional expectation components. This reparametrization can be combined with density ratio estimation techniques [40–43] to achieve efficiency augmentation without estimating conditional densities.

For any pair of interventions a ¯ and a ¯ ′ , one can identify the postintervention distribution P a ¯ , a ¯ ′ of the counterfactual nodes under random interventions ( a ¯ , Γ ¯ a ¯ ′ ) given the assumptions in Section 3.2. That is, for any data realization o that enforces a value a ¯ K in the intervention nodes, we have

For each dimension s of the statistical target parameter, one can rewrite the g-computation formulas with the following sequence of regression expressions (the dependence of Q on P is suppressed except for the last equation):

Note that each of the regression expressions can also be written as a function of conditional densities by analytically carrying out the nested conditional expectations. For example, for simplicity assuming that L k is a discrete variable, we have

The sequential regression formulation thus provides an alternative representation of the g-computation formula of Ψ s a ¯ , a ¯ ′ ( P ) as a functional of

While targeting the same g-computation formula, the plug-in estimation with sequential regression requires only estimating regression models, Q s , X a ¯ , a ¯ ′ ( P ) , of P instead of the factorized data likelihood (1), which leads to a potential gain in scalability when conditional density estimation is challenging. However, this comes at a cost of capturing less information from the data distribution P and only works for one target parameter at a time. When targeting multiple g-computation formulas, such as those defined in longitudinal settings for different outcomes Y s = ψ s ( L ¯ K ) across multiple time points and for different interventions defined by a ¯ , a ¯ ′ , Γ ¯ a ¯ , and Γ ¯ a ¯ ′ , the sequential regression approach has to define and estimate a different set of intermediate regression components Q ¯ s a ¯ , a ¯ ′ ( P ) for each dimension s of the outcome of interest and each random intervention ( a ¯ , Γ a ¯ ′ ) of interest. This potentially variational dependent representation is not intuitive, and it potentially leads to results which do not obey the parameter space; for example, when estimating a survival curve, sequential regression when independently applied to different time horizons may yield survival functions that are increasing. While the resulting estimates still respect the marginal parameter spaces for each dimension of the target parameter due to substitution estimation, sequential regression cannot guarantee the “boundedness” property [21] for the joint parameter space. The target parameter is usually multi-dimensional in longitudinal mediation analysis. For example, when the outcome of interest is the whole survival process Y ¯ = ( Y t = I T > t : t = 1 , … , K ) , the estimations of E [ Y t ( a ¯ , Γ ¯ a ¯ ′ ) ] are expected to be monotonic in t = 1 , 2 , … , K . Another example is the competing risks in Section 2.3, where the estimations also need to respect that E [ Y t ( 1 ) ( a ¯ , Γ ¯ a ¯ ′ ) ] + E [ Y t ( 2 ) ( a ¯ , Γ ¯ a ¯ ′ ) ] = E [ Y t * ( a ¯ , Γ ¯ a ¯ ′ ) ] ∈ [ 0 , 1 ] , increasing in t = 1 , … , K . Our TMLE algorithms on the other hand derive estimates of all target parameters from the same set of estimated likelihood factors and hence achieve the desired joint boundedness property. Finally, sequential regression relies on nested regression models starting at the final time point. While this first regression fit has the least data support (due to deviations from regimes of interest, right censoring, and potentially rare outcomes), challenges and potential biases in estimating the first sequential regression are inherited by subsequent regression fits.