On Efficient Inference of Causal Effects with Multiple Mediators

On Efficient Inference of Causal Effects with Multiple Mediators

Introduction

Causal inference plays a crucial role in various fields, such as epidemiology hernan2004definition, medicine hernan2000marginal, education card1999causal, and economics panizza2014public. Within this spectrum, Pearl's causal graphical models pearl2000causality,pearl2009causal have recently emerged as a powerful tool for disentangling causal structures among variables (such as confounders, exposure, mediator(s), and outcome). Causal mediation analysis, a core method for examining causal graphical models, aims to reveal the causal mechanisms underlying observed effects from exposure to outcome through the mediator(s), to evaluate the effectiveness of the intervention, and to better understand the roles of mediators pearl2012causal,pearl2014interpretation.

Existing statistical inferential tools for multiple mediators robins1992identifiability,petersen2006estimation,imai2010general,vanderweele2015explanation,chakrabortty2018inference,cai2020anoce,shi2021testing comprise the following three principal steps. Initially, causal structure learning methodologies spirtes2000constructing,chickering2002optimal,nandy2018high,li2019likelihood,yuan2019constrained,li2023inference are applied to estimate the causal graph, often presented by a directed acyclic graph (DAG), using observational data. In the absence of additional assumptions shimizu2006linear, neal2020introduction, the graph is only identified up to a Markov equivalence class (MEC), and a completed partially directed acyclic graph (CPDAG) in such a class is often used to represent the graph structure. The subsequent step is the estimation of the causal effects of mediators based on the DAG or CPDAG obtained from the initial phase. For this task, a variety of estimation techniques have been proposed, including the application of ordinary least squares (OLS) estimators vanderweele2014causal,lin2017interventional,chakrabortty2018inference, parametric models vanderweele2014mediation,vanderweele2016causal,chen2023discovery, and nonparametric methods an2022opening,brand2023recent. The final step is to conduct inferences based on the estimated effects, which often requires finding the exact (asymptotic) distributions of the estimators. As pointed out in chen2023discovery, such an inference is often regarded a separate task and has received less attention in recent causal graph literature.

Although most of the existing work on causal mediation inference is limited to scenarios with a single mediator tchetgen2012semiparametric,tchetgen2013inverse,kennedy2017non,wang2018bounded,xia2023identification, there are some studies that employ the three main standard steps to conduct mediation analysis. However, all of them fall short of comprehensive. Theoretical challenges in unknown causal structures have led to current methods for multiple mediators inference being categorized mainly into three types. One approach assumes that multivariate mediators are conditionally independent given the treatment, or a set of transformed, conditionally independent variables, significantly simplifying the analysis preacher2008asymptotic, boca2014testing, zhang2016estimating, huang2016hypothesis, guo2023statistical, yuan2023confounding. Another category, which does not impose this condition, relies on linear structural equation models (LSEMs) maathuis2009estimating,nandy2017estimating,nandy2018high,chakrabortty2018inference,zhao2022multimodal,zhao2022pathway. The last category allows for a general causal structure and correlated mediators, but uses approximations, such as assuming Gaussian conditional distributions under exposure daniel2015causal,kim2019bayesian,tai2022path, or following a Probit/logistic model for odds ratios vanderweele2014mediation,steen2017flexible,park2018causal. However, these approaches present limitations for complex applications where the causal structure may not be correctly specified.

To bridge this significant gap in addressing potential model misspecification, we consider developing a semiparametric framework to infer causal effects, adapting the general causal structure. Extensive research exists on deriving double robust and highly efficient estimates of the total causal effect of exposure when the model is misspecified scharfstein1999adjusting,bickel2001inference,bang2005doubly. Complementary to this, multiple robust estimators have been developed to quantify direct and indirect effects goetgeluk2008estimation,tchetgen2012semiparametric,chan2016globally,bhattacharya2022semiparametric,xia2023identification. A notable benefit of these multiple robust techniques is their integration of dimension reduction strategies with confounding adjustment, such that the estimators are consistent and asymptotically normal, provided that at least one of the strategies is correct van2006targeted. These methods also achieve semiparametric efficiency when all included strategies are correct van1996weak,bickel2001inference,bang2005doubly. Despite considerable progress in the field, current multiple robust estimators are limited to only a single mediator. Hence, a new semiparametric inference is on demand for inferring causal effects involving multiple interacting mediators under (potentially) unknown causal graphs.

Our Contributions

We conclude our contributions with the following three folds.

The rest of this paper is organized as follows: Section (ref) presents preliminary concepts. Section (ref) formally defines the direct and indirect causal effects of mediators. Section (ref) outlines the semiparametric efficient scores for these causal effects. Section (ref) explores the direct strategy for estimating the causal effects defined in Section (ref), along with an OLS estimation procedure for semi-linear structures. Section (ref) presents alternative estimation strategies, including the introduction of novel quadruply robust estimators. This section also provides both a general algorithm and, under specific conditions, a faster algorithm for computing these estimators. Section (ref) discusses the asymptotic properties of both the OLS and quadruply robust estimators. Section (ref) presents various simulation results that validate the theories proposed for the estimators. In Section (ref), an application of the proposed estimators is used to analyze real data collected from trauma survivors. The glossary of notations, all proofs, and additional technical materials are collected in the Appendix.

Preliminaries

Graph Terminology

Consider a graph $\mathcal{G} =({X}, E)$ with a set of nodes $X$ and a set of edges $E$. There is at most one edge between any pair of nodes. If there is an edge between $X_i$ and $X_j$, then $X_i$ and $X_j$ are adjacent. The node $X_i$ is said to be a parent of $X_j$ if there is a directed edge from $X_i$ to $X_j$. Let the set of all parents of node $X_j$ in $\mathcal{G}$ be $ \operatorname{Pa} (X_j) = \operatorname{Pa}_{\mathcal{G}} (X_j) = \operatorname{Pa}_{X_j} (\mathcal{G})$, and all adjacent nodes of $X_j$ in $\mathcal{G}$ by $\operatorname{adj}(X_j) = \operatorname{adj}_{\mathcal{G}}(X_j)$. A path from $X_i$ to $X_j$ in $\mathcal{G}$ is a sequence of distinct vertices, $\pi := \{a_0, a_1,\cdots,a_L\}\subset V$ such that $a_0 =X_i$, and $a_L=X_j$. A directed path from $X_i$ to $X_j$ is a path between $X_i$ and $X_j$ where all edges are directed towards $X_j$. A directed cycle is formed by the directed path from $X_i$ to $X_j$ together with the directed edge $X_j$ to $X_i$. A directed graph that does not contain directed cycles is called a directed acyclic graph (DAG). A directed graph is acyclic if and only if it has a topological ordering.

Causal Graph Structural Assumption

Let $A$ be a binary exposure/treatment in $\{ 0, 1\}$, $M := (M_1,M_2,\cdots,M_p)^{\top} \in \mathbb{R}^p$ be mediators with dimension $p$ in its support $\mathcal{M} = \mathcal{M}_1 \times \cdots \times \mathcal{M}_p \subseteq \mathbb{R}^p$, and $Y \in \mathbb{R}$ be the outcome of interest. Additionally, we also consider that there are $t - 1$ confounders $C: = (C_1, \ldots, C_{t - 1})^{\top} \in \mathbb{R}^{t - 1}$ in its support $\mathcal{C} \subseteq \mathbb{R}^{t - 1}$. We would just let $t = 1$ here to represent the absence of confounders, that is $C = \varnothing$. Suppose that there exists a DAG $\mathcal{G}=(X, E)$ that characterizes the causal relationship among $X=(C^{\top}, A, M^\top, Y)^\top $, where the dimension of $X$ is $d = t + p + 1$. We suppose we observe i.i.d data on $X = (C^{\top}, A, M^{\top}, Y)^{\top}$ is collected for $n$ subjects. To characterize our model, we consider the following assumptions.

In many instances, the accessible data offers an incomplete view of the inherent causal structure. To address this gap, Causal Markov Condition, Causal Faithfulness Condition, and Causal Sufficiency in the above assumption provide a sufficient condition for causal discovery in i.i.d. data contexts lee2020towards,assaad2022survey,hasan2023a. The rigorous definitions for them and related details can be found in Section 2.4 in hasan2023a. Furthermore, our structural assumptions aim at ensuring the identifiability of the causal model, which are similar to Consistency Assumption and Sequential Ignorability Assumption in tchetgen2012semiparametric, and the structure assumptions in Section 2.4 of chakrabortty2018inference.

Markov Equivalence Class

A general causal DAG, $\mathcal{G}$, may not be identifiable from the distribution of $X$. According to pearl2000causality, a DAG only encodes conditional independence relationships through the concept of $d$-separation. In general, several DAGs can encode the same conditional independence relationships, and such DAGs form a Markov equivalence class. Two DAGs belong to the same Markov equivalence class if and only if they have the same skeleton and the same v-structures kalisch2007estimating. A Markov equivalence class of DAGs can be uniquely represented by a completed partially directed acyclic graph (CPDAG) spirtes2000constructing, which is a graph that can contain both directed and undirected edges. A CPDAG satisfies the following: $X_i \leftrightarrow X_j$ in the CPDAG if the Markov equivalence class contains a DAG including $X_i \rightarrow X_j$, as well as another DAG including $X_j \rightarrow X_i$. CPDAGs can be estimated from observational data using various algorithms, such as the algorithms in kalisch2007estimating, harris2013pc, and zhang2018non. The Markov equivalence class for a fixed CPDAG $\mathcal{C}$ is denoted by $\operatorname{MEC}(\mathcal{C})$, which is a set containing all DAGs $\mathcal{G}$ that have the CPDAG structure $\mathcal{C}$. If we can obtain the true DAG from the data, we can simply treat it as a special case of the "MEC" containing only this DAG, i.e., $\operatorname{MEC}(\mathcal{G})= \{ \mathcal{G} \}$. For simplicity, we denote the corresponding causal structure for the mediators $M$ as $\mathcal{G}_M$, which can be obtained by deleting nodes $C, A, Y$, and the corresponding edges from $\mathcal{G}$. The CPDAG of mediators is similarly denoted as $\mathcal{C}_M$. For simplicity and with a minor stretch of notation, we employ $\mathcal{G}_{M}$ and $\mathcal{C}_M$ to denote the causal DAG and CPDAG of $X$, respectively, such that their corresponding mediators' causal DAG and CPDAG are represented by $\mathcal{G}_{M}$ and $\mathcal{C}_M$ exactly.

Definition of Causal Effects

In this section, we will formally give our refined definition of the causal effects of mediators. To begin with, we give the total effect $TE$, the natural direct effect that is not mediated by mediators $DE$, and the natural indirect effect that is regulated by mediators $IE$ defined in pearl2009causal.

In the above definition, $do(A=0) = do_{\mathcal{G}}(A = 0)$ is a mathematical operator to simulate physical interventions that hold $A$ constant as $0$ while keeping the rest of the model unchanged, which corresponds to remove edges into $A$ and replace $A$ by the constant $0$ in the original causal graph $\mathcal{G}$. Here, $M^{(0)}$ is the (random) value of $M$ if setting $do(A=0)$, and $M^{(1)}$ is the (random) value of $M$ if setting $do(A=1)$. One can refer to pearl2009causal for more details of `do-operator'. The expectation $\mathrm{E} [\cdot]$ is an abbreviation of $\mathrm{E}_P [\cdot]$ with $P = P_X$ is the law of $X$ under $\mathcal{G}$. Inspired by the above definition, we can give the definition of the causal effects for an individual mediator.

The Definition (ref) serves important meanings when we are concerned with different impacts of mediators. Note that $TM_j(\mathcal{G}_M)$ in our definition is an extension for the individual mediation effect in chakrabortty2018inference for LSEMs, denoted as \[ \eta_j (\mathcal{G}_M) = \frac{\partial}{\partial a} \mathrm{E} \big[ Y \mid do(A = a) \big] - \frac{\partial}{\partial a} \mathrm{E} \big[ Y \mid do_{\mathcal{G}_M}(A = a, M_j = m_j)\big], \] which can be interpreted as the change in the total causal effect of the exposure $A$ on the response $Y$ when the potential mediator $M_j$ is removed from the causal graph $\mathcal{G}$ through the intervention $d o\left(M_j=m_j\right)$. But under the non-linearity assumption with binary exposure, the above $\eta_j$ will be a function of $j$-th mediator $m_j$ (Remark 2.1 in chakrabortty2018inference). Therefore, for solving this problem, we take integral with respect to the density $f_{M_j}(m_j)$ for $M_j$, i.e. \[

\] Then by $ a \in \{ 0, 1\}$, we get the expression of $TM_j(\mathcal{G}_M)$ in Definition (ref). Our definition of $DM_j$ is a straightforward extension of Equation (6) in vansteelandt2017interventional by replacing summation with integral. The introduction of $DM_j$ and $IM_j(\mathcal{G}_M) = TM_j(\mathcal{G}_M) - DM_j$ in the above definition is driven by the need for an orthogonal decomposition of the total causal effect of mediators, a concept crucial for unraveling the intricate relationships among variables in mediation analysis, as argued in cai2020anoce. Here $DM_j$ can be interpreted as the causal effect through a particular mediator from the exposure to the outcome, i.e., $\mathrm{E} \big[ Y \mid C, A = 1, M = m\big] \big\{ f_{M_j\mid C, A} (m_j \mid C, A = 1) - f_{M_j\mid C, A} (m_j \mid C, A = 0) \big\}$, that is not regulated by any other mediators, i.e., $f_{M_{-j}\mid C, A}(m_{-j} \mid C, A = 0)$, and thus not regulated by its descendant mediators. Then $IM_j(\mathcal{G}_M) = TM_j(\mathcal{G}_M) - DM_j$ captures the indirect effect of the particular mediator $M_j$ on the outcome $Y$ regulated by descendant mediators.

Definition (ref) is generally defined for any causal graph with binary exposure $A$. In the context of particular linear causal structures, these definitions transform into concise parametric expressions, providing more intuitive `do' representations and aligning consistently with existing literature. A comprehensive discussion on this can be found in Section (ref).

Usually, we do not know the true structure of mediators $\mathcal{G}_{M}$ and we can only estimate its corresponding CPDAG $\mathcal{C}_{ M}$ maathuis2009estimating. If the number of $\operatorname{MEC}(\mathcal{C}_M)$ is larger than one, $TM_j$ and $IM_j$ based on elements $\operatorname{MEC}(\mathcal{C}_{M})$ will be not unique. As a result, we define an identifiable version of $TM_j$ based on a CPDAG as the average over $\operatorname{MEC}(\mathcal{C}_{0, M})$. Specifically,

The corresponding identifiable indirect interventional effect of mediator $M_j$ is $\overline{IM}_j := \overline{TM}_j - DM_j$ for $j \in [p]$.

Semiparametric Efficient Scores

We start with exploring the definition in Section (ref). Define propensity score $e_{a'}(x_S) := \mathrm{P} (A = a' \mid X_S = x_S)$ for $a' \in \{ 0, 1\}$, the outcome mean $\mu_{}(x_S) := \mathrm{E} \big[ Y \mid X_S = x_S\big]$, and the conditional density $\pi_{x_S}(m_T) := f_{M_T \mid X_S}(m_{T} \mid X_{S} = x_{S})$ for any subset $T \subseteq [p]$ and $S \subseteq [t + p + 1]$. Suppose all these functions belong to $\ell^2$-class. Note that in our notation, \( x_S \) and \( m_T \) can represent vectors. At times, we may abbreviate \( \mu \big((x_{S_1}^{\top}, x_{S_2}^{\top}, \ldots, x_{S_N}^{\top})^{\top} \big) \) as \( \mu (x_{S_1}, x_{S_2}, \ldots, x_{S_N}) \) and \( \pi_{(x_{S_1}^{\top}, x_{S_2}^{\top}, \ldots, x_{S_N}^{\top})^{\top}}(m_T) \) as \( \pi_{x_{S_1}, x_{S_2}, \ldots, x_{S_N}}(m_T) \) when \( x_S \) is the concatenation of these vectors, i.e., \( x_S = (x_{S_1}^{\top}, x_{S_2}^{\top}, \ldots, x_{S_N}^{\top})^{\top} \). Let $\operatorname{Pa}_j (\mathcal{G}_M) = \operatorname{Pa}_{\mathcal{G}_M}(M_j) \subseteq \{ M_1, \ldots, M_{j - 1}, \penalty 0 M_{j + 1}, \ldots, M_p\}$ denote the parent mediators of $M_j$ and $\operatorname{pa}_j (\mathcal{G}_M)$ as its realization. Denote

and

for a fixed causal graph $\mathcal{G}_M$, where $\mathcal{M}_{\operatorname{pa}_j} (\mathcal{G}_M)$ is the support of $\operatorname{Pa}_j (\mathcal{G}_M)$ and $a^* = 0$ is the reference level of the exposure. It is worth noting that all three quantities above are random due to the randomness in $C$ (and $M_j$). Given that the exposure features two levels, $0$ and $1$, for simplicity, we use the notation $\langle \cdot \rangle$ to signify the difference evaluated at these two levels. Specifically, let us define

for any function $g_{\boldsymbol{\cdot}, u} (\boldsymbol{\cdot}, v)$, where $u$ and $v$ are arbitrary parameters. Then we can have the following theorem to characterize the relationship between the interventional effects of mediators as specified in Definition (ref) and the quantities defined above.

Consider $\mathscr{M}_{\text {nonpar }}$ as the full model where the observed data likelihood is not constrained, encompassing all conventional laws $P_X$ or, equivalently, distribution $F_X$ of the observed data $X$. The aforementioned theorem establishes that the causal effects detailed in Section (ref) can be represented as regular expectations. Consequently, they function as mappings from $\mathscr{M}_{\text{nonpar}}$ to the real line $\mathbb{R}$. We assume $P_X$ satisfies the positivity assumption given below.

The efficient scores for the functionals $DE$ and $IE$ have been studied in various literature tchetgen2012semiparametric,tchetgen2013inverse,shi2020multiply. The explicit expressions and detailed analysis can be seen in Theorem 1 of tchetgen2012semiparametric. For finding the efficient scores for the functionals $DM_j$ and $IM_j$, we denote

for any $a' \in \{ 0, 1\}$ and $T, S\subseteq [p]$. Then we can derive the efficient score for $\kappa(a', C)$, $\zeta_j (a', a, C)$, and $\varrho_j (a', M_j, C \, ; \, \mathcal{G}_M)$ on $\mathscr{M}_{\text{nonpar}}$ in the following theorem.

In Theorem (ref), we retain the final two terms in $S^{\text{eff}, \text{nonpar}} \big( \mathrm{E} \varrho_j(a', M_j, C , ; , \mathcal{G}_M) \big)$, because we want to express $S^{\text{eff}, \text{nonpar}} \big( \mathrm{E} \varrho_j(a', M_j, C \, ; \, \mathcal{G}_M) \big)$ is composed of four parts. This result can be then combined with Theorem (ref) to obtain the efficient scores for $DM_j$ and $IM_j(\mathcal{G}_M)$, thus $TM_j(\mathcal{G}_M) = DM_j + IM_j(\mathcal{G}_M)$.

In Corollary (ref), the explicit formulas for \( S^{\text{eff}, \text{nonpar}} (DM_j) \), \( S^{\text{eff}, \text{nonpar}} \big( TM_j(\mathcal{G}_M)\big) \), and \( S^{\text{eff}, \text{nonpar}} \big( IM_j(\mathcal{G}_M)\big) \) can be directly derived by substituting the relevant expressions from Theorem (ref). Consequently, the semiparametric efficiency bounds for estimating \( DM_j \), \( IM_j(\mathcal{G}_M) \), and \( TM_j(\mathcal{G}_M) \) within the full nonparametric model \( \mathscr{M}_{\text {nonpar}} \) are respectively \( \mathrm{E} \big[ S^{\text{eff, nonpar}} (DM_j)\big]^2 \), \( \mathrm{E} \big[ S^{\text{eff, nonpar}} (IM_j(\mathcal{G}_M))\big]^2 \), and \( \mathrm{E} \big[ S^{\text{eff, nonpar}} (TM_j(\mathcal{G}_M))\big]^2 \) for any specified \( \mathcal{G}_M \), all with clearly delineated forms. The asymptotic variances of any regular asymptotic linear estimators in $\mathscr{M}_{\text{nonpar}}$ must be greater than or equal to these bounds. Given that $TM_j(\mathcal{G}_M) = DM_j + IM_j(\mathcal{G}_M)$, we will only focus on $DM_j$ and $IM_j(\mathcal{G}_M)$ in the subsequent sections.

Direct Strategy and Ordinary Least Squares (OLS) Estimations

An important implication of Corollary (ref) is that all regular and asymptotically linear (RAL) estimators of $DM_j$ and $IM_j(\mathcal{G}_M)$ in the model $\mathscr{M}_{\text{nonpar}}$ share the common score $S^{\text{eff, nonpar}} (DM_j)$ and $S^{\text{eff, nonpar}} (IM_j(\mathcal{G}_M))$, respectively. For illustrating this and as a motivation for multiply robust estimation when nonparametric methods are not appropriate, we provide a detailed study of different estimating strategies in this section and the next section.

Theorem (ref) gives an explicit expression for $DM_j$ and ${IM}_j(\mathcal{G}_M)$, we can correspondingly give their estimators by (i) replacing the unknown quantities $\kappa(\cdot, \cdot)$, $\zeta_j(\cdot, \cdot, \cdot)$, $\varrho_j(\cdot, \cdot, \cdot \, ; \, \mathcal{G}_M)$ with their estimators and then (ii) replacing $\mathrm{E} [\cdot]$ by $\mathbb{P}_n [\cdot] = n^{-1} \sum_{i = 1}^n [\cdot]_i$ directly. To be specific, in the step (i), we construct the following estimators: \[ \widehat\kappa^{\mathscr{M}_0} (a', c) := \widehat{\mu} (c, a'), \] \[ \widehat\zeta_j^{\mathscr{M}_0}(a', 0, c) := \int_{\mathcal{M}} \widehat{\mu}(C, 1, m) \widehat\pi_{C, a'}(m_j) \widehat\pi_{C, 0}(m_{-j}) \, \mathrm{d} m, \] and \[ \widehat\varrho_j^{\mathscr{M}_0} (a', m_j, c \, ; \, \mathcal{G}_M) := \int_{ \mathcal{M}_{\operatorname{pa}_j} (\mathcal{G}_M)} \widehat\mu \big(C, a', \operatorname{pa}_j (\mathcal{G}_M), M_j \big) \widehat\pi_{C, a'} \big(\operatorname{pa}_j(\mathcal{G}_M) \big) \, \mathrm{d} \operatorname{pa}_j (\mathcal{G}_M), \] which are consistent for $\kappa(a', C)$, $\zeta_j(a', 0, c)$ and $\varrho_j(a', m_j, c)$ for any $c \in \mathcal{C}$, $a' \in \{ 0,1 \}$, and $m \in \mathcal{M}$. Note that $\kappa(a', C)$ can be written as \[ \mathrm{E} [Y \mid C, A = a'] = \int_{\mathcal{M}} \mathrm{E} [Y \mid C, A = a', M = m] f_{M \mid C, A}(m \mid C, A = a') \, \mathrm{d} m. \] Therefore, the consistency of $\widehat\kappa^{\mathscr{M}_0} (a', c)$, $\widehat\zeta_j^{\mathscr{M}_0}(a', 0, c)$, and $\widehat\varrho_j^{\mathscr{M}_0} (a', m_j, c \, ; \, \mathcal{G}_M)$ will use the correctly specified information as follows:

Then in the step (ii), we can construct the estimators \[ \widehat{DM}_j^{\mathscr{M}_0} = \mathbb{P}_n \Big[ \big\langle \widehat\zeta_j^{\mathscr{M}_0} (\boldsymbol{\cdot}, 0, C) \big\rangle \Big], \] \[ \widehat{TM}_j^{\mathscr{M}_0} (\mathcal{G}_M) =\mathbb{P}_n \Big[ \big\langle \widehat\kappa^{\mathscr{M}_0} (\boldsymbol{\cdot}, C) - \widehat\varrho_j^{\mathscr{M}_0} (\boldsymbol{\cdot}, M_j, C \, ; \, \mathcal{G}_M) \big\rangle \Big] \] and $\widehat{IM}_j^{\mathscr{M}_0} (\mathcal{G}_M) = \widehat{TM}_j^{\mathscr{M}_0} (\mathcal{G}_M) - \widehat{DM}_j^{\mathscr{M}_0}$. Next, the estimators for the identifiable $\overline{TM}_j$ and $\overline{IM}_j$ are \[ \widehat{TM}_j^{\text{avg}, \mathscr{M}_0 } = \frac{1}{\# \operatorname{MEC}(\widehat{\mathcal{C}}_{ M})} \sum_{\mathcal{G}_{M} \in \operatorname{MEC}(\widehat{\mathcal{C}}_{ M})} \widehat{TM}_j^{\mathscr{M}_0} (\mathcal{G}_M), \] and $\widehat{IM}_j^{\text{avg}, \mathscr{M}_0} = \widehat{TM}_j^{\text{avg}, \mathscr{M}_0} - \widehat{DM}_j^{\mathscr{M}_0}$, where $\widehat{\mathcal{C}}_{M} \in \{ 0, 1\}^{p \times p}$ is the adjacency matrix of the estimated CPDAG for the mediators $M = (M_1, \ldots, M_p)^{\top} \in \mathbb{R}^p$. The consistent causal structure $\widehat{\mathcal{C}}_M$ can be obtained by after we obtain the estimated adjacency matrix $\widehat{\mathcal{C}}$ for the whole causal graph, and we then extract a subset $\widehat{\mathcal{C}}_M = \big[ \widehat{\mathcal{C}}_{kk} \big]_{k \in { t + 1, \ldots, t + p}}$ to arrive at the causal structure for the mediators. The estimated adjacency matrix $\widehat{\mathcal{C}}$ can be achieved through methods such as the PC algorithm spirtes2000constructing, greedy equivalence search (GES) chickering2002optimal, and adaptively restricted greedy equivalence search (ARGES) nandy2018high, among others.

OLS estimator under semi-linear model

The direct strategy in model $\mathscr{M}_0$ involves two unknown quantities: $\mathrm{E} [Y \mid C = \boldsymbol{\cdot}, A = \boldsymbol{\cdot}, M = \boldsymbol{\cdot}]$ and $f_{M \mid C, A}(m \mid C = \boldsymbol{\cdot}, A = \boldsymbol{\cdot})$. Specially, when we have known that the structure follows $Y \, \leftarrow \, h_Y(C, A, M) + \epsilon_{Y}$ and $M \, \leftarrow \, h_M(C, A, M) + \epsilon_{M}$\footnote{The symbol $\leftarrow$ emphasizes that the expressions should be understood as a generating mechanism rather than as a mere equation.} with given functions $h_Y(\boldsymbol{\cdot})$ and $h_M(\boldsymbol{\cdot})$, and given that the error terms $\epsilon_{Y}$ and $\epsilon_M$ belong to some classes of distributions, $\mathscr{M}_0$ will be correctly recovered. A commonly used approach for this is assuming Linear Structural Equation Models (LSEMs), i.e., $X \leftarrow B^{\top} X + \epsilon$ with mean-zero and jointly independent error vector $\epsilon$, where $B=(b_{ij})_{1\leq i\leq d,1\leq j\leq d}$ be a $d\times d$ matrix, where $b_{ij}$ is the weight of the edge $X_i\rightarrow X_j \in E$, and $b_{ij}=0$ otherwise. There are numerous rigorous theoretical findings for LSEMs, as discussed in chakrabortty2018inference,cai2020anoce,shi2022testing. However, LSEMs are not applicable when dealing with binary exposure, given that the element $A$ in $X$ is constrained to either $0$ or $1$. In lieu of LSEMs, we propose the following semi-linear structure assumption.

Through this paper, $\alpha$, $\beta$, and $B$ will always represent a scalar, vector, and matrix, respectively. Define $\theta_{MA} := \big[ (I_p - B_{MM}^{\top})^{-1} \beta_{MA} \big]$ where $I_p \in \mathbb{R}^{p \times p}$ is the identity matrix, then under the above semi-linear structural assumptions, we have the following propositions for the uniqueness of $\theta_{MA}$ under MEC, interpretation displays, and neat {parametric} expressions for the causal effects defined in Section (ref).

Proposition (ref) gives the fact that only $IM_j$ and $TM_j$ require specific DAG structure, while other quantities do not require any knowledge of the causal structure under semi-linear assumption. Meanwhile, Proposition (ref) implies that, under the semi-linear assumption, our definitions for direct/indirect individual mediation effects in Definition (ref) exactly coincides with the definitions in cai2020anoce: the first multiplier is in Proposition (ref) with the classical meaning of `natural' in the causal inference literature pearl2000causality. Thus, $DM_j$ can be interpreted as the causal effect through a particular mediator from the treatment on the outcome that is not regulated by its descendant mediators. Similarly, by the first multiplier in the $IM_j$, we know that $IM_j$ captures the indirect effect of a particular mediator on the outcome regulated by its descendant mediators.

More importantly, Proposition (ref) can imply a simply OLS estimator for the direct strategy together with Proposition (ref) as long as the sample size $n$ is larger than the dimension $d$. Indeed, we can rewrite the part of semi-linear structure (ref) as follows:

Write $\widehat{\theta}_{MA}$ as the OLS estimator of unknown parameter ${\theta}_{MA}$, similarly define the other corresponding estimated quantities as follows: \[ \left [

\right] \qquad and \qquad \left[

\right]. \] Then we will have OLS estimators for the direct strategy estimators: For $DE$ and $IE$, $\widehat{DE}^{\text{OLS}} = \widehat{\alpha}_{YA}$ and $\widehat{IE}^{\text{OLS}} = \widehat{\beta}_{YM}^{\top} \widehat\theta_{MA}$; For $DM_j$ and $IM_j$, $\widehat{DM}_j^{\text{OLS}} = \widehat{\beta}_{YM, j}\widehat\theta_{MA, j}$, $\widehat{{IM}}_j^{\text{OLS}} (\mathcal{G}_{M}) = \widehat{\theta}_{MA, j} \big\{ \widehat{\mathrm{E}}^{\text{reg}} \big[ Y \mid M_j \cup \operatorname{Pa}_j (\mathcal{G}_{M}) \cup A \cup C \big]_1 - \widehat{\beta}_{YM, j} \big\}$, and \[ \widehat{{IM}}_j^{\text{avg}, \, \text{OLS}} = \frac{1}{\# \operatorname{MEC}(\widehat{\mathcal{C}}_{ M})} \sum_{\mathcal{G}_{M} \in \operatorname{MEC}(\widehat{\mathcal{C}}_{ M})} \widehat{{IM}}_j^{\text{OLS}} (\mathcal{G}_{M}), \] where $\widehat\mathrm{E}^{\text{reg}} [X_k \mid X_l \, \cup \, X_S]_1$ is the estimated coefficient for $X_l$ obtained from the linear regression of $X_k$ on $X_l \, \cup \, X_S$, as determined from the data. Thus, when the semi-linear structure is determined, we can simplify direct strategy estimators to OLS estimators. All the these OLS estimators can be easily obtained by just applying simple regressions with nice properties, we will discuss their asymptotic properties in Section (ref).

Multiple Robust Estimators

Several Alternative Strategies

For a fixed $j \in [p]$, beyond the direct strategy above, there are alternative identification formulas for $\mathrm{E} \kappa(a', C)$, $\mathrm{E} \zeta_j (a', 0, C)$, and $\mathrm{E} \varrho_j (a', M_j, C \, ; \, \mathcal{G}_M)$. Based on these formulations, we can derive the corresponding estimators. We will discuss them one by one in the subsequent sections.

Alternative Strategy 1

The first one is using propensity score to construct the inverse probability weighting estimator. Note that we have\footnote{The calculation details are shown in (ref).} \[ \mathrm{E} \bigg[ \frac{\mathds{1}(A = a')}{e_{a'}(C)} Y\bigg] = \mathrm{E} \kappa(a', C), \]

and

Thus, corresponding estimators take the form \[ \widehat\kappa^{\mathscr{M}_1} (a', C) = \frac{\mathds{1}(A = a')}{\widehat{e}_{a'} (C)} Y, \] \[

\] and \[ \widehat\varrho_j^{\mathscr{M}_1} (a', M_j, C \, ; \, \mathcal{G}_M) = \frac{\mathds{1}(A = a')}{\widehat{e}_{a'}(C)} \widehat\pi_{C, a'} \big(\operatorname{Pa}_j (\mathcal{G}_M) \big) Y, \] respectively. Here, the propensity scores $e_{a'}(\boldsymbol{\cdot})$ and conditional densities $\pi_{\boldsymbol{\cdot}}(\boldsymbol{\cdot})$ appearing in $\widehat\kappa^{\mathscr{M}_1} (a', C)$, $\widehat\zeta_j^{\mathscr{M}_1} (a', 0, C)$, and $\widehat{\varrho}_j^{\mathscr{M}_1} (a', M_j, C \, ; \, \mathcal{G}_M)$ should be correctly estimated in the collection of quantities $\mathscr{M}_{j, \, 1}$ such that

Define

then, we can construct the estimators under $\mathscr{M}_{j, \, 1}$ is $\widehat{DM}_j^{\mathscr{M}_1} = \mathbb{P}_n \Big[ \big\langle \widehat\zeta_j^{\mathscr{M}_1} (\boldsymbol{\cdot}, 0, C) \big\rangle \Big]$, \[

\] and $\widehat{IM}_j^{\text{avg}, \mathscr{M}_1} = \widehat{TM}_j^{\text{avg}, \mathscr{M}_1} - \widehat{DM}_j^{\mathscr{M}_1}$, where $\widehat{\mathcal{C}}_M$ is the estimated adjacency matrix consistent with the true $\mathcal{C}_M$. Here, the superscript $\mathscr{M}_1$ associated with these estimators signifies that their consistency relies on the correct specification of information in $\mathscr{M}_{j, \, 1}$. For clarity and where there is no risk of confusion, we will also employ $\mathscr{M}_1$ to represent the estimation methodology behind these estimators. In the subsequent two subsections, the notations $\mathscr{M}_{2}$ and $\mathscr{M}_{3}$ bear analogous meanings.

Alternative Strategy 2

Similarly, we can verify that

and

Thus, corresponding estimators take the forms \[

\] and \[ \widehat{\varrho}_j^{\mathscr{M}_2} (a', M_j, C \, ; \, \mathcal{G}_M) = \frac{\mathds{1}(A = a')}{\widehat{e}_{a'} (C)} \int_{\mathcal{M}_j}\widehat{\mu}(C, a', \operatorname{Pa}_j(\mathcal{G}_M), m_j) \widehat{\pi}_C(m_j) \, \mathrm{d} m_j \] with estimators $\widehat{e}_{\boldsymbol{\cdot}} (\boldsymbol{\cdot})$, $\widehat\mu(\boldsymbol{\cdot})$, and $\widehat\pi_{\boldsymbol{\cdot}} (\boldsymbol{\cdot})$ appear in $\widehat\kappa^{\mathscr{M}_1} (a', C)$, $\widehat\zeta_j^{\mathscr{M}_2} (a', 0, C)$, and $\widehat{\varrho}_j^{\mathscr{M}_2} (a', M_j, C \penalty 0 \, ; \, \mathcal{G}_M)$. They use the information in $\mathscr{M}_{j, 2}$ such that

Here we notice the fact that $\pi_{c}(m_j) = \pi_{c, 0}(m_j) + \pi_{c, 1}(m_j)$, and thus, $\widehat{\kappa}^{\mathscr{M}_1}(a', C)$ will also be consistent in $\mathscr{M}_{j, \, 2}$. Then $\widehat{DM}_j^{\mathscr{M}_2} = \mathbb{P}_n \big[ \big\langle \widehat\zeta_j^{\mathscr{M}_2} (\boldsymbol{\cdot}, 0, C) \big\rangle \big]$, \[

\] and \( \widehat{IM}_j^{\text{avg}, \, \mathscr{M}_2} = \widehat{TM}_j^{\text{avg}, \, \mathscr{M}_2} - \widehat{DM}_j^{\mathscr{M}_2} \) are consistent provided that the estimated adjacency matrix \( \widehat{\mathcal{C}}_M \) is consistent to \( \mathcal{C}_{M} \). \( \overline{\mathbb{P}_n} \widehat\varrho_{j \, ; \, 1- 0}^{\mathscr{M}_2} (M_j, C \, ; \, \widehat{\mathcal{C}}_M) \) is similarly defined by substituting \( \mathscr{M}_1 \) with \( \mathscr{M}_2 \) as detailed in (ref).

Alternative Strategy 3

The last strategy is based on the third representation of the functional as follows:

and

Similarly, we can consider the estimators \[ \widehat{\zeta}_j^{\mathscr{M}_3}(a', 0, C) = \frac{\mathds{1}(A = a')}{\widehat{e}_{a'}(C)} \int_{\mathcal{M}_{-j}} \widehat{\mu} (C, 1, M_{j}, m_{-j}) \widehat\pi_{C, 0} (m_{-j})\, \mathrm{d} m_{-j}, \] and \[

\] Thus, our estimators under the third identification formulas can be written as $\widehat{DM}_j^{\mathscr{M}_3} = \mathbb{P}_n \big[ \widehat\zeta_j^{\mathscr{M}_3} (1, 0, C) - \widehat\zeta_j^{\mathscr{M}_3} (0, 0, C) \big]$, \[

\] and $\widehat{IM}_j^{\text{avg}, \, \mathscr{M}_3} = \widehat{TM}_j^{\text{avg}, \, \mathscr{M}_3} - \widehat{DM}_j^{\mathscr{M}_3}$, where the estimated adjacency matrix $\widehat{\mathcal{C}}_M$ is consistent to $\mathcal{C}_M$, and \( \overline{\mathbb{P}_n} \widehat\varrho_{j \, ; \, 1- 0}^{\mathscr{M}_3} (M_j, C \, ; \, \widehat{\mathcal{C}}_M) \) is by replacing \( \mathscr{M}_1 \) with \( \mathscr{M}_3 \) in (ref). The estimators $\widehat{e}_{\boldsymbol{\cdot}} (\boldsymbol{\cdot})$, $\widehat\mu(\boldsymbol{\cdot})$, and $\widehat\pi_{\boldsymbol{\cdot}} (\boldsymbol{\cdot})$ use the following information:

Here we note the fact that $\widehat{\kappa}^{\mathscr{M}_1}(a', C)$ will be consistent in $\mathscr{M}_{j, \, 3}$ again.

Quadruply Robust Estimator

Denote $\mathscr{M}_{j, \, \text{union}} := \mathscr{M}_0 \, \cup \, \mathscr{M}_{j, \, 1} \, \cup \, \mathscr{M}_{j, \, 2} \, \cup \, \mathscr{M}_{j, \, 3}$, then $\cup_{j = 1}^p \, \mathscr{M}_{j, \, \text{union}} \subsetneq \mathscr{M}_{\text{nonpar}}$, and $\widehat{DM}_j^{\mathscr{M}_{\ell}}$ and $\widehat{IM}_j^{\mathscr{M}_{ \ell}} (\mathcal{G}_M)$ are all mapping the estimated distribution $\widehat{F}_X$ to the true $DM_j$ and ${IM}_j (\mathcal{G}_M)$ defined in Definition (ref) for $\ell = 0, 1, 2, 3$ and any fixed $\mathcal{G}_M$, since all these representations agree on the nonparametric model $\mathscr{M}_{\text{nonpar}}$. Therefore, we may conclude that both direct strategy and alternative strategies are in fact asymptotically efficient in $\mathscr{M}_{\text{nonpar}}$ with common scores $S^{\text {eff, nonpar }}(DM_j)$ and $S^{\text {eff, nonpar }} \big(IM_j (\mathcal{G}_M)\big)$. Furthermore, from this observation, one further concludes that (asymptotic) inferences obtained using one of the four representations are identical to inferences using either of the other three representations for a fixed $\mathcal{G}_M$. However, to achieve this, each strategy need exactly correctly specified for the conditional expectation and conditional density, i.e., correctly specified for the corresponding collections $\mathscr{M}_{j, \ell}$ in Section (ref), (ref), (ref), and (ref), where we denote $\mathscr{M}_{j, \, 0} \equiv \mathscr{M}_0$ for each $j \in [p]$. In general, $\widehat{DM}_j^{\mathscr{M}_{\ell}}$, $\widehat{IM}_j^{\mathscr{M}_{\ell}} (\mathcal{G}_M)$ fail to be consistent outside of the corresponding submodel $\mathscr{M}_{j, \, \ell}$ for each $\ell \in \{ 0, 1, 2, 3\}$.

Note that the alternative strategy 1 in Section (ref) in $\mathscr{M}_{j, \, 1}$ induces Inverse Probability Weighted (IPW) estimator. A commonly-used method is combining the direct strategy estimator in Section (ref) correctly specified with the model $\mathscr{M}_0$ and IPW estimator in (ref) with the model $\mathscr{M}_0$, and getting the double robust estimator. But the double robust estimators only combine two estimation strategies, $\mathscr{M}_0$ and $\mathscr{M}_{j, \, 1}$, and ignore use other two alternative strategies. Hence, the double robust estimator may be inconsistent outside of $\mathscr{M}_0 \, \cup \, \mathscr{M}_{j, 1}$. To overcome this problem, we propose an approach that produces a quadruply robust estimator by combining the above all four strategies as follows:

where $\widehat{S}^{\text{eff, nonpar}} (\boldsymbol{\cdot})$ is equal to ${S}^{\text{eff, nonpar}} (\boldsymbol{\cdot})$ evaluated at the given consistent estimators $\widehat{e}_{a'}(\boldsymbol{\cdot})$, $\widehat{\pi}_{\boldsymbol{\cdot}}(\boldsymbol{\cdot})$, and $\widehat{\mu}(\boldsymbol{\cdot})$ for all propensity scores, the conditional densities, and the conditional expectations appearing in ${S}^{\text{eff, nonpar}} (\boldsymbol{\cdot})$. Denote \[ \widehat{\tau}_{\boldsymbol{\cdot} \, ; \, S} (C, a', M_{T}) := \int_{\mathcal{M}_{-T}} \widehat\mu (C, a', m_{-T}, M_{T}) \widehat\pi_{\boldsymbol{\cdot}} (m_{-T}) \, \mathrm{d} m_{-T}, \] as the corresponding estimator for ${\tau}_{\boldsymbol{\cdot} \, ; \, S} (C, a', M_{T})$ defined in (ref), then we have the following explicit expressions for the quadruply estimators as

and $\widehat{IM}_j^{\text{QR}} (\mathcal{G}_M) = \widehat{TM}_j^{\text{QR}} (\mathcal{G}_M) - \widehat{DM}_j^{\text{QR}}$. Then the quadruply estimator for indirect interventional mediation effect with a consistent estimated $\widehat{\mathcal{C}}_M$ is defined as $\widehat{IM}_j^{\text{avg}, \, \text{QR}} = \frac{1}{\# \operatorname{MEC}(\widehat{\mathcal{C}}_{M})} \sum_{\mathcal{G}_{M} \in \operatorname{MEC}(\widehat{\mathcal{C}}_{ M})} \penalty 0 \widehat{IM}_j^{\text{QR}} (\mathcal{G}_M)$. Compared to double robust estimators, our novel quadruply robust estimators can tolerate a higher degree of misspecification outside of $\mathscr{M}_0 \, \cup \, \mathscr{M}_{j, \, 1}$ and still achieve consistency. We will see this in Section (ref).

Subject to some mild regularity conditions, delineated in Section (ref), our quadruply estimators are asymptotic normal and efficient. Thus, based on the semiparametric efficient scores, we get the score-based variance estimators for $\widehat{DM}_j^{\text{QR}}$ and $\widehat{IM}_j^{\text{avg}, \, \text{QR}}$ as \[ \widehat\operatorname{var} ( \widehat{DM}_j^{\text{QR}} ) := \frac{1}{n^2} \sum_{i = 1}^n \left[ \widehat{S}^{\text{eff, nonpar}} (\widehat{DM}_j^{\text{QR}}) - \widehat{DM}_j^{\text{QR}}\right]^2 \] and \[

\] correspondingly. However, in practical scenarios, confidence intervals (CIs) derived using the Wald-type method, especially when grounded on score-based variance estimators, tend to be more narrow boos2013essential. This can potentially result in anti-conservatism. To achieve more concise statistical inference for our quadruply estimators, we consider utilizing the variances derived from the symmetric $t$-bootstrap approach hall1988symmetric here. A pseudocode summarizing the proposed algorithm for these quadruply estimators and their bootstrap CIs is given in Algorithm (ref). The $\log n$ truncations in Algorithm (ref) aims to achieve the numerical stability, which is a technique widely recognized in statistical literature heckman1976common,sun2020adaptive,chinot2020robust.

Practical fast implement

The formulas for the quadruply robust estimators, as shown in equations (ref) and (ref), require several numerical integrals for each $i \in [n]$, which may be computationally demanding. To address this challenge, we purpose Algorithm (ref) in the above section, in which we employ the Monte Carlo method to evaluate these integrals. However, when the data partly satisfy the semi-linear structure and both $\epsilon_M$ and $\epsilon_Y$ adhere to a mean-zero Gaussian distribution, explicit expressions for these numerical integrals can be derived, facilitating faster computation. Indeed, if we assume the linear structure in $M \, \leftarrow \, C \oplus A \oplus M$ and denote the density (or mass) function of $\epsilon_M = (\epsilon_{M, 1}, \ldots, \epsilon_{M, p})^{\top}$ as $f(x) = f(x_1, \ldots, x_p)$, then the conditional density of $M$ given $C$ and $A = 1$ is $f\big(x - {\Theta}_{MC} C - \theta_{MA} \big)$ from (ref). This allows us to compute

Similarly, we can derive explicit expressions for some other integrals in equations (ref) and (ref) as long as the linear structure in $M \, \leftarrow \, C \oplus A \oplus M$ holds. One step more, when \( \epsilon_M \) is a mean-zero Gaussian distribution, any integral in (ref) and (ref) will have an explicit expression. This leads to a more efficient implementation of (ref) and (ref). The following Algorithm (ref) and Proposition (ref) elaborates on this.

In practical scenarios where the sample size \( n \) is sufficiently large, it becomes reasonable to treat the sample means $\overline{\epsilon}_M := \frac{1}{n} \sum_{i = 1}^n \epsilon_{M, i}$ and $\overline{\epsilon}_Y := \frac{1}{n} \sum_{i = 1}^n \epsilon_{Y, i}$ as if they follow mean-zero Gaussian distributions. This permits the utilization of Proposition (ref), particularly when empirical evidence can support the linear structural relationships for \( M \, \leftarrow \, C \oplus A \oplus M \) or \( Y \, \leftarrow \, C \oplus A \oplus M \).

Asymptotic Behavior

In this section, we first give the asymptotic properties of the OLS estimators when the model satisfies Assumption (ref). Then we will establish the asymptotic normality of quadruply robust estimators, allowing the model misspecification.

Asymptotic Properties of OLS estimators

In section (ref), we highlighted that given the causal structure is appropriately specified as semi-linear according to Assumption (ref), one can employ OLS estimators by just applying two simple regressions. As we allow the number of mediators $p$ can grow with sample size $n$, some assumptions are required. The following assumptions come from portnoy1984asymptotic and portnoy1985asymptotic. They control the behavior of minimum eigenvalue will hold in probability if the observations are a sample from an appropriate distribution in $\mathbb{R}^p$. Denote the error vector as $\epsilon = (\epsilon_A^{\top}, \epsilon_M^{\top}, \epsilon_Y^{\top})^{\top}$.

We use the bold symbol $\mathbf{X}$ to represent the data matrix of any i.i.d. random observations $\{ X_i \}_{i = 1}^n $. i.e. $\mathbf{X} = [X_1, \ldots, X_n]^{\top}$. Denote the transformation $\widehat\Gamma_{\boldsymbol{\cdot}, \boldsymbol{\cdot}}: \mathbb{R}^{n \times p_1} \times \mathbb{R}^{n \times p_2} \rightarrow \mathbb{R}^{p_1 \times n}$ of two data matrix with sample size $n$ as

where $P_{\mathbf{Z}} = \mathbf{Z} (\mathbf{Z}^{\top} \mathbf{Z})^{-1} \mathbf{Z}^{\top} \in \mathbb{R}^{n \times n}$ is the projection matrix of $\mathbf{Z}$. This transformation streamlines our representation of the asymptotics for our OLS estimators.

The above theorem ensures that under mild conditions we can construct valid confidence intervals for $\widehat{DE}^{\text{OLS}}$, $\widehat{IE}^{\text{OLS}}$, and $\widehat{DM_j}^{\text{OLS}}$ when $n$ is large enough. It is worthy to note that the probabilities for $\widehat{IE}^{\text{OLS}}$ and $\widehat{DM_j}^{\text{OLS}}$ is $\geq$ instead of $=$. This distinction arises from the dual nature of the limiting distributions for these two OLS estimators: one is the standard normal, the other is not. However, as argued in chakrabortty2018inference, the non-standard asymptotic distributions here are more conservative than $\mathcal{N} (0, 1)$. Thus, we obtain $\geq$ instead of $=$. The details can be found in the proof. Notably, these asymptotic confidence intervals can be derived concurrently with the regression estimators and residuals. When applying the regression to procure these estimators, no additional steps are needed to obtain these confidence intervals.

For the estimators $\widehat{IM}_j^{\text{OLS}}$, additional assumptions are needed due to their reliance on the unknown DAG structure. This necessitates consistent CPDAG estimation, as well as more strong sparsity assumptions and restricted eigenvalue conditions, which are common in high-dimensional settings portnoy1985asymptotic, van2014asymptotically, zhang2014confidence, chakrabortty2018inference.

We now therefore obtain a valid asymptotic confidence interval for $\overline{IM}_j$ for any $j \in [p]$ alongside the regression from Theorem (ref).

Asymptotic Properties of Quadruply Robust Estimators

The quadruply robust estimators aim to obtain the robust estimators even when the model is misspeficied. The double robust estimators, which combines the direct and IPW strategies, possess commendable properties and have been the subject of extensive research as evidenced in literature such as laan2003unified,tsiatis2006semiparametric,kang2007demystifying. In this section, we will show that the proposed novel quadruply robust estimators exhibit more favorable asymptotic properties.

To present the results, we assume that the propensity score $e_{a'} (\boldsymbol{\cdot}) \in \mathcal{E}$ with some function classes $\mathcal{E}$. Similarly, for each $j \in [p]$, we assume any conditional density employed in (ref) and (ref) is \[ f_{M_T \mid X_S}(m_T \mid x_{S}) \, \in \, \mathcal{F}_{j, \, T \mid S}, \] and any conditional mean used in (ref) and (ref) adheres to \[ {\mathrm{E}} [Y \mid x_{S}] \, \in \, \mathcal{U}_{j, \, S} \] with some specific function classes $\mathcal{F}_{j, \, T \mid S}$ and $\mathcal{U}_{j, \, S}$. We propose the following assumptions concerning these function classes and the convergence rates of the estimators within these classes.

Assumption (ref) is reasonably moderate, as the function classes are user-defined. VC-type classes encompass a broad spectrum of functional categories, including but not limited to classic parametric model, neural networks and regression trees. The VC index governs the complexity of the model, typically escalating with an increase in the number of parameters within the model. We permit the VC index to diverge alongside the sample size, which serves to minimize the estimator's bias arising from model misspecification. On the other hand, an important feature of Assumption (ref) is that the required estimators' convergence rates can only be nonparametric (slower than $n^{-1 / 2}$ ) and no metric entropy condition (Donsker class for instance) is needed. In particular, $\vartheta_{j, e}^*, \vartheta_{j, \pi}^*, \vartheta_{j, \mu}^* > 1 / 4$ will perfectly admit Assumption (ref). Therefore, the estimators can be computed via standard nonparametric estimation fan2003nonlinear and supervised learning algorithms wager2018estimation,schmidt2020nonparametric. The reason both assumptions regarding the sizes of the function classes and the rate of convergence for $\widehat{DM}_j^{\text{QR}}$ and $\widehat{IM}_j^{\text{avg}, \, \text{QR}}$ are identical, which are different from conditions in Theorem (ref) and Theorem (ref), stems from the uniform convergence characteristics of our estimations for any $\mathcal{G}_M$ in $\operatorname{MEC}(\mathcal{C}_M)$.

An important result of Theorem (ref) is that: for any $j \in [p]$, the quadruply robust estimators $\widehat{DM}_j^{\text{QR}}$ and $\widehat{IM}_j^{\text{avg}, \, \text{QR}}$ are semiparametric locally efficient in the sense that they are regular and asymptotically linear under model $\mathscr{M}_{j, \, \text{union}}$, and achieve the semiparametric efficiency bound for ${DM}_j$ and $\overline{IM}_j$ under model at the intersection submodel $\mathscr{M}_0 \, \cap \, \mathscr{M}_{j, 1} \, \cap \, \mathscr{M}_{j, 2} \, \cap \, \mathscr{M}_{j, 3}$. Hence, when all models are correct, $\widehat{DM}_j^{\text{QR}}$ and $\widehat{IM}_j^{\text{avg}, \, \text{QR}}$ are semiparametric efficient in the model $\mathscr{M}_{\text{nonpar}}$ at the intersection submodel $\mathscr{M}_0 \, \cap \, \mathscr{M}_{j, 1} \, \cap \, \mathscr{M}_{j, 2} \, \cap \, \mathscr{M}_{j, 3}$ by part iv in bickel2001inference for any $j \in [p]$.

Simulation Studies

In this section, we assess the finite-sample performance of our proposed quadruply robust estimators across two simulation scenarios. The first scenario seeks to illustrate the robustness characteristics of our estimator in comparison to other estimation strategies, particularly when certain model specifications to a specific mediator are not met. In the second simulation study, we demonstrate that our method can also be superior to any other estimation strategies in estimating the both direct and indirect interventional effects across all mediators, on average, within commonly adopted model configurations.

Simulation for a single mediator

We consider the finite-sample performance of the proposed quadruply robust estimators in comparison to the estimators under direct strategy, and the alternative strategies in Section (ref), (ref), (ref), and (ref) for a single mediator. We describe the detailed setting as follows: we set $t = p = 3$, and fix the pre-specified $j$ randomly sampled from $U\{ [p] \}$. Then we design the following four data generating processes (DGPs), here $\Phi(\boldsymbol{\cdot})$ and $\operatorname{logit}(\boldsymbol{\cdot})$ are the standard normal distribution function and the inverse of the standard logistic function, and $B_{j :}$ and $B_{- j :}$ represent the $j$-th row of the matrix $B$ and the matrix $B$ with $j$-th row removed.

Here the true adjacency matrix of mediators is generated from the Erd\H{o}s-Rényi (ER) model with an expected degree as $\lfloor p / 2 \rfloor$, and the non-zero entries in $B_{MM} \in \mathbb{R}^{p \times p}$ and all the elements in $\alpha_{YA}, \beta_{MA}, \beta_{YC}, \beta_{YM}, B_{MC}$ are independently sampled from $U(-1, 1)$. In each estimation method, we consistently treat $\mathscr{M}_0$ as the underlying true model by default. We generate $n = 1000$ simulation samples, each comprising $N = 100$ independent observations, and the result for estimating the direct and indirect interventional effect of the pre-specified mediator $M_j$ is shown in Table (ref). Here, we use the PC algorithm harris2013pc to estimate the adjacency matrix of CPDAGs.

As illustrated in Table (ref), the simulation results align with the theoretical predictions made in previous sections. Specifically, when the entire distribution $F_X(\boldsymbol{\cdot})$ is correctly specified, all estimators display consistency. However, in the presence of at least one misspecified component, only the quadruply robust estimator retains consistency. In contrast, one among the other estimators, $\mathscr{M}_{\ell}$ for $\ell = 0, 1, 2, 3$, becomes inconsistent. Although we present only the continuous scenario in this part, our simulations under discrete $C$ or $M$ settings yielded similar outcomes. Importantly, under this simulation scenario, the estimator $\widehat{TM}_j^{\mathscr{M}_0} = \widehat{DM}_j^{\mathscr{M}_0} + \widehat{IM}_j^{\mathscr{M}_0}$ corresponds precisely to the estimator utilized for the individual mediation effect $\eta_j$ proposed in chakrabortty2018inference. Thus, our quadruply robust estimators outperform the estimator defined in chakrabortty2018inference.

Simulation for all mediators

Next, we consider the average performance of our quadruply robust estimators compared with other estimations under a fair model misspecification scenario in both continuous case (Section 1.4 in kang2007demystifying) and discrete case (Section 4.1 in xia2023identification). The DGPs are defined as follows:

Here $\sigma_1^2, \ldots, \sigma_p^2$ are independently drawn from the uniform distribution in $[0.5,1]$, whereas the other setting is the same as previous. In the continuous setting, instead of observing the $Z_i$ 's, we observe $C_i$ as the transformations of $Z_i$. We will always leave out the interaction and $x^{2 / 3}$ when fitting each model, and we also assume the link functions are all Probit. For computations in the continuous $M$ setting, we implement Algorithm (ref). While in the discrete $M$ context, we employ Algorithm (ref), setting the Monte Carlo sample size to $L = 100$.

As illustrated in Figure (ref), aside from the quadruply estimators (QR), other methods fail to yield consistent results. Furthermore, in most cases, our quadruply estimators exhibit a lower standard error compared to other methods. Thus, this also shows the robustness of our estimators.

Empirical Study

In this section, we illustrate our estimator in a real world application from AURORA study to explore the causal association of psychiatric disorders among trauma survivors, which is also studied in watson2023heterogeneous. In the study, our primary response of interest is the post-traumatic stress disorder (PTSD), which was assessed three months post-trauma $Y$. The focal event, in this case, is the pre-trauma insomnia $A$ that trauma survivors often experience: $A = 1$ represents survivor does have insomnia and $A = 0$ represents does not. The 4-dimensional potential mediator $M$ including Peri-traumatic PT (PTSD), stress, acute distress (ASD), and depression, gauged two weeks subsequent to the traumatic incident, are included in our analysis. This study also accounts for various confounders $C$ is a 9-dimensional vector such as age, gender, race, education level, pre-trauma physical and mental health, perceived stress level, neuroticism, and childhood trauma. The same as watson2023heterogeneous, before employing our methodology, categorical variables underwent one-hot encoding, numerical variables were centralized, and any missing data was excluded. The total number of observations is $n = 1494$ with $t = 10$ and $p = 4$. The estimated DAG of the mediators by PC algorithm harris2013pc is shown in Figure (ref). Results from the quadruply robust estimators, as obtained using Algorithm (ref) with Monte Carlo sample size $L = 100$ and a bootstrap number of $B = 500$, along with other estimation methods employing a \(\log n\) truncation and the same Monte Carlo sample size and bootstrap number, are presented in Table (ref).

As demonstrated in Table (ref), for each mediator under consideration, a substantial discrepancy is observed between the estimates of $\mathscr{M}_{\ell_1}$ and those of $\mathscr{M}_{\ell_2}$ for $\ell_2 \neq \ell_1$ when employing any of the four estimation methods $\mathscr{M}_{\ell}$ for $\ell = 0, 1, 2, 3$. Moreover, none of these estimation methods manage to identify significant direct or indirect interventional effects for any of the mediators. This highlights the pressing need for robust estimation approaches in this dataset. Notably, with the quadruply robust estimation, we discern that both the indirect interventional effect of acute distress and the direct interventional effect of peritraumatic PT are significant at the 95% confidence level, while other effects remain non-significant. These results also indicate that preventive intervention of 3-month PTSD after trauma exposure that focuses on reducing acute distress and peritraumatic PT is more likely to be effective for trauma survivors.

Discussion

The main contribution of this article is the introduction of direct and indirect interventional effects of mediators, alongside their semiparametric bounds and quadruply robust estimators. Our method accommodates continuous, categorical, and multivariate pre-treatments, mediators, and outcomes. Moreover, extending our methodology and theory to polytomous exposures is straightforward. However, extending to continuous exposure, even under LSEMs, is non-trivial in theoretical sense. A potential method is suggested in cai2021deep to replace the indicator function $\mathds{1}(A = a)$ with some kernel function $K ((A - a) / h)$ under bandwidth $h$, but as discussed in diaz2013targeted, kennedy2017non, kennedy2023semiparametric, pathwise differentiability will fail in this case, necessitating alternative estimation procedures and techniques. On the other hand, note that our framework is dimensional-free, as long as the conditional densities and expectations meet the mild convergence rate, our quadruply robust estimations will always achieve semiparametric efficiency. However, in high-dimensional cases, non-parametric estimation mentioned in this article may not achieve the rate, necessitating additional assumptions like symmetry and shape constraints, as discussed in deng2021confidence, xu2021high, rodriguez2022data. Introducing these assumptions still validates the semiparametric framework in our article under the full nonparametric model $\mathscr{M}_{\text{nonpar}}$, but our quadruply estimators may not be the most efficient under these added conditions.