Estimating Causal Mediation Effects in Multiple-Mediator Analyses With Clustered Data

Abstract

Multiple-mediator analyses with clustered data are common in educational and behavioral sciences, but limited methods exist to assess the causal mediation effects via each of multiple mediators. In this study, we extend the multiply robust method to make inferences on the causal mediation effects for two mediators with clustered data. The developed method takes into account unmeasured cluster-level confounders and can incorporate machine learning methods to nonparametrically estimate nuisance models while allowing uncertainty quantification via asymptotic standard errors and confidence intervals. We conduct simulations to evaluate the developed method for inference of both the individual-average and cluster-average causal mediation effects with clustered data. We illustrate our method using data from the Education Longitudinal Study.

Mediation analysis is an important tool for assessing causal mechanisms by which treatment affects outcomes (Baron & Kenny, 1986; VanderWeele, 2015). The causal mechanisms often involve multiple mediators, and quantifying the mediating role of each mediator is frequently of interest (Vansteelandt & Daniel, 2017).

For multiple-mediator analysis, conventional methods were based on the parametric modeling approach (Baron & Kenny, 1986; MacKinnon, 2008). Despite its popularity, the parametric modeling approach to mediation has limitations in the causal interpretation of the effect estimates. For causal interpretation in mediation analysis, the causal inference approach is increasingly popular (VanderWeele, 2015), which uses a counterfactual framework to define and identify causal estimands for mediation analysis, such as the natural indirect effects (Pearl, 2001; Robins & Greenland, 1992) and the interventional indirect effects (IIEs; VanderWeele et al., 2014; Vansteelandt & Daniel, 2017). Recently, the causal inference approach has been extended to multiple mediators (Hong, 2015; VanderWeele & Vansteelandt, 2014). For example, Qin et al. (2021) assessed the natural indirect effects for mediators that are assumed to be independent given pretreatment covariates. Daniel et al. (2015) and Imai and Yamamoto (2013) considered causal mediation for sequential mediators (i.e., one mediator affects the other), assuming no unmeasured common causes for the mediators.

In this study, we focus on a different, yet common scenario of multiple mediators, where (a) there is no substantive knowledge of the causal ordering between mediators and (b) the mediators may be dependent. To illustrate, in Figure 1, under each value of treatment $A$ , the mediators may be dependent even after conditioning pretreatment covariates $X$ . The mediators’ dependence may be due to a causal relationship (such as in [a] and [b]), sharing an unmeasured common cause (such as in [c]), or both (such as in [d]).

Figure 1.

Causal diagrams with dependent mediators. The data are unclustered. (a) M₁ affects M₂. (b) M₂ affects M₁. (c) M₁ and M₂ share an unmeasured common cause V. (d) M₁ and M₂ have a casual relationship and share an unmeasured common cause V.

For this multiple-mediator scenario (Figure 1), a practical approach to causal mediation analysis has been the interventional effects approach (Vansteelandt & Daniel, 2017). While the natural indirect effects are often the estimands with a single mediator (Imai et al., 2010), with multiple mediators the identification of natural indirect effects for each mediator requires strong assumptions (Avin et al., 2005; Loh et al., 2022). For example, under Figure 1(c) (mediators correlated due to an unmeasured common cause), the natural indirect effect for each mediator is generally unidentifiable (Vansteelandt & Daniel, 2017). To address this problem, Vansteelandt and Daniel (2017) proposed the interventional effects approach for defining causal mediation effects with multiple mediators, which requires weaker assumptions for identification; for example, the IIEs remain identifiable under Figure 1(c). As we describe in the next section, the interventional effects approach considers hypothetical interventions on mediators, which helps guide intervention development (Moreno-Betancur & Carlin, 2018) and is increasingly popular in causal mediation analyses in applied research (X. Liu & Wang, 2024; X. Liu et al., 2025; Loh & Ren, 2022; Rubinstein, Haviland, & Breslau, 2023).

The purpose of this study is to examine methods for estimating the interventional effects with clustered data (e.g., students clustered within schools), a data structure frequently occurring in multiple-mediator analyses in behavioral sciences and can add complexities to causal inference. For example, the Education Longitudinal Study (ELS) has served as a major data source in multiple-mediator analyses for the effect of extracurricular participation ( $A$ ) on adolescent academic outcome ( $Y$ ) via the adolescent’s educational expectations and beliefs ( $M_1$ and $M_2$ ), in which adolescents are clustered within schools (Ingels et al., 2007; Morris, 2016). The variables ( $A$ , $M_1$ , $M_2$ , $Y$ ) are all measured at the adolescent level, and could be affected by pretreatment covariates at not only the adolescent level (e.g., gender and race) but also at the school level (e.g., school climate).

From previous literature, clustering adds complexities to effect estimation, not only due to correlations within clusters but more importantly, due to the possibility of unmeasured cluster-level confounders (Chang & Stuart, 2022; Li et al., 2013). A key insight from the literature is that, methods adequately utilizing the clustered structure (e.g., by conditioning on clusters or cluster-mean centering) can protect against bias due to unmeasured cluster-level confounders (Suk & Kang, 2023; Zetterqvist et al., 2016). However, to our knowledge, the literature has not addressed handling clustering in the estimation of causal mediation effects in the multiple-mediator analysis (Figure 1), despite that multiple-mediator analyses with clustered data are frequently involved in the research process of educational and behavioral sciences.

Therefore, in this study, we develop methods to assess the interventional (in)direct effects of the multiple-mediator analysis (Figure 1) with clustered data. Specifically, we focus on a binary treatment and two discrete (binary) mediators. We consider that each individual’s cluster membership is an observed pretreatment characteristic. In the rest of this article, we first define the causal estimands and present the identification results. We then describe the estimation methods. Subsequently, we conduct simulations to evaluate the methods. We illustrate the application of the developed methods using data from the ELS.

Definition and Identification of Causal Mediation Effects

We start with describing the sampling process that generates the clustered sample. In the first stage, a simple random sample of $J$ clusters (e.g., schools) is selected from some study population; the selected clusters are indexed by $j=1,...,J$ . Each cluster serves a source population of individuals eligible for the study (e.g., adolescents). In the second stage, given that cluster $j$ ( $j=1,...,J$ ) is selected in the first stage, a simple random sample of $N_j$ individuals is drawn from all eligible individuals within the cluster $j$ . The cluster size $N_j$ is set to represent the cluster’s relative population size, that is, $N_j$ is either the cluster’s population size or is set proportional to the cluster’s population size (Raudenbush & Bloom, 2015). With the sampling process, we have a clustered sample of total size $\sum _{j=1}^J{N}_j$ . Let $S_j$ indicate each individual’s cluster membership, which is $S_j=j$ for individuals in cluster $j$ . Note that we have assumed simple random sampling in both stages: that is, we have equal sampling probabilities of clusters (the first stage), and given that a cluster is sampled, we have equal sampling probabilities of individuals within the cluster (the second stage); in the Supplemental Materials (Appendix A in the online version of the journal), we consider scenarios with unequal sampling probabilities in the first or second stages, and describe how our method may be modified by incorporating user-specified sampling weights in such scenarios.

After the sampling, baseline covariates are measured for each sampled cluster. The covariates have two types: (1) Measured cluster-level covariates $W_j$ on cluster $j$ (e.g., school type) that are shared by all individuals within the cluster, and (2) measured individual-level covariates $X_{\mathrm{ij}}$ that vary within the cluster (e.g., an individual’s demographics). Then, each individual selects or is assigned to either the treatment condition $A_{\mathrm{ij}}=1$ or the control condition $A_{\mathrm{ij}}=0$ . The treatment assignment may be affected by the cluster size (if $N_j$ varies), measured cluster-level covariates, unmeasured cluster-level confounder $U_j$ (e.g., school environment), and the individual’s covariates; that is, $A_{\mathrm{ij}}$ depends on $(X_{\mathrm{ij}},N_j,W_j,U_j)$ . In a follow-up, the two mediators $({M_1}_{\mathrm{ij}},{M_2}_{\mathrm{ij}})$ are measured on each individual; the mediators may be affected by the individual’s own treatment and pretreatment covariates, the cluster size, and measured/unmeasured cluster-level pretreatment covariates; that is, $({M_1}_{\mathrm{ij}},{M_2}_{\mathrm{ij}})$ depends on $(A_{\mathrm{ij}},X_{\mathrm{ij}},N_j,W_j,U_j)$ . Finally, the outcome $Y_{\mathrm{ij}}$ is measured on each individual; $Y_{\mathrm{ij}}$ may depend on the individual’s own treatment, mediators, and pretreatment covariates, as well as the cluster size and measured/unmeasured cluster-level pretreatment covariates; that is, $Y_{\mathrm{ij}}$ depends on $(A_{\mathrm{ij}},{M_1}_{\mathrm{ij}},{M_2}_{\mathrm{ij}},X_{\mathrm{ij}},N_j,W_j,U_j)$ .

To summarize, for each sampled cluster $j$ , we observed the cluster size, measured cluster-level covariates, and measured individual-level covariates, treatment, mediators, and outcome of individuals sampled within the cluster, written as $O_j=(S_j=j,N_j,W_j,\{O_{\mathrm{ij}}{\}}_{i=1}^{N_j})$ with $O_{ij}=\left(X_{ij},A_{ij},M_{1ij},M_{2ij},Y_{ij}\right)$ being the set of measured individual-level variables for each individual. We assume that the data are $J$ independently and identically distributed observations $\{O_j{\}}_{j=1}^J$ of $O_j$ . Furthermore, we assume that the cluster size $N_j$ is finite, following some distribution with finite support on positive integers (Wang et al., 2024).

Note that in our setting, the treatment variable is an individual-level variable varying within each cluster; this is to be distinguished from clustered designs with a cluster-level treatment variable (e.g., cluster randomized trials), in which all individuals in the same cluster are assigned to the same treatment value. Our setting resembles a blocked experiment that uses higher-level sampling units (e.g., schools) as blocks (often referred to as random blocks, e.g., Dong & Maynard, 2013 or structural blocks, e.g., Pashley & Miratrix, 2022) or multisite experiment (Raudenbush & Bloom, 2015), where the experimenter allocates individuals to treatment and control conditions within each cluster (within each block or within each site).

Figure 2 shows the relationships of the variables, where the mediators are dependent. The dependence structure of the mediators is unknown, analogous to Figure 1 in the unclustered setting.

Figure 2.

Causal diagrams with dependent mediators. The data are clustered. (a) ${{\mathrm{M}}_1}_{\mathrm{ij}}$ affects ${{\mathrm{M}}_2}_{\mathrm{ij}}$ . (b) ${{\mathrm{M}}_2}_{\mathrm{ij}}$ affects ${{\mathrm{M}}_1}_{\mathrm{ij}}$ . (c) ${{\mathrm{M}}_1}_{\mathrm{ij}}$ and ${{\mathrm{M}}_2}_{\mathrm{ij}}$ share an unmeasured common cause $\mathrm{V}=\{{\mathrm{V}}_{\mathrm{j}},{\mathrm{V}}_{\mathrm{ij}}\}.$ (d) ${{\mathrm{M}}_1}_{\mathrm{ij}}$ and ${{\mathrm{M}}_2}_{\mathrm{ij}}$ have a casual relationship and share an unmeasured common cause $\mathrm{V}=\{{\mathrm{V}}_{\mathrm{j}},{\mathrm{V}}_{\mathrm{ij}}\}.$

Causal Estimands

To define the effects, consider the counterfactual or potential outcomes (Rubin, 1974, 1978, 1980) for an individual. Let $Y_{\mathrm{ij}}\left(a\right)$ denote the potential outcome that an individual would have if the treatment is set to value $a$ . Furthermore, incorporating the mediators, let $Y_{\mathrm{ij}}(a,m_1,m_2)$ denote the potential outcome if the vector of treatment and mediators $(A_{\mathrm{ij}},{M_1}_{\mathrm{ij}},{M_2}_{\mathrm{ij}})$ is set to value $(a,m_1,m_2)$ . For the two mediators ${M_1}_{\mathrm{ij}}$ and ${M_2}_{\mathrm{ij}}$ , let ${M_1}_{\mathrm{ij}}\left(a\right)$ and ${M_2}_{\mathrm{ij}}\left(a\right)$ be respectively the potential mediators if the treatment is set to value $a$ . We make the stable unit treatment value assumption (SUTVA; Rubin, 1980), which assumes that there are no different versions of each treatment value and each mediator value and that there is no interference among individuals. With clustered data, the no-interference assumption requires considerations given specific substantive contexts (Hong & Raudenbush, 2006). In this study, we assume the SUTVA throughout, and relaxing this assumption (e.g., allowing for interference within clusters) is an important future research topic. The SUTVA allows us to connect the potential variables to the observed data as follows: The observed mediators ${M_1}_{\mathrm{ij}}$ and ${M_2}_{\mathrm{ij}}$ equal the potential mediators ${M_1}_{\mathrm{ij}}\left(A_{\mathrm{ij}}\right)$ and ${M_2}_{\mathrm{ij}}\left(A_{\mathrm{ij}}\right)$ under the observed treatment; the observed outcome $Y_{\mathrm{ij}}$ equals the potential outcome $Y_{\mathrm{ij}}(A_{\mathrm{ij}},{M_1}_{\mathrm{ij}},{M_2}_{\mathrm{ij}})$ under the observed treatment and mediators.

In what follows, we first define the interventional (in)direct effects (Vansteelandt & Daniel, 2017) for a specific cluster; then, aggregating these cluster-specific effects, we define the cluster-average and individual-average versions of the effects.

Cluster-Specific Causal Mediation Effects

For a specific cluster $j$ (e.g., school $j$ ), the total effect of treatment on the outcome is

$$\mathrm{TE}\left(j\right)=\frac{1}{N_j}\sum _{i=1}^{N_j}\mathrm{T}{\mathrm{E}}_{\mathrm{ij}}=\frac{1}{N_j}\sum _{i=1}^{N_j}[Y_{\mathrm{ij}}\left(1\right)-Y_{\mathrm{ij}}\left(0\right)],$$

(1)

which is the average difference between each individual’s potential outcomes under treatment versus control conditions, averaged among the $N_j$ individuals in cluster $j$ . We can equivalently write this total effect by incorporating the mediators (VanderWeele & Vansteelandt, 2009): $\mathrm{TE}\left(j\right)=\frac{1}{N_j}\sum _{i=1}^{N_j}[Y_{\mathrm{ij}}(1,{M_1}_{\mathrm{ij}}\left(1\right),{M_2}_{\mathrm{ij}}\left(1\right))-Y_{\mathrm{ij}}(0,{M_1}_{\mathrm{ij}}\left(0\right),{M_2}_{\mathrm{ij}}\left(0\right))]$ . How much do the mediators play a role in mediating the total effect? We address this next.

We first define the IIE (Vansteelandt & Daniel, 2017) via both mediators $({M_1}_{\mathrm{ij}},{M_2}_{\mathrm{ij}})$ jointly, and the interventional direct effect (IDE) not through the mediators. To ease introducing the definition, consider a hypothetical intervention in cluster $j$ , which would assign the mediators $({M_1}_{\mathrm{ij}},{M_2}_{\mathrm{ij}})$ for each individual to a random draw [denoted as ${\tilde{M}}_{\mathrm{join}{\mathrm{t}}_{\mathrm{ij}}}(a)$ ] from the mediators’ potential joint probability distribution $p({M_1}_{\mathrm{ij}}\left(a\right),{M_2}_{\mathrm{ij}}\left(a\right)\mid X_{\mathrm{ij}},S_j=j)$ under treatment value $a$ , given the individual’s pretreatment characteristics $X_{\mathrm{ij}}$ and cluster membership $S_j=j$ . Then, we define the $\mathrm{IDE}$ and joint indirect effect ( $\mathrm{II}{\mathrm{E}}_{\mathrm{joint}}$ ) for cluster $j$ as follows:

$$\mathrm{IDE}\left(j\right)=\frac{1}{N_j}\sum _{i=1}^{N_j}\mathrm{ID}{\mathrm{E}}_{\mathrm{ij}}=\frac{1}{N_j}\sum _{i=1}^{N_j}[Y_{\mathrm{ij}}(1,{\tilde{M}}_{\mathrm{join}{\mathrm{t}}_{\mathrm{ij}}}\left(0\right))-Y_{\mathrm{ij}}(0,{\tilde{M}}_{\mathrm{join}{\mathrm{t}}_{\mathrm{ij}}}\left(0\right))]$$

(2)

$$\mathrm{II}{\mathrm{E}}_{\mathrm{joint}}\left(j\right)=\frac{1}{N_j}\sum _{i=1}^{N_j}\mathrm{II}{\mathrm{E}}_{\mathrm{joint},\mathrm{ij}}=\frac{1}{N_j}\sum _{i=1}^{N_j}[Y_{\mathrm{ij}}(1,{\tilde{M}}_{\mathrm{join}{\mathrm{t}}_{\mathrm{ij}}}\left(1\right))-Y_{\mathrm{ij}}(1,{\tilde{M}}_{\mathrm{join}{\mathrm{t}}_{\mathrm{ij}}}\left(0\right))].$$

(3)

The direct effect $\mathrm{IDE}\left(j\right)$ is the cluster-specific average effect of assigning each individual in the cluster to treatment versus control, while holding the mediators’ joint probability at what it would be under control (given the individual’s covariates and cluster membership). The joint indirect effect $\mathrm{II}{\mathrm{E}}_{\mathrm{joint}}\left(j\right)$ is the cluster-specific average effect of changing the joint probability of both mediators for each individual in the cluster, from what it would be under control to under treatment (given the individual’s covariates and cluster membership) while assigning the individual to the treatment $a=1$ . These two effects decompose the conditional expectation of the total effect as follows: Given $X_{\mathrm{ij}}$ and $S_j=j$ , the expected total effect for an individual in cluster $j$ is $\mathrm{E}[\mathrm{T}{\mathrm{E}}_{\mathrm{ij}}\mid X_{\mathrm{ij}},S_j=j]=\mathrm{E}[\mathrm{ID}{\mathrm{E}}_{\mathrm{ij}}+\mathrm{II}{\mathrm{E}}_{\mathrm{joint},\mathrm{ij}}\mid X_{\mathrm{ij}},S_j=j]$ ,where the expectation is over the mediators’ distribution $p({M_1}_{\mathrm{ij}}\left(a\right),{M_2}_{\mathrm{ij}}\left(a\right)\mid X_{\mathrm{ij}},S_j=j)$ .

Next, we define the IIEs for each mediator. To do so, consider a hypothetical intervention in cluster $j$ , which would assign the two mediators ${M_1}_{\mathrm{ij}}$ and ${M_2}_{\mathrm{ij}}$ for an individual as two independent random draws, denoted by ${\tilde{M}}_{1\mathrm{ij}}\left(a\right)$ and ${\tilde{M}}_{2\mathrm{ij}}\left(a\right)$ , respectively. ${\tilde{M}}_{1\mathrm{ij}}\left(a\right)$ is a random draw from $p({M_1}_{\mathrm{ij}}\left(a\right)\mid X_{\mathrm{ij}},S_j=j)$ , the potential probability distribution of mediator ${M_1}_{\mathrm{ij}}$ under treatment value $a$ given the covariates and cluster membership $S_j=j$ . Analogously, ${\tilde{M}}_{2_{\mathrm{ij}}}\left(a\right)$ is a random draw from $p({M_2}_{\mathrm{ij}}\left(a\right)\mid X_{\mathrm{ij}},S_j=j)$ . Then, we define the IIE of treatment on outcome via $M_{1\mathrm{ij}}$ for cluster $j$ as

$$\mathrm{II}{\mathrm{E}}_{M_1}\left(j\right)=\sum _{i=1}^{N_j}\frac{\mathrm{II}{\mathrm{E}}_{M_1,\mathrm{ij}}}{N_j}=\frac{1}{N_j}\sum _{i=1}^{N_j}[Y_{\mathrm{ij}}(1,{\tilde{M}}_{1\mathrm{ij}}\left(1\right),{\tilde{M}}_{2\mathrm{ij}}\left(0\right))-Y_{\mathrm{ij}}(1,{\tilde{M}}_{1\mathrm{ij}}\left(0\right),{\tilde{M}}_{2\mathrm{ij}}\left(0\right))].$$

(4)

This effect via $M_1$ (4) is the cluster-specific average effect of changing the probability of mediator ${M_1}_{\mathrm{ij}}$ for each individual in the cluster, from what it would be under control to what it would be under treatment (given the individual’s covariates and cluster membership) while assigning the individual to the treatment $a=1$ and holding the probability of the mediator ${M_2}_{\mathrm{ij}}$ at what it would be under control. This effect (4) captures the part of the treatment effect mediated by mediator $M_1$ alone, but not by the mediator’s causal descendants (i.e., variables affected by the mediator). To illustrate, consider the diagrams in Figure 2. The IIE via $M_1$ captures the pathway $A\to M_1\to Y$ in Figure 2(a), (c), and (d); while in Figure 2(b), the IIE via $M_1$ captures the combination of the pathways $A\to M_1\to Y$ and $A\to M_2\to M_1\to Y$ .

Corresponding to ${\mathrm{IIE}}_{M_1}\left(j\right)$ , we define the IIE of treatment on outcome via mediator ${M_2}_{\mathrm{ij}}$ for cluster $j$ as

$$\mathrm{II}{\mathrm{E}}_{M_2}\left(j\right)=\sum _{i=1}^{N_j}\frac{\mathrm{II}{\mathrm{E}}_{M_2,\mathrm{ij}}}{N_j}=\frac{1}{N_j}\sum _{i=1}^{N_j}[Y_{\mathrm{ij}}(1,{\tilde{M}}_{1\mathrm{ij}}\left(1\right),{\tilde{M}}_{2\mathrm{ij}}\left(1\right))-Y_{\mathrm{ij}}(1,{\tilde{M}}_{1\mathrm{ij}}\left(1\right),{\tilde{M}}_{2\mathrm{ij}}\left(0\right))].$$

(5)

This effect via $M_2$ (5) is the effect of changing the probability of ${M_2}_{\mathrm{ij}}$ for each individual in cluster $j$ from what it would be under control to treatment while assigning the individual to the treatment and fixing the probability of ${M_1}_{\mathrm{ij}}$ at the distribution under treatment (given the individual’s covariates and cluster membership). In Figure 2, the IIE via $M_2$ captures the pathway $A\to M_2\to Y$ in the diagrams of (b) and (c), and captures the combination of the pathways $A\to M_2\to Y$ and $A\to M_1\to M_2\to Y$ in the diagrams of (a) and (d).

Finally, we define the IIE due to the mediators’ mutual dependence (Vansteelandt & Daniel, 2017), as the difference between the joint indirect effect (3) and the two separate indirect effects through each mediator (4) and (5):

$$\begin{array}{c}\multicolumn{1}{c}{\mathrm{II}{\mathrm{E}}_{\mathrm{mu}}\left(j\right)=\frac{1}{N_j}\sum _{i=1}^{N_j}\mathrm{II}{\mathrm{E}}_{\mathrm{joint},\mathrm{ij}}-\mathrm{II}{\mathrm{E}}_{M_1,\mathrm{ij}}-\mathrm{II}{\mathrm{E}}_{M_2,\mathrm{ij}}}\\ =\frac{1}{N_j}\sum _{i=1}^{N_j}[\{Y_{\mathrm{ij}}(1,{\tilde{M}}_{{\mathrm{joint}}_{\mathrm{ij}}}\left(1\right))-Y_{\mathrm{ij}}(1,{\tilde{M}}_{1\mathrm{ij}}\left(1\right),{\tilde{M}}_{2\mathrm{ij}}\left(1\right))\}\hfill \\ -\{Y_{\mathrm{ij}}(1,{\tilde{M}}_{{\mathrm{joint}}_{\mathrm{ij}}}\left(0\right))-Y_{\mathrm{ij}}(1,{\tilde{M}}_{1\mathrm{ij}}\left(0\right),{\tilde{M}}_{2\mathrm{ij}}\left(0\right))\}].\hfill \end{array}$$

(6)

This (6) expresses the effect of changing the strength of dependence between the mediators from under control to treatment, given the covariates $X_{\mathrm{ij}}$ and cluster membership $S_j=j$ . $\mathrm{II}{\mathrm{E}}_{\mathrm{mu}}\left(j\right)$ would be zero in expectation, when the potential mediators are independent given $X_{\mathrm{ij}}$ and cluster membership $S_j=j$ . Nonetheless, when the potential mediators are dependent, $\mathrm{II}{\mathrm{E}}_{\mathrm{mu}}\left(j\right)$ can be nonzero in expectation and capture part of the mediation pathways. In what follows, we collectively refer to the IDE (2), indirect effects via each mediator (4) and (5), and indirect effect due to mediators’ mutual dependence (6), as “causal mediation effects.”

Individual- and Cluster-Average Effects

With the clustered data, an average effect can be defined for two types of average units (Miratrix et al., 2021; Raudenbush & Bloom, 2015): (i) The average effect for an average cluster (i.e., weighting clusters equally), or the cluster-average effect and (ii) the average effect for an average individual (i.e., weighting individuals equally), or the individual-average effect. Specifically, let $\mathrm{Effect}\left(j\right)$ be one of the cluster-specific causal mediation effects, $\mathrm{Effect}\left(j\right)\in \{\mathrm{IDE}\left(j\right),\mathrm{II}{\mathrm{E}}_{M_1}\left(j\right),\mathrm{II}{\mathrm{E}}_{M_2}\left(j\right),\mathrm{II}{\mathrm{E}}_{\mathrm{mu}}\left(j\right)\}$ . Then, we can express the cluster-average effect, $\mathrm{Effect}\left(\mathrm{c}\mathrm{-}\mathrm{avg}\right)$ , and the individual-average effect, $\mathrm{Effect}\left(\mathrm{i}\mathrm{-}\mathrm{avg}\right)$ , as the following (weighted) averages of the cluster-specific effects:

$$\mathrm{Effect}\left(\mathrm{c}\mathrm{-}\mathrm{avg}\right)=\frac{1}{J}\sum _{j=1}^J\mathrm{Effect}\left(j\right),\mathrm{Effect}\left(\mathrm{i}\mathrm{-}\mathrm{avg}\right)=\frac{1}{J}\sum _{j=1}^J\frac{N_j}{\overline{N}}\mathrm{Effect}\left(j\right),$$

(7)

where $\overline{N}=\sum _{j=1}^J{N}_j/J$ is the mean cluster size. The definition (7) is a sample version of the average effect definition, as the average is across the $J$ clusters in the sample.

Alternatively, we have the superpopulation version of the average effect definition (Miratrix et al., 2021), which will be the focus of our study. This definition considers the characteristics of the sampled clusters, including the cluster size and cluster-specific effect $\{N_j,\mathrm{Effect}\left(j\right)\}$ , are representative of some much broader population of clusters (i.e., the superpopulation with infinite clusters). Following the literature (Wang et al., 2024), we define the superpopulation version of the cluster-average effect, $\mathrm{Effect}\left(\mathrm{C}\mathrm{-}\mathrm{avg}\right)$ , and individual-average effect, $\mathrm{Effect}\left(\mathrm{I}\mathrm{-}\mathrm{avg}\right)$ , using the following expectations of the cluster-specific effect:

$$\mathrm{Effect}\left(\mathrm{C}\mathrm{-}\mathrm{avg}\right)=\mathrm{E}\left[\frac{N_j\mathrm{Effect}\left(j\right)}{N_j}\right]=\mathrm{E}\left[\mathrm{Effect}\left(j\right)\right],\mathrm{Effect}\left(\mathrm{I}\mathrm{-}\mathrm{avg}\right)=\frac{\mathrm{E}[N_j·\mathrm{Effect}\left(j\right)]}{\mathrm{E}[N_j]},$$

(8)

where the expectation “ $\mathrm{E}[]$ ” is over the distribution of $\{N_j,\mathrm{Effect}\left(j\right)\}$ across repeated simple random sampling of a cluster, and then simple random sampling of $N_j$ individuals in the sampled cluster (with cluster size $N_j$ equal to or proportional to the cluster’s population size).

Identification

To estimate the causal mediation effects (defined above), the first step is to express these counterfactual-based estimands as functions only of the observed data distribution (i.e., identification). Extending Vansteelandt and Daniel (2017) to the clustered setting, we make the following assumptions for identification:

• A1. The effect of treatment on outcome is conditionally unconfounded, $Y_{\mathrm{ij}}(a,m_1,m_2)\perp \perp A_{\mathrm{ij}}\mid \{X_{\mathrm{ij}},S_j=j\}$ ;

• A2. The effect of both mediators on outcome is conditionally unconfounded, $Y_{\mathrm{ij}}(a,m_1,m_2)\perp \perp \{{M_1}_{\mathrm{ij}},{M_2}_{\mathrm{ij}}\}\mid \{A_{\mathrm{ij}}=a,X_{\mathrm{ij}},S_j=j\}$

• A3. The effect of treatment on both mediators is conditionally unconfounded, $\{{M_1}_{\mathrm{ij}}\left(a\right),{M_2}_{\mathrm{ij}}\left(a\right)\}\perp \perp A_{\mathrm{ij}}\mid \{X_{\mathrm{ij}},S_j=j\}$ .

Furthermore, we make a within-cluster positivity assumption as follows. Within a given cluster $j$ , there is a nonzero probability for each individual to have any value of the treatment and any value of the mediator vector so that $p(A_{\mathrm{ij}}=a\mid X_{\mathrm{ij}},S_j=j)>0$ and $p({M_1}_{\mathrm{ij}}=m_1,{M_2}_{\mathrm{ij}}=m_2\mid A_{\mathrm{ij}}=a,X_{\mathrm{ij}},S_j=j)>0$ for all $a$ , $m_1$ , and $m_2$ values.

The above identification assumptions allow for unmeasured cluster-level confounder $U_j$ (Figure 2), by conditioning on the cluster membership $S_j$ and correspondingly making a stronger, within-cluster positivity assumption. Under A1 to A3 and the positivity assumption, we identify the cluster- and individual-average causal mediation effects as follows. For each cluster $j$ , we use $\mathrm{\mu }(a,m_1,m_2,X_{\mathrm{ij}},j)$ to denote the conditional mean outcome $\mathrm{E}(Y_{\mathrm{ij}}\mid A_{\mathrm{ij}}=a,{M_1}_{\mathrm{ij}}=m_1,{M_2}_{\mathrm{ij}}=m_2,X_{\mathrm{ij}},S_j=j)$ , use $p(m_1,m_2\mid a,X_{\mathrm{ij}},j)$ to denote the joint probability distribution of both mediators $p({M_1}_{\mathrm{ij}}=m_1,{M_2}_{\mathrm{ij}}=m_2\mid A_{\mathrm{ij}}=a,X_{\mathrm{ij}},S_j=j)$ , and use $p(m_1\mid a,X_{\mathrm{ij}},j)$ and $p(m_2\mid a,X_{\mathrm{ij}},j)$ to respectively denote the probability distribution of each mediator $p({M_1}_{\mathrm{ij}}=m_1\mid A_{\mathrm{ij}}=a,X_{\mathrm{ij}},S_j=j)$ and $p({M_2}_{\mathrm{ij}}=m_2\mid A_{\mathrm{ij}}=a,X_{\mathrm{ij}},S_j=j)$ .

Then, for the direct effect (IDE in Equation 2), the cluster-average and individual-average potential outcomes under $a,a^{*}\in \{0,1\}$ are identified by

$$\begin{array}{c}\mathrm{E}\left[\frac{\sum _{i=1}^{N_j}{Y}_{\mathrm{ij}}(a,{\tilde{M}}_{{\mathrm{joint}}_{\mathrm{ij}}}\left(a^{*}\right))}{N_j}\right]={\mathrm{\theta }}_{\mathrm{joint}}(a,a^{*},\mathrm{C}\mathrm{-}\mathrm{avg})=\mathrm{E}\left[\frac{\sum _{i=1}^{N_j}{\theta }_{\mathrm{joint}}(a,a^{*},j)}{N_j}\right],\hfill \\ \frac{\mathrm{E}\left[\sum _{i=1}^{N_j}{Y}_{\mathrm{ij}}(a,{\tilde{M}}_{{\mathrm{joint}}_{\mathrm{ij}}}\left(a^{*}\right))\right]}{E[N_j]}={\mathrm{\theta }}_{\mathrm{joint}}(a,a^{*},\mathrm{I}\mathrm{-}\mathrm{avg})=\frac{\mathrm{E}\left[\sum _{i=1}^{N_j}{\mathrm{\theta }}_{\mathrm{joint}}(a,a^{*},j)\right]}{\mathrm{E}[N_j]},\hfill \\ \mathrm{where}{\mathrm{\theta }}_{\mathrm{joint}}(a,a^{*},j)=\sum _{x,m_1,m_2}\mathrm{\mu }(a,m_1,m_2,x,j)p(m_1,m_2\mid a^{*},x,j)p(X_{\mathrm{ij}}=x\mid S_j=j).\hfill \end{array}$$

(9)

In the indirect effects (IIEs) via each mediator, the cluster-average and individual-average potential outcomes for $a,a^{*},a^{**}\in \{0,1\}$ are identified by

$$\begin{array}{c}\mathrm{E}\left[\frac{\sum _{i=1}^{N_j}{Y}_{\mathrm{ij}}(a,{\tilde{M}}_{1\mathrm{ij}}\left(a^{*}\right),{\tilde{M}}_{2\mathrm{ij}}\left(a^{**}\right))}{N_j}\right]={\mathrm{\theta }}_{\mathrm{each}}(a,a^{*},a^{**},\mathrm{C}\mathrm{-}\mathrm{avg})=\mathrm{E}\left[\frac{\sum _{i=1}^{N_j}{\mathrm{\theta }}_{\mathrm{each}}(a,a^{*},a^{**},j)}{N_j}\right],\hfill \\ \frac{E\left[\sum _{i=1}^{N_j}{Y}_{\mathrm{ij}}(a,{\tilde{M}}_{1\mathrm{ij}}\left(a^{*}\right),{\tilde{M}}_{2\mathrm{ij}}\left(a^{**}\right))\right]}{\mathrm{E}[N_j]}={\mathrm{\theta }}_{\mathrm{each}}(a,a^{*},a^{**},\mathrm{I}\mathrm{-}\mathrm{avg})=\frac{\mathrm{E}\left[\sum _{i=1}^{N_j}{\mathrm{\theta }}_{\mathrm{each}}(a,a^{*},a^{**},j)\right]}{\mathrm{E}[N_j]},\hfill \\ \mathrm{where}{\mathrm{\theta }}_{\mathrm{each}}(a,a^{*},a^{**},j)=\hfill \\ \sum _{x,m_1,m_2}\mathrm{\mu }(a,m_1,m_2,x,j)p(m_1\mid a^{*},x,j)p(m_2\mid a^{**},x,j)p(X_{\mathrm{ij}}=x\mid S_j=j).\hfill \end{array}$$

(10)

We discuss several components of the identification results (9) and (10). First, conditioning on the cluster membership $S_j=j$ implies conditioning on all cluster-level variables, including the cluster size and measured/unmeasured cluster-level covariates $(N_j,W_j,U_j)$ . Second, the expectation is over the distribution of the cluster size, cluster-level covariates, and individual-level covariates, $(N_j,W_j,U_j,\{X_{\mathrm{ij}}{\}}_{i=1}^{N_j})$ , across the repeated sampling of clusters and then individuals within the clusters. Third, in the effect definition and identification, we have considered the cluster sizes $N_j,j=1,...,J$ to represent the relative population sizes of the clusters $j=1,...,J$ . By considering $N_j$ as equal to the cluster’s population size, this assumes that the $N_j$ individuals in the data at hand constitute the cluster’s population (e.g., all eligible children in a school). Alternatively, by considering $N_j$ as proportional to the cluster’s population size, this assumes that the $N_j$ individuals in the data at hand are a random sample representative of the cluster’s population (e.g., a random sample of all eligible children in a school). Fourth, the identification result (9) and (10) extends the previous causal mediation literature (Vansteelandt & Daniel, 2017) in several ways: (a) it accounts for unmeasured cluster-level confounders, by utilizing the clustered data structure; (b) it distinguishes between the cluster-average and individual-average effect types, both of which can be of substantive interest in the clustered setting; and (c) it explicitly considers varying cluster sizes $N_j$ , and via conditioning on $S_j=j$ , it allows for the influence of $N_j$ on the treatment, mediators, outcome, and their relationships.

With (9) and (10) together, we express the cluster-average and individual-average causal mediation effects (IIEs and IDE) as contrasts of the average potential outcomes following the effect definitions. For example, for the IIEs, the cluster-average version is identified by $\mathrm{II}{\mathrm{E}}_{M_1}\left(\mathrm{C}\mathrm{-}\mathrm{avg}\right)={\mathrm{\theta }}_{\mathrm{each}}(1,1,0,\mathrm{C}\mathrm{-}\mathrm{avg})-{\mathrm{\theta }}_{\mathrm{each}}(1,0,0,\mathrm{C}\mathrm{-}\mathrm{avg})$ and $\mathrm{II}{\mathrm{E}}_{M_2}\left(\mathrm{C}\mathrm{-}\mathrm{avg}\right)={\mathrm{\theta }}_{\mathrm{each}}(1,1,1,\mathrm{C}\mathrm{-}\mathrm{avg})-{\mathrm{\theta }}_{\mathrm{each}}(1,1,0,\mathrm{C}\mathrm{-}\mathrm{avg})$ ; for the IDE, the individual-average version is identified by $\mathrm{IDE}\left(\mathrm{I}\mathrm{-}\mathrm{avg}\right)={\mathrm{\theta }}_{\mathrm{joint}}(1,0,\mathrm{I}\mathrm{-}\mathrm{avg})-{\mathrm{\theta }}_{\mathrm{joint}}(0,0,\mathrm{I}\mathrm{-}\mathrm{avg})$ .

Estimation

In the previous section, we have defined and identified the causal mediation effects. We now examine how to construct an effect estimator and make inferences. One approach could be a regression-based approach, which requires consistently estimating models for the mediators and outcome in the identification result (Loh et al., 2022; Vansteelandt & Daniel, 2017); another approach could be extending the weighting-based approach, which requires consistently estimating models for the mediators and treatment (Hong, 2015). These approaches often estimate the models parametrically, which could induce bias due to overly simplified modeling assumptions; furthermore, if the model estimation is data-adaptive (e.g., via model selection or machine learning methods), the sampling distribution of the resulting effect estimate is generally unknown, creating difficulties for uncertainty quantification and statistical inference. This motivates us to take a different estimation approach that helps overcome this issue.

Specifically, we develop an estimation method based on the multiply robust estimation approach (Benkeser & Ran, 2021). Multiply robust approaches (Kennedy, 2023; Tchetgen Tchetgen & Shpitser, 2012) often yield estimators with desirable properties (e.g., multiple robustness so that the estimator remains consistent even when one of the nuisance models is misspecified); furthermore, they allow the model estimation to be data-adaptive (e.g., via machine learning algorithms), while providing uncertainty quantification via asymptotic variance and Wald-type confidence intervals. For the multiple-mediator analysis, the multiply robust method with unclustered data has been developed (Benkeser & Ran, 2021; Rubinstein, Haviland, & Breslau, 2023). Nonetheless, it has not been extended to the clustered setting where unmeasured cluster-level confounders may exist (Figure 2). Therefore, in what follows, we first describe the multiply robust estimator. Then, we develop methods for (a) estimating the nuisance models allowing for unmeasured cluster-level confounders and (b) making inferences (variance estimation and confidence intervals) considering within-cluster correlations.

The Multiply Robust Estimator

The multiply robust estimator combines estimates of nuisance models for the outcome $Y_{\mathrm{ij}}$ , mediators $({M_1}_{\mathrm{ij}},{M_2}_{\mathrm{ij}})$ , and treatment $A_{\mathrm{ij}}$ , in a way that allows for errors in estimating any one model. To ease describing the estimator, we introduce some additional notation for the nuisance models. We denote the “fully marginalized mean outcomes” in the identification results as

$$\begin{array}{c}{\mathrm{\mu }}_{\mathrm{each}}(a,a^{*},a^{**},X_{\mathrm{ij}},j)=\sum _{m_1,m_2}\mathrm{\mu }(a,m_1,m_2,X_{\mathrm{ij}},j)p(m_1\mid a^{*},X_{\mathrm{ij}},j)p(m_2\mid a^{**},X_{\mathrm{ij}},j)\hfill \\ {\mathrm{\mu }}_{\mathrm{joint}}(a,a^{*},X_{\mathrm{ij}},j)=\sum _{m_1,m_2}\mathrm{\mu }(a,m_1,m_2,X_{\mathrm{ij}},j)p(m_1,m_2\mid a^{*},X_{\mathrm{ij}},j).\hfill \end{array}$$

These mean outcomes are “fully marginalized,” in the sense that they marginalize over both mediators: ${\mathrm{\mu }}_{\mathrm{each}}$ marginalizes the distributions of each mediator, and it will enter the estimator for the indirect effects via each mediator (IIEs); ${\mathrm{\mu }}_{\mathrm{joint}}$ marginalizes the joint distribution of both mediators, and it will be part of the estimator for the direct effect (IDE). We also have the following “partially marginalized mean outcomes”:

$$\begin{array}{c}{\mathrm{\mu }}_{m_1^{*}}(a,{M_2}_{\mathrm{ij}},X_{\mathrm{ij}},j)=\sum _{m_1}\mathrm{\mu }(a,m_1,{M_2}_{\mathrm{ij}},X_{\mathrm{ij}},j)p(m_1\mid a^{*},X_{\mathrm{ij}},j),\hfill \\ {\mathrm{\mu }}_{m_2^{**}}(a,{M_1}_{\mathrm{ij}},X_{\mathrm{ij}},j)=\sum _{m_2}\mathrm{\mu }(a,{M_1}_{\mathrm{ij}},m_2,X_{\mathrm{ij}},j)p(m_2\mid a^{**},X_{\mathrm{ij}},j),\hfill \end{array}$$

which marginalize only one of the mediators, while still being a function of the other mediator.

Furthermore, we introduce the following mediator weights

$$w_{\mathrm{each}}({M_1}_{\mathrm{ij}},{M_2}_{\mathrm{ij}},A_{\mathrm{ij}},X_{\mathrm{ij}},j)=\frac{p({M_1}_{\mathrm{ij}}\mid a^{*},X_{\mathrm{ij}},j)p({M_2}_{\mathrm{ij}}\mid a^{**},X_{\mathrm{ij}},j)}{p({M_1}_{\mathrm{ij}},{M_2}_{\mathrm{ij}}\mid A_{\mathrm{ij}},X_{\mathrm{ij}},j)}$$

(11)

$$w_{\mathrm{joint}}({M_1}_{\mathrm{ij}},{M_2}_{\mathrm{ij}},A_{\mathrm{ij}},X_{\mathrm{ij}},j)=\frac{p({M_1}_{\mathrm{ij}},{M_2}_{\mathrm{ij}}\mid a^{*},X_{\mathrm{ij}},j)}{p({M_1}_{\mathrm{ij}},{M_2}_{\mathrm{ij}}\mid A_{\mathrm{ij}},X_{\mathrm{ij}},j)}.$$

(12)

In these mediator weights, the numerators are the mediator distributions in the identification results (10) and (9), and the denominator is the distribution of the mediators given the observed treatment. Intuitively, these weights adjust the observational distribution of the mediators to be the desired one (i.e., that in the identification result). Finally, we use $p(a\mid X_{\mathrm{ij}},j)$ to denote $p(A_{\mathrm{ij}}=a\mid X_{\mathrm{ij}},S_j=j)$ , the conditional probability of treatment, or the “propensity score” (Rosenbaum & Rubin, 1983).

The multiply robust estimator (Benkeser & Ran, 2021) combines estimates of the above nuisance models with the observed data as follows. For the identification result ${\mathrm{\theta }}_{\mathrm{each}}(a,a^{*},a^{**},j)$ in (10), it has the following cluster-specific estimating function, which is a function of the observed data $O_j$ in cluster $j$ :

$$\begin{array}{c}\multicolumn{1}{c}{{\hat{\mathrm{\phi }}}_{\mathrm{each}}(a,a^{*},a^{**},O_j)=\frac{1}{N_j}\sum _{i=1}^{N_j}{\hat{\mathrm{\phi }}}_{\mathrm{each}}(a,a^{*},a^{**},O_{\mathrm{ij}}),\mathrm{where}{\hat{\mathrm{\phi }}}_{\mathrm{each}}(a,a^{*},a^{**},O_{\mathrm{ij}})}\\ =\frac{1(A_{\mathrm{ij}}=a)}{\hat{p}(a\mid X_{\mathrm{ij}},j)}{\hat{w}}_{\mathrm{each}}({M_1}_{\mathrm{ij}},{M_2}_{\mathrm{ij}},a,X_{\mathrm{ij}},j)\{Y_{\mathrm{ij}}-\hat{\mathrm{\mu }}(a,{M_1}_{\mathrm{ij}},{M_2}_{\mathrm{ij}},X_{\mathrm{ij}},j)\}\hfill \end{array}$$

(13)

$$+\frac{1(A_{\mathrm{ij}}=a^{*})}{\hat{p}(a^{*}\mid X_{\mathrm{ij}},j)}\{{\hat{\mathrm{\mu }}}_{m_2^{**}}(a,{M_1}_{\mathrm{ij}},X_{\mathrm{ij}},j)-{\hat{\mathrm{\mu }}}_{\mathrm{each}}(a,a^{*},a^{**},X_{\mathrm{ij}},j)\}$$

(14)

$$\begin{array}{c}+\frac{1(A_{\mathrm{ij}}=a^{**})}{\hat{p}(a^{**}\mid X_{\mathrm{ij}},j)}\{{\hat{\mathrm{\mu }}}_{m_1^{*}}(a,{M_2}_{\mathrm{ij}},X_{\mathrm{ij}},j)-{\hat{\mathrm{\mu }}}_{\mathrm{each}}(a,a^{*},a^{**},X_{\mathrm{ij}},j)\}\hfill \\ +{\hat{\mathrm{\mu }}}_{\mathrm{each}}(a,a^{*},a^{**},X_{\mathrm{ij}},j).\hfill \end{array}$$

(15)

For the identification result ${\mathrm{\theta }}_{\mathrm{joint}}(a,a^{*},j)$ in (9), the cluster-specific estimating function is

$$\begin{array}{c}\multicolumn{1}{c}{{\hat{\mathrm{\phi }}}_{\mathrm{joint}}(a,a^{*},O_j)=\frac{1}{N_j}\sum _{i=1}^{N_j}{\hat{\mathrm{\phi }}}_{\mathrm{joint}}(a,a^{*},O_{\mathrm{ij}}),\mathrm{where}{\hat{\mathrm{\phi }}}_{\mathrm{joint}}(a,a^{*},O_{\mathrm{ij}})}\\ =\frac{1(A_{\mathrm{ij}}=a)}{\hat{p}(a\mid X_{\mathrm{ij}},j)}{\hat{w}}_{\mathrm{joint}}({M_1}_{\mathrm{ij}},{M_2}_{\mathrm{ij}},a,X_{\mathrm{ij}},j)\{Y_{\mathrm{ij}}-\hat{\mathrm{\mu }}(a,{M_1}_{\mathrm{ij}},{M_2}_{\mathrm{ij}},X_{\mathrm{ij}},j)\}\hfill \end{array}$$

(16)

$$\begin{array}{c}+\frac{1(A_{\mathrm{ij}}=a^{*})}{\hat{p}(a^{*}\mid X_{\mathrm{ij}},j)}\{\hat{\mathrm{\mu }}(a,{M_1}_{\mathrm{ij}},{M_2}_{\mathrm{ij}},X_{\mathrm{ij}},j)-{\hat{\mathrm{\mu }}}_{\mathrm{joint}}(a,a^{*},X_{\mathrm{ij}},j)\}\hfill \\ +{\hat{\mathrm{\mu }}}_{\mathrm{joint}}(a,a^{*},X_{\mathrm{ij}},j).\hfill \end{array}$$

(17)

Further, we perform weight normalization for the above cluster-specific estimating functions. Weight normalization helps alleviate the potential adverse influence of unstable weights, which has been suggested in previous literature on multiply robust estimation (Farbmacher et al., 2022). Specifically, within each cluster, we divide each weighting term by the cluster-mean weight. That is, for ${\hat{\mathrm{\phi }}}_{\mathrm{each}}(a,a^{*},a^{**},O_j)$ , we divide the first line (13) by the cluster-mean weight $\frac{1}{N_j}\sum _{i=1}^{N_j}1(A_{\mathrm{ij}}=a)/\hat{p}(a\mid X_{\mathrm{ij}},j)·{\hat{w}}_{\mathrm{each}}({M_1}_{\mathrm{ij}},{M_2}_{\mathrm{ij}},X_{\mathrm{ij}},j)$ , divide the second line (14) by $\frac{1}{N_j}\sum _{i=1}^{N_j}1(A_{\mathrm{ij}}=a^{*})/\hat{p}(a^{*}\mid X_{\mathrm{ij}},j)$ , and divide the third line (15) by $\frac{1}{N_j}\sum _{i=1}^{N_j}1(A_{\mathrm{ij}}=a^{**})/\hat{p}(a^{**}\mid X_{\mathrm{ij}},j)$ ; for ${\hat{\mathrm{\phi }}}_{\mathrm{joint}}(a,a^{*},O_j)$ , we divide the first line (16) by $\frac{1}{N_j}\sum _{i=1}^{N_j}1(A_{\mathrm{ij}}=a)/\hat{p}(a\mid X_{\mathrm{ij}},j)·{\hat{w}}_{\mathrm{joint}}({M_1}_{\mathrm{ij}},{M_2}_{\mathrm{ij}},X_{\mathrm{ij}},j)$ , and divide the second line (17) by $\frac{1}{N_j}\sum _{i=1}^{N_j}1(A_{\mathrm{ij}}=a^{*})/\hat{p}(a^{*}\mid X_{\mathrm{ij}},j)$ .

After the weight normalization for Equations 13 to 17, we contrast the cluster-specific estimating functions to obtain the causal mediation effects:

$$\begin{array}{c}\widehat{\mathrm{IDE}}\left(O_j\right)={\hat{\mathrm{\phi }}}_{\mathrm{joint}}(1,0,O_j)-{\hat{\mathrm{\phi }}}_{\mathrm{joint}}(0,0,O_j),\hfill \\ {\widehat{\mathrm{IIE}}}_{M_1}\left(O_j\right)={\hat{\mathrm{\phi }}}_{\mathrm{each}}(1,1,0,O_j)-{\hat{\mathrm{\phi }}}_{\mathrm{each}}(1,0,0,O_j),{\widehat{\mathrm{IIE}}}_{M_2}\left(O_j\right)={\hat{\mathrm{\phi }}}_{\mathrm{each}}(1,1,1,O_j)-{\hat{\mathrm{\phi }}}_{\mathrm{each}}(1,1,0,O_j)\hfill \\ {\widehat{\mathrm{IIE}}}_{\mathrm{mu}}\left(O_j\right)={\hat{\mathrm{\phi }}}_{\mathrm{joint}}(1,1,O_j)-{\hat{\mathrm{\phi }}}_{\mathrm{each}}(1,1,1,O_j)-[{\hat{\mathrm{\phi }}}_{\mathrm{joint}}(1,0,O_j)-{\hat{\mathrm{\phi }}}_{\mathrm{each}}(1,0,0,O_j)].\hfill \end{array}$$

(18)

Then, let $\widehat{\mathrm{Effect}}\left(O_j\right)\in \{\widehat{\mathrm{IDE}}\left(O_j\right),{\widehat{\mathrm{IIE}}}_{M_1}\left(O_j\right),{\widehat{\mathrm{IIE}}}_{M_2}\left(O_j\right),{\widehat{\mathrm{IIE}}}_{\mathrm{mu}}\left(O_j\right)\}$ be each of these cluster-specific estimating functions. We obtain multiply robust estimators for the cluster-average and individual-average effects as the following (weighted) averages:

$$\widehat{\mathrm{Effect}}\left(\mathrm{C}\mathrm{-}\mathrm{avg}\right)=\frac{1}{J}\sum _{j=1}^J\widehat{\mathrm{Effect}}\left(O_j\right),\widehat{\mathrm{Effect}}\left(\mathrm{I}\mathrm{-}\mathrm{avg}\right)=\frac{1}{J}\sum _{j=1}^J\frac{N_j}{\overline{N}}·\widehat{\mathrm{Effect}}\left(O_j\right),$$

(19)

where $\overline{N}=\sum _{j=1}^J{N}_j/J$ is the average cluster size (as denoted earlier). Analogous to the effect definitions (8), the effect estimator for the cluster-average effect weights clusters equally, and the effect estimator for the individual-average effect weights individuals equally (i.e., by weighting clusters proportionally to cluster sizes).

The multiply robust estimator (19) provides multiple robustness for consistent estimation of the causal mediation effects. Specifically, it would remain consistent under consistent estimation for at least one of the following sets of nuisance models (Benkeser & Ran, 2021): (a) the conditional mean outcome $\hat{\mathrm{\mu }}(a,m_1,m_2,X_{\mathrm{ij}},j)$ and the joint mediator distribution $\hat{p}(m_1,m_2\mid a,X_{\mathrm{ij}},j)$ are consistently estimated; or (b) the propensity score $\hat{p}(a\mid X_{\mathrm{ij}},j)$ , the conditional mean outcome, and either mediator distribution $\hat{p}(m_1\mid a,X_{\mathrm{ij}},j)$ or $\hat{p}(m_2\mid a,X_{\mathrm{ij}},j)$ , are consistently estimated; or (c) the propensity score and the joint mediator distribution are consistently estimated.

Estimating the Nuisance Models Under Unmeasured Cluster-Level Confounders

To control for the unmeasured cluster-level confounder $U_j$ , all the nuisance models are conditioned on the cluster membership $S_j=j$ . Conditioning on cluster membership $S_j=j$ implies conditioning on (a) the means of individual-level variables within the cluster and (b) a full set of dummy indicators of the cluster membership. Therefore, we implement the conditioning on the cluster membership $S_j=j$ in the nuisance model estimation via the following two approaches.

First, we control for cluster means and perform cluster-mean centering for individual-level variables in all the nuisance models. Specifically, cluster-mean centered variables are orthogonal to any cluster-level variable, including unmeasured cluster-level confounders; thus, cluster-mean centering helps remove the influence of unmeasured cluster-level confounders (Raudenbush & Bryk, 2002; Suk & Kang, 2023; Zetterqvist et al., 2016). In addition, while conditioning on cluster means of individual-level variables does not necessarily imply conditioning on the cluster membership ( $S_j$ ), variation in cluster means can result from the influence of unmeasured cluster-level variables; therefore, conditioning on the cluster means can also help control for unmeasured cluster-level confounders.

Second, we also control for dummy indicators of an individual’s cluster membership in all the nuisance models. Conditioning on the full set of cluster dummies necessarily implies conditioning on the cluster membership ( $S_j$ ), and has been a common approach to handling unmeasured cluster-level confounders in treatment effect estimation (Chang & Stuart, 2022; Li et al., 2013). A potential issue is that with many clusters, incorporating a large number of cluster dummies into the nuisance models may result in unstable model estimates, which in turn may reduce the precision of effect estimation (Li et al., 2013). Therefore, in the next section, we will explore omitting the cluster dummies (as an alternative to including the cluster dummies) in the nuisance model estimation through simulations.

Specifically, we implement the above two approaches to controlling for the unmeasured cluster-level confounder as follows. We first calculate the cluster means of individual-level covariates, treatment, mediators, and outcome: ${\overline{X}}_j=\sum _{i=1}^{N_j}{X}_{\mathrm{ij}}/N_j$ , ${\overline{\mathrm{A}}}_j=\sum _{i=1}^{N_j}{A}_{\mathrm{ij}}/N_j$ , ${\overline{M_1}}_j=\sum _{i=1}^{N_j}{M_1}_{\mathrm{ij}}/N_j$ , ${\overline{M_2}}_j=\sum _{i=1}^{N_j}{M_2}_{\mathrm{ij}}/N_j$ , and ${\overline{Y}}_j=\sum _{i=1}^{N_j}{Y}_{\mathrm{ij}}/N_j$ . We also create a full set of cluster dummies, $\{S_j^{\mathrm{dumm}}{\}}_{j=1}^J$ , where $S_j^{\mathrm{dumm}}=1$ for individuals with cluster membership $j$ (i.e., $S_j=j$ ) and $S_j^{\mathrm{dumm}}=0$ otherwise. Then, we incorporate these cluster means and cluster dummies into estimating the nuisance models for the treatment, mediators, and outcome.

For the treatment, we estimate its conditional probability (propensity score) using a regression model for $A_{\mathrm{ij}}$ , regressing it on the cluster-mean centered covariates $X_{\mathrm{ij}}-{\overline{X}}_j$ , its cluster mean ${\overline{A}}_j$ , measured cluster-level covariates $W_j$ , the cluster size $N_j$ , and the set of cluster dummies $\{S_j^{\mathrm{dumm}}{\}}_{j=1}^J$ . The cluster-mean treatment serves as a cluster-level summary statistic to account for the influence of unmeasured cluster-level confounder (as described above), which was also used in propensity score estimation with unmeasured cluster-level confounders (Zetterqvist et al., 2016). Using the fitted model, we obtain the estimate (i.e., predicted value) of the conditional probability of treatment $a\in \{0,1\}$ given an individual’s covariates and cluster membership: $\hat{p}(a\mid X_{\mathrm{ij}},j)=\hat{p}(A_{\mathrm{ij}}=a\mid X_{\mathrm{ij}}-{\overline{X}}_j,{\overline{A}}_j,W_j,N_j,\{S_j^{\mathrm{dumm}}{\}}_{j=1}^J)$ .

For the mediators, we estimate the joint probability distribution by factorizing it as the distribution of one mediator and the conditional distribution of the other mediator, $p({M_1}_{\mathrm{ij}},{M_2}_{\mathrm{ij}}\mid a,X_{\mathrm{ij}},j)=p({M_1}_{\mathrm{ij}}\mid a,X_{\mathrm{ij}},j)p({M_2}_{\mathrm{ij}}\mid {M_1}_{\mathrm{ij}},a,X_{\mathrm{ij}},j)$ . Once the joint distribution is estimated, we obtain the distribution of the other mediator as $p({M_2}_{\mathrm{ij}}\mid a,X_{\mathrm{ij}},j)=\sum _{m_1}p(m_1,{M_2}_{\mathrm{ij}}\mid a,X_{\mathrm{ij}},j)$ .

Specifically, to estimate the factorized joint mediator distribution, we fit a regression model for ${M_1}_{\mathrm{ij}}$ , regressing it on $A_{\mathrm{ij}}-{\overline{A}}_j$ , $X_{\mathrm{ij}}-{\overline{X}}_j$ , ${\overline{M_1}}_j$ , $W_j$ , $N_j$ , and $\{S_j^{\mathrm{dumm}}{\}}_{j=1}^J$ (i.e., the cluster-mean centered treatment and covariates, its cluster mean, measured cluster-level covariates, cluster size, cluster dummies). Using the fitted mediator model, we obtain the estimate for the probability of mediator ${M_1}_{\mathrm{ij}}=m_1$ given the treatment value $a\in \{0,1\}$ , conditional on an individual’s covariates and cluster membership: $\hat{p}({M_1}_{\mathrm{ij}}=m_1\mid a,X_{\mathrm{ij}},j)=\hat{p}({M_1}_{\mathrm{ij}}=m_1\mid a-{\overline{A}}_j,X_{\mathrm{ij}}-{\overline{X}}_j,{\overline{M_1}}_j,W_j,N_j,\{S_j^{\mathrm{dumm}}{\}}_{j=1}^J)$ . Then, to estimate the conditional mediator distribution, we fit a regression model for the other mediator ${M_2}_{\mathrm{ij}}$ , regressing it on ${M_1}_{\mathrm{ij}}-{\overline{M_1}}_j$ , $A_{\mathrm{ij}}-{\overline{A}}_j$ , $A_{\mathrm{ij}}{M_1}_{\mathrm{ij}}-{\overline{A{M}_1}}_j$ , $X_{\mathrm{ij}}-{\overline{X}}_j$ , ${\overline{M_2}}_j$ , $W_j$ , $N_j$ , and $\{S_j^{\mathrm{dumm}}{\}}_{j=1}^J$ . In this conditional mediator model, the interaction term is included for the treatment and the other mediator $A_{\mathrm{ij}}{M_1}_{\mathrm{ij}}$ , as they are variables of interest; the interaction term is also centered at its cluster mean ${\overline{A{M}_1}}_j=\sum _{i=1}^{N_j}{A}_{\mathrm{ij}}{M_1}_{\mathrm{ij}}/N_j$ . From this fitted model, we obtain the estimate for the probability of mediator ${M_2}_{\mathrm{ij}}=m_2$ for given values of the treatment and the other mediator $a,m_1\in \{0,1\}$ , conditional on an individual’s covariates and cluster membership: $\hat{p}({M_2}_{\mathrm{ij}}=m_2\mid a,m_1,X_{\mathrm{ij}},j)=\hat{p}({M_2}_{\mathrm{ij}}=m_2\mid a-{\overline{A}}_j,m_1-{\overline{M_1}}_j,\mathrm{am}-{\overline{A{M}_1}}_j,X_{\mathrm{ij}}-{\overline{X}}_j,{\overline{M_2}}_j,W_j,N_j,\{S_j^{\mathrm{dumm}}{\}}_{j=1}^J)$ .

For the outcome, we estimate the conditional mean using a regression model for $Y_{\mathrm{ij}}$ , regressing it on the cluster-mean centered covariates $X_{\mathrm{ij}}-{\overline{X}}_j$ , the cluster-mean centered treatment and mediators and their interactions $A_{\mathrm{ij}}-{\overline{A}}_j$ , ${M_1}_{\mathrm{ij}}-{\overline{M_1}}_j$ , ${M_2}_{\mathrm{ij}}-{\overline{M_2}}_j$ , $A_{\mathrm{ij}}{M_1}_{\mathrm{ij}}-{\overline{A{M}_1}}_j$ , $A_{\mathrm{ij}}{M_2}_{\mathrm{ij}}-{\overline{A{M}_2}}_j$ , ${M_1}_{\mathrm{ij}}{M_2}_{\mathrm{ij}}-{\overline{M_1{M}_2}}_j$ , $A_{\mathrm{ij}}{M_1}_{\mathrm{ij}}{M_2}_{\mathrm{ij}}-{\overline{A{M}_1{M}_2}}_j$ , cluster-mean outcome ${\overline{Y}}_j$ , measured cluster-level covariates $W_j$ , cluster size $N_j$ , and cluster dummies $\{S_j^{\mathrm{dumm}}{\}}_{j=1}^J$ . In the outcome model, all interactions among the treatment and mediators are considered; these interaction terms are cluster-mean centered, where the cluster mean is the mean of the interaction term within the cluster, for example, ${\overline{A{M}_1{M}_2}}_j=\sum _{i=1}^{N_j}{A}_{\mathrm{ij}}{M_1}_{\mathrm{ij}}{M_2}_{\mathrm{ij}}/N_j$ . Using the fitted outcome model, we obtain the estimate of the conditional mean outcome for given values of the treatment and mediators $a,m_1,m_2\in \{0,1\}$ , conditional on an individual’s covariates and cluster membership: $\hat{\mathrm{\mu }}(a,m_1,m_2,X_{\mathrm{ij}},j)=\hat{\mathrm{E}}[Y_{\mathrm{ij}}\mid a-{\overline{A}}_j,m_1-{\overline{M_1}}_j,m_2-{\overline{M_2}}_j,a{m}_1-{\overline{A{M}_1}}_j,a{m}_2-{\overline{A{M}_2}}_j,m_1{m}_2-{\overline{M_1{M}_2}}_j,a{m}_1{m}_2-{\overline{A{M}_1{M}_2}}_j,X_{\mathrm{ij}}-{\overline{X}}_j,{\overline{Y}}_j,W_j,N_j,\{S_j^{\mathrm{dumm}}{\}}_{j=1}^J]$ .

Using Machine Learning to Estimate the Nuisance Models

To fit the above regressions and obtain nuisance model estimates, parametric modeling (e.g., generalized linear modeling) may be applicable when we have prior knowledge of correct model specifications; however, this knowledge is generally unavailable in practice. Therefore, we consider the use of machine learning methods to obtain the nuisance model estimates.

Based on the literature (Kennedy, 2023; Rubinstein, Branson, & Kennedy, 2023), although machine learning methods often yield nuisance model estimates with small errors (e.g., regularization errors), multiply robust estimators (including those presented above) are robust to small errors in the nuisance model estimates. Specifically, unlike regression-based estimation (i.e., directly plug-in nuisance model estimates into the identification results), with which the asymptotic bias contains first-order terms of errors of the nuisance model estimates, the multiply robust estimator has no first-order bias asymptotically. The asymptotic bias of the multiply robust estimator contains only products of errors of two estimated nuisance models, which is of second order (see Rubinstein, Branson, & Kennedy, 2023 for detailed expressions). Thus, the bias of the multiply robust estimator will be asymptotically negligible, when the product of errors in estimating any two nuisance models goes to zero at a rate faster than $1/\sqrt{J}$ , which is satisfied when the nuisance models are estimated by correct parametric models or by machine learning methods with sufficiently fast convergence rate (achievable by commonly available machine learning methods such as lasso, random forest, and neural network; Díaz, 2020). Therefore, even when machine learning methods are used, the multiply robust estimator (Equation 19) minus the true effect would asymptotically behave as a sample average of $J$ independent, mean-zero random variables, which converges to a normal distribution at $1/\sqrt{J}$ rate; this asymptotic normality permits the use of the variance estimator and Wald-type confidence intervals (described in the next subsection) for uncertainty quantification and inference.

With machine learning methods to estimate nuisance models, cross-fitting is necessary to obtain estimates (i.e., predicted values) of the nuisance models for calculating the multiply robust estimator (Chernozhukov et al., 2018; Zheng & van der Laan, 2011). Cross-fitting helps avoid bias due to overfitting issues that may occur with data-adaptive machine learning regressions (Chernozhukov et al., 2018). To perform cross-fitting while considering the presence of the unmeasured cluster-level confounder, we randomly split the dataset within clusters to obtain $V$ folds (e.g., five folds) of approximately equal number of individuals (so that each fold contains individuals from all clusters). For each fold, we let it be the validation dataset and use the remaining $V-1$ folds as the training dataset to fit the nuisance models. Then, using the fitted models (fitted with the training dataset), we obtain estimates (i.e., predicted values) for nuisance models of individuals in the validation dataset. After obtaining the nuisance model estimates for all folds, we use these estimates to calculate the multiply robust estimators (19).

Inference

As mentioned in the previous subsection, we obtain an asymptotic variance estimator and the Wald-type confidence interval for inference with the multiply robust estimator (19). We estimate the asymptotic variance as the empirical variance of the cluster-specific estimating function $[s_j\widehat{\mathrm{Effect}}\left(O_j\right)]$ divided by the number of clusters $J$ , where $s_j=1$ for the cluster-average effect and $s_j=N_j/\overline{N}$ for the individual-average effect; that is, the variance estimator is

$$\hat{\mathrm{V}}\mathrm{ar}[\widehat{\mathrm{Effect}}\left(\mathrm{C}\mathrm{-}\mathrm{avg}\right)]=\frac{1}{J}\mathrm{var}[\widehat{\mathrm{Effect}}\left(O_j\right)],\hat{\mathrm{V}}\mathrm{ar}[\widehat{\mathrm{Effect}}\left(\mathrm{I}\mathrm{-}\mathrm{avg}\right)]=\frac{1}{J}\mathrm{var}\left[\frac{N_j}{\overline{N}}\widehat{\mathrm{Effect}}\left(O_j\right)\right],$$

(20)

where $\mathrm{var}[·]$ is taking the empirical variance. This variance estimator considers clusters as independent units while accounting for correlations of observations within clusters.

Then, using the variance estimator, we construct a Wald-type confidence interval using the $t$ -distribution. The use of $t$ -distribution (instead of the normal distribution) is to improve the confidence interval coverage when there are few clusters (Wang et al., 2024). Specifically, for the effects contrasting two potential outcomes, namely the indirect effects via each mediator ( $\mathrm{II}{\mathrm{E}}_{M_1}$ , $\mathrm{II}{\mathrm{E}}_{M_2}$ ) and the direct effect ( $\mathrm{IDE}$ ), we use the $t$ -distribution with $J-2$ degrees of freedom, and correspondingly, we multiply the variance estimator by $\frac{J}{J-2}$ . That is, the 95% interval for $\mathrm{II}{\mathrm{E}}_{M_1}$ is: ${\widehat{\mathrm{IIE}}}_{M_1}\pm 0_{.975,\mathrm{df}=J-2}\sqrt{\frac{J}{J-2}\hat{\mathrm{V}}\mathrm{ar}[{\widehat{\mathrm{IIE}}}_{M_1}]}$ ; the 95% intervals for $\mathrm{II}{\mathrm{E}}_{M_2}$ , and $\mathrm{IDE}$ are obtained in the same way. For the indirect effect due to mutual dependence $\mathrm{II}{\mathrm{E}}_{\mathrm{mu}}$ , which contrasts four potential outcomes, we use the degrees of freedom $J-4$ , and correspondingly multiply its variance estimator by $\frac{J}{J-4}$ ; thus, the 95% interval for $\mathrm{II}{\mathrm{E}}_{\mathrm{mu}}$ is ${\widehat{\mathrm{IIE}}}_{\mathrm{mu}}\pm 0_{.975,\mathrm{df}=J-4}\sqrt{\frac{J}{J-4}\hat{\mathrm{V}}\mathrm{ar}[{\widehat{\mathrm{IIE}}}_{\mathrm{mu}}]}$ .

We provide R package “MediationClustered” for implementing our method (available at https://github.com/xliu12/MediationClustered).

Simulation Study

We conduct a simulation study to evaluate the proposed method. The R code for the simulations is available at https://osf.io/skcba/?view_only=None. We examined two types of sample size conditions: (1) a smaller number of larger clusters and (2) a larger number of smaller clusters. For (1), the number of clusters was $J=$ J = 10, 20, or 40, with the cluster size $N_j$ uniformly distributed between 50 and 100, $N_j~\mathrm{Uniform}[N_{\min }=50,N_{\max }=100]$ . For (2), the number of clusters was $J=$ J = 40, 70, or 100, with $N_j$ uniformly distributed between 5 and 20, $N_j~\mathrm{Uniform}[N_{\min }=5,N_{\max }=20]$ . Larger numbers of clusters often occur when large-scale surveys are data sources for analyses (e.g., 70–100 schools); smaller numbers of clusters can occur in smaller-scale multisite studies (e.g., 10 sites).

Data Generation

For each cluster $j$ , an unmeasured cluster-level confounder $U_j$ was generated from the standard normal distribution. For each individual in the cluster, three individual-level covariates were generated $X_{\mathrm{ij}}=\{X_{\left(p\right)\mathrm{ij}},p=1,2,3\}$ . Each was $X_{\left(p\right)\mathrm{ij}}=X_{\left(p\right)j}^b+X_{\left(p\right)\mathrm{ij}}^w$ , a summation of its between-cluster component $X_{\left(p\right)j}^b$ and within-cluster component $X_{\left(p\right)\mathrm{ij}}^w$ ; $X_{\left(p\right)j}^b$ and $X_{\left(p\right)\mathrm{ij}}^w$ were each generated by normal distributions with means of zero and variances of 0.2 and 0.8, respectively. The unmeasured $U_j$ and measured $X_{\mathrm{ij}}$ , as well as cluster size $N_j$ , affected the treatment, mediators, and outcome. Specifically, let $N_j^{\mathrm{scaled}}=\frac{N_j-N_{\min }}{N_{\max }-N_{\min }}$ be the cluster size scaled to range [0, 1]. The data-generation models for $A_{\mathrm{ij}}$ , ${M_1}_{\mathrm{ij}}$ , ${M_2}_{\mathrm{ij}}$ , and $Y_{\mathrm{ij}}$ were as follows (with $1\{·\}$ being the indicator function):

$$\begin{array}{c}\multicolumn{1}{c}{A_{\mathrm{ij}}=1\{X_{\mathrm{ij}}^{\top }{\mathrm{\lambda }}_1+{\mathrm{\lambda }}_2{U}_j+s_j^a+{ϵ}_{\mathrm{ij}}^a>0\}}\\ \multicolumn{1}{c}{{M_1}_{\mathrm{ij}}=1\{({\mathrm{\alpha }}_1+{\mathrm{\alpha }}_1{U}_j+{\mathrm{\alpha }}_1{N_j}^{\mathrm{scaled}})A_{\mathrm{ij}}+X_{\mathrm{ij}}^{\top }{\mathrm{\alpha }}_2+{\mathrm{\alpha }}_3{U}_j+s_j^{m_1}+{ϵ}_{\mathrm{ij}}^{m_1}>0\}}\\ \multicolumn{1}{c}{{M_2}_{\mathrm{ij}}=1\{({\mathrm{\gamma }}_1+{\mathrm{\gamma }}_1{U}_j+{\mathrm{\gamma }}_1{N_j}^{\mathrm{scaled}})A_{\mathrm{ij}}+{\mathrm{\gamma }}_2{M_1}_{\mathrm{ij}}+{\mathrm{\gamma }}_3{A}_{\mathrm{ij}}{M_1}_{\mathrm{ij}}+X_{\mathrm{ij}}^{\top }{\mathrm{\gamma }}_4+{\mathrm{\gamma }}_5{U}_j+s_j^{m_2}+{ϵ}_{\mathrm{ij}}^{m_2}>0\}}\\ \multicolumn{1}{c}{Y_{\mathrm{ij}}={\mathrm{\beta }}_0+({\mathrm{\beta }}_1+{\mathrm{\beta }}_1{U}_j+{\mathrm{\beta }}_1{N_j}^{\mathrm{scaled}})A_{\mathrm{ij}}+({\mathrm{\beta }}_2+{\mathrm{\beta }}_2{U}_j){M_1}_{\mathrm{ij}}+({\mathrm{\beta }}_3+{\mathrm{\beta }}_3{U}_j){M_2}_{\mathrm{ij}}+}\\ {\mathrm{\beta }}_4{M_1}_{\mathrm{ij}}{M_2}_{\mathrm{ij}}+{\mathrm{\beta }}_5{A}_{\mathrm{ij}}{M_1}_{\mathrm{ij}}{M_2}_{\mathrm{ij}}+X_{\mathrm{ij}}^{\top }{\mathrm{\beta }}_6+{\mathrm{\beta }}_7{U}_j+s_j^y+{ϵ}_{\mathrm{ij}}^y,\hfill \end{array}$$

(21)

where the treatment effects can be heterogeneous across clusters due to cluster size $N_j$ . Coefficients in the data-generation models are described in the Supplemental Materials (Appendix B in the online version of the journal).

Simulation Conditions

In addition to the sample size conditions ( $J$ and $N_j$ ), we varied the following data-generation conditions: (1) Null versus non-null effects and (2) linear versus nonlinear terms of covariates. Specifically, for (1), we simulated the null condition to evaluate the Type I error rates. In the null condition, all terms involving the treatment and mediators were set to zero in the data-generation model for $Y_{\mathrm{ij}}$ . Thus, all the individual-average and cluster-average causal mediation effects were zero in the null condition. For (2), we manipulated the functional forms of the covariates to evaluate whether machine learning estimation for nuisance models can yield reduced bias compared to parametric estimation when the data generation has nonlinearities. In the “linear” scenario, we used the data-generation models as shown above (Equation 21), in which the covariates $X_{\mathrm{ij}}$ only had linear terms. In the “nonlinear” scenario, we replaced each covariate $X_{\left(p\right)\mathrm{ij}}$ ( $p=1,2,3$ ) with its quadratic term $\frac{X_{\left(p\right)\mathrm{ij}}^2-1}{2}$ in all data-generation models (Equation 21) for $A_{\mathrm{ij}}$ , ${M_1}_{\mathrm{ij}}$ , ${M_2}_{\mathrm{ij}}$ , and $Y_{\mathrm{ij}}$ .

Analysis Methods

The unmeasured cluster-level confounder $U_j$ was unavailable (i.e., omitted) in the observed data for estimation. Then, to control for the unmeasured $U_j$ , we examined the following options for including clusters in the nuisance model estimation. Option (A), labeled as “cluster means + cluster dummies,” is the method described above (in section “Estimation”). Option (B), labeled “cluster means,” omits the cluster dummies $\{S_j^{\mathrm{dumm}}{\}}_{j=1}^J$ from Option (A). In both Option (A) and Option (B), the cluster size $N_j^{\mathrm{scaled}}$ was controlled for and $W_j$ (which denotes a set of measured cluster-level covariates) was an empty set in all the nuisance models.

We also examined the following factors for estimating nuisance models. First, we compared (a) parametric estimation, labeled as “parametric,” versus (b) machine learning estimation, labeled as “nonparametric”. For the parametric estimation, we used generalized linear regressions via R function glm(). For the nonparametric estimation, we used fivefold cross-fitting to obtain the nuisance model estimates; each nuisance model was estimated by the super learner ensemble procedure (Polley & van der Laan, 2017; van der Laan et al., 2007) with generalized additive model and neural network via R package “SuperLearner” (i.e., “SL.gam” and “SL.nnet” were specified in R function “SuperLearner”). Second, we examined different orders of estimating the mediator distributions: (a) “M1 before M2” (i.e., the order described in section “Estimation”), where the joint mediator distribution is factorized as $p({M_1}_{\mathrm{ij}},{M_2}_{\mathrm{ij}}\mid a,X_{\mathrm{ij}},j)=p({M_1}_{\mathrm{ij}}\mid a,X_{\mathrm{ij}},j)p({M_2}_{\mathrm{ij}}\mid {M_1}_{\mathrm{ij}},a,X_{\mathrm{ij}},j)$ and we first estimated the distribution for ${M_1}_{\mathrm{ij}}$ and then for ${M_2}_{\mathrm{ij}}$ ; (b) “M2 before M1”, where the joint mediator distribution is factorized as $p({M_1}_{\mathrm{ij}},{M_2}_{\mathrm{ij}}\mid a,X_{\mathrm{ij}},j)=p({M_2}_{\mathrm{ij}}\mid a,X_{\mathrm{ij}},j)p({M_1}_{\mathrm{ij}}\mid {M_2}_{\mathrm{ij}},a,X_{\mathrm{ij}},j)$ and we first estimated the distribution for ${M_2}_{\mathrm{ij}}$ and then for ${M_1}_{\mathrm{ij}}$ . This simulation factor allows us to check that using the different factorization of the joint mediator distribution would not appreciably influence the simulation results.

With the nuisance model estimates yielded (i.e., from a combination of the above implementation factors), we obtained the estimate and 95% Wald-type confidence interval for each of the cluster-average and individual-average causal mediation effects, namely $\mathrm{Effect}\left(\mathrm{C}\mathrm{-}\mathrm{avg}\right)$ and $\mathrm{Effect}\left(\mathrm{I}\mathrm{-}\mathrm{avg}\right)$ with $\mathrm{Effect}\in \{\mathrm{II}{\mathrm{E}}_{M_1},\mathrm{II}{\mathrm{E}}_{M_2},\mathrm{II}{\mathrm{E}}_{\mathrm{mu}},\mathrm{IDE}\}$ .

Evaluation Criteria

Under each condition, we run 1,000 simulation replications in the R software. We evaluated the performance for each of the cluster-average and individual-average causal mediation effects (i.e., $\mathrm{Effect}\left(\mathrm{C}\mathrm{-}\mathrm{avg}\right)$ and $\mathrm{Effect}\left(\mathrm{I}\mathrm{-}\mathrm{avg}\right)$ with $\mathrm{Effect}\in \{\mathrm{II}{\mathrm{E}}_{M_1},\mathrm{II}{\mathrm{E}}_{M_2},\mathrm{II}{\mathrm{E}}_{\mathrm{mu}},\mathrm{IDE}\}$ ). For each effect, we calculated the true value by using the data-generation models to generate the potential outcomes for a large number of clusters ( $2\times {10}^5$ clusters) and combining them to calculate the effect definition. Table B1 in the Supplemental Materials (Appendix B in the online version of the journal) lists the true effect values in the non-null condition; in the null condition, the true effect values were all zero.

To evaluate the performance, we used the following criteria: (1) bias, (2) mean squared error (MSE), (3) coverage rate of 95% confidence intervals, and (4) Type I error rate with 5% nominal significance level or statistical power. Specifically, for each effect, Bias =  $\frac{\sum _{i=1}^{1,000}{\hat{\mathrm{\tau }}}_i-\mathrm{\tau }}{1,000}$ , where $\mathrm{\tau }$ denotes the true value and ${\hat{\mathrm{\tau }}}_i$ denotes the estimate from the $i$ -th simulation replication. The MSE evaluates the overall accuracy, which was MSE =  $\frac{\sum _{i=1}^{1,000}{({\hat{\tau }}_i-\tau )}^2}{1,000}$ . The coverage rate was the proportion of replications where the confidence interval covered the true value. The type I error rate (under the null condition) or power (under the non-null condition) was the proportion of replications where the 95% confidence interval excluded zero.

Simulation Results

Figures B1 to B8 in the Supplemental Materials (Appendix B in the online version of the journal) display the simulation results for the individual-average causal mediation effects (IDE and IIEs) when they were non-null or null. The results for the corresponding cluster-average effects (Figures B9–B16 in the online version of the journal) show similar patterns. We observe several patterns in the simulation results.

Under the “linear” scenario, the multiply robust method had near zero bias, confidence interval coverage rates close to the nominal 95% (ranged from 91% to 99%), and Type I error rate close to the nominal 5% (ranged from 1% to 9%), across all combinations of implementation factors. The power of testing each effect increased with increasing the number of clusters $J$ , holding constant the cluster size condition. The MSE decreased with increasing $J$ . This is so for all the causal mediation effects, including both the cluster-average and individual-average versions. Comparing the nonparametric estimators with different methods to incorporate clusters, we observe that the method “cluster means” (incorporating the cluster means while omitting the cluster dummies) led to better power and MSE than the method “cluster means + cluster dummies” (incorporating both the cluster means and cluster dummies). Nonetheless, in Figures B8 and B16 (available in the online version of this article) under the null condition, with a relatively large number of clusters ( $J=100$ ), the other estimators (the nonparametric estimator with “cluster means” and the parametric estimators) led to more anti-conservative Type I error rates (e.g., 9%), while the nonparametric estimator with “cluster means + cluster dummies” yielded more satisfactory Type I error rates (4%–6%).

Under the “nonlinear” scenario, the bias was lower for the multiply robust estimators with nonparametric nuisance estimation than with parametric nuisance estimation. Comparing the nonparametric estimators, the MSE and power were better for the nonparametric estimator incorporating the cluster means while omitting the cluster dummies (“cluster means”), than for the nonparametric estimator incorporating both the cluster means and cluster dummies (“cluster means + cluster dummies”). Nonetheless, in the presence of the unmeasured cluster-level confounder $U_j$ , the most satisfactory confidence interval coverage and Type I error rates were yielded by the nonparametric estimator incorporating both the cluster means and cluster dummies (“cluster means + cluster dummies”). In addition, we observe that in the “nonlinear” scenario, the coverage rate and bias with the nonparametric estimator could become less satisfactory as the number of clusters increased, which is unlike the pattern in the “linear” scenario. This implies that in the “nonlinear” scenario, the errors of the estimated nuisance models did not go to zero at a sufficiently fast rate as the number of independent units (i.e., clusters) increased (Rudolph et al., 2024). Despite this, in most cases, the coverage rates were above 90% with the best-performing multiply robust estimator (i.e., the nonparametric estimator with “cluster means + cluster dummies”). The order of estimating mediator distributions (i.e., estimating $M_2$ before $M_1$ , or estimating $M_1$ before $M_2$ ) did not appreciably affect the performance in general. In the Supplemental Materials, we ran additional simulations (Appendix C in the online version of the journal) where we manipulated the nuisance model estimation in the “nonlinear” scenario to satisfy certain convergence rates, following the simulation approach in Rubinstein, Branson, and Kennedy (2023); results from these simulations (Figures C1–C6 in the online version of the journal) confirmed that when the convergence rates of the nuisance model estimation are sufficiently fast to satisfy the theoretical requirement for the multiply robust estimator (i.e., the product of errors in estimates of any two nuisance models goes to zero at rate $o_p\left(J^{-1/2}\right)$ ), the bias goes down to zero as $J$ increases and the confidence interval coverage maintains around 95%.

Empirical Example

To illustrate the developed method for multiple-mediator analysis with clustered data, we assess mediators that were hypothesized to play a role in mediating the impact of adolescent extracurricular activity participation (treatment $A$ ) on math outcome ( $Y$ ) using data from the ELS of 2002 (Ingels et al., 2007). Based on substantive research (Haghighat & Knifsend, 2019), we consider two mediators related to the adolescent’s educational expectation ( $M_1$ ) and perceived importance of good grades ( $M_2$ ). (The dataset and R code used in this illustration are available at https://osf.io/skcba/?view_only=None.)

The analysis sample was 1,079 adolescents in $J=76$ public non-vocational high schools in rural areas. The eligible schools ( $J=76$ ) were schools in which there were 10 or more sampled adolescents (i.e., cluster size $N_j\ge$ 10) and the sampled adolescents had different treatment status (described below) on the extracurricular activity participation (i.e., the treatment proportion $\sum _{i=1}^{N_j}1(A_{\mathrm{ij}}=1)/N_j$ was above 0 and below 1). Across the schools, the number of sampled adolescents $N_j$ ranged from 10 to 20 with a median of 14. The schools in the analysis sample were considered to be a simple random sample from some underlying population of schools. Within each school, the sampled adolescents were considered as a simple random sample from some underlying population of adolescents served by the school, with the cluster size $N_j$ representing the relative population size of the school.

The treatment ( $A_{\mathrm{ij}}$ ) was an adolescent’s reported participation in school-sponsored extracurricular activities (e.g., school band and school academic clubs) in the 2001 to 2002 school year (1 = yes, 0 = no). The outcome ( $Y_{\mathrm{ij}}$ ) was the math standardized test score in the first follow-up (i.e., 2004). The mediators ( $M_{1\mathrm{ij}}$ and $M_{2\mathrm{ij}}$ ) were based on the adolescent’s reports in the base-year assessment in spring of 2002. The mediator of educational expectation ( $M_{1\mathrm{ij}}$ ) was whether the adolescent expected to obtain a master’s degree or above (1 = yes, 0 = no). The mediator of the perceived importance of good grades ( $M_{2\mathrm{ij}}$ ) was whether the adolescent perceived good grades as very important (1 = yes, 0 = no). In the analysis, the controlled individual-level characteristics ( $X_{\mathrm{ij}}$ ) were adolescent pretreatment characteristics, which were: The adolescent’s sex, race (Hispanic, Black, White, or other), English fluency, generational status (i.e., whether the student and mother born in the United States), Head Start program attendance, family socioeconomic status, family composition (mother and father, or other), number of siblings, parents’ highest level of education, total family income (above 50k or not), and home literacy resources. The controlled cluster-level covariate was the cluster size ( $N_j$ ).

We conducted the analysis using the developed multiply robust method. For the nuisance model estimation, we implemented the method described in the section on estimation (i.e., the method labeled as “nonparametric” and “cluster means + cluster dummies” in the simulations; the order of estimating mediator distributions was ${M_1}_{\mathrm{ij}}$ before ${M_2}_{\mathrm{ij}}$ ). In the implementation, we used fivefold cross-fitting to obtain estimates of the nuisance models. The estimation of each nuisance model was done by the super learner ensemble with generalized additive model, neural network, and generalized linear model (i.e., “SL.glm,”“SL.gam,” and “SL.nnet” via R package “SuperLearner”; van der Laan et al., 2007). Using the nuisance model estimates, we calculated the multiply robust estimator and the Wald-type 95% confidence interval for each of the estimands. Table 1 shows the results.

Table 1.

Results of the Multiple-Mediator Analysis in the Illustrative Example

Interventional Effect	Individual-Average		Cluster-Average
	Estimate	95% CI	Estimate	95% CI
Direct ( $\mathrm{IDE}$ )	1.561	[0.151, 2.971]	1.511	[0.087, 2.934]
Indirect via educational expectation ( $\mathrm{II}{\mathrm{E}}_{M_1}$ )	0.563	[0.236, 0.891]	0.565	[0.235, 0.895]
Indirect via importance of good grades ( $\mathrm{II}{\mathrm{E}}_{M_2}$ )	0.276	[0.006, 0.546]	0.271	[0.000, 0.542]
Indirect via the mediators’ mutual dependence ( $\mathrm{II}{\mathrm{E}}_{\mathrm{mu}}$ )	0.037	[−0.047, 0.121]	0.034	[−0.058, 0.126]

Note. CI = confidence interval; IDE = interventional direct effect; $\mathrm{IIE}$  = interventional indirect effects.

For the individual-average effects, there were positive estimates for the IIEs of extracurricular activity participation ( $A$ ) on math achievement ( $Y$ ) via educational expectation and via the importance of good grades (i.e., $\mathrm{II}{\mathrm{E}}_{M_1}$ and $\mathrm{II}{\mathrm{E}}_{M_2}$ ). This suggests that for an average adolescent, an increase in the math achievement outcome would be achieved, by changing the probability distribution of each of the mediators (e.g., educational expectation, $M_1$ ) from what it would be if the adolescent had no-participation (i.e., treatment status $A=0$ ) to what it would be if the adolescent had participated (i.e., $A=1$ ) (given the adolescent’s pretreatment characteristics and school membership) while fixing the treatment status to no-participation and fixing the probability distribution of the other mediator (i.e., the importance of good grades, $M_2$ ). For the IDE of extracurricular activity participation on math achievement (i.e., $\mathrm{IDE}$ ), the positive estimate indicates that, for an average adolescent, an increase in math achievement would be achieved by changing the adolescent’s treatment status from no-participation ( $A=0$ ) to participation ( $A=1$ ), while fixing the mediators’ joint probability distribution to what it would be under no-participation (given the adolescent’s pretreatment characteristics and school membership).

For the cluster-average effects, the results were similar to those for the individual-average effects. Nonetheless, we would interpret these cluster-average effects as the effects of an average school. For example, consider the cluster-average indirect effect via educational expectation ( $\mathrm{II}{\mathrm{E}}_{M_1}$ ). The positive estimate suggests that, for an average school, an increase in the school mean math achievement outcome would be achieved, by changing the probability distribution of educational expectation ( $M_1$ ) for each adolescent in the school, from what it would be if the adolescent had no-participation ( $A=0$ ) to what it would be if the adolescent had participated (i.e., $A=1$ ) (given the adolescent’s pretreatment characteristics and school membership) while fixing the treatment status at no-participation and fixing the probability distribution of the other mediator ( $M_2$ ) at what it would be under no-participation (given the adolescent’s pretreatment characteristics and school membership).

For comparison, we also obtained estimation results ignoring the clustering in the nuisance model estimation (Table D1 in Appendix D of the Supplemental Materials in the online version of the journal). That is, we omitted all the cluster dummies and cluster means and did not perform cluster-mean centering in estimating the nuisance models. In these results omitting the clusters (Table D1), the effect estimates were smaller (in magnitudes) than the corresponding estimates in Table 1 controlling for the clusters (e.g., the estimate of the cluster-average direct effect $\mathrm{IDE}$ was 1.511 in Table 1, whereas it was 0.823 in Table D1; similar patterns also apply to the estimates of the other effects). By omitting the clusters (i.e., school membership of an adolescent), the validity of these estimates (Table D1) is threatened by the possibility that certain school-level variables (e.g., school environment) could influence the adolescent’s participation in extracurricular activities (treatment), educational expectation and perceived importance of good grades (mediators), and the achievement outcome, even after conditioning on measured baseline covariates. Compared to Table 1 (which accounts for unmeasured school-level confounders), the results in Table D1 underestimated the interventional (in)direct effects of extracurricular activity participation on the achievement outcome, which could be resulting from bias due to uncontrolled school-level confounders in obtaining Table D1.

Discussion

For multiple-mediator analysis with clustered data, a common scenario in behavioral sciences, there has been limited methodological research on assessing the causal mediation effects. This study reduces this methodological gap by developing and evaluating the method to assess the cluster-average and individual-average causal mediation effects with preexisting clusters where the individual-level treatment, mediators, and outcome may be affected by unmeasured cluster-level confounders. Through simulations, we recommend using the multiply robust estimator with both the cluster means and cluster dummies included in the nuisance models, to adequately control for the unmeasured cluster-level confounder. In addition, when there is no prior knowledge of correct model specification, we recommend using data-adaptive methods (e.g., machine learning algorithms) to estimate the nuisance models. Below, we discuss considerations and extensions of the methods presented in this study.

Choosing Between the Estimands

We have provided methods for both the cluster-average and individual-average versions of the effects. In practice, the choice between the different versions would depend on which version is more relevant to the research question in the substantive context (Kahan et al., 2023; Miratrix et al., 2021). Suppose in a substantive application, the investigator is interested in the effect (one of the causal mediation effects) for all individuals (i.e., considering individuals equally), the individual-average effects would be more relevant. For instance, when the treatment is a policy affecting all individuals, a policy investigator may be interested in the policy’s impact on a typical individual (i.e., individual-average effects). In comparison, when the substantive investigator is interested in the effect of a cluster (e.g., how the school mean outcome would be affected by a school’s decision about intervention assignments for individuals in the school), the cluster-average effects would be more relevant, as they capture the effect for a typical cluster.

Limitations and Future Directions

The current study has limitations, and several topics await future research. First, for estimating the mediator distributions, it is possible to circumvent the estimation by extending recent methods (e.g., alternative representations of multiply robust estimators) in the literature on nonparametric mediation estimation (R. Liu et al., 2024; Rudolph et al., 2024). Second, future research can further evaluate and improve upon the multiply robust method with clustering presented here. For instance, alternative options for handling unmeasured cluster-level confounders await investigation, such as grouping similar clusters (Lee et al., 2021). For the variance estimation using correct parametric nuisance models, the finite sample performance could be potentially improved by incorporating uncertainty in estimating the parametric nuisance models (e.g., by stacking estimating equations of the model parameters); on this topic, recent progress has been made with weighting-based estimation (Bein et al., 2018), and future research may explore extensions to the multiply robust estimator. For the cross-fitting procedure that combines the multiply robust estimator with machine learning estimation of nuisance models, the impact of the number of folds could be worthwhile for future evaluation; while recent research found that the impact tends to be minor (Ellul et al., 2024), it is unclear whether this is generalizable to the multiple-mediator analysis with clustered data.

Third, with clustered data, the no-interference assumption would be violated when individuals sharing the same cluster membership may interact and influence each other’s outcomes (Hong & Raudenbush, 2006). In this case, we may make a cluster-level version of the SUTVA that allows for interference within the same clusters while assuming no interference between different clusters. Under the cluster-level SUTVA, it can be interesting to examine a cluster-level version of the estimands, such as causal mediation effects of the cluster-mean treatment via cluster means mediators (VanderWeele et al., 2013). Furthermore, the handling of missing data in clustered designs (Enders et al., 2020; X. Liu, 2024b) for multiple-mediator analyses could await more research.

Fourth, sensitivity analysis could be an important future research area in multiple-mediator analysis with clustered data. For unclustered data, fruitful developments have been made in quantifying the sensitivity to unmeasured confounders in mediation analyses (X. Liu et al., 2024; Park & Esterling, 2021; Rubinstein, Branson, & Kennedy, 2023; Zhang & Ding, 2023); for example, Zhang and Ding (2023) derived omitted-variable bias formulas using partial correlations, and extended the robustness value approach (Cinelli & Hazlett, 2020) to mediation analysis. It could be useful to extend the sensitivity analysis literature to the interventional (in)direct effects with clustered data.

Fifth, exploring the inference for the sample version of the effects could be an interesting direction, such as in the context of randomized experiments (Ding et al., 2017; Miratrix et al., 2021). For the sample effects, we make inferences viewing the sample at hand as the inferential population. Correspondingly, we would consider the clusters in the sample as fixed, and the uncertainty quantification would no longer incorporate the randomness due to cluster sampling. To quantify uncertainty in inference on the sample effects, we could take the randomization-based perspective (e.g., Part II in Ding, 2024). For instance, a (hypothetical) randomized experiment for assessing the interventional (in)direct effects (Moreno-Betancur & Carlin, 2018) could entail randomizing the treatment and then randomizing the mediators given the treatment, covariates, and cluster membership. Then, in making inferences for the sample effects, we would consider quantifying uncertainty in the observed outcomes with respect to randomness due to random assignments of the treatment and mediators. A related direction is to examine the inference of mediation for a finite population; for example, Abadie et al. (2023) proposed an alternative inferential framework for clustered data, allowing the sampled clusters to constitute a possibly large fraction of a finite population at hand (e.g., states in the United States). Furthermore, we have focused on two-level clustering (i.e., assuming the clusters as independent units). Nonetheless, designs with three or more levels of clustering and partially clustered designs also occur frequently in educational and behavioral research (Hedges & Borenstein, 2014; X. Liu, 2024a; X. Liu et al., 2023; Lohr et al., 2014). With more complex clustering structures, developing and evaluating methods for causal inference in multiple-mediator analyses await future studies.

Supplemental Material

sj-pdf-1-jeb-10.3102_10769986251318093 – Supplemental material for Estimating Causal Mediation Effects in Multiple-Mediator Analyses With Clustered Data

Supplemental material, sj-pdf-1-jeb-10.3102_10769986251318093 for Estimating Causal Mediation Effects in Multiple-Mediator Analyses With Clustered Data by Xiao Liu in Journal of Educational and Behavioral Statistics