Modeling Continuous Mediator Distributions

Why the Mediator Distribution Matters

An interventional mediation calculation combines two ingredients: an outcome model and a distribution from which mediator values are drawn. A good outcome model is not enough if the modeled mediator distribution gives the wrong probabilities to values that matter for the outcome. This is especially important when the outcome changes nonlinearly with the mediator: two distributions can have the same mean and variance but give quite different probabilities of crossing a clinically meaningful threshold.

The task here is to model the distribution of a continuous mediator conditional on exposure and baseline covariates. The distribution model is then used to average the outcome model over mediator values, as in g-computation. The causal identification conditions for this calculation are distinct from the distributional modeling choices: a flexible or well-fitting mediator model does not, by itself, remove confounding or establish a causal interpretation.

Definition: Let $A$ denote exposure, $M$ a continuous mediator, $C$ measured baseline covariates, and $Y$ the outcome. Write $f(m\mid A=a,C=c)$ for the conditional density of the mediator, and $\mu(a,m,c)=E[Y\mid A=a,M=m,C=c]$ for the conditional outcome mean. The interventional mean under exposure level $a$ and mediator reference exposure level $a^*$ is $$ \psi(a,a^*)=E_C\left[\int \mu(a,m,C)f(m\mid A=a^*,C)\,dm\right]. $$ It averages the outcome model at exposure $a$ over mediator values drawn from the mediator distribution associated with $a^*$, within covariate strata, and then averages over the target covariate distribution.

In applications, this formula is identified only under the relevant causal assumptions, including consistency, appropriate exchangeability conditions, and positivity. Those assumptions concern the causal interpretation and available data support. Modeling adds another requirement: the estimated density and outcome model must adequately represent the quantities being integrated. A numerical integration can be exact for the fitted models and still estimate the wrong causal quantity if those models are misspecified.

Three Ways to Model a Continuous Mediator

A Gaussian conditional model is a natural starting point when the mediator is approximately symmetric within exposure and covariate strata. A common specification is $M\mid A,C\sim N\{\eta(A,C),\sigma^2\}$, where $\eta(A,C)$ is modeled with a regression and the conditional variance is constant. The model specifies the entire distribution, not just its mean: it implies symmetric tails and, under constant variance, equal residual spread across covariate strata.

A transformed model is useful when the mediator has a restricted range or strong right skew. For a positive mediator, one option is to model $\log(M)$ as Gaussian conditional on $A,C$, then transform simulated values back to the original scale. This ensures generated mediator values are positive. It also implies a particular right-skewed distribution on the original scale. A transformation does not automatically solve distributional problems: the transformed variable still needs a suitable conditional model, and the resulting original-scale distribution may still have the wrong tails or shape.

Flexible models aim to represent features such as skewness, heavy tails, or multiple modes without forcing a single Gaussian shape. Possibilities include conditional mixture distributions, distributional regression, or empirical resampling within suitably defined covariate and exposure strata. Each approach has its own assumptions and data requirements. Flexibility is not a guarantee of accuracy: sparse data in parts of the mediator range can make the estimated distribution unstable, and a model still needs to condition on the relevant covariates.

Model	Useful feature	Key restriction to consider
Gaussian conditional model	Simple estimation and simulation; represents symmetric distributions compactly	May poorly represent skewness, heavy tails, or multiple modes
Transformed Gaussian model	Can respect a boundary such as $M>0$ and represent skewness on the original scale	Results depend on the transformation and Gaussian assumption on the transformed scale
Flexible conditional distribution	Can represent more varied shapes, including mixtures and nonlinear changes in shape	May require more data and careful control of instability or overfitting

Whichever model is used, the target is $f(m\mid A,C)$, not merely $E[M\mid A,C]$. A regression of the mediator mean is sufficient for some outcome models, but not in general. If $\mu(a,m,c)$ is nonlinear in $m$, the integral depends on more than the mediator mean.

Integrating the Outcome Model

For a fitted model, the interventional mean can be computed by drawing mediator values from the fitted conditional distribution and averaging predicted outcomes. With covariate values $C_i$ representing the target population, and $K$ simulated mediator draws for each person, a Monte Carlo approximation is

\widehat{\psi}(a,a^*)\approx \frac{1}{nK}\sum_{i=1}^n\sum_{k=1}^K \widehat{\mu}(a,M_{ik},C_i), \qquad M_{ik}\sim \widehat{f}(m\mid A=a^*,C_i).

This averages over both the modeled mediator distribution and the empirical covariate distribution. Increasing $K$ can reduce simulation noise in this calculation, but it cannot correct a misspecified density or outcome model. In simple examples, the integral can be evaluated directly; simulation is particularly useful when the outcome model or mediator distribution makes direct integration difficult.

Worked Example: A Gaussian Mediator and a Quadratic Outcome Mean

Suppose a hypothetical study has no baseline covariates. The mediator model is $M\mid A=a\sim N(2+a,1)$, so the conditional means are 2 when $A=0$ and 3 when $A=1$, with variance 1 at each exposure level. Let the outcome be continuous, with conditional mean $\mu(a,m)=2+0.4a+0.3m+0.05m^2$. Calculate $\psi(1,0)$, $\psi(0,0)$, and $\psi(1,1)$.

For a normal random variable with mean $r$ and variance 1, $E[M]=r$ and $E[M^2]=r^2+1$. Thus, under the mediator reference distribution for $A=0$, $E[M]=2$ and $E[M^2]=5$. At exposure $a=1$,

\psi(1,0)=2+0.4(1)+0.3(2)+0.05(5)=3.25.

At exposure $a=0$ using the same mediator reference distribution,

\psi(0,0)=2+0.4(0)+0.3(2)+0.05(5)=2.85.

For $\psi(1,1)$, the mediator has mean 3 and second moment $3^2+1=10$:

\psi(1,1)=2+0.4(1)+0.3(3)+0.05(10)=3.80.

The interventional direct contrast using the $A=0$ mediator distribution is $\psi(1,0)-\psi(0,0)=3.25-2.85=0.40$ outcome units. The interventional indirect contrast at $A=1$ is $\psi(1,1)-\psi(1,0)=3.80-3.25=0.55$ outcome units. The mediator variance contributes because the outcome mean contains $M^2$; simply substituting the mediator mean into the outcome model would omit that contribution.

When a Gaussian Approximation Changes the Effect

A common way to simplify a skewed mediator model is to match a Gaussian distribution to the mediator's conditional mean and variance. This preserves two summaries but not the full distribution. If the outcome response emphasizes a tail or threshold, the resulting interventional mean can differ materially from the one obtained under the skewed distribution.

Worked Example: A Lognormal Mediator and a Threshold Outcome

Consider hypothetical data with no baseline covariates. The true mediator distributions are $\log(M)\mid A=0\sim N(0,1)$ and $\log(M)\mid A=1\sim N(1,1)$. The outcome is binary, with $\Pr(Y=1\mid A=a,M=m)=0.10+0.60I(m>e)$. The risk is therefore 0.10 below the threshold and 0.70 above it. Under the true mediator model, the threshold probabilities are

P(M>e\mid A=0)=P\{\log(M)>1\mid A=0\}=1-\Phi(1)=0.1587,

and

P(M>e\mid A=1)=P\{\log(M)>1\mid A=1\}=0.5000.

Consequently, $\psi(1,0)=0.10+0.60(0.1587)=0.1952$, while $\psi(1,1)=0.10+0.60(0.5000)=0.4000$. The interventional indirect contrast at $A=1$ is $0.4000-0.1952=0.2048$, or about 20.48 percentage points in outcome risk.

Now replace each lognormal distribution with a Gaussian distribution having the same original-scale mean and variance. If $\log(M)\sim N(\nu,1)$, then $E[M]=e^{\nu+1/2}$ and $\operatorname{Var}(M)=(e-1)e^{2\nu+1}$. For $A=0$, these are 1.6487 and 4.6708, giving a moment-matched Gaussian standard deviation of 2.161. Its threshold probability is approximately $1-\Phi\{(e-1.6487)/2.161\}=0.310$. For $A=1$, the mean is 4.4817 and variance is about 34.5126, giving standard deviation 5.875 and Gaussian threshold probability approximately $1-\Phi\{(e-4.4817)/5.875\}=0.618$.

The moment-matched Gaussian model therefore gives $\psi(1,0)\approx0.10+0.60(0.310)=0.286$ and $\psi(1,1)\approx0.10+0.60(0.618)=0.471$. Its indirect contrast is approximately $0.471-0.286=0.185$, rather than the true-model value of 0.2048. These are approximations based on the unrounded tail probabilities; subtracting rounded intermediate values can change the last reported digit. The discrepancy arises because matching the first two moments does not match the probability above the outcome threshold.

Flexible Distributions and Multimodality

A flexible model can be important even when the mediator is not strongly skewed. For example, a mixture distribution may represent two common mediator ranges. A single Gaussian model can smooth over the low-density region between those ranges and assign very different probabilities to a threshold. The next example isolates that issue with an explicitly specified mixture.

Worked Example: A Mixture Mediator Near an Outcome Threshold

Suppose the hypothetical mediator distribution is a two-component Gaussian mixture. At $A=0$, half the mediator values follow $N(-1,0.2^2)$ and half follow $N(1,0.2^2)$. At $A=1$, the components are $N(-0.5,0.2^2)$ and $N(1.5,0.2^2)$, again with equal weights. Let $Y=I(M>0)$. This outcome is binary, and its conditional mean is the probability that the mediator exceeds zero.

At $A=0$, the mixture probability above zero is $0.5\{1-\Phi(5)\}+0.5\{1-\Phi(-5)\}=0.5000$, to four decimal places. At $A=1$, it is $0.5\{1-\Phi(2.5)\}+0.5\{1-\Phi(-7.5)\}\approx0.5031$. Thus the interventional contrast $\psi(1,1)-\psi(1,0)$ is about 0.0031, or 0.31 percentage points in risk.

The mixture means are 0 and 0.5. Both variances are $1+0.2^2=1.04$: the between-component variance is 1, and the within-component variance is 0.04. If each mixture is replaced with a Gaussian having the same mean and variance, the threshold probabilities become 0.5000 at $A=0$ and approximately $\Phi(0.5/\sqrt{1.04})=0.688$ at $A=1$. The Gaussian approximation therefore suggests a contrast of about 0.188, far larger than the mixture-model contrast. The two models agree on the means and variances; they differ in the shape of the distribution around the outcome threshold.

Choosing and Checking a Distribution Model

Model choice should be guided by the mediator's scale, the observed conditional distributions, and the way the outcome model uses mediator values. A Gaussian model can be adequate when conditional distributions are roughly symmetric and its spread is reasonable across exposure and covariate levels. A transformation is attractive when it respects a meaningful boundary or reduces skewness, but the transformed-scale model still needs scrutiny. Flexible models can represent complex shapes, but they need enough observations to estimate those shapes reliably.

The practical target is not distributional fit in the abstract. It is accurate estimation of the integrals used in $\psi(a,a^*)$. In particular, assess whether a model represents the mediator values where $\mu(a,m,c)$ changes sharply, as at a threshold. Comparing estimated conditional means and variances is not enough when tail probabilities or multimodality affect the outcome mean. It is also useful to compare plausible distributional models and examine whether the resulting interventional contrasts are sensitive to that choice.

Keep the modeling layers distinct. The outcome regression estimates how the mean outcome varies with exposure, mediator, and covariates; the mediator model supplies the distribution used for averaging. Both can be misspecified. Correctly specifying one does not compensate automatically for misspecifying the other, and neither modeling choice replaces the causal assumptions required for identification.

Common mistake: “The Gaussian mediator model matches the observed mean and variance, so it must give the same interventional effect.” This is wrong whenever the outcome model depends on other features of the mediator distribution. The correction is to evaluate the outcome integral using the full modeled distribution and to pay particular attention to tails, boundaries, and thresholds that influence the outcome.

Check Your Understanding

Use these questions to check whether you can connect mediator-distribution modeling to interventional-effect calculations.

In the definition of $\psi(a,a^*)$, which exposure level determines the mediator distribution, and which exposure level is used in the outcome model?
Why can two mediator distributions with the same mean and variance yield different interventional means?
What restriction does a Gaussian conditional mediator model place on the shape of the distribution?
For a positive mediator, what feature does modeling $\log(M)$ as Gaussian automatically preserve when simulated values are transformed back?
In the mixture example, why did the moment-matched Gaussian model give a different threshold probability even though its mean and variance matched those of the mixture?