A Causal Inference Approach for Mediated Moderation with Multiple Mediators

Abstract

Understanding how effect heterogeneity across subgroups is mediated by multiple mediators is important yet under-studied. Despite the growing popularity of causal mediation analysis, existing methods rarely address the mediation of moderated treatment effects, particularly when multiple mediators are of interest. This study develops a causal inference approach for mediated moderation analysis with multiple mediators, decomposing a moderated treatment effect into mediated moderation effects and remaining moderation. We present causal estimands and extend a multiply robust estimator that can incorporate machine learning techniques to relax modeling assumptions. Simulations were conducted to evaluate the method’s performance. An empirical example about adolescent mental health illustrates the application. We hope this study provides a novel causal inference-based approach to understanding multiple mediating mechanisms underlying subgroup heterogeneity in treatment effects.

Mediated moderation analyses are widely used to examine the extent to which a between-subgroup difference in the effect of a treatment on an outcome can be accounted for by between-subgroup differences in effects of the treatment on mediators (Fairchild & MacKinnon, 2009; Morgan-Lopez et al., 2003; Muller et al., 2005). It helps researchers address questions such as “why do treatment effects vary across subgroups” and thereby deepen understandings about effect heterogeneity. For example, mediated moderation analysis was used to investigate how the socioeconomic status (SES) differences in the effect of a stressful exposure on mental health could be accounted by mediators related to psychosocial resources (Miller et al., 2013).

To conduct mediated moderation analyses, most existing methods rely on conventional parametric modeling approaches, such as parametric structural equation modeling (Baron & Kenny, 1986; Fairchild & MacKinnon, 2009; Hayes, 2013; Muller et al., 2005). Despite popularity of the parametric modeling-based approaches, they lack formal frameworks to clarify the assumptions under which the estimates have a causal interpretation. Over the past decades, causal inference methods have been fast-growing, which provide frameworks (e.g., potential outcomes; Rubin, 1974, 1978) that enable researchers to define causal effects and explicate the assumptions for identifying causal effects from data. More recently, causal inference methods have been extended to mediation analysis, which contributes various tools to define and identify mediation effects (e.g., the natural indirect effect; Pearl, 2001; Robins & Greenland, 1992) without requiring any parametric models, as well as robust and flexible estimation approaches (e.g., machine learning-based estimators; Benkeser & Ran, 2021) for the defined effects (Hong, 2015; Imai et al., 2010; VanderWeele, 2015).

However, such causal inference-based methods for mediated moderation analysis are underdeveloped. An exception is Liu et al. (2025), which developed causal estimands and robust estimation methods for a single mediator. Despite the contribution of their work, in practice, mediated moderation analyses in behavioral sciences often involve multiple mediators, in which the research questions concern the extent to which each mediator contributes to the effect moderation between different subgroups. To date, no studies have developed methods for addressing multiple mediators in mediated-moderation analysis using causal inference approaches.

To reduce the methodological gap, this study aims to develop causal estimands, identification results, and robust estimation methods that can incorporate machine learning approaches for a common two-mediators scenario in mediated moderation analysis. Specifically, we consider the subgroup status to be an observed pretreatment variable, which is common in studies about effect moderation (Fairchild & MacKinnon, 2009; Qin & Wang, 2023). Furthermore, following the distinction between “effect moderation” (or effect modification) and “interaction” in the causal inference literature (VanderWeele, 2009), our study focuses on the former, that is, we consider that manipulating the subgroup status is of no substantive interest, and the substantive interest concerns how the effect of treatment varies between the subgroups. This focus naturally allows non-manipulable subgroup status (e.g., sex, race/ethnicity) to be used in our methods.

In the rest of this manuscript, we first define causal estimands for mediated moderation with two mediators, based on the interventional effects approach (Liu et al., 2025; Vansteelandt & Daniel, 2017) reviewed in the next section. Then, we present identification results for the defined causal estimands. Subsequently, we develop methods for estimation and inference by extending nonparametric estimators from the causal mediation literature (Benkeser & Ran, 2021; Rudolph et al., 2024). Under achievable conditions, the developed methods provide multiple robustness in consistency and asymptotically proper statistical inferences even with the use of machine learning models. A simulation study is presented in Supplementary Appendix I (available in the online version of this article). To illustrate the application, we provide an empirical example on adolescent mental health. We conclude by discussing limitations of the current study and future research directions.

Causal Estimands of Mediated Moderation

In this section, we define causal estimands of mediated moderation within the potential outcome framework (Rubin, 1974, 1978). To ease presentation, we introduce a toy example (Figure 1) about the moderation of the effect of a home visiting program (treatment $T$ : 1 = program, 0 = no program) on child outcome $Y$ between avoidant and non-avoidant subgroups (Cluxton-Keller et al., 2014). Specifically, in the avoidant subgroup (coded as $R={\mathrm{r}}_{\mathrm{foc}}$ , our focal subgroup), the mothers showed an avoidant characteristic at baseline (i.e., pretreatment), an attachment characteristic reflecting discomforts with trust; in the non-avoidant subgroup (coded as $R={\mathrm{r}}_{\mathrm{ref}}$ , our reference subgroup), the mothers had no such characteristic at baseline. While the program effect was beneficial for the non-avoidant subgroup, the avoidant subgroup experienced no such benefits. To investigate potential mediating pathways for extending the program benefits to the avoidant subgroup, two mediators were hypothesized: maternal depression and parenting stress ( $M_1$ and $M_2$ , respectively), which were concurrently measured during follow-up. The hypothesis was that, perhaps, compared with the non-avoidant subgroup, the avoidant subgroup experienced less beneficial (or even adverse) program effects on maternal depression and/or parenting stress, which, in turn, contributed to the less beneficial program effect on the child’s outcome.

Figure 1.

Conceptual illustration of the mediated moderation scenario in the toy example.

Let $M_1(t)$ and $M_2(t)$ denote, respectively, the potential mediators for $M_1$ and $M_2$ under treatment value $t$ , for example, the maternal depression and parenting stress levels that a mother would experience under the home visiting program ( $t=1$ ) or under no program ( $t=0$ ). Let $Y(t,m_1,m_2)$ denote the potential outcome if the treatment $T$ , and mediators $M_1$ and $M_2$ were set to values $t$ , $m_1$ , and $m_2$ , respectively. These definitions require the stable unit treatment value assumption (Rubin, 1980, 1986), which assumes that (a) there are no multiple versions of each treatment value and (b) there is no interference among individuals. Under this assumption, the observed mediators equal the potential mediators under the observed value of the treatment variable, $M_1=M_1(T)$ and $M_2=M_2(T)$ ; the observed outcome equals the potential outcome under the observed values of the treatment and mediators,

$$Y=Y(T,M_1,M_2)=Y(T,M_1(T),M_2(T)).$$

A Brief Review of the Interventional Indirect Effects with Multiple Mediatorsand the Decomposition of Total Moderation

To ease introducing our approach, which builds on the interventional effects defined in Vansteelandt and Daniel (2017) and Liu et al. (2025), we briefly review the effect definitions in these prior studies. Then, we define our estimands in the next subsection.

Interventional Indirect Effects with Multiple Mediators

Considering two mediators, Vansteelandt and Daniel (2017) decomposed the total effect of treatment $T$ on $Y$ into four components: the Interventional Indirect Effects (IIEs) via $M_1$ and via $M_2$ , the IIE due to mutual dependence between $M_1$ and $M_2$ , and the interventional direct effect of $T$ on $Y$ not through the mediators. Below, we review the definition of the IIE via $M_1$ , which allows us to illustrate how we build on the definitions in Vansteelandt and Daniel (2017) to define mediated moderation effects with two mediators.

Let $X$ denote a set of pretreatment variables (e.g., in Figure 1, $X=\{R,C\}$ , where $C$ contains the covariates), and let $\mathrm{p}(x)$ denote the distribution of $X$ in the population. The IIE of treatment $T$ on outcome $Y$ via mediator $M_1$ defined in Vansteelandt and Daniel (2017) can be written as ${\mathrm{IIE}}_{M_1}=\sum _x{\mathrm{IIE}}_{M_1}(x)\mathrm{p}(x)$ , with

$$\begin{array}{c}{\mathrm{IIE}}_{M_1}(x)=\mathsf{E}[Y(1,{\mathcal{M}}_1(1\mid x),{\mathcal{M}}_2(0\mid x))-Y(1,{\mathcal{M}}_1(0\mid x),{\mathcal{M}}_2(0\mid x))\mid x],\hfill \end{array}$$

(1)

where $Y(1,{\mathcal{M}}_1(1\mid x),{\mathcal{M}}_2(0\mid x))$ represents the potential outcome of an individual if the treatment is assigned to value $1$ and mediators $M_1$ and $M_2$ are assigned as two independent random draws from, respectively, the potential marginal distribution of $M_1$ under treatment value $1$ given covariate value $x$ and the potential marginal distribution of $M_2$ under treatment value $0$ given covariate value $x$ . Similarly, $Y(1,{\mathcal{M}}_1(0\mid x),{\mathcal{M}}_2(0\mid x))$ represents the potential outcome if $T$ is assigned to value $t$ , and $M_1$ and $M_2$ are assigned as two independent random draws from, respectively, the potential marginal distribution of $M_1$ and that of $M_2$ under treatment value $0$ given covariate value $x$ . Supplementary Appendix A (available in the online version of this article) provides the mathematical definition of ${\mathrm{IIE}}_{M_1}$ that explicitly expresses it in terms of the mediator distributions of interest. Throughout the main text, for notational simplicity, we use summation notation such as $\sum _{x}f(x)\mathrm{p}(x)$ (where $f$ is a function) to represent the corresponding integral $\int f(x)\mathrm{p}(x)\mathrm{d}v(x)$ (where $v$ is a dominating measure for the distribution of $X$ ).

The effect defined above (Equation 1) corresponds to an indirect effect via $M_1$ alone, in the sense that it measures the expected difference in the outcome due to shifting the distribution of $M_1$ from its potential distribution under treatment value $0$ to that under treatment value $1$ , while the treatment is fixed and the distribution of the other mediator is also fixed. Thus, this effect captures, within strata of covariates, the part of the potential outcome difference between two treatment values ( $0$ vs. $1$ ) that is attributable to the difference between the potential distributions of $M_1$ under these two different treatment values.

Building on this definition (Equation 1), we define the mediated moderation via each mediator (e.g., $M_1$ ) by instead contrasting the potential distributions of $M_1$ between individuals in different subgroups under each treatment value, while the subgroup status is fixed and the potential distribution of the other mediator is also fixed. Similarly, we extend the other definitions in Vansteelandt and Daniel (2017) and define the mediated moderation via $M_2$ and due to the mediators’ mutual dependence.

As described later, these separate mediated moderation effects via $M_1$ , via $M_2$ , and due to their mutual dependence decompose the mediated moderation defined in Liu et al. (2025) when treating $M=(M_1,M_2)$ as a single, joint mediator vector; the definition in Liu et al. (2025) is reviewed next.

Decomposition of Total Moderation

Liu et al. (2025) decomposed the total moderation into the mediated moderation via mediator $M$ and the remaining moderation. Specifically, the total moderation, $\mathrm{TotMod}$ , is defined by contrasting the treatment effects for subgroups $R={\mathrm{r}}_{\mathrm{foc}}$ vs. $R={\mathrm{r}}_{\mathrm{ref}}$ with balanced distributions of pretreatment covariates $C$ : $\mathrm{TotMod}=\sum _c\mathrm{TotMod}(c)\mathsf{p}(c)$ , where $\mathsf{p}(c)=\mathsf{p}(C=c)$ represents the distribution of the covariates in the population and, letting $Y(t)$ denote the potential outcome under treatment value $t$ ,

$$\begin{array}{c}\multicolumn{1}{c}{\mathrm{TotMod}(c)=\mathsf{E}[Y(1)-Y(0)\mid R={\mathrm{r}}_{\mathrm{foc}},c]}\end{array}$$

(2)

$$\begin{array}{c}-\mathsf{E}[Y(1)-Y(0)\mid R={\mathrm{r}}_{\mathrm{ref}},c].\hfill \end{array}$$

(3)

In our example (Figure 1), $\mathrm{TotMod}$ quantifies the moderation in the home visiting treatment effect between the avoidant versus non-avoidant subgroups ( ${\mathrm{r}}_{\mathrm{foc}}$ vs. ${\mathrm{r}}_{\mathrm{ref}}$ ), controlling for subgroup differences in pretreatment covariates.

To decompose the total moderation, consider the potential outcomes for the avoidant subgroup ${\mathrm{r}}_{\mathrm{foc}}$ if their mediator levels were each set as a random draw from the distribution of potential mediator levels of the non-avoidant subgroup ${\mathrm{r}}_{\mathrm{ref}}$ , within covariate stratum $c$ . Denote this random draw as $\mathsf{M}(t\mid {\mathrm{r}}_{\mathrm{ref}},c)$ for $t\in \{0,1\}$ , which is drawn from the distribution $\mathrm{p}(\mathsf{M}(t)\mid R={\mathrm{r}}_{\mathrm{ref}},c)$ with $\mathsf{M}(t)=(M_1(t),M_2(t))$ .

The mediated moderation via $M=(M_1,M_2)$ , ${\mathrm{MedMod}}_{\mathrm{joint}}$ , is defined as ${\mathrm{MedMod}}_{\mathrm{joint}}=\sum _c{\mathrm{MedMod}}_{\mathrm{joint}}(c)\mathrm{p}(c)$ , where

$$\begin{array}{c}{\mathrm{MedMod}}_{\mathrm{joint}}(c)=\mathsf{E}[Y(1,\mathsf{M}(1))-Y(0,\mathsf{M}(0))\mid R={\mathrm{r}}_{\mathrm{foc}},c]\hfill \end{array}$$

(4)

$$\begin{array}{c}-\mathsf{E}[Y(1,\mathsf{M}(1\mid {\mathrm{r}}_{\mathrm{ref}},c))-Y(0,\mathsf{M}(0\mid {\mathrm{r}}_{\mathrm{ref}},c))\mid R={\mathrm{r}}_{\mathrm{foc}},c].\hfill \end{array}$$

(5)

${\mathrm{MedMod}}_{\mathrm{joint}}$ measures the change in the treatment effect for the focal subgroup ${\mathrm{r}}_{\mathrm{foc}}$ if their mediator levels were set to the levels drawn from the reference subgroup ${\mathrm{r}}_{\mathrm{ref}}$ , within strata of the covariates. Then, the remaining moderation, $\mathrm{RemainMod}$ , considers the moderation that would remain under this counterfactual scenario, that is: $\mathrm{RemainMod}=\sum _c\mathrm{RemainMod}(c)\mathrm{p}(c)$ , where

$$\begin{array}{c}\mathrm{RemainMod}\left(c\right)=\mathsf{E}[Y(1,\mathsf{M}(1\mid {\mathrm{r}}_{\mathrm{ref}},c))-Y(0,\mathsf{M}(0\mid {\mathrm{r}}_{\mathrm{ref}},c))\mid R={\mathrm{r}}_{\mathrm{foc}},c]\hfill \end{array}$$

(6)

$$\begin{array}{c}-\mathsf{E}[Y(1,M(1))-Y(0,M(0))\mid R={\mathrm{r}}_{\mathrm{ref}},c].\hfill \end{array}$$

(7)

$\mathrm{RemainMod}$ measures the moderation in the treatment effect that would still remain if the potential mediator distributions among the focal subgroup ${\mathrm{r}}_{\mathrm{foc}}$ were shifted to the distributions among the reference subgroup ${\mathrm{r}}_{\mathrm{ref}}$ , within strata of the covariates. Supplementary Appendix A (available in the online version of this article) provides the mathematical definitions of the above two estimands that explicitly express them in terms of the mediator distributions of interest.

The mediated and remaining moderation decompose the total moderation, $\mathrm{TotMod}=\mathrm{MedMo}{\mathrm{d}}_{\mathrm{joint}}+\mathrm{RemainMod}$ . This follows by noticing the equality $Y(t)=Y(t,M(t))$ . Inspecting this decomposition, ${\mathrm{MedMod}}_{\mathrm{joint}}$ and $\mathrm{RemainMod}$ captures, respectively, the treatment effect moderation that would reduce and that would remain, if the subgroups were equalized in the potential distributions of mediator $M$ within strata of the covariates.

This decomposition in Liu et al. (2025), however, provides little insights into the role of each mediator in mediating the heterogeneity (i.e., moderation).

Mediated Moderation via Each Mediator

To disentangle the role of each mediator, we extend Vansteelandt and Daniel (2017) and decompose the jointly mediated moderation, ${\mathrm{MedMod}}_{\mathrm{joint}}$ , into three components: the mediated moderation via $M_1$ , via $M_2$ , and due to the mediators’ mutual dependence.

As described earlier, we define the mediated moderation via $M_1$ (or $M_2$ ) by contrasting the potential distributions of mediator $M_1$ (or $M_2$ ) between the two subgroups. Specifically, for $t\in \{0,1\},r\in \{{\mathrm{r}}_{\mathrm{foc}},{\mathrm{r}}_{\mathrm{ref}}\}$ , we use ${\mathrm{M}}_1(t\mid r,c)$ to denote a random draw from $\mathrm{p}(M_1(t)\mid r,c)$ , the potential distribution of mediator $M_1$ among individuals in subgroup $r$ within the covariate stratum $c$ if they were assigned to treatment value $t$ . Analogously, for mediator $M_2$ , we use ${\mathrm{M}}_2(t\mid r,c)$ to denote a random draw from $\mathrm{p}(M_2(t)\mid r,c)$ .

We define the mediated moderation via $M_1$ as ${\mathrm{MedMod}}_{M_1}=\sum _c{\mathrm{MedMod}}_{M_1}(c)\mathrm{p}(c)$ , where

$$\begin{array}{c}\mathsf{E}[Y(1,{\mathrm{M}}_1(1\mid {\mathbf{r}}_{\mathbf{foc}},c),{\mathrm{M}}_2(1\mid {\mathrm{r}}_{\mathrm{ref}},c))-Y(0,{\mathrm{M}}_1(0\mid {\mathbf{r}}_{\mathbf{foc}},c),M_2(0\mid {\mathrm{r}}_{\mathrm{ref}},c))\mid R={\mathrm{r}}_{\mathrm{foc}},c]\hfill \end{array}$$

(8)

$$\begin{array}{c}-\mathsf{E}[Y(1,{\mathrm{M}}_1(1\mid {\mathbf{r}}_{\mathbf{ref}},c),{\mathrm{M}}_2(1\mid {\mathrm{r}}_{\mathrm{ref}},c))-Y(0,{\mathrm{M}}_1(0\mid {\mathbf{r}}_{\mathbf{ref}},c),{\mathrm{M}}_2(0\mid {\mathrm{r}}_{\mathrm{ref}},c))\mid R={\mathrm{r}}_{\mathrm{foc}},c],\hfill \end{array}$$

(9)

where the highlighted parts (in bold) are contrasted between the first and the second lines. This effect expresses the change in the treatment effect for the focal subgroup ${\mathrm{r}}_{\mathrm{foc}}$ , if, under each treatment value, the distribution of mediator $M_1$ were shifted from the mediator’s potential distribution among them to that among the reference subgroup ${\mathrm{r}}_{\mathrm{ref}}$ (given covariates), while holding the distribution of $M_2$ at its potential distribution among ${\mathrm{r}}_{\mathrm{ref}}$ (given covariates).

${\mathrm{MedMod}}_{M_1}$ captures the moderated treatment effect for the focal (e.g., avoidant) subgroup due to the between-subgroup difference in the potential distribution of mediator $M_1$ (e.g., maternal depression) under each treatment value, while the potential distribution of the other mediator ( $M_2$ ) under each treatment value is fixed. Specifically, in Equation (8) of the definition, the distribution of $M_1$ is kept at its potential distribution among subgroup ${\mathrm{r}}_{\mathrm{foc}}$ , by assigning $M_1$ as a random draw from this distribution (i.e., ${\mathrm{M}}_1(t\mid {\mathrm{r}}_{\mathrm{foc}},c)$ , for $t\in \{0,1\}$ }); this is contrasted with Equation (9), where the distribution of $M_1$ is set to its potential distribution among subgroup ${\mathrm{r}}_{\mathrm{ref}}$ . In both lines, the potential distribution of $M_2$ under each treatment value is fixed to that among a specific subgroup; similarly, the subgroup among which different potential outcomes are compared is always the focal subgroup.

Corresponding to the definition for $M_1$ , we define the mediated moderation via mediator $M_2$ as ${\mathrm{MedMod}}_{M_2}=\sum _c{\mathrm{MedMod}}_{M_2}(c)\mathrm{p}(c)$ , where

$$\begin{array}{c}\mathsf{E}[Y(1,{\mathrm{M}}_1(1\mid {\mathrm{r}}_{\mathrm{foc}},c),{\mathrm{M}}_2(1\mid {\mathbf{r}}_{\mathbf{foc}},c))-Y(0,{\mathrm{M}}_1(0\mid {\mathrm{r}}_{\mathrm{foc}},c),{\mathrm{M}}_2(0\mid {\mathbf{r}}_{\mathbf{foc}},c))\mid R={\mathrm{r}}_{\mathrm{foc}},c]\hfill \end{array}$$

(10)

$$\begin{array}{c}-\mathsf{E}[Y(1,{\mathrm{M}}_1(1\mid {\mathrm{r}}_{\mathrm{foc}},c),{\mathrm{M}}_2(1\mid {\mathbf{r}}_{\mathbf{ref}},c))-Y(0,{\mathrm{M}}_1(0\mid {\mathrm{r}}_{\mathrm{foc}},c),{\mathrm{M}}_2(0\mid {\mathbf{r}}_{\mathbf{ref}},c))\mid R={\mathrm{r}}_{\mathrm{foc}},c],\hfill \end{array}$$

(11)

where the highlighted parts (in bold) are contrasted between the first and the second lines. Similar to ${\mathrm{MedMod}}_{M_1}$ , ${\mathrm{MedMod}}_{M_2}$ expresses the change in the treatment effect for the focal subgroup ${\mathrm{r}}_{\mathrm{foc}}$ , if the distribution of $M_2$ under each treatment value were shifted from its potential distribution among ${\mathrm{r}}_{\mathrm{foc}}$ to that among ${\mathrm{r}}_{\mathrm{ref}}$ , while fixing the distribution of $M_1$ under each treatment value at its potential distribution among ${\mathrm{r}}_{\mathrm{foc}}$ . Supplementary Appendix A (available in the online version of this article) provides the mathematical definitions of ${\mathrm{MedMod}}_{M_1}$ and ${\mathrm{MedMod}}_{M_2}$ , which explicitly express them in terms of the mediator distributions of interest.

The above definitions of ${\mathrm{MedMod}}_{M_1}$ and ${\mathrm{MedMod}}_{M2}$ consider the two mediators independently; they omit how the dependence between mediators, which may differ between subgroups, contributes to the effect moderation. This contribution is captured by the last estimand, the mediated moderation due to mediators’ mutual dependence, which is defined as ${\mathrm{MedMod}}_{\mathrm{mu}}=\sum _c{\mathrm{MedMod}}_{\mathrm{mu}}(c)\mathrm{p}(c)$ , where

$$\begin{array}{c}\multicolumn{1}{c}{{\mathrm{MedMod}}_{\mathrm{mu}}(c)={\mathrm{MedMod}}_{\mathrm{joint}}(c)-{\mathrm{MedMod}}_{M_1}(c)-{\mathrm{MedMod}}_{M_2}(c)=}\\ \multicolumn{1}{c}{\begin{array}{c}=\mathsf{E}[Y(1,M(1))-Y(0,M(0))\mid R={\mathrm{r}}_{\mathrm{foc}},c]\hfill \end{array}}\end{array}$$

(12)

$$\begin{array}{c}-\mathsf{E}[Y(1,{\mathrm{M}}_1(1\mid {\mathrm{r}}_{\mathrm{foc}},c),{\mathrm{M}}_2(1\mid {\mathrm{r}}_{\mathrm{foc}},c))-Y(0,{\mathrm{M}}_1(0\mid {\mathrm{r}}_{\mathrm{foc}},c),M_2(0\mid {\mathrm{r}}_{\mathrm{foc}},c))\mid R={\mathrm{r}}_{\mathrm{foc}},c]\hfill \end{array}$$

(13)

$$\begin{array}{c}+\mathsf{E}[Y(1,{\mathrm{M}}_1(1\mid {\mathrm{r}}_{\mathrm{ref}},c),{\mathrm{M}}_2(1\mid {\mathrm{r}}_{\mathrm{ref}},c))-Y(0,{\mathrm{M}}_1(0\mid {\mathrm{r}}_{\mathrm{ref}},c),{\mathrm{M}}_2(0\mid {\mathrm{r}}_{\mathrm{ref}},c))\mid R={\mathrm{r}}_{\mathrm{foc}},c]\hfill \end{array}$$

(14)

$$\begin{array}{c}-\mathsf{E}[Y(1,\mathrm{M}(1\mid {\mathrm{r}}_{\mathrm{ref}},c))-Y(0,\mathrm{M}(0\mid {\mathrm{r}}_{\mathrm{ref}},c))\mid R={\mathrm{r}}_{\mathrm{foc}},c].\hfill \end{array}$$

(15)

${\mathrm{MedMod}}_{\mathrm{mu}}$ expresses the change in the treatment effect for the focal subgroup ${\mathrm{r}}_{\mathrm{foc}}$ , if under each treatment value, the dependence of potential levels of $M_1$ and $M_2$ were shifted from the dependence strength among them (as captured by Equations 12–13) to the dependence strength among the reference subgroup ${\mathrm{r}}_{\mathrm{ref}}$ (as captured by Equations 14–15), within each stratum of covariates.

The three mediated moderation effects defined here, plus the remaining moderation, achieve a decomposition of the total moderation: $\mathrm{TotMod}=\mathrm{Med}\mathrm{Mo}{\mathrm{d}}_{M_1}+\mathrm{MedMo}{\mathrm{d}}_{M_2}+\mathrm{MedMo}{\mathrm{d}}_{\mathrm{mu}}+\mathrm{RemainMod}$ . In what follows, we collectively refer to these causal estimands $\{\mathrm{MedMo}{\mathrm{d}}_{M_1},{\mathrm{MedMod}}_{M_2},{\mathrm{MedMod}}_{\mathrm{mu}},\mathrm{RemainMod}\}$ as the “mediated moderation estimands.”

Features of the Proposed Causal Estimands

The proposed estimands are defined using the interventional effect approach (for review of different approaches to effect definition, see Nguyen et al., 2021; VanderWeele, 2015). That is, our estimands consider potential outcomes if the distributions of mediators are shifted (as if by means of an intervention for the population). This gives our estimands several features, and whether a feature is advantageous or disadvantageous depends on the substantive research context.

First, the proposed approach does not consider subgroup status $R$ as a “second” treatment variable for which manipulation is of substantive interest. This feature is misaligned with studies where the interest is not about moderation but rather about the “interaction” effects (VanderWeele, 2009, p. 864) of manipulating both $T$ and $R$ . This feature, however, would be desirable in studies investigating effect moderation by $R$ , especially when subgroup status $R$ is considered not manipulable, which is the focus of the current study. Indeed, a motivation for our use of the interventional effect approach is that non-manipulable subgroup variables, such as sex and race/ethnicity, are commonly used moderating variables in empirical studies. Nonetheless, in Supplementary Appendix G (available in the online version of this article), we present analogues of the mediated moderation estimands where $R$ is manipulated at an individual level; as described there, these analogues can be identified by the same quantities and estimated using the same methods as the mediated moderation estimands in the main text.

Second, as with the interventional effects in Vansteelandt and Daniel (2017), the identification of our causal estimands (described in the next section) is achievable without imposing any assumptions on the mediators’ causal structure (e.g., causal ordering or the absence of unmeasured common causes of mediators). Concurrent mediators with unknown causal ordering occur frequently in empirical studies. Thus, this feature of our approach offers a practical advantage, in comparison to the natural indirect effect approach in standard causal mediation analysis (Daniel et al., 2015; Valeri & VanderWeele, 2013). The latter approach generally requires making assumptions on the mediators’ causal structure that may or may not be plausible, especially with concurrent mediators (Loh et al., 2022a).

Third, irrespective of the mediators’ causal structure, ${\mathrm{MedMod}}_{M1}$ captures the effect moderation attributable to the subgroup differences in $M_1$ , but not descendants of $M_1$ ; this analogously applies to ${\mathrm{MedMod}}_{M2}$ . To illustrate, suppose the unknown causal structure between mediators is $M_1\to M_2$ . The subgroup difference in $M_1$ can contribute to the subgroup difference in $M_2$ . ${\mathrm{MedMod}}_{M2}$ captures the effect moderation due to the subgroup differences in the potential distributions of $M_2$ (i.e., $\mathrm{p}(M_2(t)\mid R={\mathrm{r}}_{\mathrm{foc}},c)$ vs. $\mathrm{p}(M_2(t)\mid R={\mathrm{r}}_{\mathrm{ref}},c)$ )), without distinguishing which part of this difference is due to $M_1$ and which is not. ${\mathrm{MedMod}}_{M_1}$ captures the effect moderation due to the subgroup differences in the potential distributions of $M_1$ , and this difference cannot be due to $M_2$ under the assumed causal ordering (i.e., $M_1\to M_2$ ). Additionally, the subgroups can differ in the strengths of dependence (e.g., covariance) between potential levels of $M_1$ and $M_2$ , which, in turn, can contribute to the effect moderation; this contribution is captured by ${\mathrm{MedMod}}_{\mathrm{mu}}$ . This feature is disadvantageous when researchers know the causal ordering between mediators and would like this knowledge to be utilized in the definition of causal estimands. Nonetheless, this feature is advantageous when no knowledge is available for the mediators’ causal structure (i.e., the scenario considered here; see Figure 2) because it allows the interpretations of the causal estimands to hold, regardless of the underlying causal structure between mediators.

Figure 2.

Causal diagram of study variables.

Fourth, the proposed causal estimands involve covariates $C$ , a feature shared by many interventional effect-based estimands (Nguyen et al., 2021; Vansteelandt & Daniel, 2017). This feature would be disadvantageous if it is undesirable to compare subgroups after controlling for the subgroup differences in the covariates. However, this feature is useful when, as often occurs in mediated moderation analysis, the goal is to decompose a differential (i.e., moderated) treatment effect between subgroups after controlling for the subgroups’ pretreatment differences (Liu et al., 2025). We further discuss this feature and its potential limitations in the “Discussion” section.

Identification

The causal estimands defined above involve potential outcomes that are unobserved. Therefore, assumptions are needed to identify them using the observed data distribution. Identification assumptions have been given for the interventional effects in causal mediation analysis (Vansteelandt & Daniel, 2017) and in mediated moderation analysis with a single mediator (Liu et al., 2025).

For identifying the mediated moderation estimands defined here, we make the following assumptions under the causal diagram in Figure 2: (a) positivity and (b) unconfoundedness. Positivity assumes that within stratum $C=c$ , there are nonzero probabilities for each treatment value $t$ , subgroup status $r$ , and mediator value $(m_1,m_2)$ , that is, $\mathrm{p}(T=t,R=r\mid C=c)>0$ and $\mathrm{p}(M_1=m_1,M_2=m_2\mid T=t,R=r,C=c)>0$ for all $t$ , $r$ , $(m_1,m_2)$ , and $c$ values.

Unconfoundedness assumes that, given observed pretreatment variables $\{R,C\}$ : (A1) the treatment–outcome relationship is unconfounded, $Y(t,m_1,m_2)\perp \perp \backslash \mathrm{doubleperp}T\mid \{R,C\}$ ; (A2) the mediators–outcome relationship is unconfounded, $Y(t,m_1,m_2)\perp \perp \{M_1,M_2\}\mid \{T=t,R,C\}$ ; and (A3) the treatment–mediator relationship is unconfounded, $\{M_1(t),M_2(t)\}\perp \perp T\mid \{R,C\}$ .

The above assumptions have several implications. First, the positivity assumption requires that the subgroups being compared have sufficient overlap in values of covariates $C$ . If the overlap is insufficient in the data at hand, the positivity assumption will likely have practical violations, which, in turn, can result in unreliable inferences on any comparisons between the subgroups controlling for $C$ (Li et al., 2013), including the mediated moderation estimands. Furthermore, the positivity assumption requires that each value of the treatment and mediators is possible to occur for each subgroup within each covariate stratum.

Second, there is no unconfoundedness assumption for the relations of subgroup status $R$ with other variables. This is because no counterfactual manipulation of $R$ is considered in the potential outcomes defined above. Rather, the subgroups are compared from a controlled descriptive perspective (J. W. Jackson, 2017; Li et al., 2013), whereby subgroups who are balanced in pretreatment covariates are compared in the potential outcomes and potential mediators. This perspective aligns with substantive interests in mediated moderation analyses, which often concerns why, after controlling for the subgroups’ pretreatment differences, the treatment effects still meaningfully vary between the subgroups (Liu et al., 2025).

Third, the unconfoundedness assumptions (A1) and (A3) require no unmeasured confounders of the treatment–outcome and treatment–mediators relations, and they are satisfied when the treatment is randomized. (A2) is a strong assumption, which requires that conditional on treatment, subgroup status, and covariates, there are no unmeasured confounders of the relation between $(M_1,M_2)$ and $Y$ . However, (A2) is still achievable by randomizing mediators within strata of $(T,R,C)$ , which is weaker than the assumption required in standard causal mediation analyses that assess the nature indirect effects; the latter requires a stringent assumption (referred to as the cross-world counterfactual independence assumption; Moreno-Betancur et al., 2021; Pearl, 2001 that cannot be satisfied even by randomization (Richardson & Robins, 2013).

For $t\in \{0,1\}$ and $r,r^{*},r^{**}\in \{{\mathrm{r}}_{\mathrm{foc}},{\mathrm{r}}_{\mathrm{ref}}\}$ , we define:

$$\begin{array}{c}{\mathrm{\theta }}_{\mathrm{joint}}(t,r,r^{*})=\sum _{c,m_1,m_2}\mathsf{E}[Y\mid m_1,m_2,t,r,c]\mathrm{p}(m_1,m_2\mid t,r^{*},c)\mathrm{p}(c),\hfill \end{array}$$

(16)

$$\begin{array}{c}{\mathrm{\theta }}_{\mathrm{each}}(t,r,r^{*},r^{**})=\sum _{c,m_1,m_2}\mathsf{E}[Y\mid m_1,m_2,t,r,c]\mathrm{p}(m_1\mid t,r^{*},c)\mathrm{p}(m_2\mid t,r^{**},c)\mathrm{p}(c).\hfill \end{array}$$

(17)

Under the identification assumptions, the mean potential outcomes in the causal estimands are identified as $\sum _c\mathsf{E}[Y(t,\mathrm{M}(t\mid r^{*},c))\mid R=r,c]\mathrm{p}(c)={\mathrm{\theta }}_{\mathrm{joint}}(t,r,r^{*})$ , $\sum _c\mathsf{E}[Y(t,M(t))\mid R=r,c]\mathrm{p}(c)={\mathrm{\theta }}_{\mathrm{joint}}(t,r,r)$ , and $\sum _c\mathsf{E}[Y(t,{\mathrm{M}}_1(t\mid r^{*},c),{\mathrm{M}}_2(t\mid r^{**},c))\mid R=r,c]\mathrm{p}(c)={\mathrm{\theta }}_{\mathrm{each}}(t,r,r^{*},r^{**})$ . Proofs can be found in Supplementary Appendix B (available in the online version of this article).

Thus, by contrasting ${\mathrm{\theta }}_{\mathrm{each}}$ and/or ${\mathrm{\theta }}_{\mathrm{joint}}$ , we identify the mediated and remaining moderation as ${\mathrm{MedMod}}_{M_1}={\mathrm{\theta }}_{\mathrm{each}}(1,{\mathrm{r}}_{\mathrm{foc}},{\mathrm{r}}_{\mathrm{foc}},{\mathrm{r}}_{\mathrm{ref}})-{\mathrm{\theta }}_{\mathrm{each}}(1,{\mathrm{r}}_{\mathrm{foc}},{\mathrm{r}}_{\mathrm{ref}},{\mathrm{r}}_{\mathrm{ref}})-\{{\mathrm{\theta }}_{\mathrm{each}}(0,{\mathrm{r}}_{\mathrm{foc}},{\mathrm{r}}_{\mathrm{foc}},{\mathrm{r}}_{\mathrm{ref}})-{\mathrm{\theta }}_{\mathrm{each}}(0,{\mathrm{r}}_{\mathrm{foc}},{\mathrm{r}}_{\mathrm{ref}},{\mathrm{r}}_{\mathrm{ref}})\}$ , ${\mathrm{MedMod}}_{M_2}={\mathrm{\theta }}_{\mathrm{each}}(1,{\mathrm{r}}_{\mathrm{foc}},{\mathrm{r}}_{\mathrm{foc}},{\mathrm{r}}_{\mathrm{foc}})-{\mathrm{\theta }}_{\mathrm{each}}(1,{\mathrm{r}}_{\mathrm{foc}},{\mathrm{r}}_{\mathrm{foc}},{\mathrm{r}}_{\mathrm{ref}})-\{{\mathrm{\theta }}_{\mathrm{each}}(0,{\mathrm{r}}_{\mathrm{foc}},{\mathrm{r}}_{\mathrm{foc}},{\mathrm{r}}_{\mathrm{foc}})-{\mathrm{\theta }}_{\mathrm{each}}(0,{\mathrm{r}}_{\mathrm{foc}},{\mathrm{r}}_{\mathrm{foc}},{\mathrm{r}}_{\mathrm{ref}})\}$ , $\mathrm{RemainMod}={\mathrm{\theta }}_{\mathrm{joint}}(1,{\mathrm{r}}_{\mathrm{foc}},{\mathrm{r}}_{\mathrm{ref}})-{\mathrm{\theta }}_{\mathrm{joint}}(0,{\mathrm{r}}_{\mathrm{foc}},{\mathrm{r}}_{\mathrm{ref}})-\{{\mathrm{\theta }}_{\mathrm{joint}}(1,{\mathrm{r}}_{\mathrm{ref}},{\mathrm{r}}_{\mathrm{ref}})-{\mathrm{\theta }}_{\mathrm{joint}}(0,{\mathrm{r}}_{\mathrm{ref}},{\mathrm{r}}_{\mathrm{ref}})\}$ , ${\mathrm{MedMod}}_{\mathrm{joint}}={\mathrm{\theta }}_{\mathrm{joint}}(1,{\mathrm{r}}_{\mathrm{foc}},{\mathrm{r}}_{\mathrm{foc}})-{\mathrm{\theta }}_{\mathrm{joint}}(0,{\mathrm{r}}_{\mathrm{foc}},{\mathrm{r}}_{\mathrm{foc}})-\{{\mathrm{\theta }}_{\mathrm{joint}}(1,{\mathrm{r}}_{\mathrm{foc}},{\mathrm{r}}_{\mathrm{ref}})-{\mathrm{\theta }}_{\mathrm{joint}}(0,{\mathrm{r}}_{\mathrm{foc}},{\mathrm{r}}_{\mathrm{ref}})\}$ , and ${\mathrm{MedMod}}_{\mathrm{mu}}={\mathrm{MedMod}}_{\mathrm{joint}}-{\mathrm{MedMod}}_{M_1}-{\mathrm{MedMod}}_{M_2}$ .

In the next section, we develop estimation methods. To simplify notation, we write ${\mathrm{\theta }}_{\mathrm{each}}(t,r,r^{*},r^{**})$ (Equation 17) and ${\mathrm{\theta }}_{\mathrm{joint}}(t,r,r^{*})$ (Equation 16) as ${\mathrm{\theta }}_{\mathrm{each}}$ , and ${\mathrm{\theta }}_{\mathrm{joint}}$ , respectively, when there is no confusion.

Estimation

In this section, we develop methods for estimation and inference of the mediated moderation estimands. We begin by examining the product-of-coefficient approach conventionally used in the parametric modeling framework (Baron & Kenny, 1986). Then, to reduce reliance on modeling assumptions that may be not grounded in prior knowledge, we develop nonparametric methods based on multiply robust estimators in causal mediation literature (Benkeser & Ran, 2021; Rudolph et al., 2024), incorporating the use of machine learning methods.

Connections with the Product-of-Coefficient Approach

The product-of-coefficient approach is typically used with continuous mediators and outcome. It involves specifying models for the mediators’ distribution $\mathrm{p}(M_1,M_2\mid T,R,C)$ and the conditional mean outcome $\mathsf{E}[Y\mid M_1,M_2,T,R,C]$ . A common specification is the following:

$$\begin{array}{c}{M}_1={\mathrm{\alpha }}_{1,0}+{\mathrm{\alpha }}_{1,r}R+{\mathrm{\alpha }}_{1,t}T+{\mathrm{\alpha }}_{1,\mathrm{rt}}\mathrm{RT}+C^{\top }{\mathrm{\alpha }}_{1,c}+{ϵ}^{m1}\hfill \end{array}$$

(18)

$$\begin{array}{c}{M}_2={\mathrm{\alpha }}_{2,0}+{\mathrm{\alpha }}_{2,r}R+{\mathrm{\alpha }}_{2,t}T+{\mathrm{\alpha }}_{2,\mathrm{rt}}\mathrm{RT}+C^{\top }{\mathrm{\alpha }}_{2,c}+{ϵ}^{m2}\hfill \end{array}$$

(19)

$$\begin{array}{c}Y={\mathrm{\beta }}_0+{\mathrm{\beta }}_{r}R+{\mathrm{\beta }}_{t}T+{\mathrm{\beta }}_{\mathrm{rt}}\mathrm{RT}+{\mathrm{\beta }}_{m1}{M}_1+{\mathrm{\beta }}_{m2}{M}_2+C^{\top }{\mathrm{\beta }}_c+{ϵ}^y,\hfill \end{array}$$

(20)

where the mediator residuals $({ϵ}^{m1},{ϵ}^{m2})$ are assumed as jointly normally distributed with means of zero and a constant covariance matrix, and the outcome residual ${ϵ}^y$ is assumed to be independent of the mediator residuals $({ϵ}^{m1},{ϵ}^{m2})$ and follows a normal distribution with mean zero and a constant variance.

Plugging the model-implied mediator distributions and conditional mean outcome models (Equations 18–20) into the identification formulas (i.e., ${\mathrm{\theta }}_{\mathrm{each}}$ and ${\mathrm{\theta }}_{\mathrm{joint}}$ in Equations 17–16), it yields the following estimators for the mediated moderation effects and remaining moderation:

$$\begin{array}{c}\begin{array}{c}{\widehat{\mathrm{MedMod}}}_{M_1}={\mathrm{\alpha }}_{1,\mathrm{rt}}{\mathrm{\beta }}_{m1},{\widehat{\mathrm{MedMod}}}_{M_2}={\mathrm{\alpha }}_{2,\mathrm{rt}}{\mathrm{\beta }}_{m2}\hfill \end{array}\hfill \\ \begin{array}{c}{\widehat{\mathrm{MedMod}}}_{\mathrm{mu}}=0,\widehat{\mathrm{RemainMod}}={\mathrm{\beta }}_{\mathrm{rt}}.\hfill \end{array}\hfill \end{array}$$

(21)

The estimators for ${\mathrm{MedMod}}_{M_1}$ and ${\mathrm{MedMod}}_{M_2}$ take the form of the product-of-coefficient estimators (Baron & Kenny, 1986; Fairchild & MacKinnon, 2009). This implies that, under the identification assumptions and assuming correct model specification, an estimate of ${\mathrm{\alpha }}_{1,\mathrm{rt}}{\mathrm{\beta }}_{m1}$ (or ${\mathrm{\alpha }}_{2,\mathrm{rt}}{\mathrm{\beta }}_{m2}$ ) from the product-of-coefficient approach has a causal interpretation as the mediated moderation via $M_1$ (or $M_2$ ) defined above in the potential outcome framework, and an estimate of ${\mathrm{\beta }}_{\mathrm{rt}}$ can be interpreted as the remaining moderation ( $\mathrm{RemainMod}$ ), which captures the treatment effect moderation that would still remain after equalizing the potential distributions of both mediators across the subgroups. Lastly, ${\widehat{\mathrm{MedMod}}}_{\mathrm{mu}}=0$ , indicating that the estimator of ${\mathrm{MedMod}}_{\mathrm{mu}}$ is restricted to zero ${\mathrm{MedMod}}_{\mathrm{mu}}$ . This is because the outcome model (20) assumes no interactions between the mediators; see Loh et al. (2022a) for discussions in the context of causal mediation analysis.

The above regression-based estimation is relatively simple to conduct. However, it has important limitations. First, if the models are misspecified, the resulting estimators for the mediated moderation estimands can be biased, even when the identification assumptions hold. Second, if data-adaptive modeling approaches (e.g., model selection, machine learning) are used to estimate components of the identification formulas (i.e., $\mathrm{p}(M_1,M_2\mid T,R,C)$ or $\mathsf{E}[Y\mid M_1,M_2,T,R,C]$ ), the resulting estimators’ sampling distributions generally have no established forms and standard nonparametric bootstrapping may not be appropriate (Miles et al., 2020; van der Laan & Rose, 2011), making it difficult to quantify the estimators’ uncertainty and obtain inferences.

Multiply Robust Estimator

To address limitations of the regression-based estimation in the previous subsection, we develop nonparametric estimators based on multiply robust methods in causal mediation literature (Benkeser & Ran, 2021; Liu, 2025; Rudolph et al., 2024). These methods often yield estimators that allow for data-adaptive modeling approaches, while maintaining the ability to provide asymptotically proper inferences (e.g., confidence intervals [CIs]).

Background of Influence Function-Based Estimation and Cross-Fitting

A central tool to develop multiply robust estimators is the influence function (IF) in semiparametric estimation (for review, see Kennedy, 2024; van der Vaart, 2002), which has contributed to many machine learning-based methods for causal inference such as targeted learning (van der Laan & Rose, 2018) and debiased machine learning (Chernozhukov et al., 2018), especially with the use of cross-fitting. To ease introducing our methods, below, we briefly review these techniques (recent tutorials can be found in Kennedy, 2024; Renson et al., 2025; van der Laan & Rose, 2011, 2018).

Intuitively, the IF acts as a first-order derivative of a functional such as ${\mathrm{\theta }}_{\mathrm{each}}$ (the identification formula, Equation 17) with respect to a model of the data distribution $\mathrm{P}$ . As a first-order derivative, the IF allows us to express how large a difference between values of ${\mathrm{\theta }}_{\mathrm{each}}$ (e.g., ${\mathrm{\theta }}_{\mathrm{each}}(\hat{\mathrm{P}})-{\mathrm{\theta }}_{\mathrm{each}}(\mathrm{P})$ )) will be produced by small, first-order deviations between two models of the data distribution (e.g., a misspecified model $\hat{\mathrm{P}}$ and a correct model $\mathrm{P}$ ). Thus, using the IF, we can construct estimators by correcting for the part of the estimation error ( ${\mathrm{\theta }}_{\mathrm{each}}(\hat{\mathrm{P}})-{\mathrm{\theta }}_{\mathrm{each}}(\mathrm{P})$ ) produced by first-order deviations between the misspecified and correct models.

An important advantage of IF-based estimators is that they can provide asymptotically proper inferences (e.g., CIs), even when components of the data distribution are estimated using data-adaptive machine learning methods. Specifically, for making inferences under the use of machine learning, a key problem is how to handle an overfitting issue related to a term of the estimator’s error (referred to as the empirical process term; Kennedy, 2024). This term would be non-negligible and result in anti-conservative CIs (Naimi et al., 2023), if we use the same data twice, first to estimate components of the data distribution via very complex machine learning models, and then to compute the IF-based estimator. An effective solution, without restricting complexity of the machine learning models, is called “cross-fitting” (or sample splitting), whereby we estimate (or train) the models using a subset of the data, and then compute the estimator using the independent remaining data (Chernozhukov et al., 2018; Zheng & van der Laan, 2011). It has been found that under common conditions, cross-fitted, IF-based estimators provide asymptotically proper inferences even under the use of machine learning (Kennedy, 2024).

Influence Functions

From the above motivation, we construct estimators using the IFs for ${\mathrm{\theta }}_{\mathrm{each}}$ and ${\mathrm{\theta }}_{\mathrm{joint}}$ . To present them, we introduce some additional notation that represents components of the IFs.

For subgroup status $r$ and treatment value $t$ , define the following mean outcome functions:

$$\begin{array}{c}\mathrm{\mu }(M_1,M_2,t,r,C)=E(Y\mid M_1,M_2,t,r,C)\hfill \end{array}$$

(22)

$$\begin{array}{c}{\mathrm{\mu }}_1(M_1,t,r,C)=\sum _{m_2}\mathrm{\mu }(M_1,m_2,t,r,C)\mathrm{p}(m_2\mid r^{**},t,C),\hfill \end{array}$$

(23)

$$\begin{array}{c}{\mathrm{\mu }}_2(M_2,t,r,C)=\sum _{m_1}\mathrm{\mu }(m_1,M_2,t,r,C)\mathrm{p}(m_1\mid r^{*},t,C),\hfill \end{array}$$

(24)

where $r$ , $r^{*},r^{**}\in \{{\mathrm{r}}_{\mathrm{foc}},{\mathrm{r}}_{\mathrm{ref}}\}$ . Equations (23) and (24) are partially marginalized mean outcome functions, which average the conditional mean outcome (Equation 22) over one mediator. For example, ${\mathrm{\mu }}_2(M_2,t,r,C)$ averages over the “marginal” distribution of $M_1$ given $R=r^{*}$ , $T=t$ , and covariates $C$ ; this distribution is “marginal” in the sense that it does not depend on the other mediator. In addition, define:

$$\begin{array}{c}{\mathrm{\mu }}_{\mathrm{each}}(t,r,C)=\sum _{m_1,m_2}\mathrm{\mu }(m_1,m_2,t,r,C)\mathrm{p}(m_1\mid r^{*},t,C)\mathrm{p}(m_2\mid r^{**},t,C)\hfill \end{array}$$

(25)

$$\begin{array}{c}{\mathrm{\mu }}_{\mathrm{joint}}(t,r,C)=\sum _{m_1,m_2}\mathrm{\mu }(m_1,m_2,t,r,C)\mathrm{p}(m_1,m_2\mid r^{*},t,C).\hfill \end{array}$$

(26)

These functions (Equations 25 and 26) are “fully marginalized” mean outcome functions, which average the conditional mean outcome (Equation 22) over both mediators.

Furthermore, define the following functions representing ratios of mediator probabilities:

$${\mathrm{w}}_{\mathrm{each}}(M_1,M_2,t,r,C)=\frac{\mathsf{p}(M_1\mid r^{*},t,C)\mathsf{p}(M_2\mid r^{**},t,C)}{\mathsf{p}(M_1,M_2\mid r,t,C)}$$

(27)

$${\mathrm{w}}_{\mathrm{joint}}(M_1,M_2,t,r,C)=\frac{\mathsf{p}(M_1,M_2\mid r^{*},t,C)}{\mathsf{p}(M_1,M_2\mid r,t,C)}.$$

(28)

${\mathrm{w}}_{\mathrm{each}}$ and ${\mathrm{w}}_{\mathrm{joint}}$ provide weights that adjust the observational distribution of mediators to the mediator distributions of interest in the respective estimands ${\mathrm{\theta }}_{\mathrm{each}}$ and ${\mathrm{\theta }}_{\mathrm{joint}}$ . Similar weights have been used in weighting-based estimators in causal mediation analysis (Hong, 2010; Loh et al., 2020). Lastly, we introduce the “propensity score” (Rosenbaum & Rubin, 1983), $\mathsf{p}(T,R\mid C)$ , representing the conditional distribution of the treatment and subgroup variables given the covariates.

Let $\eta =\{\mathrm{\mu },{\mathrm{\mu }}_1,{\mathrm{\mu }}_2,{\mathrm{\mu }}_{\mathrm{each}},{\mathrm{\mu }}_{\mathrm{joint}},{\mathrm{w}}_{\mathrm{each}},{\mathrm{w}}_{\mathrm{joint}},\mathrm{p}(t,r\mid c)\}$ , which comprises all “nuisance” functions of the data distribution that can influence ${\mathrm{\theta }}_{\mathrm{each}}$ and ${\mathrm{\theta }}_{\mathrm{joint}}$ . Let $O=(C,R,T,M_1,M_2,Y)$ collect the observed variables of an individual. For ${\mathrm{\theta }}_{\mathrm{each}}$ and ${\mathrm{\theta }}_{\mathrm{joint}}$ , the IFs are respectively ${\mathrm{\phi }}_{\mathrm{each}}(O;\mathrm{\eta })-{\mathrm{\theta }}_{\mathrm{each}}$ and ${\mathrm{\phi }}_{\mathrm{joint}}(O;\mathrm{\eta })-{\mathrm{\theta }}_{\mathrm{joint}}$ , where the uncentered IFs ${\mathrm{\phi }}_{\mathrm{each}}(O;\mathrm{\eta })$ and ${\mathrm{\phi }}_{\mathrm{joint}}(O;\mathrm{\eta })$ take the following forms:

$$\begin{array}{c}\multicolumn{1}{c}{{\mathrm{\phi }}_{\mathrm{each}}(O;\mathrm{\eta })=\frac{1\{T=t,R=r\}}{\mathsf{p}(t,r\mid C)}{\mathrm{w}}_{\mathrm{each}}(M_1,M_2,t,r,C)[Y-\mathrm{\mu }(M_1,M_2,t,r,C)]}\\ \multicolumn{1}{c}{+\frac{1\{T=t,R=r^{*}\}}{\mathsf{p}(t,r^{*}\mid C)}[{\mathrm{\mu }}_1(M_1,t,r,C)-{\mathrm{\mu }}_{\mathrm{each}}(t,r,C)]}\\ \multicolumn{1}{c}{+\frac{1\{T=t,R=r^{**}\}}{\mathsf{p}(t,r^{**}\mid C)}[{\mathrm{\mu }}_2(M_2,t,r,C)-{\mathrm{\mu }}_{\mathrm{each}}(t,r,C)]+{\mathrm{\mu }}_{\mathrm{each}}(t,r,C),}\end{array}$$

(29)

and

$$\begin{array}{l}{\mathrm{\phi }}_{\mathrm{joint}}(O;\eta )=\frac{\mathbf{1}\{T=t,R=r\}}{\mathsf{p}(t,r\mid C)}{\mathrm{w}}_{\mathrm{joint}}(M_1,M_2,t,r,C)\left[Y-\mathrm{\mu }(M_1,M_2,t,r,C)\right]\\ +\frac{\mathbf{1}\{T=t,R=r^{\ast }\}}{\mathsf{p}(t,r^{\ast }\mid C)}[\mathrm{\mu }(M_1,M_2,t,r,C)-{\mathrm{\mu }}_{\mathrm{joint}}(t,r,C)]\\ +{\mathrm{\mu }}_{\mathrm{joint}}(t,r,C).\end{array}$$

(30)

Derivations of the IFs can be found in Supplementary Appendix F (available in the online version of this article). Intuitively, these IFs account for first-order errors of the estimation of nuisance functions. For ${\mathrm{\phi }}_{\mathrm{joint}}$ , in the third line, ${\mathrm{\mu }}_{\mathrm{joint}}$ depends on both the conditional mean outcome and the joint mediator distribution; thus, an estimation error of ${\mathrm{\mu }}_{\mathrm{joint}}$ can arise due to estimation errors in any of these components. Then, the first line corrects for the estimation error of the conditional mean outcome $\mathrm{\mu }$ , by adding to ${\mathrm{\mu }}_{\mathrm{joint}}$ a “weighted residual bias correction” (Schafer & Kang, 2008, p. 301), where the weight is the product of the mediator probability ratio ${\mathrm{w}}_{\mathrm{joint}}$ (Equation 28) and the inverse propensity score; the second line corrects for the estimation error of the joint mediator distribution, by adding a weighted residual bias correction, with the weight given by the inverse propensity score. Similarly, for ${\mathrm{\phi }}_{\mathrm{each}}$ , in the fourth line ${\mathrm{\mu }}_{\mathrm{each}}$ (Equation 25, the expectation of which is the estimand ${\mathrm{\theta }}_{\mathrm{each}}$ of interest) depends on the conditional mean outcome $\mathrm{\mu }$ and the marginal distributions of $M_1$ and $M_2$ . The first line corrects for the estimation error of $\mathrm{\mu }$ by adding a weighted residual bias correction, with the weight given by the product of the mediator probability ratio ${\mathrm{w}}_{\mathrm{each}}$ and the inverse propensity score; the second and third lines correct for the estimation errors of the marginal distributions of $M_1$ and $M_2$ , respectively, by adding weighted residual bias corrections with the weights being the inverse of corresponding propensity scores.

Cross-Fitting and Machine Learning for Nuisance Function Estimation

To estimate the nuisance functions, parametric models could be used assuming the correct model specification is available. However, this assumption is often strong in practice. Alternatively, the nuisance functions can be data-adaptively estimated using machine learning methods.

As reviewed above, cross-fitting (Chernozhukov et al., 2018; Zheng & van der Laan, 2012) allows quantifying the estimation uncertainty even when complex machine learning methods are used. Therefore, to avoid restricting the complexity of estimated nuisance functions $\hat{\mathrm{\eta }}$ from machine learning methods, we use cross-fitting to obtain $\hat{\mathrm{\eta }}$ . Specifically, we randomly split the sample index set $\{1,\dots ,n\}$ into $K$ (e.g., $K=5$ ) prediction folds of approximately the same size, denoted by ${\mathcal{V}}_1,\dots ,{\mathcal{V}}_K$ ; for each prediction fold $k$ , the remaining $K-1$ folds form the training data, denoted as ${\mathcal{T}}_k=\{1,\dots ,n\}\setminus {\mathcal{V}}_k$ , which is used to estimate (or train) machine learning models for nuisance functions $\mathrm{\eta }$ . Let ${\hat{\mathrm{\eta }}}_k$ represent the estimated nuisance functions produced using training data ${\mathcal{T}}_k$ . For individual $i$ , let $k(i)$ indicate the index of the prediction fold that contains $i$ . Cross-fitting allows us to obtain the estimated $\hat{\mathrm{\eta }}$ for individual $i$ as $\hat{\mathrm{\eta }}={\hat{\mathrm{\eta }}}_{k(i)}$ , which is estimated using only the training data that exclude individual $i$ . Then, for the observed data $O_i$ of individual $i$ , we compute the value of ${\mathrm{\phi }}_{\mathrm{each}}(O_i;{\hat{\mathrm{\eta }}}_{k(i)})$ and ${\mathrm{\phi }}_{\mathrm{joint}}(O_i;{\hat{\mathrm{\eta }}}_{k(i)})$ .

With $\hat{\mathrm{\eta }}={\hat{\mathrm{\eta }}}_{k(i)}$ obtained via cross-fitting, we compute estimators of ${\mathrm{\theta }}_{\mathrm{each}}$ and ${\mathrm{\theta }}_{\mathrm{joint}}$ as:

$${\widehat{\mathrm{\theta }}}_{\mathrm{each}}=\frac{1}{n}\sum _{i=1}^n{\mathrm{\phi }}_{\mathrm{each}}(O_i;{\widehat{\eta }}_{k(i)}),\mathrm{and}{\widehat{\mathrm{\theta }}}_{\mathrm{joint}}=\frac{1}{n}\sum _{i=1}^n{\mathrm{\phi }}_{\mathrm{joint}}(O_i;{\widehat{\eta }}_{k(i)}).$$

(31)

Inference

Examining the estimators’ asymptotic behaviors reveals that they are asymptotically normal and multiply robust in consistency, under achievable assumptions about the estimated nuisance functions and boundedness of the weights. Below, we state these results, which then allows us to obtain inferences on ${\mathrm{\theta }}_{\mathrm{each}}$ and ${\mathrm{\theta }}_{\mathrm{joint}}$ .

Asymptotic Sampling Distributions

Assume (1) the weights $\frac{1\{T=t,R=r\}}{\mathrm{p}(t,r\mid C)}$ , ${\mathrm{w}}_{\mathrm{each}}(M_1,M_2,t,r,C)$ , and ${\mathrm{w}}_{\mathrm{joint}}(M_1,M_2,t,r,C)$ and their estimates are bounded; and (2) nuisance functions are estimated consistently with a sufficient convergence rate, such that for each estimator ${\hat{\mathrm{\theta }}}_{\mathrm{each}}$ (or ${\hat{\mathrm{\theta }}}_{\mathrm{joint}}$ ), the estimator’s second-order error term $R_2(\hat{\mathrm{\eta }},\mathrm{\eta })$ , where $\hat{\mathrm{\eta }}={\hat{\mathrm{\eta }}}_{k(i)}$ , satisfy the following: $R_2(\hat{\mathrm{\eta }},\mathrm{\eta })=o_P(n^{-1/2})$ . Supplementary Appendix C (available in the online version of this article) presents the expression of the second-order error term, which depends on only products of errors of estimated nuisance functions.

Then, letting convergence in distribution be denoted by ⇝, the estimators are $\sqrt{n}$ -consistent and asymptotically normal:

$$\begin{array}{c}\sqrt{n}({\hat{\mathrm{\theta }}}_{\mathrm{each}}-{\mathrm{\theta }}_{\mathrm{each}})⇝N(0,{\mathrm{\sigma }}_{\mathrm{each}}^2),\sqrt{n}({\hat{\mathrm{\theta }}}_{\mathrm{joint}}-{\mathrm{\theta }}_{\mathrm{joint}})⇝N(0,{\mathrm{\sigma }}_{\mathrm{joint}}^2),\hfill \end{array}$$

(32)

where ${\mathrm{\sigma }}_{\mathrm{each}}^2=\mathrm{Var}[{\mathrm{\phi }}_{\mathrm{each}}(O;\mathrm{\eta })]$ and ${\sigma }_{\mathrm{joint}}^2=\mathrm{Var}[{\mathrm{\phi }}_{\mathrm{joint}}(O;\eta )]$ are variances of the IFs. Proofs can be found in Supplementary Appendix D (available in the online version of this article).

In practice, the assumption in (1) requires the positivity assumption to hold with the data at hand. The assumption in (2) can be satisfied, for example, when the rate of estimating each nuisance function is faster than $n^{-1/4}$ , a rate achievable by many machine learning methods (e.g., regression trees, random forest, neural networks; Chernozhukov et al., 2018).

The above result enables asymptotically proper uncertainty quantification (standard errors, CIs) for the estimators even when data-adaptive approaches (e.g., model selection, machine learning) are used. It also allows us to obtain asymptotic sampling distributions of the mediated moderation estimands (Equations 8–12) via the delta method and make inferences on them. For example, for ${\mathrm{MedMod}}_{M_1}$ , the estimator can be computed as ${\mathrm{MedMod}}_{M_1}({\hat{\mathrm{\theta }}}_{\mathrm{each}})=\frac{1}{n}\sum _{i=1}^n{\mathrm{\phi }}_{{\mathrm{MedMod}}_{M_1}}(O_i;{\hat{\mathrm{\eta }}}_{k(i)})$ , where ${\mathrm{\phi }}_{{\mathrm{MedMod}}_{M_1}}(O_i;{\widehat{\eta }}_{k(i)})=[{\mathrm{\phi }}_{\mathrm{each}(1,{\mathrm{r}}_{\mathrm{foc}},{\mathrm{r}}_{\mathrm{foc}},{\mathrm{r}}_{\mathrm{ref}})}$ $(O_i;{\widehat{\eta }}_{k(i)})-{\mathrm{\phi }}_{\mathrm{each}(1,{\mathrm{r}}_{\mathrm{foc}}{\mathrm{r}}_{\mathrm{foc}},{\mathrm{r}}_{\mathrm{ref}},{\mathrm{r}}_{\mathrm{ref}})}(O_i;{\widehat{\eta }}_{k(i)})]-[{\mathrm{\phi }}_{\mathrm{each}(0,{\mathrm{r}}_{\mathrm{foc}},{\mathrm{r}}_{\mathrm{foc}},{\mathrm{r}}_{\mathrm{ref}})}(O_i;{\widehat{\eta }}_{k(i)})-{\mathrm{\phi }}_{\mathrm{each}(0,{\mathrm{r}}_{\mathrm{foc}},{\mathrm{r}}_{\mathrm{ref}},{\mathrm{r}}_{\mathrm{ref}})}(O_i;{\widehat{\eta }}_{k(i)})]$ is the contrast of function ${\mathrm{\phi }}_{\mathrm{each}}$ (Equation 29) evaluated with the different values of $(t,r,r^{*},r^{**})$ following the definition of ${\mathrm{MedMod}}_{M_1}$ . The variance of the estimator can be estimated as the sample variance of ${\mathrm{\phi }}_{{\mathrm{MedMod}}_{M_1}}(O_i;{\hat{\mathrm{\eta }}}_{k(i)})$ divided by the sample size, represented as ${\hat{\mathrm{\sigma }}}_{{\mathrm{MedMod}}_{M_1}}^2=\frac{1}{n}\mathrm{var}[{\mathrm{\phi }}_{{\mathrm{MedMod}}_{M_1}}(O_i;{\hat{\mathrm{\eta }}}_{k(i)})]$ . A 95% Wald-type CI for ${\mathrm{MedMod}}_{M_1}$ can then be constructed as ${\mathrm{MedMod}}_{M_1}({\hat{\mathrm{\theta }}}_{\mathrm{each}})\pm 1.96\sqrt{{\hat{\mathrm{\sigma }}}_{{\mathrm{MedMod}}_{M_1}}^2}$ .

The above asymptotic property, which justifies the use of the Wald-type CIs, requires that all nuisance functions can be estimated consistently at sufficiently fast rates. However, this is not required for the consistency of the estimators. As shown next, the estimators remain consistent under weaker assumptions that allow inconsistent estimation of some nuisance function.

Multiple Robustness

The estimators ${\hat{\mathrm{\theta }}}_{\mathrm{each}}$ and ${\hat{\mathrm{\theta }}}_{\mathrm{joint}}$ are multiply robust in the sense that they are consistent so long as one of multiple subsets of nuisance functions is estimated consistently. Specifically, ${\hat{\mathrm{\theta }}}_{\mathrm{each}}$ is consistent when we achieve at least one of the following: (a) the mean outcome $\mathrm{\mu }$ (Equation 22), the mediators’ joint distribution $\mathrm{p}(m_1,m_2\mid r,t,c)$ , and the marginal mediator distributions $\mathrm{p}(m_1\mid r,t,c)$ and $\mathrm{p}(m_2\mid r,t,c)$ are consistently estimated; or (b) $\mathrm{\mu }$ , the propensity score $\mathrm{p}(t,r\mid c)$ and $M_1$ ’s marginal distribution $\mathrm{p}(m_1\mid r,t,c)$ are consistently estimated; or (c) $\mathrm{\mu }$ , $\mathrm{p}(t,r\mid c)$ , and $M_2$ ’s marginal distribution $\mathrm{p}(m_2\mid r,t,c)$ are consistently estimated; or (d) $\mathrm{p}(t,r\mid c)$ , and the joint and marginal mediator distributions $\mathrm{p}(m_1,m_2\mid r,t,c)$ , $\mathrm{p}(m_1\mid r,t,c)$ , and $\mathrm{p}(m_2\mid r,t,c)$ are consistently estimated.

For ${\hat{\mathrm{\theta }}}_{\mathrm{joint}}$ , the consistency is achieved when at least one of the following is satisfied: (a′) the mean outcome $\mathrm{\mu }$ and the mediators’ joint distribution $\mathrm{p}(m_1,m_2\mid r,t,c)$ are consistently estimated; or (b′) $\mathrm{\mu }$ and the propensity score $\mathrm{p}(t,r\mid c)$ are consistently estimated; or (c′) $\mathrm{p}(t,r\mid c)$ and the joint mediator distribution $\mathrm{p}(m_1,m_2\mid r,t,c)$ are consistently estimated. Proofs for the above robustness properties of the estimators ${\hat{\mathrm{\theta }}}_{\mathrm{each}}$ and ${\hat{\mathrm{\theta }}}_{\mathrm{joint}}$ can be found in Supplementary Appendix E (available in the online version of this article).

Together, these robustness properties imply that, for consistent estimation of all the mediated moderation estimands, consistent estimation of the mediators’ joint distribution or either mediator’s marginal distribution is important. Nonetheless, consistent estimation of all the joint or marginal mediator distributions is unnecessary, when either (b) or (c) is achieved. Furthermore, consistent estimation of the mean outcome (or propensity score) is unnecessary, when (d) (or [a]) is achieved.

We provide an R package “multiMedMod,” which implements the developed method (available at https://osf.io/n78kq/?view_only=618e75fab3da4528807c6f743467c17a). In the R package, the nuisance functions are estimated using an ensemble machine learning procedure known as the “super learner” (van der Laan et al., 2007). This procedure allows a user to specify multiple machine learning prediction methods (e.g., parametric regressions, generalized additive models, random forest, neural network) and combines them to optimize cross-validated predictive performance. Tutorials of the super learner ensemble method can be found in, for example, Naimi and Balzer (2018), Polley et al. (2024), and van der Laan and Rose (2011). Our R package accommodates both binary or continuous outcome and mediators. For continuous mediators, we implement the multiply robust estimators by re-parameterizing the nuisance functions such that the estimation of mediator densities is avoided. This re-parameterization follows Rudolph et al. (2024) on causal mediation analysis. Implementation details for handling continuous mediators are provided in Supplementary Appendix H (available in the online version of this article).

Simulation Study

A simulation study is presented in Supplementary Appendix I (available in the online version of this article), where we evaluated the performance of the estimation method with binary and continuous mediators and outcomes. We examined sample sizes at n = 200, 300, 500, 800, and 1,000. Key findings are the following. First, the performance was generally better with binary mediators than with continuous mediators. This may be due to the instability of weighting terms when the mediators were continuous. Previous research also found that, owing to instability of the weights, mediation effect estimators involving weights can be prone to finite-sample bias, especially with continuous mediators and not-large sample sizes (Loh et al., 2020). Second, the bias was generally larger when the sample size was small (e.g., 200) and became closer to zero with larger sample sizes (e.g., 800 or above). The mean squared error decreased as the sample size increased. Third, the CI coverage rates ranged from 91% to 97%, which fell within the coverage rate criteria (91%–98%) suggested in Muthén and Muthén (2002). These findings suggest that the multiply robust estimators can work for making inferences on the mediated moderation estimands, but larger sample sizes would be needed, especially when working with continuous mediators.

Illustrative Example

This example is based on studies highlighting the adverse effects of perceived discrimination on adolescent mental health (Maene et al., 2022; Smith & Nicholson, 2022). Studies indicate that these effects could vary by gender, often being more severe for females than males. Additionally, psychosocial resources, such as mastery and perceived family emotional support, could be mediating mechanisms that contribute to this effect heterogeneity by gender (Erving et al., 2022). To illustrate the proposed approach, we examine how these psychosocial resources mediate gender differences in the effects of perceived discrimination on depressive symptoms among adolescents.

Example Data and Analysis

The dataset used in this example consists of $n=$  = 1,115 African American adolescents in the National Survey of American Life–Adolescent Supplement (J. S. Jackson et al., 2016; Seaton et al., 2008). The analysis includes the following variables. The subgroup status $R$ is the adolescent-reported sex, with male as the reference subgroup ( $R={\mathrm{r}}_{\mathrm{ref}}$ ) and female as the focal subgroup ( $R={\mathrm{r}}_{\mathrm{foc}}$ ). The treatment $T$ is the adolescent’s perceived discrimination from teacher, coded based on the count of the adolescent-reported perceived teacher discrimination events in day-to-day life (e.g., “your teachers treat you with less respect than other students”; $T=0$ for a zero count, and $T=1$ for one or more counts). The outcome $Y$ is the adolescent’s depressive symptoms, measured as a composite score of the Epidemiologic Studies Depression Scale (Radloff, 1977). The mediator $M_1$ is a binary indicator of mastery, coded based on the adolescent’s agreement with “having little control over the things that happen to me” (coded as 1 for somewhat or strongly agree, and 0 otherwise). The mediator $M_2$ is a binary indicator of perceived family emotional support, coded based on the adolescent’s response to “how often do your family members make you feel loved and cared for” (coded as 1 for very often, and 0 otherwise). The covariates $C$ include background information about the adolescent’s age, family income, school program type (general, academic, vocational, or other), urbanity, whether the adolescent ever repeated a grade, whether ever suspended from school for a day or longer, whether felt family has enough money, number of places ever lived, family composition (numbers of children and adults in the household, whether raised by mother, whether raised by father), and region of the country. The variables were selected based on related substantive literature (Miller et al., 2013).

We used the developed method to estimate the mediated moderation for the two mediators. Specifically, we used four-fold cross-fitting. The nuisance functions were estimated by a super learner ensemble of methods in the R package SuperLearner (Polley et al., 2024; van der Laan et al., 2007), with an ensemble of an intercept-only model, generalized linear model, lasso regression, multivariate adaptive regression splines (Friedman, 1991), random forest (Wright & Ziegler, 2017), and neural network (Ripley et al., 2016). The implementation was done using the R package accompanying this manuscript. The dataset and R code used in this example are available at https://osf.io/n78kq/?view_only=618e75fab3da4528807c6f743467c17a. The results are shown in Table 1 and Figure 3.

Table 1.

Mediated Moderation Analysis Results in the Example

Effect	Interpretation	Estimate	95% CI
$\mathrm{TotMod}=\mathrm{RemainMod}+{\mathrm{MedMod}}_{\mathrm{joint}}$
$\mathrm{TotMod}$	Total Moderation	1.51	[0.14, 2.88]
$\mathrm{RemainMod}$	Remaining moderation (not via $\{M_1,M_2\}$ )	0.95 (63% of $\mathrm{TotMod}$ )	[−0.38, 2.29]
${\mathrm{MedMod}}_{\mathrm{joint}}$	Mediated moderation via both mediators $\{M_1,M_2\}$ jointly	0.56 (37% of $\mathrm{TotMod}$ )	[0.08, 1.04]
${\mathrm{MedMod}}_{\mathrm{joint}}={\mathrm{MedMod}}_{M_1}+{\mathrm{MedMod}}_{M_2}+{\mathrm{MedMod}}_{\mathrm{mu}}$
${\mathrm{MedMod}}_{M_1}$	Mediated moderation via mastery (mediator $M_1$ )	0.38	[0.04, 0.72]
${\mathrm{MedMod}}_{M_2}$	Mediated moderation via perceived family emotional support (mediator $M_2$ )	0.27	[−0.02, 0.57]
${\mathrm{MedMod}}_{\mathrm{mu}}$	Mediated moderation due to mutual dependence of the mediators	−0.09	[−0.24, 0.06]

Note. CI: confidence interval.

Figure 3.

Estimates of expected potential outcomes between the subgroups in the example.

Results of the Illustrative Example

The total moderation (TotMod), quantifying the difference in the impact of perceived teacher discrimination on depressive symptoms between males and females, is estimated at 1.51 (95% CI [0.14, 2.88]). The total moderation is decomposed into the mediated moderation and the remaining moderation, $\mathrm{TotMod}=\mathrm{RemainMod}+\mathrm{MedMo}{\mathrm{d}}_{\mathrm{joint}}$ . The remaining moderation ( $\mathrm{RemainMod}$ ), corresponding to the part of this gender-based moderation not mediated by the two psychosocial resources (i.e., mastery $M_1$ and perceived family emotional support $M_2$ ) is estimated at 0.95 (95% CI [−0.38, 2.29]). This indicates that the gender difference in the impact would become less pronounced (with the CI covering zero), if, within strata of the covariates, the potential levels of both psychosocial resources (i.e., the mediators) were equalized across gender. The jointly mediated moderation ( ${\mathrm{MedMod}}_{\mathrm{joint}}$ ), corresponding to the part of the gender-based moderation attributable both mediators jointly, is estimated at 0.56 (95% CI [0.08, 1.04]), or 37% of the total moderation. This suggests that the adverse impact for females would reduce by 0.56, if, within strata of the covariates, the potential levels of both psychosocial resources among them were jointly made to match those among males.

Next, we further decompose ${\mathrm{MedMod}}_{\mathrm{joint}}$ to examine the role of each of the two psychosocial resources (i.e., mastery $M_1$ and perceived family emotional support $M_2$ ) in contributing to the effect moderation by gender. First, we estimate that the mediated moderation via mastery ( ${\mathrm{MedMod}}_{{\mathrm{M}}_1}$ ) is 0.38 (95% CI [0.04, 0.72]), accounting for 25% of the total moderation. This suggests that among females, the adverse impact of perceived teacher discrimination on depressive symptoms would decrease by approximately 0.38 if, within strata of the covariates, the potential levels of mastery ( $M_1$ ) among them were made equivalent to those among males, while keeping the potential levels of perceived family emotional support ( $M_2$ ) held at the levels among the males. Second, we estimate that the mediated moderation via perceived family emotional support ( ${\mathrm{MedMod}}_{M_2}$ ) is 0.27 (95% CI [−0.02, 0.57]), or 18% of the total moderation. This indicates that compared to mastery ( $M_1$ ), perceived family emotional support ( $M_2$ ) appears to play a less significant role in contributing to the gender-based moderation. Lastly, we found that, at 0.05 significance level, there is no statistically significant mediated moderation due to mutual dependence of these two psychosocial resources ( ${\mathrm{MedMod}}_{\mathrm{mu}}$ : estimate = −0.09, 95% CI [−0.24, 0.06]).

Figure 3 displays the expected potential outcomes estimated in the mediated moderation analysis. The estimates show that, under non-exposure to perceived teacher discrimination ( $T=0$ ; left panel), there is little change in the expected potential depressive symptom outcome among the females, regardless of whether or not shifting the potential levels of the two mediators among them. In contrast, under exposure to perceived teacher discrimination ( $T=1$ ; right panel), there would be an appreciable reduction in the expected potential depressive symptom outcome among the females, if the potential levels of the psychosocial resources (mastery $M_1$ and perceived family emotional support $M_2$ ) among them were changed to equal those among the males, within strata of the covariates.

Taken together, these results suggest that interventions aimed at improving psychosocial resources (mastery and perceived family emotional support) could help mitigate the adverse impact of perceived teacher discrimination among female adolescents. If the potential levels of both psychosocial resources (mediators) among the females were jointly made to match those among the males (within strata of the covariates), the gender-based difference in the effect of perceived teacher discrimination on depressive symptoms would reduce (by an estimated 37%).

Discussion

Assessing mediated moderation allows addressing important questions about the mechanisms that lead to effect moderation between different subgroups (Fairchild & MacKinnon, 2009; Liu et al., 2025; Muller et al., 2005). However, methods that allow clarifying causal interpretations of mediated moderation estimates formally have been underdeveloped, especially when more than one mediator is of interest. This study reduces this gap by developing causal estimands and identification results for assessing the roles of two mediators in the effect moderation, and by developing robust estimation methods allowing the use of machine learning approaches.

Particularly, this study builds upon and extends three major streams of causal mediation literature. First, the causal estimands defined here extend the counterfactual approach to quantifying mediation (VanderWeele, 2015; Vansteelandt & Daniel, 2017). In doing so, our estimands for mediated moderation do not rely on specific parametric models. This improves upon conventional parametric modeling approaches to defining mediated moderation, as the latter requires modeling assumptions that may or may not be plausible in practice. Second, our decomposition approach builds on recent work on interventional effects (J. W. Jackson & VanderWeele, 2018; Liu et al., 2025; Vansteelandt & Daniel, 2017), and extends this literature by decomposing the effect moderation with respect to two mediators, yielding new interventional effect-based estimands for mediated moderation. The resulting estimands make it possible to test the role of each mediator in generating the subgroup differences in treatment effects. Third, our estimation method extends a fast-growing literature on doubly robust machine learning approaches in causal inference (Chernozhukov et al., 2018; van der Laan & Rose, 2018) to the mediated moderation analysis. The developed method allows relaxing modeling assumptions and improves robustness to the estimation errors from machine learning methods.

Several limitations of the current study present opportunities for future research. First, by the unconfoundedness assumption (A2), we have assumed no posttreatment confounders of the mediators–outcome relation (or known as treatment-induced confounders; VanderWeele et al., 2014), within strata of $(T,R,C)$ . To handle posttreatment confounders in the identification, (A2) could be modified by including them in the conditioning set. The corresponding estimation methods would require future developments, which may be built on recent estimation methods considering both posttreatment confounders and multiple mediators (Mittinty et al., 2019). Furthermore, it is important to develop methods for assessing the sensitivity of the mediated moderation results to unmeasured pretreatment or posttreatment confounders. Previous literature has provided parametric and nonparametric sensitivity analysis methods for mediation analysis (Ding & Vanderweele, 2016; Qin & Yang, 2022; Zhang & Ding, 2023), and future studies could examine the feasibility and potential extensions of these methods for mediated moderation.

Second, the interventional effects approach in the current study could have limitations for defining mediation moderation effects. In our current approach, all covariates $C$ are used for both the definition and identification of mediation moderation estimands. This essentially treats all pretreatment covariates that serve for confounding control as allowable to be conditioned on for defining the between-subgroup difference in the treatment effect. However, this covariate allowability judgment may not necessarily align with the judgment in all substantive contexts. For example, a healthcare researcher examining effect moderation by race may judge some covariates in $C$ (e.g., SES) as non-allowable for defining racial differences (J. W. Jackson, 2021). Therefore, an important future direction could be to distinguish between allowable and non-allowable covariates for defining causal estimands of mediation moderation, such as by incorporating a meaningful causal decomposition perspective (J. W. Jackson, 2021).

Third, in the developed multiply robust method, we have computed the estimators’ variances using empirical values of ${\mathrm{\phi }}_{\mathrm{each}}(O;\hat{\mathrm{\eta }})$ and ${\mathrm{\phi }}_{\mathrm{joint}}(O;\hat{\mathrm{\eta }})$ , which do not consider uncertainty of estimated nuisance functions $\hat{\mathrm{\eta }}$ . While this is a standard variance estimation approach for IF-based estimators (Benkeser & Ran, 2021; Kennedy, 2024), it could yield anti-conservative variance estimates as sample sizes became smaller or violations of positivity assumptions became severer (Tran et al., 2023; van der Laan & Gruber, 2011). Future research could examine methods (e.g., modified bootstrap-based methods in Tran et al., 2023) to improve the uncertainty quantification (variance estimation, CI construction) for the estimators of the mediated moderation estimands.

Fourth, extending the analysis to more than two mediators could be a valuable direction. For example, Loh et al. (2022b) defined IIEs with high-dimensional mediators, presenting an effect definition approach that may be explored for mediated moderation analyses. Furthermore, accommodating continuous or multivariate moderator variables in the causal estimand definition, identification, and estimation for the mediated moderation analysis awaits future studies. Toward this end, studies could explore methods for continuous or multivariate exposures (Hejazi et al., 2023) and moderators (Valente et al., 2023), as well as approaches for estimating heterogeneous treatment effects given high-dimensional moderators (Athey & Imbens, 2019; Kennedy, 2023). Additionally, we have focused on comparing pretreatment subgroups. Nonetheless, the controlled descriptive perspective (which we adopted for subgroup comparisons) could be extended to the latent subgroups defined using potential status of a posttreatment variable in the principal stratification framework (Frangakis & Rubin, 2002; Page et al., 2015); in that case, a controlled descriptive comparison could entail comparing those latent subgroups, while conditioning on pretreatment covariates.

Conclusion

In conclusion, this study provides a new tool for researchers investigating mediators underlying effect heterogeneity between subgroups. We hope that by enabling the investigation from a causal inference approach, our method can help researchers gain insights on the development or adaptations of interventions for specific subgroups.

Supplemental Material

sj-pdf-1-jeb-10.3102_10769986251410929 – Supplemental material for A Causal Inference Approach for Mediated Moderation with Multiple Mediators

Supplemental material, sj-pdf-1-jeb-10.3102_10769986251410929 for A Causal Inference Approach for Mediated Moderation with Multiple Mediators by Xiao Liu and Monica E. Romero in Journal of Educational and Behavioral Statistics