Computing Natural Effects From Nonlinear Outcome Models

From Logistic Models to Natural Effects

A logistic-regression coefficient describes a change in log odds conditional on the variables in the model. It is not, by itself, a natural direct or indirect effect. Natural effects instead compare outcome risks predicted under specific combinations of treatment and mediator values, while accounting for how the mediator is distributed under each treatment level. Because logistic models predict probabilities nonlinearly, those predictions must be averaged on the risk scale before the effects are contrasted.

This tutorial develops that computation for a binary treatment $A$, a binary mediator $M$, and a binary outcome $Y$. Baseline covariates $L$ can be included by calculating within covariate strata and then averaging over a stated target population. The examples use hypothetical fitted-model quantities; they are calculations, not claims about real studies.

Definition: Let $M(a)$ be the mediator value that would be observed under treatment $a$, and let $Y(a,m)$ be the outcome that would be observed under treatment $a$ if the mediator were set to $m$. Define the natural-effect risk $$ R(a,a^*)=E\{Y(a,M(a^*))\}. $$ It is the mean outcome under treatment $a$, with the mediator taking the value it would naturally have taken under treatment $a^*$. With binary $A$, the natural direct effect at mediator reference level $a^*$ is $R(1,a^*)-R(0,a^*)$; the natural indirect effect at outcome-treatment level $a$ is $R(a,1)-R(a,0)$. These are risk differences when $Y$ is binary.

A commonly used additive decomposition is $R(1,1)-R(0,0)=\{R(1,0)-R(0,0)\}+\{R(1,1)-R(1,0)\}$. The first term is a natural direct effect using the mediator distribution under $A=0$; the second is a natural indirect effect evaluated at $A=1$. With treatment–mediator interaction, the effects depend on these reference choices, so state the particular contrasts being reported.

The Plug-In Mediation Calculation

For baseline covariates $L$, the mediation formula expresses the natural-effect risk as an average of outcome risks over the mediator distribution under $a^*$. For a binary mediator, its conditional version is

R(a,a^*)= E_L\left[ \sum_{m=0}^{1} E(Y\mid A=a,M=m,L)\, P(M=m\mid A=a^*,L) \right].

The outer expectation averages over the baseline covariate distribution in the target population. In estimation, replace the conditional outcome mean and mediator probability with fitted-model predictions, calculate the sum for each target-population covariate value, and average those sums. If the target is the study population, its observed baseline covariate distribution is often used for this final average.

This formula is identified as a causal natural-effect risk only under appropriate assumptions. These include consistency; positivity for the treatment and mediator values needed in the calculation; no unmeasured confounding of treatment and outcome, treatment and mediator, or mediator and outcome conditional on baseline covariates; and the cross-world independence condition required for natural effects. The simple formula also presumes there is no exposure-induced mediator–outcome confounder that must be handled. Model-based estimates additionally require adequate specification of the conditional outcome and mediator models, or other valid ways to estimate those quantities. The computation does not repair a failure of causal identification.

For the examples, an outcome logistic model has conditional risk $p(a,m,l)=\operatorname{expit}\{\eta(a,m,l)\}$, where $\operatorname{expit}(x)=1/(1+\exp(-x))$. A mediator logistic model can provide $q(a,l)=P(M=1\mid A=a,L=l)$. The contribution at covariate value $l$ is then

r(a,a^*\mid l) = p(a,1,l)q(a^*,l)+p(a,0,l)\{1-q(a^*,l)\}.

The order matters: first transform the outcome-model linear predictor to a probability, then average those probabilities using the mediator distribution. Averaging linear predictors or odds and treating the result as a risk is generally incorrect.

Worked Example: Binary Mediator Without an Interaction

Suppose hypothetical fitted models for a population without covariates give $P(M=1\mid A=0)=0.25$ and $P(M=1\mid A=1)=0.60$. The fitted outcome model is $\operatorname{logit}P(Y=1\mid A,M)=-2+0.8A+1.1M$. Assume the causal identification conditions just stated are appropriate for this illustration, and use the fitted models as the estimation models.

The four outcome risks are obtained by applying the expit function to each linear predictor:

\begin{aligned} p(0,0)&=\operatorname{expit}(-2)=0.1192,\\ p(1,0)&=\operatorname{expit}(-1.2)=0.2315,\\ p(0,1)&=\operatorname{expit}(-0.9)=0.2891,\\ p(1,1)&=\operatorname{expit}(-0.1)=0.4750. \end{aligned}

For $R(0,0)$, predict outcomes at $A=0$ and average over the mediator distribution under $A=0$: $R(0,0)=0.75(0.1192)+0.25(0.2891)=0.1617$. For $R(1,0)$, predict at $A=1$ but retain that same mediator distribution: $R(1,0)=0.75(0.2315)+0.25(0.4750)=0.2924$.

Using the mediator distribution under $A=1$, the other two risks are $R(0,1)=0.40(0.1192)+0.60(0.2891)=0.2211$ and $R(1,1)=0.40(0.2315)+0.60(0.4750)=0.3776$. Thus the natural direct effect using the $A=0$ mediator distribution is $R(1,0)-R(0,0)=0.2924-0.1617=0.1307$, or about 13.1 percentage points. The natural indirect effect at $A=1$ is $R(1,1)-R(1,0)=0.3776-0.2924=0.0852$, or about 8.52 percentage points. Their sum, 0.2159, equals the total natural-effect contrast $R(1,1)-R(0,0)=0.3776-0.1617=0.2159$, up to rounding.

Why the Nonlinear Model Requires Prediction and Averaging

In a logistic model, the coefficient for treatment is a conditional log-odds contrast when the other modeled variables are held fixed. The natural-effect calculation asks a different question: what risk results after combining model predictions across mediator values and, where relevant, baseline covariates? The expit transformation is nonlinear, so transforming an averaged linear predictor will generally not equal averaging the transformed predictions.

When the model includes an $A\times M$ interaction, outcome predictions must be made for all treatment–mediator combinations separately. A single treatment coefficient cannot represent the treatment contrast at every mediator value, and a single mediator coefficient cannot represent the mediator contrast at every treatment level.

Worked Example: Treatment–Mediator Interaction

Consider another hypothetical population with $P(M=1\mid A=0)=0.20$ and $P(M=1\mid A=1)=0.70$. Suppose the fitted outcome model is $\operatorname{logit}P(Y=1\mid A,M)=-2+0.6A+1.0M+0.9AM$. The causal and model assumptions for plug-in estimation are taken to hold for this example.

The interaction changes the linear predictor for the joint treatment–mediator condition. The four predicted risks are $p(0,0)=\operatorname{expit}(-2)=0.1192$, $p(1,0)=\operatorname{expit}(-1.4)=0.1978$, $p(0,1)=\operatorname{expit}(-1)=0.2689$, and $p(1,1)=\operatorname{expit}(0.5)=0.6225$. For example, at $A=1,M=1$, the predictor is $-2+0.6+1.0+0.9=0.5$.

Averaging over the mediator distribution under $A=0$ gives $R(0,0)=0.80(0.1192)+0.20(0.2689)=0.1492$ and $R(1,0)=0.80(0.1978)+0.20(0.6225)=0.2827$. Averaging over the distribution under $A=1$ gives $R(0,1)=0.30(0.1192)+0.70(0.2689)=0.2240$ and $R(1,1)=0.30(0.1978)+0.70(0.6225)=0.4951$.

The direct effect using the $A=0$ mediator distribution is $0.2827-0.1492=0.1336$, or 13.4 percentage points. The indirect effect at $A=1$ is $0.4951-0.2827=0.2123$, or 21.2 percentage points. These sum to the total contrast $0.4951-0.1492=0.3459$, or 34.6 percentage points.

The other reference choices give different effects: the direct effect using the $A=1$ mediator distribution is $R(1,1)-R(0,1)=0.2710$, while the indirect effect evaluated at $A=0$ is $R(0,1)-R(0,0)=0.0749$. This does not contradict the first decomposition. Each pair sums to the same total contrast, but interaction makes the component effects depend on which treatment level is held fixed.

Standardizing Over Baseline Covariates

When baseline covariates predict the mediator or outcome, compute the mediator probabilities and outcome risks conditional on those covariates. Then combine them within each covariate value and average using the target population's covariate distribution. This standardization is not the same as reading a coefficient from a regression table. It makes explicit which population's average risk is being reported.

Worked Example: Averaging Over Two Baseline Strata

Suppose a hypothetical target population has equal numbers in two baseline strata, $L=0$ and $L=1$. Fitted mediator models give $P(M=1\mid A=0,L=0)=0.20$, $P(M=1\mid A=1,L=0)=0.50$, $P(M=1\mid A=0,L=1)=0.40$, and $P(M=1\mid A=1,L=1)=0.80$. The outcome model is $\operatorname{logit}P(Y=1\mid A,M,L)=-2+0.5A+0.8M+0.4L+0.6AM$. Assume the identification conditions and model specification are adequate.

For $L=0$, the four risks in order $(a,m)=(0,0),(1,0),(0,1),(1,1)$ are $0.1192, 0.1824, 0.2315,$ and $0.4750$. For $L=1$, they are $0.1680, 0.2497, 0.3100,$ and $0.5744$. Each value is the expit of the model predictor at that combination; for example, at $A=1,M=1,L=1$, the predictor is $-2+0.5+0.8+0.4+0.6=0.3$, giving risk $0.5744$.

For $L=0$, averaging over the mediator distribution under $A=0$ gives $r(0,0\mid 0)=0.80(0.1192)+0.20(0.2315)=0.1417$ and $r(1,0\mid 0)=0.80(0.1824)+0.20(0.4750)=0.2409$. Using the distribution under $A=1$ gives $r(0,1\mid 0)=0.50(0.1192)+0.50(0.2315)=0.1753$ and $r(1,1\mid 0)=0.50(0.1824)+0.50(0.4750)=0.3287$.

For $L=1$, the corresponding calculations are $r(0,0\mid 1)=0.60(0.1680)+0.40(0.3100)=0.2248$, $r(1,0\mid 1)=0.60(0.2497)+0.40(0.5744)=0.3796$, $r(0,1\mid 1)=0.20(0.1680)+0.80(0.3100)=0.2816$, and $r(1,1\mid 1)=0.20(0.2497)+0.80(0.5744)=0.5095$.

The target population gives each stratum weight $0.50$. Therefore, $R(0,0)=0.50(0.1417)+0.50(0.2248)=0.1832$, $R(1,0)=0.50(0.2409)+0.50(0.3796)=0.3103$, $R(0,1)=0.50(0.1753)+0.50(0.2816)=0.2285$, and $R(1,1)=0.50(0.3287)+0.50(0.5095)=0.4191$. The direct effect using the $A=0$ mediator distribution is $0.3103-0.1832=0.1271$, or 12.7 percentage points. The indirect effect at $A=1$ is $0.4191-0.3103=0.1088$, or 10.9 percentage points. Their sum is the total contrast $0.4191-0.1832=0.2359$, or 23.6 percentage points. The reported risks describe the specified equally weighted target population; a different target distribution could yield different marginal effects.

Calculation Logic and Common Pitfalls

The risk calculation follows directly by conditioning on the mediator. For binary $M(a^*)$, its probability of being one is $q(a^*,l)$, and its probability of being zero is $1-q(a^*,l)$. Thus, within $L=l$, the expectation of the outcome prediction is the sum of its two possible values weighted by those probabilities. Averaging this conditional quantity over $L$ yields $R(a,a^*)$. Subtracting two such risks gives the stated risk-difference effect.

The decomposition is a telescoping identity, not an additional modeling assumption:

\begin{aligned} \{R(1,0)-R(0,0)\}+\{R(1,1)-R(1,0)\} &=R(1,1)-R(0,0). \end{aligned}

The intermediate term $R(1,0)$ cancels. The algebra guarantees the decomposition for the chosen risks, but it does not guarantee that those risks have a causal interpretation; that still depends on identification and estimation assumptions.

A frequent mistake is to report the logistic treatment coefficient as the natural direct effect, or the mediator coefficient as the natural indirect effect. For example, a treatment coefficient of $0.6$ in a logistic model is a conditional log-odds coefficient, not a risk difference of 0.6 and not a natural direct effect. The correction is to predict risks for each needed $(A,M,L)$ combination, average them over the appropriate mediator distribution and target covariate distribution, and then subtract the relevant marginal risks.

Do not average on the wrong scale. Apply the expit function to each predictor before averaging outcome risks.
Do not omit the mediator distribution. A predicted risk at one fixed mediator value is not $R(a,a^*)$.
Do not ignore interaction terms. Use the full model predictor for every treatment–mediator combination.
Do not leave the reference choices implicit. Report whether the direct effect uses $M(0)$ or $M(1)$, and whether the indirect effect is evaluated at $A=0$ or $A=1$.
Do not mistake computation for identification. Correct arithmetic cannot compensate for unmeasured confounding, unsupported treatment–mediator combinations, or an inadequate model.

Key takeaway: Natural effects from a nonlinear outcome model are calculated by predicting outcome risks under each treatment–mediator combination, averaging those risks over the mediator distribution specified by the natural effect, standardizing over the target baseline population, and finally taking the required risk differences.

Check Your Understanding

Use the definitions and calculation sequence above to answer these questions.

In $R(1,0)$, which treatment value is used for the outcome prediction, and which treatment value determines the mediator distribution?
Why must a logistic linear predictor be transformed to a probability before averaging over mediator values?
With treatment–mediator interaction, why can the natural direct effect using $M(0)$ differ from the natural direct effect using $M(1)$?
What does averaging stratum-specific natural-effect risks over a target covariate distribution determine?
If fitted models give valid predictions but mediator–outcome confounding remains unmeasured, does the plug-in calculation alone identify a causal natural effect? Explain.