Worked Multiple-Mediator Analysis With Dependence

Why the Relationship Between Mediators Matters

When an exposure may affect two mediators, analyzing each mediator separately can miss how they occur together. For example, two post-exposure health measures may be positively associated under one exposure level and negatively associated under another. If the outcome responds to their combination, an intervention that preserves that association can produce a different mean outcome from one that draws the mediators independently.

This tutorial works through both interventions. The example distinguishes changes in the mediators’ marginal distributions from changes in their dependence, then shows how those differences enter interventional direct and indirect effects. All data and fitted quantities below are hypothetical.

Estimands and Assumptions

Definition: Let $A$ be a binary exposure taking values 0 and 1, let $M_1$ and $M_2$ be two mediators measured after exposure, and let $C$ be a vector of pre-exposure covariates. Write $q_b(m_1,m_2\mid c)$ for the joint conditional distribution of the mediator vector under exposure level $b$. A joint mediator intervention draws the pair $(M_1,M_2)$ together from $q_b(\cdot,\cdot\mid c)$, preserving its dependence. An independent-draw intervention draws $M_1$ and $M_2$ separately from their respective marginals under $b$; its distribution is $q_b^I(m_1,m_2\mid c)=q_{b1}(m_1\mid c)q_{b2}(m_2\mid c)$.

Let $\mu(a,m_1,m_2,c)=E[Y\mid A=a,M_1=m_1,M_2=m_2,C=c]$ be the outcome regression. For a mediator distribution $q$, define the standardized mean $\Psi(a;q)$ as the mean outcome if exposure were set to $a$, while the mediator pair were drawn from $q$, averaging over the target population’s distribution of $C$:

\Psi(a;q)=\int \sum_{m_1,m_2}\mu(a,m_1,m_2,c)q(m_1,m_2\mid c)\,dF_C(c).

The sum is over mediator values when they are discrete; for continuous mediators it is an integral. The notation $q_b^J$ will denote the joint distribution $q_b$, and $q_b^I$ the product of its marginal distributions. The interventional direct effect at mediator reference distribution $q_b^K$ is $\Psi(1;q_b^K)-\Psi(0;q_b^K)$; the interventional indirect effect at exposure level $a$ is $\Psi(a;q_1^K)-\Psi(a;q_0^K)$. Here $K$ is either $J$ or $I$. Each contrast compares means under specified exposure and mediator interventions.

A causal interpretation requires, among other conditions, consistency, a well-defined exposure and mediator intervention, and no interference between individuals. It also requires adequate positivity: the exposure and mediator values needed for each contrast must have support in the relevant covariate strata. For identification from observed data, assume no unmeasured confounding of exposure and outcome or of exposure and the mediator vector given $C$, and no unmeasured mediator-outcome confounding given $A,C$. In this simple mediation setup, assume there is no exposure-induced mediator-outcome confounder omitted from the model. These assumptions allow the standardized mean above to be identified using the observed outcome regression and mediator distributions; they are not verified merely by fitting models.

Estimating the Standardized Means

A g-computation analysis estimates the outcome regression and the mediator distributions conditional on baseline covariates, then averages predicted outcomes over the desired intervention distribution. For the joint intervention, estimate the joint mediator distribution, including dependence. For the independent-draw intervention, estimate each marginal distribution and multiply them within each covariate stratum. In either case, use the same target distribution of $C$ when comparing means.

Estimate the required models.
Fit $\mu(a,m_1,m_2,c)$ and the exposure-specific joint mediator distribution. Check that model forms and support are adequate for the values used in standardization.

Specify the intervention.
Use the joint distribution to preserve dependence, or multiply the two fitted marginal distributions to define an independent-draw intervention.

Average predicted outcomes.
For each target individual or covariate profile, average the outcome predictions over the selected mediator distribution, then average over the target population.

Form the contrasts.
Subtract standardized means to obtain the direct effect at a chosen mediator distribution and the indirect effect at a chosen exposure level. Keep the intervention type and reference distribution explicit.

Worked Example: Standardizing Over a Joint Mediator Distribution

Suppose a hypothetical analysis has binary mediators $M_1,M_2$, no baseline covariates, and an outcome measured in points. Its assumed outcome regression is

\mu(a,m_1,m_2)=10+2a+3m_1+4m_2+6m_1m_2+2am_1+4am_1m_2.

The following hypothetical exposure-specific mediator probabilities are estimated. Each row gives the probability of the mediator pair in the order shown.

Exposure level	$P(M_1=0,M_2=0)$	$P(M_1=1,M_2=0)$	$P(M_1=0,M_2=1)$	$P(M_1=1,M_2=1)$
0	0.50	0.20	0.20	0.10
1	0.05	0.25	0.25	0.45

Both rows sum to 1, and all four mediator pairs have positive probability at each exposure level. Under $q_0$, each mediator has marginal probability $0.30$ of being 1, and $P(M_1=1,M_2=1)=0.10$. Under $q_1$, each has marginal probability $0.70$, and the joint probability of both being 1 is $0.45$.

For binary mediators, let $p=E[M_1]=E[M_2]$ in this example and $r=E[M_1M_2]=P(M_1=1,M_2=1)$. Averaging the outcome regression gives

\Psi(a;q)=10+2a+7p+6r+2ap+4ar.

For example, under exposure level 0 and the joint distribution $q_0$, substitute $a=0$, $p=0.30$, and $r=0.10$:

\Psi(0;q_0^J)=10+7(0.30)+6(0.10)=12.70\text{ points}.

Applying the same formula to all four combinations gives:

\begin{aligned} \Psi(1;q_0^J)&=10+2+7(0.30)+6(0.10)+2(0.30)+4(0.10)=15.70,\\ \Psi(0;q_1^J)&=10+7(0.70)+6(0.45)=17.60,\\ \Psi(1;q_1^J)&=10+2+7(0.70)+6(0.45)+2(0.70)+4(0.45)=22.80. \end{aligned}

Thus, for the joint interventions, the direct effect at $q_0^J$ is $15.70-12.70=3.00$ points, and the indirect effect at exposure level 1 is $22.80-15.70=7.10$ points. The total contrast between the two exposure-specific joint mediator regimes is $22.80-12.70=10.10$ points.

What Independent Draws Change

The independent-draw intervention retains the exposure-specific marginal distributions but replaces the joint distribution with the product of those marginals. It therefore changes the probability of mediator combinations unless the mediators are conditionally independent given $A,C$. The intervention is not simply another way to estimate the same joint-intervention mean: it specifies a different distribution of mediator pairs.

Worked Example: Calculating Independent-Draw Effects

Use the same hypothetical outcome regression and mediator marginals. Under $q_0$, each mediator has probability $0.30$ of being 1, so the independent probabilities for $(0,0),(1,0),(0,1),(1,1)$ are $0.49,0.21,0.21,0.09$. Under $q_1$, each marginal probability is $0.70$, giving independent probabilities $0.09,0.21,0.21,0.49$. These probabilities are the products of the corresponding Bernoulli marginals.

For $q_0^I$, $p=0.30$ and $r=0.30\times0.30=0.09$. For $q_1^I$, $p=0.70$ and $r=0.70\times0.70=0.49$. Substitution into the standardized-mean formula yields:

\begin{aligned} \Psi(0;q_0^I)&=10+7(0.30)+6(0.09)=12.64\text{ points},\\ \Psi(1;q_0^I)&=10+2+7(0.30)+6(0.09)+2(0.30)+4(0.09)=15.60\text{ points},\\ \Psi(0;q_1^I)&=10+7(0.70)+6(0.49)=17.84\text{ points},\\ \Psi(1;q_1^I)&=10+2+7(0.70)+6(0.49)+2(0.70)+4(0.49)=23.20\text{ points}. \end{aligned}

The independent-draw direct effect at $q_0^I$ is $15.60-12.64=2.96$ points, while the independent-draw indirect effect at exposure level 1 is $23.20-15.60=7.60$ points. The exposure-specific independent-draw regime contrast is $23.20-12.64=10.56$ points. Each calculation uses the same outcome scale and target population as the joint analysis, but a different mediator intervention.

Deriving the Dependence Contribution

The example’s interaction terms make the role of dependence explicit. For binary mediators, their covariance under a distribution $q$ is $E[M_1M_2]-E[M_1]E[M_2]$. The independent-draw distribution has the same marginal means as the joint distribution but has covariance zero. In the example, the outcome regression includes the term $(6+4a)M_1M_2$, so only this term distinguishes the joint and independent standardized means when their marginals match.

Derivation: Write $p_1=E[M_1]$, $p_2=E[M_2]$, and $r=E[M_1M_2]$ under the joint distribution. The product distribution has $E[M_1M_2]=p_1p_2$. The additive terms in the outcome regression have identical expectations under both distributions, while the interaction term contributes $(6+4a)r$ under the joint distribution and $(6+4a)p_1p_2$ under independent draws. Therefore, $$ \Psi(a;q^J)-\Psi(a;q^I)=(6+4a)(r-p_1p_2)=(6+4a)\operatorname{Cov}(M_1,M_2). $$

Under $q_0$, the covariance is $0.10-(0.30)(0.30)=0.01$. The joint mean therefore exceeds the independent mean by $(6+4a)(0.01)$: $0.06$ points at $a=0$ and $0.10$ points at $a=1$. Under $q_1$, the covariance is $0.45-(0.70)(0.70)=-0.04$. The joint mean is lower than the independent mean by $0.24$ points at $a=0$ and $0.40$ points at $a=1$. These differences match the four pairs of standardized means calculated above.

The change in dependence also changes the effects, not just the component means. At exposure level 1, the joint indirect effect is $7.10$ points, compared with $7.60$ points under independent draws. The joint direct effect at the $q_0$ reference is $3.00$ points, compared with $2.96$ points under independent draws. The direct-effect difference arises here because the outcome’s mediator interaction varies with exposure through the $4aM_1M_2$ term. With no exposure interaction in the outcome’s dependence-sensitive terms, the direct-effect contrast at a fixed reference distribution would not change in this particular model.

For either intervention type, the direct and indirect effects add to its corresponding regime contrast when they use compatible reference choices. For example, at the joint intervention, $3.00+7.10=10.10$ points; at the independent intervention, $2.96+7.60=10.56$ points. This is an algebraic decomposition of the specified contrasts, not a claim that the two intervention types represent the same total effect.

Worked Example: Checking the Effect Decomposition

Use the joint means from the first calculation to verify the alternative decomposition, which uses the $q_1^J$ mediator distribution for the direct effect and exposure level 0 for the indirect effect.

The direct effect at $q_1^J$ is $22.80-17.60=5.20$ points. The indirect effect at exposure level 0 is $17.60-12.70=4.90$ points. Their sum is $5.20+4.90=10.10$ points, equal to the joint regime contrast $22.80-12.70=10.10$ points. The direct and indirect components differ from the earlier decomposition because the reference mediator distribution and exposure level differ, even though their sum is unchanged.

For independent draws, the analogous check gives a direct effect at $q_1^I$ of $23.20-17.84=5.36$ points and an indirect effect at exposure level 0 of $17.84-12.64=5.20$ points. Their sum, $10.56$ points, equals the independent-draw regime contrast. The arithmetic confirms the decomposition but does not make its components reference-free.

Interpretation and Common Pitfalls

The joint intervention answers a question about shifting the mediator vector while retaining how its components co-occur under each exposure level. The independent-draw intervention answers a different question: what would the mean be if each mediator were shifted to its exposure-specific marginal distribution, but the pair were generated independently within covariate strata? Such a regime may be scientifically useful when separate components can be manipulated separately, but it should not be described as the natural joint mediator distribution.

A common mistake is to calculate a product of marginal probabilities and call it the joint mediator distribution. That is correct only under conditional independence of the mediators given exposure and baseline covariates. If the mediators are dependent, the product changes the probability of combinations and can change outcome means, as the example demonstrates. Another mistake is to interpret a difference between joint and independent estimates as proof that dependence itself is causal. The contrast depends on the specified intervention and outcome model, and causal interpretation still rests on the identification assumptions.

In applied work, report the intervention distribution, the reference exposure or mediator distribution used for each component effect, and the outcome units. Use a joint mediator model capable of representing dependence when estimating joint interventions; multiplying separately fitted marginals cannot recover the joint distribution. Uncertainty calculations should reflect estimation of both mediator and outcome models, and diagnostics should examine fit and support for mediator combinations used in standardization.

Key takeaway: Joint and independent-draw interventions can share mediator marginals while assigning different probabilities to mediator combinations. When outcomes respond to those combinations, the standardized means and interventional effects can differ. Name the intervention and reference distribution whenever reporting a multiple-mediator effect.

Check Your Understanding

Use the definitions and calculations above to answer these questions.

What feature of the mediator distribution is preserved by a joint intervention but generally changed by an independent-draw intervention?
If each binary mediator has marginal probability $0.40$ of being 1, what is the probability of $(M_1,M_2)=(1,1)$ under independent draws?
In the worked model, why does the covariance between the mediators affect $\Psi(a;q)$? Identify the relevant outcome-regression terms.
For the joint example, verify the indirect effect at exposure level 0 using the two standardized means at $a=0$.
Why is the independent-draw regime contrast not automatically the natural total effect of exposure?