Causal Mediation: Why a x b Stops Working

อ่านฉบับภาษาไทย (Thai version)

Abstract

Mediation analysis asks how much of an exposure's effect on an outcome runs through a mediator, a variable on the causal path between them. The familiar recipe multiplies two slopes. The product a x b equals the natural indirect effect only when both models are linear and there is no exposure-mediator interaction, and it has a causal meaning only if there is no unmeasured exposure-outcome, mediator-outcome or exposure-mediator confounding and no mediator-outcome confounder that is itself affected by the exposure. In a simulated cohort of 2,000 adults, exercise lowers blood pressure partly through weight change, and each kilogram matters more with more exercise. For 2.5 versus 0 hours of weekly exercise, the true natural indirect effect in the simulated population is -4.95 mmHg, while the product method targets -4.06 mmHg. This article defines the effects with potential outcomes, gives closed forms that keep the interaction, shows Stata and R code with delta-method intervals, covers a common binary outcome and sets out the cross-world assumption.

🔍 Click to enlarge Click anywhere to close

A product of two slopes, and a reviewer's question

A community health team has run an exercise programme for a year. Among 2,000 adults, those who exercised more had lower systolic blood pressure (SBP) at the end of the year, and they had also lost more weight. The team wants to know how much of the fall in pressure works through weight loss. The programme is fictional and its data are simulated.

The analyst reaches for a familiar recipe. Regress weight change on exercise, regress blood pressure on exercise and weight change, and multiply the two slopes. The product is reported as the indirect effect, and the exercise slope of the second model as the direct effect.

A reviewer reads the draft and points to a scatter plot. Among people who exercise more, each kilogram of weight lost goes with a larger fall in pressure. Exercise and weight change interact on the mmHg scale of blood pressure, and the reviewer asks which slope for weight belongs in the product.

This article answers that question: it pins down when a product of two coefficients is the indirect effect, and how to compute the effects when it is not. The effects are first defined as comparisons between hypothetical worlds, then computed from models that keep the interaction. That route also shows which assumptions the answer rests on.

The product method and where it holds

Write $A$ for the exposure, hours of exercise per week, and $M$ for the mediator, weight change in kg over the year, where a negative value is a loss. $Y$ is the outcome, SBP in mmHg, and $C$ stands for baseline covariates, here age and sex, that may confound these relations. A mediator is a variable on the causal path from the exposure to the outcome: exercise changes weight, and weight changes pressure.

The classic regression approach to mediation fits linear models for the mediator and for the outcome [1]. The mediator model regresses $M$ on $A$ and $C$; its exercise slope, $\beta_1$, is the change in mean weight change per extra weekly hour. The outcome model regresses $Y$ on $A$, $M$ and $C$ with no product term. Its weight slope, $b$, is the change in mean SBP per kg, and its exercise slope, $c'$, is read as the direct effect per weekly hour.

Here $a^*$ and $a$ are the two exercise levels being compared, $a^* = 0$ and $a = 2.5$ hours per week, so the contrast width $a - a^*$ is 2.5 hours. The direct effect over the same contrast is $c' \times (a - a^*)$. In the familiar a x b recipe, the exposure-to-mediator slope is called $a$; this article keeps $a$ for an exposure level and writes that slope as $\beta_1$.

The product a x b equals the natural indirect effect only when both models are linear and there is no exposure-mediator interaction, and it has a causal meaning only if there is no unmeasured exposure-outcome, mediator-outcome or exposure-mediator confounding and no mediator-outcome confounder that is itself affected by the exposure [2]. The natural indirect effect gets its formal definition below. In the linear case with no interaction on the mmHg difference scale, exercise shifts mean weight by $\beta_1$ per weekly hour, every kilogram shifts mean pressure by the same $b$, and the two shifts multiply. The signs matter as much as the sizes, as the hand example shows.

The product method in the simulated cohort

The simulated cohort used throughout this article holds 2,000 adults aged 30 to 80. For each person it records weekly exercise, weight change over one year, SBP at one year, age and sex. A hat marks an estimate: fitting the two linear models with age and sex as covariates gives $\hat\beta_1 = -1.19$ kg per weekly hour and $\hat b = 1.37$ mmHg per kg. The 95% intervals printed in the two outputs below cover the values these slopes take in the simulated population.

Over the 2.5-hour contrast, the product of the two slopes and the width is -4.09 mmHg, with a 95% confidence interval of -4.81 to -3.37 mmHg from the delta method. The delta method gives a standard error for a function of estimated coefficients, such as a product, by treating the function as nearly linear close to the estimates. For a product of two slopes it is also known as the Sobel standard error. Every interval in the two result tables is a delta-method 95% interval; the outputs of Stata's mediate command and of the R mediation package use other interval methods, as their captions say.

Stata: the mediator model

* Mediator model: weight change on exercise and the confounders (age, sex)
regress wtchg exercise age male

. * Mediator model: weight change on exercise and the confounders (age, sex)
. regress wtchg exercise age male

      Source |       SS           df       MS      Number of obs   =     2,000
-------------+----------------------------------   F(3, 1996)      =    412.53
       Model |  7937.34882         3  2645.78294   Prob > F        =    0.0000
    Residual |  12801.4031     1,996  6.41352861   R-squared       =    0.3827
-------------+----------------------------------   Adj R-squared   =    0.3818
       Total |  20738.7519     1,999  10.3745632   Root MSE        =    2.5325

------------------------------------------------------------------------------
       wtchg | Coefficient  Std. err.      t    P>|t|     [95% conf. interval]
-------------+----------------------------------------------------------------
    exercise |  -1.194437   .0435007   -27.46   0.000    -1.279749   -1.109126
         age |   .0236362   .0044603     5.30   0.000     .0148888    .0323836
        male |   .2156961   .1139233     1.89   0.058     -.007725    .4391173
       _cons |  -.1197085   .3035733    -0.39   0.693    -.7150622    .4756452
------------------------------------------------------------------------------

Stata: the outcome model without the interaction term

* Product method: outcome model without the interaction gives b, then a x b over the 2.5-hour contrast
regress sbp exercise wtchg age male

. * Product method: outcome model without the interaction gives b, then a x b over the 2.5-hour contrast
. regress sbp exercise wtchg age male

      Source |       SS           df       MS      Number of obs   =     2,000
-------------+----------------------------------   F(4, 1995)      =    403.71
       Model |  263039.698         4  65759.9245   Prob > F        =    0.0000
    Residual |  324960.139     1,995  162.887288   R-squared       =    0.4473
-------------+----------------------------------   Adj R-squared   =    0.4462
       Total |  587999.837     1,999  294.146992   Root MSE        =    12.763

------------------------------------------------------------------------------
         sbp | Coefficient  Std. err.      t    P>|t|     [95% conf. interval]
-------------+----------------------------------------------------------------
    exercise |  -1.979195   .2573194    -7.69   0.000    -2.483838   -1.474552
       wtchg |   1.369749   .1128015    12.14   0.000     1.148528     1.59097
         age |   .4485983   .0226358    19.82   0.000      .404206    .4929905
        male |   3.890667   .5746421     6.77   0.000     2.763706    5.017629
       _cons |   103.2818   1.529944    67.51   0.000     100.2814    106.2823
------------------------------------------------------------------------------

Three ways the product method breaks

The recipe assumes one number for the effect of weight on pressure. To let that effect change with exercise, add a product term to the outcome model:

Here $E[Y \mid a, m, c]$ is the mean SBP at exercise level $a$, weight change $m$ and covariates $c$. The coefficient $\theta_1$ is the exercise slope at zero weight change, $\theta_2$ is the weight slope at zero exercise, and $\theta_4$ holds the covariate coefficients. The coefficient $\theta_3$ is the exposure-mediator interaction: how much the weight slope changes per extra weekly hour of exercise. The weight slope at exercise level $a$ is therefore $\theta_2 + \theta_3 a$.

The simulated cohort was generated with $\theta_2 = 0.4$ and $\theta_3 = 0.5$. Each kilogram of weight change therefore moves SBP by 0.4 mmHg at no exercise and by $0.4 + 0.5 \times 2.5 = 1.65$ mmHg at 2.5 hours a week. Neither number is the $b$ of the product method. A fit without the product term blends the weight slopes across the exercise levels in the cohort, and in a very large simulated sample $b$ settles at 1.35 mmHg per kg.

That blend is the first failure. The product method's target in the simulated population is -4.06 mmHg, while the natural indirect effect, defined in the next section, is -4.95 mmHg. The gap comes from using one averaged weight slope where the indirect effect needs the slope at the exercise level being compared.

The second failure is a non-linear outcome model. For a yes-or-no outcome such as hypertension, the outcome model is usually a logistic regression, whose coefficients are changes in log odds. A product of $\beta_1$ and a log-odds coefficient is not an effect on the risk-difference scale. Multiplied by the contrast width and exponentiated, that product approximates the natural indirect effect on the odds-ratio scale, but only when the outcome is rare and the logistic model has no exposure-mediator product term, that is, no interaction on the log-odds scale.

The third failure is a binary mediator, such as whether a person joins a salt-reduction plan during the programme. Its model is then logistic too, and $\beta_1$ becomes a change in the log odds of the mediator, not a shift in its mean. The indirect effect then depends on how far the probability of the mediator moves. That distance depends on each person's covariates, because the same change in log odds moves the probability by different amounts at different starting probabilities.

Fitted to the simulated cohort, the outcome model with the product term gives $\hat\theta_3 = 0.55$ mmHg per kg per weekly hour, with $\hat\theta_1 = -0.84$ mmHg per weekly hour and $\hat\theta_2 = 0.29$ mmHg per kg. The 95% intervals in the output below cover the values the data were generated from, 0.5, -1 and 0.4, so each estimate is within sampling error of its generating value.

Stata: the outcome model with the exposure-mediator interaction

* Outcome model with the exposure-mediator interaction (theta1 A + theta2 M + theta3 A M)
regress sbp exercise wtchg exw age male

. * Outcome model with the exposure-mediator interaction (theta1 A + theta2 M + theta3 A M)
. regress sbp exercise wtchg exw age male

      Source |       SS           df       MS      Number of obs   =     2,000
-------------+----------------------------------   F(5, 1994)      =    358.74
       Model |  278452.465         5  55690.4931   Prob > F        =    0.0000
    Residual |  309547.372     1,994  155.239404   R-squared       =    0.4736
-------------+----------------------------------   Adj R-squared   =    0.4722
       Total |  587999.837     1,999  294.146992   Root MSE        =     12.46

------------------------------------------------------------------------------
         sbp | Coefficient  Std. err.      t    P>|t|     [95% conf. interval]
-------------+----------------------------------------------------------------
    exercise |  -.8351658   .2762007    -3.02   0.003    -1.376838   -.2934935
       wtchg |   .2898677   .1545066     1.88   0.061    -.0131437    .5928791
         exw |   .5457527   .0547717     9.96   0.000     .4383369    .6531686
         age |   .4668235   .0221736    21.05   0.000     .4233377    .5103092
        male |   3.880244   .5609905     6.92   0.000     2.780054    4.980433
       _cons |   101.6947   1.502064    67.70   0.000     98.74896    104.6405
------------------------------------------------------------------------------

Naming the effects with potential outcomes

A potential outcome is the value an outcome would take if the exposure, and possibly the mediator, were set to given levels [3, 4]. Write $M^{a}$ for the weight change a person would have with exercise set to $a$. Write $Y^{a,m}$ for the SBP they would have with exercise set to $a$ and weight change set to $m$. For any person, at most one of these values is observed, so the definitions below compare averages, written $E[\cdot]$.

Nesting the two gives $Y^{a,M^{a^*}}$: SBP with exercise at $a$, but with weight change at the level it would have reached with exercise at $a^*$. With exercise at $a$, weight change takes its natural value $M^{a}$, so $Y^{a} = Y^{a,M^{a}}$. Each natural effect below changes one thing at a time.

The controlled direct effect fixes weight change at the same value $m$ for everyone and changes only exercise. It answers what the extra exercise would do if weight change could be held at $m$ by intervention.

The natural direct effect changes exercise from $a^*$ to $a$ but leaves each person's weight change where it would be without the extra exercise. Only the paths that bypass weight are switched on.

The natural indirect effect keeps exercise at $a$ and moves each person's weight change from its level under $a^*$ to its level under $a$. Only the path through weight is switched on.

The total effect compares the two exercise levels and lets weight respond freely. The two natural effects add up to it, because the middle term $Y^{a,M^{a^*}}$ appears in both and cancels.

The proportion mediated is $\mathrm{NIE}/\mathrm{TE}$, the indirect effect as a share of the total on the same scale. The section 'The proportion mediated, read with care' shows why it needs care. In a linear model with no interaction on the mmHg scale, the controlled direct effect is the same at every $m$ and equals the natural direct effect; with an interaction, they generally differ.

The quantity $Y^{a,M^{a^*}}$ mixes two worlds: exercise as in one, weight change as in the other. It is never observed for anyone. That is what makes natural effects informative about mechanism. It is also what makes them hard to identify, that is, to compute from the observed data under stated assumptions (see 'What it takes to read these effects causally').

Closed forms when exercise changes the weight effect

With a linear mediator model and a linear outcome model that keeps the product term, each effect has a closed form, a formula in the model coefficients [2, 5]. The mediator model is

where $\beta_0$ is the intercept, $\beta_1$ the exercise slope and $\beta_2$ the covariate coefficients. With the outcome model above, the effects for the contrast from $a^*$ to $a$, at covariate values $c$, are as follows.

The controlled direct effect depends on the level $m$ at which weight change is fixed, because the exercise slope at that level is $\theta_1 + \theta_3 m$.

The natural direct effect is the controlled direct effect evaluated at the mean weight change under $a^*$, which is $\beta_0 + \beta_1 a^* + \beta_2' c$.

The natural indirect effect multiplies the shift in mean weight change, $\beta_1 (a - a^*)$, by the weight slope at exercise level $a$, which is $\theta_2 + \theta_3 a$.

Set $\theta_3 = 0$ and the formulas collapse to the product method. The natural direct effect becomes $\theta_1 (a - a^*)$, and the natural indirect effect becomes $\theta_2 \beta_1 (a - a^*)$, with $\theta_2$ then equal to $b$. With $\theta_3 \ne 0$, the indirect effect uses the weight slope at $a$, and the direct effect uses the mean weight change at $a^*$. Because both models are linear, plugging in the covariate means gives the average effect over the cohort.

The simulated cohort: estimates beside the simulated truth

The same closed forms, applied to the models fitted to the 2,000 simulated adults, give the estimates in the table. Both models include age and sex, which enter the formulas at their sample means. Each interval combines the two models' estimated variances through the delta method.

Every delta-method interval in the table covers its simulated truth. In this sample the natural indirect effect lands close to its simulated truth, -4.94 mmHg against -4.95. The natural direct effect lands further away, -0.32 mmHg against -1.00. It is a small difference between two larger estimated terms, so its delta-method interval, -1.87 to 1.23 mmHg, is wide.

The product method, fitted to the same data, gives -4.09 mmHg with a delta-method interval of -4.81 to -3.37. That interval covers the product method's own target, -4.06, and excludes the simulated true natural indirect effect, -4.95. A larger sample would only tighten the interval around the wrong target.

The panes below show the closed forms in code. Stata's nlcom command, short for nonlinear combinations of estimates, applies the delta method once the two fits are stacked; the R version takes numerical derivatives of the same formulas. Stata's built-in mediate command, the user-written paramed command and the R mediation package, which implements the simulation approach of Imai and colleagues [6], give the same or very close estimates.

Continuous outcome: simulated truth, estimate and delta-method interval

Stata: the closed forms and their delta-method intervals with nlcom

* Stack both models (block-diagonal variance) so nlcom can apply the delta method
matrix coleq bm = med
matrix coleq by = out
matrix bb = bm, by
local names : colfullnames bb
matrix VV = (Vm, J(4, 6, 0) \ J(6, 4, 0), Vy)
matrix colnames VV = `names'
matrix rownames VV = `names'
ereturn post bb VV

* Closed forms (VanderWeele), covariates at their means; m0 = E[M | a* = 0, mean covariates]
local m0 "(_b[med:_cons] + _b[med:age]*(`agebar') + _b[med:male]*(`malebar'))"
local nde "((_b[out:exercise] + _b[out:exw]*`m0')*2.5)"
local nie "((_b[out:wtchg] + _b[out:exw]*2.5)*_b[med:exercise]*2.5)"
local te "(`nde' + `nie')"
nlcom (nde: `nde') (nie: `nie') (te: `te') (pm: `nie'/`te') (cde0: _b[out:exercise]*2.5) (pnie: _b[out:wtchg]*_b[med:exercise]*2.5) (tnde: (_b[out:exercise] + _b[out:exw]*(`m0' + _b[med:exercise]*2.5))*2.5)

. * Stack both models (block-diagonal variance) so nlcom can apply the delta method
. matrix coleq bm = med

. matrix coleq by = out

. matrix bb = bm, by

. local names : colfullnames bb

. matrix VV = (Vm, J(4, 6, 0) \ J(6, 4, 0), Vy)

. matrix colnames VV = `names'

. matrix rownames VV = `names'

. ereturn post bb VV

.
. * Closed forms (VanderWeele), covariates at their means; m0 = E[M | a* = 0, mean covariates]
. local m0 "(_b[med:_cons] + _b[med:age]*(`agebar') + _b[med:male]*(`malebar'))"

. local nde "((_b[out:exercise] + _b[out:exw]*`m0')*2.5)"

. local nie "((_b[out:wtchg] + _b[out:exw]*2.5)*_b[med:exercise]*2.5)"

. local te "(`nde' + `nie')"

. nlcom (nde: `nde') (nie: `nie') (te: `te') (pm: `nie'/`te') (cde0: _b[out:exercise]*2.5) (pnie: _b[out:wtchg]*_b[med:e
> xercise]*2.5) (tnde: (_b[out:exercise] + _b[out:exw]*(`m0' + _b[med:exercise]*2.5))*2.5)

         nde: ((_b[out:exercise] + _b[out:exw]*(_b[med:_cons] + _b[med:age]*(55.313) + _b[med:male]*(.504)))*2.5)
         nie: ((_b[out:wtchg] + _b[out:exw]*2.5)*_b[med:exercise]*2.5)
          te: (((_b[out:exercise] + _b[out:exw]*(_b[med:_cons] + _b[med:age]*(55.313) + _b[med:male]*(.504)))*2.5) + ((_
> b[out:wtchg] + _b[out:exw]*2.5)*_b[med:exercise]*2.5))
          pm: ((_b[out:wtchg] + _b[out:exw]*2.5)*_b[med:exercise]*2.5)/(((_b[out:exercise] + _b[out:exw]*(_b[med:_cons]
> + _b[med:age]*(55.313) + _b[med:male]*(.504)))*2.5) + ((_b[out:wtchg] + _b[out:exw]*2.5)*_b[med:exercise]*2.5))
        cde0: _b[out:exercise]*2.5
        pnie: _b[out:wtchg]*_b[med:exercise]*2.5
        tnde: (_b[out:exercise] + _b[out:exw]*((_b[med:_cons] + _b[med:age]*(55.313) + _b[med:male]*(.504)) + _b[med:exe
> rcise]*2.5))*2.5

------------------------------------------------------------------------------
             | Coefficient  Std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
         nde |  -.3191431    .792894    -0.40   0.687    -1.873187    1.234901
         nie |  -4.939743   .3844027   -12.85   0.000    -5.693158   -4.186327
          te |  -5.258886   .6638817    -7.92   0.000     -6.56007   -3.957701
          pm |   .9393135   .1441129     6.52   0.000     .6568574     1.22177
        cde0 |  -2.087915   .6905018    -3.02   0.002    -3.441273   -.7345559
        pnie |  -.8655719   .4624469    -1.87   0.061    -1.771951    .0408073
        tnde |  -4.393314   .6362454    -6.91   0.000    -5.640332   -3.146296
------------------------------------------------------------------------------

R: the same closed forms with a numerical delta method

# Closed forms (VanderWeele), covariates at their means; m0 = E[M | a* = 0, mean covariates]
effects <- function(t) {
  m0 <- t[["med_(Intercept)"]] + t[["med_age"]] * agebar + t[["med_male"]] * malebar
  b1 <- t[["med_exercise"]]
  nde <- (t[["out_exercise"]] + t[["out_exw"]] * m0) * 2.5
  nie <- (t[["out_wtchg"]] + t[["out_exw"]] * 2.5) * b1 * 2.5
  c(nde = nde, nie = nie, te = nde + nie, pm = nie / (nde + nie),
    cde0 = t[["out_exercise"]] * 2.5,
    pnie = t[["out_wtchg"]] * b1 * 2.5,
    tnde = (t[["out_exercise"]] + t[["out_exw"]] * (m0 + b1 * 2.5)) * 2.5)
}
# Delta method with a central-difference Jacobian
est <- effects(th)
J <- sapply(seq_along(th), function(j) {
  h <- 1e-6 * max(1, abs(th[j]))
  up <- th; dn <- th
  up[j] <- up[j] + h
  dn[j] <- dn[j] - h
  (effects(up) - effects(dn)) / (2 * h)
})
se <- sqrt(diag(J %*% V %*% t(J)))
for (k in c("nde", "nie", "te", "pm", "cde0")) {
  canon(k, est[[k]])
  canon(paste0(k, ".lb"), est[[k]] - z975 * se[[k]])
  canon(paste0(k, ".ub"), est[[k]] + z975 * se[[k]])
}
canon("pnie", est[["pnie"]])
canon("tnde", est[["tnde"]])

> effects <- function(t) {
+     m0 <- t[["med_(Intercept)"]] + t[["med_age"]] * agebar +
+         t[["med_male"]] * malebar
+     b1 <- t[["med_exercise"]]
+     nde <- (t[["out_exercise"]] + t[["out_exw"]] * m0) * 2.5
+     nie <- (t[["out_wtchg"]] + t[["out_exw"]] * 2.5) * b1 * 2.5
+     c(nde = nde, nie = nie, te = nde + nie, pm = nie/(nde + nie),
+         cde0 = t[["out_exercise"]] * 2.5, pnie = t[["out_wtchg"]] *
+             b1 * 2.5, tnde = (t[["out_exercise"]] + t[["out_exw"]] *
+             (m0 + b1 * 2.5)) * 2.5)
+ }

> est <- effects(th)

> J <- sapply(seq_along(th), function(j) {
+     h <- 1e-06 * max(1, abs(th[j]))
+     up <- th
+     dn <- th
+     up[j] <- up[j] + h
+     dn[j] <- dn[j] - h
+     (effects(up) - effects(dn))/(2 * h)
+ })

> se <- sqrt(diag(J %*% V %*% t(J)))

> for (k in c("nde", "nie", "te", "pm", "cde0")) {
+     canon(k, est[[k]])
+     canon(paste0(k, ".lb"), est[[k]] - z975 * se[[k]])
+     canon(paste0(k, ".ub"), est[[k]] + z975 * se[[k]])
+ }
CANON w6.nde -0.3191
CANON w6.nde.lb -1.8732
CANON w6.nde.ub 1.2349
CANON w6.nie -4.9397
CANON w6.nie.lb -5.6932
CANON w6.nie.ub -4.1863
CANON w6.te -5.2589
CANON w6.te.lb -6.5601
CANON w6.te.ub -3.9577
CANON w6.pm 0.9393
CANON w6.pm.lb 0.6569
CANON w6.pm.ub 1.2218
CANON w6.cde0 -2.0879
CANON w6.cde0.lb -3.4413
CANON w6.cde0.ub -0.7346

> canon("pnie", est[["pnie"]])
CANON w6.pnie -0.8656

> canon("tnde", est[["tnde"]])
CANON w6.tnde -4.3933

Stata: the built-in mediate command

capture noisily mediate (sbp age male) (wtchg age male) (exercise, continuous(0 2.5)), nie nde pnie tnde te

. capture noisily mediate (sbp age male) (wtchg age male) (exercise, continuous(0 2.5)), nie nde pnie tnde te

Iteration 0:  EE criterion =  2.075e-26
Iteration 1:  EE criterion =  9.065e-28

Causal mediation analysis                                Number of obs = 2,000

Outcome model:     Linear
Mediator model:    Linear
Mediator variable: wtchg
Treatment type:    Continuous

Continuous treatment levels:
  0: exercise =   0 (control)
  1: exercise = 2.5

------------------------------------------------------------------------------------
                   |               Robust
               sbp | Coefficient  std. err.      z    P>|z|     [95% conf. interval]
-------------------+----------------------------------------------------------------
NIE                |
          exercise |
         (1 vs 0)  |  -4.939743   .3481939   -14.19   0.000     -5.62219   -4.257295
-------------------+----------------------------------------------------------------
NDE                |
          exercise |
         (1 vs 0)  |  -.3191431   .7607765    -0.42   0.675    -1.810238    1.171951
-------------------+----------------------------------------------------------------
PNIE               |
          exercise |
         (1 vs 0)  |  -.8655719    .481407    -1.80   0.072    -1.809112    .0779685
-------------------+----------------------------------------------------------------
TNDE               |
          exercise |
         (1 vs 0)  |  -4.393314   .5850942    -7.51   0.000    -5.540077    -3.24655
-------------------+----------------------------------------------------------------
TE                 |
          exercise |
         (1 vs 0)  |  -5.258886   .6896579    -7.63   0.000     -6.61059   -3.907181
------------------------------------------------------------------------------------

Stata: the user-written paramed command

    capture noisily paramed sbp, avar(exercise) mvar(wtchg) cvars(age male) a0(0) a1(2.5) m(0) yreg(linear) mreg(linear)

.     capture noisily paramed sbp, avar(exercise) mvar(wtchg) cvars(age male) a0(0) a1(2.5) m(0) yreg(linear) mreg(linea
> r)

      Source |       SS           df       MS      Number of obs   =     2,000
-------------+----------------------------------   F(5, 1994)      =    358.74
       Model |  278452.465         5  55690.4931   Prob > F        =    0.0000
    Residual |  309547.372     1,994  155.239404   R-squared       =    0.4736
-------------+----------------------------------   Adj R-squared   =    0.4722
       Total |  587999.837     1,999  294.146992   Root MSE        =     12.46

-----------------------------------------------------------------------------------
              sbp | Coefficient  Std. err.      t    P>|t|     [95% conf. interval]
------------------+----------------------------------------------------------------
         exercise |  -.8351658   .2762007    -3.02   0.003    -1.376838   -.2934936
            wtchg |   .2898677   .1545066     1.88   0.061    -.0131437    .5928791
_exercise_X_wtchg |   .5457527   .0547717     9.96   0.000     .4383369    .6531685
              age |   .4668235   .0221736    21.05   0.000     .4233377    .5103092
             male |   3.880244   .5609906     6.92   0.000     2.780054    4.980433
            _cons |   101.6947   1.502064    67.70   0.000     98.74896    104.6405
-----------------------------------------------------------------------------------

      Source |       SS           df       MS      Number of obs   =     2,000
-------------+----------------------------------   F(3, 1996)      =    412.53
       Model |  7937.34882         3  2645.78294   Prob > F        =    0.0000
    Residual |  12801.4031     1,996  6.41352861   R-squared       =    0.3827
-------------+----------------------------------   Adj R-squared   =    0.3818
       Total |  20738.7519     1,999  10.3745632   Root MSE        =    2.5325

------------------------------------------------------------------------------
       wtchg | Coefficient  Std. err.      t    P>|t|     [95% conf. interval]
-------------+----------------------------------------------------------------
    exercise |  -1.194437   .0435007   -27.46   0.000    -1.279749   -1.109126
         age |   .0236362   .0044603     5.30   0.000     .0148888    .0323836
        male |   .2156961   .1139233     1.89   0.058     -.007725    .4391173
       _cons |  -.1197085   .3035733    -0.39   0.693    -.7150622    .4756452
------------------------------------------------------------------------------

             |   Estimate     Std Err   P>|z|   [95% Conf   Interval]
-------------+-------------------------------------------------------
         cde | -2.0879146   .69050183   0.002  -3.4412982    -.734531
         nde | -.31914317   .79289403   0.687  -1.8732155   1.2349291
         nie | -4.9397425   .38440274   0.000  -5.6931719  -4.1863132
         mte | -5.2588857    .6638817   0.000  -6.5600938  -3.9576776

cde:controlled direct effect, nde:natural direct effect, nie:natural indirect effect, mte:marginal total effect

R: the mediation package

# Cross-check with the mediation package (quasi-Bayesian simulation, 1000 draws)
set.seed(202664)
fy_int <- lm(sbp ~ exercise * wtchg + age + male, data = d)
med <- mediation::mediate(fm, fy_int, treat = "exercise", mediator = "wtchg",
                          control.value = 0, treat.value = 2.5, sims = 1000)
# d1 = indirect effect with exposure at a (NIE); z0 = direct effect with M at M(a*) (NDE)
canon("mediation.nie.r", med$d1)
canon("mediation.nie.lb.r", med$d1.ci[1])
canon("mediation.nie.ub.r", med$d1.ci[2])
canon("mediation.nde.r", med$z0)
canon("mediation.nde.lb.r", med$z0.ci[1])
canon("mediation.nde.ub.r", med$z0.ci[2])
canon("mediation.pnie.r", med$d0)
canon("mediation.tnde.r", med$z1)
canon("mediation.te.r", med$tau.coef)
canon("mediation.te.lb.r", med$tau.ci[1])
canon("mediation.te.ub.r", med$tau.ci[2])

> set.seed(202664)

> fy_int <- lm(sbp ~ exercise * wtchg + age + male,
+     data = d)

> med <- mediation::mediate(fm, fy_int, treat = "exercise",
+     mediator = "wtchg", control.value = 0, treat.value = 2.5,
+     sims = 1000)

> canon("mediation.nie.r", med$d1)
CANON w6.mediation.nie.r -4.9542

> canon("mediation.nie.lb.r", med$d1.ci[1])
CANON w6.mediation.nie.lb.r -5.7222

> canon("mediation.nie.ub.r", med$d1.ci[2])
CANON w6.mediation.nie.ub.r -4.1960

> canon("mediation.nde.r", med$z0)
CANON w6.mediation.nde.r -0.3032

> canon("mediation.nde.lb.r", med$z0.ci[1])
CANON w6.mediation.nde.lb.r -1.8452

> canon("mediation.nde.ub.r", med$z0.ci[2])
CANON w6.mediation.nde.ub.r 1.3059

> canon("mediation.pnie.r", med$d0)
CANON w6.mediation.pnie.r -0.8804

> canon("mediation.tnde.r", med$z1)
CANON w6.mediation.tnde.r -4.3771

> canon("mediation.te.r", med$tau.coef)
CANON w6.mediation.te.r -5.2574

> canon("mediation.te.lb.r", med$tau.ci[1])
CANON w6.mediation.te.lb.r -6.5889

> canon("mediation.te.ub.r", med$tau.ci[2])
CANON w6.mediation.te.ub.r -3.9108

A yes-or-no outcome: the scale question returns

The simulated cohort also records hypertension at one year, a yes-or-no outcome, for the same 2,000 people with the same exercise and weight change. Hypertension is common, with a prevalence of 0.312 in the sample and 0.316 in the simulated population. The outcome model is now a logistic regression with the exercise-by-weight product term, while the mediator model stays linear.

On the risk-difference scale, each natural effect is a difference between two average risks. With a logistic outcome model there is no simple closed form, so the effects are computed with the mediational g-formula. For each person, the fitted outcome model predicts the risk of hypertension with exercise set to $a^*$ or $a$ and weight change spread over its fitted distribution under $a^*$ or $a$; these risks are then averaged over people. The three effects compare the averaged risks in pairs, exactly as in the definitions above.

On the odds-ratio scale, closed forms exist for a logistic outcome with a linear mediator [2]. They rest on the rare-outcome approximation: they treat odds ratios as if they were risk ratios, which is close only when the outcome is uncommon at every level of the exposure, the mediator and the covariates. Hypertension in this cohort is not rare.

Hypertension: simulated truth, estimate and delta-method interval

Reading the two scales

On the risk-difference scale, the delta-method intervals for the natural direct, natural indirect and total effects each cover their simulated truth. The estimated natural direct effect, -0.068, sits far from -0.015, again because the direct part is the noisier one. The delta-method interval for the proportion mediated, 0.31 to 0.87, misses the simulated truth of 0.89. A 95% interval is expected to miss now and then, and a ratio of two noisy differences is where the delta method's straight-line approximation is weakest.

On the odds-ratio scale, the rare-outcome closed forms give 1.03 for the natural direct effect and 0.51 for the natural indirect effect in the simulated population. The exact odds ratios there, at the mean age and sex mix, are 0.93 and 0.53. With this common outcome, the approximation even puts the direct effect on the wrong side of 1.

The sample estimates, 0.76 (delta-method interval 0.53 to 1.08) and 0.58 (delta-method interval 0.49 to 0.68), estimate the approximate quantities, not the exact ones. Each interval covers its rare-outcome target, 1.03 and 0.51. At this sample size they also cover the exact values, so the data alone would not show that the formulas aim at the wrong quantity.

Stata: the logistic outcome model with the interaction

* Outcome model: logistic regression with the exposure-mediator interaction
logit htn exercise wtchg exw age male

. * Outcome model: logistic regression with the exposure-mediator interaction
. logit htn exercise wtchg exw age male

Iteration 0:  Log likelihood = -1241.3831
Iteration 1:  Log likelihood = -1015.5316
Iteration 2:  Log likelihood = -1000.2137
Iteration 3:  Log likelihood = -999.57991
Iteration 4:  Log likelihood = -999.57787
Iteration 5:  Log likelihood = -999.57787

Logistic regression                                     Number of obs =  2,000
                                                        LR chi2(5)    = 483.61
                                                        Prob > chi2   = 0.0000
Log likelihood = -999.57787                             Pseudo R2     = 0.1948

------------------------------------------------------------------------------
         htn | Coefficient  Std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
    exercise |  -.1819008   .0538403    -3.38   0.001    -.2874259   -.0763757
       wtchg |   .1060281   .0314984     3.37   0.001     .0442925    .1677638
         exw |   .0314067   .0156946     2.00   0.045     .0006459    .0621675
         age |   .0514611   .0045331    11.35   0.000     .0425763    .0603459
        male |   .1396875   .1104594     1.26   0.206    -.0768089    .3561839
       _cons |  -3.431735   .3109377   -11.04   0.000    -4.041162   -2.822308
------------------------------------------------------------------------------

R: the risk-difference effects by the mediational g-formula

# Mediational g-formula on the risk-difference scale: average over people and 20 normal quantiles of M
zq <- qnorm((seq_len(20) - 0.5) / 20)
i <- rep(seq_len(nrow(db)), each = 20)
z <- rep(zq, times = nrow(db))
m0 <- cm[["(Intercept)"]] + cm[["age"]] * db$age[i] + cm[["male"]] * db$male[i] + sig * z
m1 <- m0 + cm[["exercise"]] * 2.5
xc <- g[["(Intercept)"]] + g[["age"]] * db$age[i] + g[["male"]] * db$male[i]
r00 <- mean(plogis(xc + g[["wtchg"]] * m0))
r10 <- mean(plogis(xc + g[["exercise"]] * 2.5 + g[["wtchg"]] * m0 + g[["exw"]] * 2.5 * m0))
r11 <- mean(plogis(xc + g[["exercise"]] * 2.5 + g[["wtchg"]] * m1 + g[["exw"]] * 2.5 * m1))
canon("bin.rd_nde", r10 - r00)
canon("bin.rd_nie", r11 - r10)
canon("bin.rd_te", r11 - r00)
canon("bin.rd_pm", (r11 - r10) / (r11 - r00))
# delta-method 95% CIs on the risk-difference scale (NDE, NIE, TE, proportion mediated)
for (k in 1:4) {
  e <- c("rd_nde", "rd_nie", "rd_te", "rd_pm")[k]
  canon(paste0("bin.", e, ".lb"), f0b[k] - z975 * seb[k])
  canon(paste0("bin.", e, ".ub"), f0b[k] + z975 * seb[k])
}

> zq <- qnorm((seq_len(20) - 0.5)/20)

> i <- rep(seq_len(nrow(db)), each = 20)

> z <- rep(zq, times = nrow(db))

> m0 <- cm[["(Intercept)"]] + cm[["age"]] * db$age[i] +
+     cm[["male"]] * db$male[i] + sig * z

> m1 <- m0 + cm[["exercise"]] * 2.5

> xc <- g[["(Intercept)"]] + g[["age"]] * db$age[i] +
+     g[["male"]] * db$male[i]

> r00 <- mean(plogis(xc + g[["wtchg"]] * m0))

> r10 <- mean(plogis(xc + g[["exercise"]] * 2.5 + g[["wtchg"]] *
+     m0 + g[["exw"]] * 2.5 * m0))

> r11 <- mean(plogis(xc + g[["exercise"]] * 2.5 + g[["wtchg"]] *
+     m1 + g[["exw"]] * 2.5 * m1))

> canon("bin.rd_nde", r10 - r00)
CANON w6.bin.rd_nde -0.0679

> canon("bin.rd_nie", r11 - r10)
CANON w6.bin.rd_nie -0.0988

> canon("bin.rd_te", r11 - r00)
CANON w6.bin.rd_te -0.1667

> canon("bin.rd_pm", (r11 - r10)/(r11 - r00))
CANON w6.bin.rd_pm 0.5925

> for (k in 1:4) {
+     e <- c("rd_nde", "rd_nie", "rd_te", "rd_pm")[k]
+     canon(paste0("bin.", e, ".lb"), f0b[k] - z975 * seb[k])
+     canon(paste0("bin.", e, ".ub"), f0b[k] + z975 * seb[k])
+ }
CANON w6.bin.rd_nde.lb -0.1317
CANON w6.bin.rd_nde.ub -0.0042
CANON w6.bin.rd_nie.lb -0.1300
CANON w6.bin.rd_nie.ub -0.0676
CANON w6.bin.rd_te.lb -0.2170
CANON w6.bin.rd_te.ub -0.1164
CANON w6.bin.rd_pm.lb 0.3116
CANON w6.bin.rd_pm.ub 0.8735

What it takes to read these effects causally

The definitions above are causal, and the regressions are not. Moving from one to the other needs four assumptions, with $C$ the measured covariates [2, 5]:

1. No unmeasured confounding of exposure and outcome, given $C$.
2. No unmeasured confounding of mediator and outcome, given $A$ and $C$.
3. No unmeasured confounding of exposure and mediator, given $C$.
4. The cross-world assumption: given $C$, the potential outcome $Y^{a,m}$ is independent of the potential mediator $M^{a^*}$ under the other exposure level, written $Y^{a,m} \perp M^{a^*} \mid C$.

The controlled direct effect needs only the first two. The natural effects need all four, because they contain $Y^{a,M^{a^*}}$, which joins exercise at $a$ to weight change from a world where exercise was $a^*$.

The fourth assumption fails whenever a mediator-outcome confounder is itself affected by the exposure, even a measured one. Suppose exercise improves diet quality, and diet affects both weight and blood pressure. Adjusting for diet then blocks part of the effect of exercise, while leaving it out leaves weight and pressure confounded. Natural effects are then not identified, though related interventional effects, which set the mediator to a random draw from its distribution under each exposure level, can be [7].

Why the cross-world assumption cannot be checked

Even when no such variable exists, the cross-world assumption cannot be checked from data. No person can be observed with exercise at $a$ and, at the same time, with the weight change they would have had at $a^*$. No experiment can produce that pairing either. Randomising exercise secures the first and third assumptions, but not the second or the fourth.

In the simulated cohort all four assumptions hold by construction. Age and sex are the only common causes and both are measured. No variable affected by exercise confounds weight change and pressure, and the random parts of exercise, weight change and blood pressure were drawn independently. Real data offer no such guarantee, which is why the sentence below has two halves.

Sensitivity analysis for unmeasured mediator-outcome confounding

Randomising the exposure secures the first and third assumptions but leaves the second and the fourth open. People who lose more weight may differ in other ways that also lower their pressure, such as sleep or alcohol intake. A sensitivity analysis asks how strong such an unmeasured confounder, $U$, would have to be to change the conclusion.

Two methods are in common use. Bias formulas [5] ask the analyst to state how strongly $U$ affects blood pressure and how its distribution differs between exercise levels among people with the same weight change. They then return natural direct and indirect effects corrected for that much confounding. The approach of Imai and colleagues [6] uses the correlation between the error terms of the mediator and outcome models, which is zero when no such confounder exists, and reports how large that correlation would have to be for the indirect effect to vanish.

Working through either method is beyond this article. Neither method shows that no such confounder exists. Each turns an untestable assumption into a statement a reader can judge, and each is best planned before the data are analysed.

A workflow from question to report

The steps below put the definitions, models and assumptions in the order a protocol needs them.

1. Draw the DAG. Place exposure, mediator and outcome, add the confounders of each of the three relations, and mark any mediator-outcome confounder that the exposure affects.
2. Name the estimand, the exact quantity to be estimated: the effect (controlled at a stated $m$, or natural), the two exposure levels $a^*$ and $a$, the scale and the population.
3. Fit the models with the interaction. Keep the exposure-mediator product term in the outcome model unless a reason to drop it is stated in advance.
4. Compute the effects. Use closed forms when both models are linear, and the odds-ratio forms only for a rare outcome; otherwise use the mediational g-formula or simulation.
5. Attach intervals. Delta-method intervals are quick. Bootstrap intervals, which refit the models on many resampled copies of the data, are a common alternative, especially for ratios such as the proportion mediated.
6. Run a sensitivity analysis for unmeasured mediator-outcome confounding.
7. Report the DAG, the estimand, both models, the effects with their intervals and the assumptions, following the AGReMA reporting guideline for mediation analyses [8].

This order mirrors causal inference in general: the question and the estimand come first, and the models serve them.

The proportion mediated, read with care

The proportion mediated is $\mathrm{NIE}/\mathrm{TE}$ on one scale. In the hand example it is 0.83: of the 5.95 mmHg total fall in SBP, 4.95 mmHg runs through weight in the natural-effects split. Three properties limit how far that number can be read.

First, with an exposure-mediator interaction the split is not unique [3]. The total effect also equals the pure indirect effect (PNIE), with exercise held at $a^*$ while weight change shifts, plus the total direct effect (TNDE), with weight change held at its level under $a$. In this language, the natural direct and indirect effects used so far are the pure direct and the total indirect effect.

Here, with $a^* = 0$, $\mathrm{PNIE} = \theta_2 \beta_1 (a - a^*)$, which is $0.4 \times (-3.0) = -1.20$ mmHg, uses the weight slope at no exercise. The total direct effect uses the mean weight change under $a$, $1.2 + (-1.2)(2.5) = -1.8$ kg. Its exercise slope is then $-1 + 0.5 \times (-1.8) = -1.9$ and $\mathrm{TNDE} = -1.9 \times 2.5 = -4.75$ mmHg. Measured with the pure indirect effect, the proportion mediated is $-1.20/-5.95 = 0.20$, not 0.83, from the same model and the same people.

Second, the proportion mediated is a ratio, and a ratio is unstable when its denominator is small or uncertain. In the simulated cohort it is estimated at 0.94, with a delta-method interval of 0.66 to 1.22 that covers the simulated truth of 0.83 and runs past 1. A share above 1 would mean the direct effect points the other way, so the interval spans quite different stories.

Third, the proportion mediated does not say what would happen if the mediator were removed. That question fixes weight change by intervention, and the controlled direct effect answers it. With weight change held at 0 kg for everyone, the extra exercise would lower SBP by 2.50 mmHg in the simulated population, against a total effect of 5.95 mmHg.

A proportion mediated is best reported, if at all, beside the natural direct, natural indirect and total effects with their intervals, with the split it uses named.

What to do in your own analysis

• Write the contrast ($a^*$ and $a$), the effect (controlled or natural) and the scale into the analysis plan before fitting any model.
• Keep the exposure-mediator product term in the outcome model unless there is a stated prior reason to drop it, and report its estimate.
• Report the natural direct, natural indirect and total effects together, each with its interval and the method behind the interval.
• For a common binary outcome, consider the risk-difference scale with the mediational g-formula rather than the rare-outcome odds-ratio formulas.
• List the mediator-outcome confounders you adjusted for, name any that the exposure affects, and plan a sensitivity analysis for those you could not measure.
• Consider reporting the controlled direct effect as well when the mediator is something that could be set by intervention, such as a drug dose.

Related in the wiki: [[causal-interaction-in-clinical-epidemiology-concepts-measurement-and-interpretation]]