23  Mediation Analysis: Opening the Causal Black Box

A trial shows your new drug lowers mortality. Good news — but a colleague asks a deeper question: does the drug save lives directly, or only because it lowers blood pressure? The distinction is not academic. If the benefit runs entirely through blood pressure, then any cheaper drug that lowers blood pressure as much should help equally, and a patient whose blood pressure is already controlled may gain nothing. If instead the drug has a direct protective effect of its own, the mechanism is different and the clinical story changes. Mediation analysis is how we prise open that black box — splitting a treatment’s total effect into the part that travels through an intermediate variable and the part that does not. In the previous chapters we learned to estimate the total effect of a treatment while fighting off confounding. Here we go one level deeper and ask how that effect comes about.

23.1 Introduction

A mediator is a variable that sits on the causal path between a treatment and an outcome: the treatment changes the mediator, and the mediator in turn changes the outcome. Drug lowers blood pressure; lower blood pressure lowers mortality. Mediation analysis decomposes the treatment’s overall effect into:

• the indirect (mediated) effect — the part that flows through the mediator, and
• the direct effect — everything else, the part that does not go through that mediator.

Their sum is the total effect we already know how to estimate:

\[ \text{Total effect} = \text{Direct effect} + \text{Indirect effect} \]

Why does a clinician care? Three reasons. First, mechanism: understanding how a treatment works guides drug development and tells us which patients will benefit. Second, targeting: if the effect is fully mediated by blood pressure, monitoring and managing blood pressure becomes the lever. Third, honesty about surrogates: regulators increasingly ask whether an intermediate marker truly carries a treatment’s benefit, which is precisely a mediation question.

This chapter introduces the modern, causal way to do this — distinguishing it from an older statistical recipe that, while historically dominant, can give the wrong answer.

23.2 The Mediation Triangle

Every mediation problem starts with the same picture: an exposure (treatment) that affects an outcome both directly and indirectly through a mediator.

graph LR
    A["Exposure<br/>(drug)"] --> M["Mediator<br/>(blood pressure)"]
    M --> Y["Outcome<br/>(mortality)"]
    A --> Y
    style A fill:#eef3fb,color:#111
    style M fill:#eef3fb,color:#111
    style Y fill:#eef3fb,color:#111

Figure 23.1: The mediation triangle. The exposure affects the outcome through two routes: indirectly via the mediator (exposure -> mediator -> outcome) and directly (exposure -> outcome). Mediation analysis estimates the size of each route.

The arrow exposure -> mediator -> outcome is the indirect path. The arrow exposure -> outcome is the direct path. The art of mediation analysis is measuring the size of each without confusing them — and, crucially, without being fooled by confounding that can lurk on any of the three arrows.

23.3 From Baron and Kenny to the Causal Framework

23.3.1 The historical approach

For decades, mediation was done with the Baron and Kenny procedure (Baron and Kenny 1986): fit a regression of the outcome on the exposure (the total effect), then add the mediator; the amount by which the exposure’s coefficient shrinks was called the “mediated” effect, and the product of two regression coefficients estimated the indirect effect.

This is intuitive and still widely taught, but it has serious limitations that matter clinically, as reviewed for epidemiologists by Richiardi et al. (2013):

1. It breaks down when exposure and mediator interact. If the drug’s direct protection is stronger in patients whose blood pressure dropped a lot, the simple product-of-coefficients no longer corresponds to any well-defined effect.
2. It does not work cleanly for binary or time-to-event outcomes — exactly the outcomes (death, MACE, readmission) clinicians care about most — because effects on the odds-ratio or hazard scale do not add up the way the method assumes.
3. It hides its causal assumptions, so users rarely check whether they hold.

23.3.2 The modern causal framework

VanderWeele and others rebuilt mediation on the potential-outcomes foundation from the earlier chapters, giving the effects a precise causal meaning that survives interactions and non-linear outcomes (VanderWeele 2016). Two definitions do the work.

Recall that a potential outcome is what would happen under a specified scenario. In mediation we need nested potential outcomes — what the outcome would be if we set the exposure to one value and set the mediator to the value it would have taken under some (possibly different) exposure.

• The natural direct effect (NDE) is the change in outcome when we switch the exposure on but hold the mediator at the value it would have had with the exposure off. In words: the drug’s effect if it could not act through blood pressure at all — blood pressure frozen where it would have sat without the drug.
• The natural indirect effect (NIE) is the change in outcome from letting the mediator move as it naturally would when the exposure switches on, while holding the exposure itself fixed. In words: the effect of just the blood-pressure change the drug causes, with the drug’s other actions held constant.

Their sum is the total effect: \(\text{Total} = \text{NDE} + \text{NIE}\). There is also the controlled direct effect (CDE) — the exposure’s effect when we fix the mediator at one specific value for everyone (e.g. “what if we forced every patient’s blood pressure to 120?”). The CDE answers an intervention question (“what if we clamped the mediator?”), whereas the natural effects answer a decomposition question (“how much of the effect runs through the mediator naturally?”).

A widely reported summary is the proportion mediated — the indirect effect as a fraction of the total — which tells a clinician, in one number, how much of the benefit travels through the mediator.

NoteNatural vs controlled, in one breath

• Natural direct/indirect effects: let the mediator behave as it naturally would. Used to decompose a total effect into mechanism. Sum to the total effect.
• Controlled direct effect: pin the mediator to a chosen value. Answers “what if we intervened on the mediator?” Does not, in general, add up with an indirect effect.

23.3.3 The assumptions: the part everyone underestimates

Mediation asks more of the data than estimating a total effect, because we are now estimating effects on two arrows (exposure -> mediator and mediator -> outcome) and the direct arrow. Identifying natural direct and indirect effects requires no unmeasured confounding of four relationships:

1. exposure -> outcome,
2. exposure -> mediator,
3. mediator -> outcome, and
4. critically, no mediator-outcome confounder that is itself affected by the exposure.

graph LR
    A["Exposure"] --> M["Mediator"]
    M --> Y["Outcome"]
    A --> Y
    C1["Exposure-mediator<br/>confounder"] --> A
    C1 --> M
    C2["Mediator-outcome<br/>confounder"] --> M
    C2 --> Y
    style A fill:#eef3fb,color:#111
    style M fill:#eef3fb,color:#111
    style Y fill:#eef3fb,color:#111
    style C1 fill:#fbf3e8,color:#111
    style C2 fill:#fbf3e8,color:#111

Figure 23.2: Mediation needs more than a clean total effect. Confounders may sit on the exposure-mediator path (C1) and on the mediator-outcome path (C2). The mediator-outcome confounders especially must be measured and adjusted, or the indirect effect is biased.

A subtlety the older method ignored: the exposure-mediator interaction. If the drug’s direct effect depends on how much blood pressure moved, a correct analysis must include an interaction term between exposure and mediator in the outcome model. The modern regression-based estimators do this automatically and report effects that remain meaningful in its presence. Because these no-unmeasured-confounding assumptions are strong and untestable, a sensitivity analysis — such as the E-value extended to mediation — is not optional; it is part of the analysis.

23.4 A Worked Example

We simulate the clinical story directly: a drug that lowers a continuous mediator (think blood pressure), which in turn affects a continuous outcome, plus a direct effect of the drug on the outcome. Because we build the decomposition, we know the right answer and can check the methods recover it.

R

ImportantInstall this one first: regmedint is on CRAN, CMAverse is not

The two leading R packages for causal mediation install in completely different ways, and this trips up nearly everyone.

regmedint is on CRAN, so the ordinary route works:

install.packages("regmedint")

CMAverse has never been on CRAN. It lives on GitHub only, so install.packages("CMAverse") fails — and it fails with a message that misleadingly suggests your R version is the problem:

> install.packages("CMAverse")
Warning message:
package 'CMAverse' is not available for this version of R

Your R version is fine. The package simply is not in the repository. Install it from GitHub instead:

install.packages("remotes")                      # once
remotes::install_github("BS1125/CMAverse")

If you use this book’s pixi environment, pixi run -e full post_install already does both for you.

Because of this, the worked example below uses regmedint as its primary tool: it is on CRAN, it implements the same regression-based estimator with exposure-mediator interaction, and it covers linear, logistic, Poisson, Cox, and accelerated-failure-time outcomes. CMAverse follows as the second example, for its built-in sensitivity analysis.

Option 1: regmedint (on CRAN)

# install.packages("regmedint")
library(regmedint)   # regression-based causal mediation with EM interaction

# --- Simulate exposure -> mediator -> outcome with a KNOWN decomposition ---
set.seed(2024)
n <- 2000
conf <- rnorm(n)                              # a measured confounder
exposure <- rbinom(n, 1, plogis(0.2 * conf))  # drug, depends on confounder

# Mediator (e.g. change in blood pressure): the drug moves it by 2.0
mediator <- 2.0 * exposure + 0.5 * conf + rnorm(n)

# Outcome: a DIRECT drug effect (0.8) PLUS an effect via the mediator (0.6),
# so the indirect effect is 2.0 * 0.6 = 1.2 and the total effect is 2.0
outcome <- 0.8 * exposure + 0.6 * mediator + 0.4 * conf + rnorm(n)

dat <- data.frame(exposure, mediator, outcome, conf)

fit <- regmedint(
  data = dat,
  yvar = "outcome",                 # outcome
  avar = "exposure",                # exposure
  mvar = "mediator",                # mediator
  cvar = "conf",                    # baseline confounders
  a0 = 0, a1 = 1,                   # contrast: exposure 1 vs 0
  m_cde = 0,                        # value to fix the mediator at, for the CDE
  c_cond = mean(dat$conf),          # covariate values to report effects at
  mreg = "linear",                  # mediator model
  yreg = "linear",                  # outcome model
  interaction = TRUE,               # allow exposure-mediator interaction
  casecontrol = FALSE
)

summary(fit)

What the code shows. We build a cohort where the drug has two routes to the outcome: a direct effect (coefficient 0.8) and an indirect one through the mediator (the drug moves the mediator by 2.0, and each mediator unit moves the outcome by 0.6, so the indirect effect is \(2.0 \times 0.6 = 1.2\) and the total is 2.0). regmedint() fits both the mediator model and the outcome model — the latter including the exposure-mediator interaction, because interaction = TRUE — and then reports the decomposition with delta-method standard errors.

Running it recovers what we built in:

Row in the output	Meaning	Estimate	What we built
pnde	Pure natural direct effect — the drug’s effect with the mediator held where it would have sat without the drug	0.72 (95% CI 0.58 to 0.87)	0.8
tnie	Total natural indirect effect — the effect of just the mediator change the drug causes	1.29 (1.16 to 1.42)	1.2
te	Total effect — and note that it equals pnde + tnie	2.01 (1.91 to 2.12)	2.0
pm	Proportion mediated	0.64 (0.57 to 0.71)	0.60
cde	Controlled direct effect — the drug’s effect if we forced every patient’s mediator to m_cde = 0	0.72 (0.58 to 0.87)	—

Read clinically: about 64% of the drug’s benefit flows through its effect on blood pressure; the remaining third is a direct effect by some other mechanism. The confidence intervals contain the values we built in, which is the check that matters.

NoteWhy the output has four direct/indirect effects, not two

regmedint prints pnde, tnde, pnie, and tnie. The prefixes are “pure” and “total”, and the distinction only exists when exposure and mediator interact.

• pnde and tnie decompose the total effect one way: \(\text{TE} = \text{pnde} + \text{tnie}\). Here the direct effect is evaluated with the mediator behaving as it would without exposure, and the interaction is bundled into the indirect part.
• tnde and pnie decompose it the other way: \(\text{TE} = \text{tnde} + \text{pnie}\), with the interaction bundled into the direct part.

Both are valid; they just assign the interaction to different sides of the ledger. When there is no interaction, pnde = tnde and pnie = tnie, and the distinction vanishes. Report which decomposition you used, because the proportion mediated differs between them. When in doubt, the pnde / tnie pair is the conventional default and the one most papers mean by “natural direct and indirect effect”.

Option 2: CMAverse (GitHub only)

CMAverse implements several causal mediation estimators behind one interface and — its main advantage — includes the E-value sensitivity analysis built in.

# NOT on CRAN. Install with:
#   remotes::install_github("BS1125/CMAverse")
library(CMAverse)   # unified causal mediation analysis + E-value

# --- Same simulated cohort as above ---------------------------------------
set.seed(2024)
n <- 2000
conf <- rnorm(n)
exposure <- rbinom(n, 1, plogis(0.2 * conf))
mediator <- 2.0 * exposure + 0.5 * conf + rnorm(n)
outcome <- 0.8 * exposure + 0.6 * mediator + 0.4 * conf + rnorm(n)
dat <- data.frame(exposure, mediator, outcome, conf)

# --- Estimate natural direct/indirect effects, regression-based -----------
res <- cmest(
  data = dat, model = "rb",            # rb = regression-based
  outcome = "outcome", exposure = "exposure",
  mediator = "mediator", basec = "conf",   # baseline confounders
  EMint = TRUE,                        # allow exposure-mediator interaction
  mreg = list("linear"), yreg = "linear",
  astar = 0, a = 1, mval = list(0),    # contrast: exposure 1 vs 0
  estimation = "imputation", inference = "bootstrap", nboot = 500
)
summary(res)

# --- Sensitivity analysis for unmeasured mediator-outcome confounding ----
cmsens(object = res, sens = "uc")

What the code shows. cmest() with model = "rb" fits the same regression-based estimator as regmedint, so on this data the two packages agree to three decimal places: pnde 0.723 (regmedint: 0.724), tnie 1.290 (1.290), te 2.012 (2.014), proportion mediated 0.641 (0.641). That is a reassuring cross-check and worth doing once in real work — if two independent implementations of the same estimator disagree, you have specified something differently in one of them.

cmest() also prints a fuller decomposition than regmedint: alongside the natural effects it reports intref and intmed, which split out the part of the effect attributable to the exposure-mediator interaction specifically. Here they are near zero with intervals spanning zero, which is correct — we simulated no interaction.

The cmsens() call is the reason to bother with the GitHub install. It performs the E-value sensitivity analysis for each effect in the decomposition, on the risk-ratio scale. For the indirect effect it returns an E-value of about 3.4 (and 3.1 for the confidence limit nearest the null), meaning an unmeasured mediator-outcome confounder would have to be associated with both the mediator and the outcome by a risk ratio of at least 3.4 to explain away the mechanism we found. That is a strong confounder — so this particular result is reasonably robust.

The classic mediation package (Imai and colleagues) is also on CRAN and still widely used and well documented, though regmedint and CMAverse are the current first choices when you need causal effect definitions that survive an exposure-mediator interaction.

Python

Python’s causal-mediation ecosystem is genuinely thinner than R’s — there is no single package with the breadth of CMAverse. The honest options are to compute the regression-based decomposition manually with statsmodels (entirely feasible for the common linear case shown here) or to use a general causal library. Below we do it by hand so the mechanics are transparent.

import numpy as np
import pandas as pd
import statsmodels.formula.api as smf

# --- Simulate the SAME exposure -> mediator -> outcome structure ---
rng = np.random.default_rng(2024)
n = 2000
conf     = rng.normal(0, 1, n)
exposure = rng.binomial(1, 1/(1+np.exp(-0.2*conf)))
mediator = 2.0*exposure + 0.5*conf + rng.normal(0, 1, n)
outcome  = 0.8*exposure + 0.6*mediator + 0.4*conf + rng.normal(0, 1, n)
df = pd.DataFrame(dict(exposure=exposure, mediator=mediator,
                       outcome=outcome, conf=conf))

# --- Regression-based decomposition (linear case, no EM interaction) ---
# Mediator model: how much does exposure move the mediator?
m_model = smf.ols("mediator ~ exposure + conf", data=df).fit()
beta_exp_on_med = m_model.params["exposure"]          # ~ 2.0

# Outcome model: direct effect of exposure, and effect of mediator
y_model = smf.ols("outcome ~ exposure + mediator + conf", data=df).fit()
direct  = y_model.params["exposure"]                  # NDE ~ 0.8
med_on_y = y_model.params["mediator"]                 # ~ 0.6

indirect = beta_exp_on_med * med_on_y                 # NIE ~ 1.2
total    = direct + indirect                          # ~ 2.0
prop_med = indirect / total                           # ~ 0.6

print(f"Natural direct effect (NDE):   {direct:.3f}")
print(f"Natural indirect effect (NIE): {indirect:.3f}")
print(f"Total effect:                  {total:.3f}")
print(f"Proportion mediated:           {prop_med:.1%}")

# --- Bootstrap confidence interval for the indirect effect ---
def nie(data):
    bm = smf.ols("mediator ~ exposure + conf", data=data).fit().params["exposure"]
    by = smf.ols("outcome ~ exposure + mediator + conf", data=data).fit().params["mediator"]
    return bm * by

boot = np.array([nie(df.sample(len(df), replace=True, random_state=i))
                 for i in range(500)])
lo, hi = np.percentile(boot, [2.5, 97.5])
print(f"NIE 95% bootstrap CI: ({lo:.3f}, {hi:.3f})")

What the code shows. We rebuild the identical data and decompose the effect by hand. The mediator model estimates how strongly the drug moves blood pressure (the exposure coefficient, near 2.0). The outcome model gives two numbers: the drug’s direct effect holding the mediator fixed (near 0.8, the natural direct effect) and how much each unit of mediator moves the outcome (near 0.6). Multiplying the first by the last gives the indirect effect (near 1.2, the natural indirect effect); adding the two gives the total (near 2.0); and their ratio is the proportion mediated (about 60%). We then bootstrap the indirect effect to get an honest confidence interval. The clinical reading matches the R output here — because we simulated no exposure–mediator interaction, the simple product-of-coefficients and the CMAverse estimate agree; with a real interaction or a binary/survival outcome they would not, which is the limitation flagged below. Either way: most of the benefit runs through blood pressure, with a meaningful direct component remaining.

WarningThe Python product-of-coefficients shortcut has limits

The manual calculation above is the classic product-of-coefficients estimator and is correct for a linear outcome with no exposure-mediator interaction. The moment your outcome is binary or time-to-event, or you suspect the drug’s direct effect depends on the mediator, this shortcut breaks — and Python lacks a polished package to rescue you. For those cases, the practical recommendation is to do the mediation analysis in R with CMAverse or regmedint. We say this plainly rather than pretend the Python tooling is equivalent.

23.5 Sensitivity Analysis: the E-value for Mediation

Because mediation leans on the strong no-unmeasured-confounding assumptions listed above — especially no unmeasured mediator-outcome confounding, which is often the weakest link — every mediation result should be accompanied by a sensitivity analysis. The E-value, which we met for total effects in the causal-inference chapter, has been extended by VanderWeele and colleagues to direct and indirect effects. It answers: how strongly would an unmeasured confounder have to be associated with both the mediator and the outcome to fully explain away the indirect (or direct) effect we reported? A large E-value means only an implausibly strong hidden confounder could erase the mechanism you found; an E-value near 1 means a mild, common confounder could. In CMAverse this comes free via cmsens(), as shown above (remembering that CMAverse installs from GitHub, not CRAN); reporting it turns “the drug works through blood pressure” into a claim a sceptical reviewer can weigh.

If you would rather stay on CRAN, the EValue package computes the same quantity from a reported effect estimate and its confidence interval — you supply the indirect effect on a ratio scale and read off the E-value, exactly as in Section 22.8. It is less convenient than cmsens(), which knows about your fitted model, but it does the job.

23.6 Exercises

TipExercise 1: Decompose a known effect (R or Python)

Simulate exposure -> mediator -> outcome with a continuous mediator, a direct effect of 1.0, a mediator effect of 0.5, and an exposure effect on the mediator of 1.5.

1. What are the true NDE, NIE, total effect, and proportion mediated?
2. Estimate them and confirm your method recovers the truth.
3. Interpret the proportion mediated for a clinician in one sentence.

NoteSolution

R

# =============================================================================
# Chapter 17c, Exercise 1: Decompose a Known Mediation Effect
# Simulate exposure -> mediator -> outcome and recover NDE, NIE, total, prop. med.
# =============================================================================
# We first implement the regression-based / product-of-coefficients estimator by
# hand with base R lm(), so that every step is visible. For a continuous mediator
# and continuous outcome with NO exposure-mediator interaction this is exactly
# what the packages compute. A bootstrap gives the CI for the indirect effect.
# Part (d) then cross-checks the answer against the regmedint package.
#
# Libraries -------------------------------------------------------------------
# Base R only for parts (a)-(c). Part (d) is optional and needs:
#   install.packages("regmedint")
# (Note: install.packages("CMAverse") does NOT work -- that package is on GitHub
#  only. See the chapter callout, or use remotes::install_github("BS1125/CMAverse").)

set.seed(42)

# --- Simulate exposure -> mediator -> outcome with a KNOWN decomposition ---
# Data-generating coefficients (the "truth"):
#   exposure -> mediator (a):        1.5
#   mediator -> outcome (b):         0.5
#   direct exposure -> outcome (c'): 1.0
a_true  <- 1.5   # effect of exposure on mediator
b_true  <- 0.5   # effect of mediator on outcome
cp_true <- 1.0   # direct effect of exposure on outcome

n <- 5000
exposure <- rbinom(n, 1, 0.5)                         # randomized-like exposure
mediator <- a_true * exposure + rnorm(n)              # continuous mediator
outcome  <- cp_true * exposure + b_true * mediator + rnorm(n)  # continuous outcome

dat <- data.frame(exposure, mediator, outcome)

# -----------------------------------------------------------------------------
# (a) TRUE natural direct/indirect effects (by construction)
# -----------------------------------------------------------------------------
# For a linear model with no exposure-mediator interaction:
#   NIE  = a * b       (path exposure -> mediator -> outcome)
#   NDE  = c'          (direct path)
#   Total = NDE + NIE
#   Proportion mediated = NIE / Total
nie_true   <- a_true * b_true          # 1.5 * 0.5 = 0.75
nde_true   <- cp_true                  # 1.0
total_true <- nde_true + nie_true      # 1.75
prop_true  <- nie_true / total_true    # 0.75 / 1.75 = 0.4286

# -----------------------------------------------------------------------------
# (b) ESTIMATE the decomposition from the data
# -----------------------------------------------------------------------------
# Mediator model: M ~ X  -> coefficient on exposure is the 'a' path
m_model <- lm(mediator ~ exposure, data = dat)
a_hat   <- coef(m_model)["exposure"]

# Outcome model: Y ~ X + M -> exposure coef is NDE (c'), mediator coef is 'b'
y_model <- lm(outcome ~ exposure + mediator, data = dat)
nde_hat <- coef(y_model)["exposure"]   # natural direct effect
b_hat   <- coef(y_model)["mediator"]   # mediator -> outcome

nie_hat   <- a_hat * b_hat             # natural indirect effect (product of coefs)
total_hat <- nde_hat + nie_hat
prop_hat  <- nie_hat / total_hat

# Bootstrap 95% CI for the indirect effect (NIE)
n_boot <- 1000
boot_nie <- numeric(n_boot)
for (i in seq_len(n_boot)) {
  idx <- sample(seq_len(n), n, replace = TRUE)
  d   <- dat[idx, ]
  bm  <- coef(lm(mediator ~ exposure, data = d))["exposure"]
  by  <- coef(lm(outcome ~ exposure + mediator, data = d))["mediator"]
  boot_nie[i] <- bm * by
}
nie_ci <- quantile(boot_nie, c(0.025, 0.975))

# -----------------------------------------------------------------------------
# Print true vs estimated
# -----------------------------------------------------------------------------
cat("=== Exercise 1: Mediation decomposition (true vs estimated) ===\n\n")
res <- data.frame(
  Quantity   = c("NDE (direct)", "NIE (indirect)", "Total effect", "Prop. mediated"),
  True       = c(nde_true, nie_true, total_true, prop_true),
  Estimated  = c(nde_hat,  nie_hat,  total_hat,  prop_hat)
)
res$True      <- round(res$True, 4)
res$Estimated <- round(res$Estimated, 4)
print(res, row.names = FALSE)

cat(sprintf("\nNIE 95%% bootstrap CI: (%.3f, %.3f)\n", nie_ci[1], nie_ci[2]))
cat(sprintf("Truth NIE = 0.75 lies inside CI: %s\n",
            nie_ci[1] <= 0.75 && 0.75 <= nie_ci[2]))

# -----------------------------------------------------------------------------
# (c) Clinician interpretation of the proportion mediated
# -----------------------------------------------------------------------------
cat("\n=== (c) One-sentence interpretation for a clinician ===\n")
cat(sprintf(
  "\"About %.0f%% of the treatment's total benefit travels through the mediator,\n",
  100 * prop_hat
))
cat("so a cheaper intervention that moved the mediator by the same amount would\n")
cat("capture roughly that share of the benefit -- but the majority of the effect\n")
cat("works by some other route, and would be lost.\"\n")

# -----------------------------------------------------------------------------
# (d) OPTIONAL cross-check against regmedint
# -----------------------------------------------------------------------------
# Never trust a hand-rolled estimator you have not checked against a package.
# regmedint is on CRAN, so this really does install and run.
if (requireNamespace("regmedint", quietly = TRUE)) {
  library(regmedint)
  fit <- regmedint(
    data = dat,
    yvar = "outcome", avar = "exposure", mvar = "mediator",
    cvar = NULL,
    a0 = 0, a1 = 1, m_cde = 0, c_cond = NULL,
    mreg = "linear", yreg = "linear",
    interaction = FALSE, casecontrol = FALSE
  )
  # summary_myreg holds the mediation decomposition: rows cde, pnde, tnie,
  # tnde, pnie, te, pm; columns est, se, Z, p, lower, upper.
  decomp <- summary(fit)$summary_myreg
  cat("\n=== (d) Cross-check with regmedint ===\n")
  print(round(decomp[c("pnde", "tnie", "te", "pm"), c("est", "lower", "upper")], 4))
  cat("\nCompare with the hand calculation above:\n")
  cat(sprintf("  pnde (= NDE) %.4f | tnie (= NIE) %.4f | te %.4f | pm %.4f\n",
              nde_hat, nie_hat, total_hat, prop_hat))
  cat("They agree, as they must: with no exposure-mediator interaction, the\n")
  cat("product-of-coefficients estimator IS the regression-based causal estimator.\n")
} else {
  cat("\n(d) Skipped: install.packages(\"regmedint\") to run the cross-check.\n")
}

Python

# =============================================================================
# Chapter 17c, Exercise 1: Decompose a Known Mediation Effect
# Simulate exposure -> mediator -> outcome and recover NDE, NIE, total, prop. med.
# =============================================================================
# NOTE: The chapter demonstrates this with R's CMAverse::cmest(). Here we
# implement the regression-based / product-of-coefficients estimator manually
# with statsmodels, which is exactly equivalent to CMAverse's "rb" model for a
# continuous mediator and continuous outcome WITHOUT an exposure-mediator
# interaction. A bootstrap gives the CI for the indirect effect.

import numpy as np
import pandas as pd
import statsmodels.formula.api as smf

np.random.seed(42)

# --- Simulate exposure -> mediator -> outcome with a KNOWN decomposition ---
# Data-generating coefficients (the "truth"):
#   exposure -> mediator (a):        1.5
#   mediator -> outcome (b):         0.5
#   direct exposure -> outcome (c'): 1.0
a_true, b_true, cp_true = 1.5, 0.5, 1.0

n = 5000
exposure = np.random.binomial(1, 0.5, n)                       # randomized-like
mediator = a_true * exposure + np.random.normal(0, 1, n)       # continuous mediator
outcome = cp_true * exposure + b_true * mediator + np.random.normal(0, 1, n)

df = pd.DataFrame(dict(exposure=exposure, mediator=mediator, outcome=outcome))

# -----------------------------------------------------------------------------
# (a) TRUE natural direct/indirect effects (by construction)
# -----------------------------------------------------------------------------
# Linear model, no exposure-mediator interaction:
#   NIE  = a * b, NDE = c', Total = NDE + NIE, Prop. mediated = NIE / Total
nie_true = a_true * b_true            # 1.5 * 0.5 = 0.75
nde_true = cp_true                    # 1.0
total_true = nde_true + nie_true      # 1.75
prop_true = nie_true / total_true     # 0.4286

# -----------------------------------------------------------------------------
# (b) ESTIMATE the decomposition from the data
# -----------------------------------------------------------------------------
# Mediator model M ~ X: exposure coefficient is the 'a' path
m_model = smf.ols("mediator ~ exposure", data=df).fit()
a_hat = m_model.params["exposure"]

# Outcome model Y ~ X + M: exposure coef = NDE (c'), mediator coef = 'b'
y_model = smf.ols("outcome ~ exposure + mediator", data=df).fit()
nde_hat = y_model.params["exposure"]  # natural direct effect
b_hat = y_model.params["mediator"]    # mediator -> outcome

nie_hat = a_hat * b_hat               # natural indirect effect (product of coefs)
total_hat = nde_hat + nie_hat
prop_hat = nie_hat / total_hat


# Bootstrap 95% CI for the indirect effect (NIE)
def nie_boot(data):
    bm = smf.ols("mediator ~ exposure", data=data).fit().params["exposure"]
    by = smf.ols("outcome ~ exposure + mediator", data=data).fit().params["mediator"]
    return bm * by


boot = np.array([nie_boot(df.sample(len(df), replace=True, random_state=i))
                 for i in range(1000)])
lo, hi = np.percentile(boot, [2.5, 97.5])

# -----------------------------------------------------------------------------
# Print true vs estimated
# -----------------------------------------------------------------------------
print("=== Exercise 1: Mediation decomposition (true vs estimated) ===\n")
res = pd.DataFrame({
    "Quantity": ["NDE (direct)", "NIE (indirect)", "Total effect", "Prop. mediated"],
    "True": [nde_true, nie_true, total_true, prop_true],
    "Estimated": [nde_hat, nie_hat, total_hat, prop_hat],
}).round(4)
print(res.to_string(index=False))

print(f"\nNIE 95% bootstrap CI: ({lo:.3f}, {hi:.3f})")
print(f"Truth NIE = 0.75 lies inside CI: {lo <= 0.75 <= hi}")

# -----------------------------------------------------------------------------
# (c) Clinician interpretation of the proportion mediated
# -----------------------------------------------------------------------------
# About 43% of the exposure's total effect on the outcome travels through the
# mediator, so roughly two-fifths of the benefit could in principle be captured
# by acting on the mediator alone, while the majority is a direct effect that a
# mediator-targeting intervention would miss.

TipExercise 2: Why Baron and Kenny can mislead (Conceptual)

1. Explain, in plain words, why adding the mediator to a regression and watching the exposure coefficient shrink can give the wrong answer when exposure and mediator interact.
2. For a binary outcome analysed by logistic regression, why does the simple product-of-coefficients indirect effect become unreliable?

NoteSolution

R

# =============================================================================
# Chapter 17c, Exercise 2: Why Baron and Kenny Can Mislead (Conceptual)
# Limits of the coefficient-shrinkage / product-of-coefficients approach.
# =============================================================================
# This exercise is conceptual; the answers are written as structured comments.
# A tiny simulation illustrates part (a) so the point is concrete rather than
# asserted. (The chapter shows the CMAverse equivalent that handles these cases.)

set.seed(42)

# -----------------------------------------------------------------------------
# (a) Coefficient shrinkage misleads under exposure-mediator INTERACTION
# -----------------------------------------------------------------------------
# The Baron-Kenny logic: fit Y ~ X (total effect c), then Y ~ X + M and read the
# shrunk exposure coefficient (c') as the "direct" effect, with c - c' the
# "mediated" effect. This assumes ONE number describes the direct effect for
# everyone.
#
# When exposure and mediator INTERACT, the exposure's direct effect is not a
# single number: it depends on the value of the mediator. The true outcome model
# is  Y = c'*X + b*M + d*(X*M) + e , so the effect of switching X on is c' + d*M,
# which changes from patient to patient. A model that omits the X:M term forces a
# single average coefficient, so:
#   - the reported "direct effect" is an ill-defined blend that matches no
#     specific patient, and
#   - the natural direct and indirect effects (which are properly defined even
#     WITH interaction) are not recovered by the simple shrinkage.
# The modern causal estimators keep the X:M term and define NDE/NIE as contrasts
# of nested potential outcomes, which stay meaningful under interaction.
#
# Illustration: data WITH an exposure-mediator interaction.

n <- 5000
X <- rbinom(n, 1, 0.5)
M <- 1.5 * X + rnorm(n)                       # exposure moves the mediator
# Outcome with a strong X:M interaction (d = 0.8)
Y <- 1.0 * X + 0.5 * M + 0.8 * (X * M) + rnorm(n)

# Naive Baron-Kenny: no interaction term
naive <- lm(Y ~ X + M)
cat("Naive Y ~ X + M (ignores interaction):\n")
cat(sprintf("  'direct' exposure coefficient = %.3f  (a single blended number)\n",
            coef(naive)["X"]))

# Correct model: include the interaction
correct <- lm(Y ~ X * M)
cat("\nCorrect Y ~ X * M (models interaction):\n")
cat(sprintf("  exposure main effect  = %.3f\n", coef(correct)["X"]))
cat(sprintf("  X:M interaction       = %.3f  (direct effect depends on M!)\n",
            coef(correct)["X:M"]))
cat("\n=> The naive single coefficient hides that the direct effect grows with M.\n")

# -----------------------------------------------------------------------------
# (b) Product-of-coefficients is unreliable for a BINARY outcome (logistic)
# -----------------------------------------------------------------------------
# For a linear outcome, effects add on the same (natural) scale, so
# NIE = a * b and NDE = c' decompose the total effect cleanly.
#
# For a binary outcome fitted by logistic regression, the coefficients live on
# the LOG-ODDS scale, and log-odds are NON-COLLAPSIBLE and NON-LINEAR:
#   1. Odds ratios do not add or multiply to reproduce the total-effect odds
#      ratio, so "a * b" and "c'" no longer sum to the total effect.
#   2. The exposure coefficient changes when the mediator is added even with NO
#      mediation and NO confounding, purely because of non-collapsibility of the
#      odds ratio - so the shrinkage is not a mediated effect at all.
#   3. The estimand implied by the product depends on the (arbitrary) outcome
#      scale, so the "proportion mediated" is not stable or interpretable.
# The fix is the causal, counterfactual definition of NDE/NIE (e.g. via the
# Valeri-VanderWeele regression formulas or simulation-based estimators), which
# are defined on the probability/risk scale and remain valid for logistic,
# Poisson, Cox, and AFT outcome models. This is exactly what CMAverse and
# regmedint implement.

Python

# =============================================================================
# Chapter 17c, Exercise 2: Why Baron and Kenny Can Mislead (Conceptual)
# Limits of the coefficient-shrinkage / product-of-coefficients approach.
# =============================================================================
# This exercise is conceptual; the answers are written as structured comments.
# A tiny simulation illustrates part (a) so the point is concrete rather than
# asserted. (The chapter shows the CMAverse equivalent that handles these cases.)

import numpy as np
import pandas as pd
import statsmodels.formula.api as smf

np.random.seed(42)

# -----------------------------------------------------------------------------
# (a) Coefficient shrinkage misleads under exposure-mediator INTERACTION
# -----------------------------------------------------------------------------
# The Baron-Kenny logic: fit Y ~ X (total effect c), then Y ~ X + M and read the
# shrunk exposure coefficient (c') as the "direct" effect, with c - c' the
# "mediated" effect. This assumes ONE number describes the direct effect for
# everyone.
#
# When exposure and mediator INTERACT, the exposure's direct effect is not a
# single number: it depends on the value of the mediator. The true outcome model
# is  Y = c'*X + b*M + d*(X*M) + e , so the effect of switching X on is c' + d*M,
# which changes from patient to patient. A model that omits the X:M term forces a
# single average coefficient, so:
#   - the reported "direct effect" is an ill-defined blend that matches no
#     specific patient, and
#   - the natural direct and indirect effects (which are properly defined even
#     WITH interaction) are not recovered by the simple shrinkage.
# The modern causal estimators keep the X:M term and define NDE/NIE as contrasts
# of nested potential outcomes, which stay meaningful under interaction.
#
# Illustration: data WITH an exposure-mediator interaction.

n = 5000
X = np.random.binomial(1, 0.5, n)
M = 1.5 * X + np.random.normal(0, 1, n)                 # exposure moves the mediator
# Outcome with a strong X:M interaction (d = 0.8)
Y = 1.0 * X + 0.5 * M + 0.8 * (X * M) + np.random.normal(0, 1, n)
df = pd.DataFrame(dict(X=X, M=M, Y=Y))

naive = smf.ols("Y ~ X + M", data=df).fit()             # ignores interaction
print("Naive Y ~ X + M (ignores interaction):")
print(f"  'direct' exposure coefficient = {naive.params['X']:.3f}  "
      f"(a single blended number)")

correct = smf.ols("Y ~ X * M", data=df).fit()           # models interaction
print("\nCorrect Y ~ X * M (models interaction):")
print(f"  exposure main effect  = {correct.params['X']:.3f}")
print(f"  X:M interaction       = {correct.params['X:M']:.3f}  "
      f"(direct effect depends on M!)")
print("\n=> The naive single coefficient hides that the direct effect grows with M.")

# -----------------------------------------------------------------------------
# (b) Product-of-coefficients is unreliable for a BINARY outcome (logistic)
# -----------------------------------------------------------------------------
# For a linear outcome, effects add on the same (natural) scale, so
# NIE = a * b and NDE = c' decompose the total effect cleanly.
#
# For a binary outcome fitted by logistic regression, the coefficients live on
# the LOG-ODDS scale, and log-odds are NON-COLLAPSIBLE and NON-LINEAR:
#   1. Odds ratios do not add or multiply to reproduce the total-effect odds
#      ratio, so "a * b" and "c'" no longer sum to the total effect.
#   2. The exposure coefficient changes when the mediator is added even with NO
#      mediation and NO confounding, purely because of non-collapsibility of the
#      odds ratio - so the shrinkage is not a mediated effect at all.
#   3. The estimand implied by the product depends on the (arbitrary) outcome
#      scale, so the "proportion mediated" is not stable or interpretable.
# The fix is the causal, counterfactual definition of NDE/NIE (e.g. via the
# Valeri-VanderWeele regression formulas or simulation-based estimators), which
# are defined on the probability/risk scale and remain valid for logistic,
# Poisson, Cox, and AFT outcome models. This is exactly what CMAverse and
# regmedint implement.

TipExercise 3: Assumptions and a DAG (Conceptual)

Draw the mediation triangle for “exercise -> weight loss -> blood pressure”.

1. Name one plausible exposure-mediator confounder and one plausible mediator-outcome confounder.
2. Which of the four no-unmeasured-confounding assumptions do you think is hardest to satisfy here, and why?
3. Suppose a mediator-outcome confounder is itself caused by exercise. Why does this specifically threaten the natural indirect effect?

NoteSolution

R

# =============================================================================
# Chapter 17c, Exercise 3: Assumptions and a DAG (Conceptual)
# Mediation triangle for exercise -> weight loss -> blood pressure.
# =============================================================================
# Conceptual exercise; answers as structured comments. The script prints the
# mediation triangle as text so the reasoning is anchored to a concrete DAG.

cat("Mediation triangle for the clinical question:\n\n")
cat("      Exercise (X) ------------------> Blood pressure (Y)   [direct path]\n")
cat("           |                               ^\n")
cat("           |                               |\n")
cat("           v                               |\n")
cat("      Weight loss (M) -------------------- +               [indirect path]\n\n")
cat("  Indirect: Exercise -> Weight loss -> Blood pressure\n")
cat("  Direct:   Exercise -> Blood pressure (via mechanisms other than weight)\n\n")

# -----------------------------------------------------------------------------
# (a) One exposure-mediator confounder and one mediator-outcome confounder
# -----------------------------------------------------------------------------
# Exposure(exercise)-mediator(weight loss) confounder:
#   BASELINE DIET / caloric intake. People who eat a healthier, lower-calorie
#   diet tend to both exercise more AND lose more weight, creating a spurious
#   exercise-weight-loss association not due to exercise itself.
#   (Other valid answers: baseline motivation, socioeconomic status, age.)
#
# Mediator(weight loss)-outcome(blood pressure) confounder:
#   DIETARY SODIUM INTAKE. High salt intake both impedes weight loss (fluid
#   retention, dietary pattern) AND directly raises blood pressure, confounding
#   the weight-loss -> blood-pressure arrow.
#   (Other valid answers: alcohol intake, antihypertensive medication use.)

# -----------------------------------------------------------------------------
# (b) Hardest no-unmeasured-confounding assumption, and why
# -----------------------------------------------------------------------------
# The four assumptions are: no unmeasured confounding of (1) X->Y, (2) X->M,
# (3) M->Y, and (4) no exposure-affected M->Y confounder.
#
# HARDEST here: assumption (3), no unmeasured mediator(weight loss)-outcome
# (blood pressure) confounding. Even in a randomized exercise trial, weight loss
# is NOT randomized - it arises naturally within each arm. Whatever drives how
# much weight a person loses (diet, sodium, alcohol, adherence, metabolic
# health) also tends to affect blood pressure directly. Randomizing exercise
# fixes assumptions (1) and (2) but does nothing for the M->Y arrow, so the
# indirect effect leans on adjusting for confounders we can rarely measure fully.

# -----------------------------------------------------------------------------
# (c) Why an exposure-CAUSED mediator-outcome confounder threatens the NIE
# -----------------------------------------------------------------------------
# Suppose exercise causes better sleep, and sleep both further influences weight
# loss and directly lowers blood pressure. Then sleep is a mediator-outcome
# confounder that is ITSELF caused by the exposure (assumption 4 violated).
#   - It sits on the causal path, so we cannot simply "adjust" for it: adjusting
#     for a variable on the exercise -> ... -> blood pressure path BLOCKS part of
#     exercise's real effect (over-adjustment / mediator-of-a-mediator problem).
#   - Yet NOT adjusting leaves the weight-loss -> blood-pressure arrow confounded.
# There is no way to condition our way out: the natural indirect effect is no
# longer identified by standard regression, because any choice biases it in one
# direction or the other. Handling it requires more advanced methods
# (interventional/randomized-interventional effects, or g-methods), not the
# simple product-of-coefficients.

Python

# =============================================================================
# Chapter 17c, Exercise 3: Assumptions and a DAG (Conceptual)
# Mediation triangle for exercise -> weight loss -> blood pressure.
# =============================================================================
# Conceptual exercise; answers as structured comments. The script prints the
# mediation triangle as text so the reasoning is anchored to a concrete DAG.

print("Mediation triangle for the clinical question:\n")
print("      Exercise (X) ------------------> Blood pressure (Y)   [direct path]")
print("           |                               ^")
print("           |                               |")
print("           v                               |")
print("      Weight loss (M) -------------------- +               [indirect path]\n")
print("  Indirect: Exercise -> Weight loss -> Blood pressure")
print("  Direct:   Exercise -> Blood pressure (via mechanisms other than weight)\n")

# -----------------------------------------------------------------------------
# (a) One exposure-mediator confounder and one mediator-outcome confounder
# -----------------------------------------------------------------------------
# Exposure(exercise)-mediator(weight loss) confounder:
#   BASELINE DIET / caloric intake. People who eat a healthier, lower-calorie
#   diet tend to both exercise more AND lose more weight, creating a spurious
#   exercise-weight-loss association not due to exercise itself.
#   (Other valid answers: baseline motivation, socioeconomic status, age.)
#
# Mediator(weight loss)-outcome(blood pressure) confounder:
#   DIETARY SODIUM INTAKE. High salt intake both impedes weight loss (fluid
#   retention, dietary pattern) AND directly raises blood pressure, confounding
#   the weight-loss -> blood-pressure arrow.
#   (Other valid answers: alcohol intake, antihypertensive medication use.)

# -----------------------------------------------------------------------------
# (b) Hardest no-unmeasured-confounding assumption, and why
# -----------------------------------------------------------------------------
# The four assumptions are: no unmeasured confounding of (1) X->Y, (2) X->M,
# (3) M->Y, and (4) no exposure-affected M->Y confounder.
#
# HARDEST here: assumption (3), no unmeasured mediator(weight loss)-outcome
# (blood pressure) confounding. Even in a randomized exercise trial, weight loss
# is NOT randomized - it arises naturally within each arm. Whatever drives how
# much weight a person loses (diet, sodium, alcohol, adherence, metabolic
# health) also tends to affect blood pressure directly. Randomizing exercise
# fixes assumptions (1) and (2) but does nothing for the M->Y arrow, so the
# indirect effect leans on adjusting for confounders we can rarely measure fully.

# -----------------------------------------------------------------------------
# (c) Why an exposure-CAUSED mediator-outcome confounder threatens the NIE
# -----------------------------------------------------------------------------
# Suppose exercise causes better sleep, and sleep both further influences weight
# loss and directly lowers blood pressure. Then sleep is a mediator-outcome
# confounder that is ITSELF caused by the exposure (assumption 4 violated).
#   - It sits on the causal path, so we cannot simply "adjust" for it: adjusting
#     for a variable on the exercise -> ... -> blood pressure path BLOCKS part of
#     exercise's real effect (over-adjustment / mediator-of-a-mediator problem).
#   - Yet NOT adjusting leaves the weight-loss -> blood-pressure arrow confounded.
# There is no way to condition our way out: the natural indirect effect is no
# longer identified by standard regression, because any choice biases it in one
# direction or the other. Handling it requires more advanced methods
# (interventional/randomized-interventional effects, or g-methods), not the
# simple product-of-coefficients.

TipExercise 4: Sensitivity analysis (Conceptual + code)

Take your Exercise 1 result.

1. Work out how strong an unmeasured mediator-outcome confounder would need to be to explain away the indirect effect. Compute the E-value from the closed-form formula (which is what the packages report), and if you have CMAverse installed from GitHub, cross-check it with cmsens(object, sens = "uc").
2. If the E-value for the indirect effect were 1.3, how would you temper your conclusion about mechanism?

NoteSolution

R

# =============================================================================
# Chapter 17c, Exercise 4: Sensitivity Analysis - the E-value (Conceptual + code)
# How strong must unmeasured M-Y confounding be to explain away the indirect effect?
# =============================================================================
# The chapter obtains this "free" via CMAverse::cmsens(object, sens = "uc").
# CMAverse is not installed here, so we compute the E-value with the closed-form
# formula (VanderWeele & Ding 2017), which is exactly what cmsens() reports for
# an effect expressed on the risk-ratio scale.

set.seed(42)

# -----------------------------------------------------------------------------
# E-value formula (risk-ratio scale)
# -----------------------------------------------------------------------------
# For a point estimate RR (with RR >= 1):
#     E-value = RR + sqrt(RR * (RR - 1))
# It is the minimum strength of association (on the RR scale) that an unmeasured
# confounder would need with BOTH the mediator and the outcome, above and beyond
# measured covariates, to fully explain away the observed indirect effect.
# For a protective effect (RR < 1), first transform: RR <- 1 / RR.

evalue <- function(rr) {
  if (rr < 1) rr <- 1 / rr          # put protective effects on the >=1 scale
  rr + sqrt(rr * (rr - 1))
}

# -----------------------------------------------------------------------------
# (a) Reason about the strength needed, with a worked calculation
# -----------------------------------------------------------------------------
# Suppose the indirect (mediated) effect from Exercise 1, re-expressed on a
# risk-ratio scale for a binary version of the outcome, corresponds to an
# indirect-effect RR of about 1.50. We compute its E-value:
rr_indirect <- 1.50
ev <- evalue(rr_indirect)

cat("=== Exercise 4: E-value for the indirect effect ===\n\n")
cat(sprintf("Indirect-effect risk ratio (example): RR = %.2f\n", rr_indirect))
cat(sprintf("E-value = RR + sqrt(RR*(RR-1)) = %.3f\n\n", ev))
cat(sprintf(
  "Interpretation: an unmeasured mediator-outcome confounder would need to be\n"))
cat(sprintf(
  "associated with BOTH the mediator and the outcome by a risk ratio of at\n"))
cat(sprintf(
  "least %.2f each (beyond measured covariates) to reduce the indirect effect\n", ev))
cat(sprintf(
  "to the null. Weaker confounding could shift but not erase it.\n\n"))

# For reference, the E-value at a range of indirect-effect RRs:
cat("E-value as the indirect effect grows:\n")
for (rr in c(1.1, 1.3, 1.5, 2.0, 3.0)) {
  cat(sprintf("  RR = %.1f  ->  E-value = %.2f\n", rr, evalue(rr)))
}

# -----------------------------------------------------------------------------
# (b) How to temper the conclusion if the E-value were 1.3
# -----------------------------------------------------------------------------
# An E-value of 1.3 is SMALL. It says a fairly modest unmeasured mediator-outcome
# confounder - one associated with both the mediator and the outcome by only
# about a 1.3-fold risk ratio each - would be enough to explain away the entire
# indirect effect. Confounders of that magnitude are common and plausible in
# observational health data (e.g. an unmeasured lifestyle or comorbidity factor).
# So we would report the mechanistic claim cautiously: "the data are CONSISTENT
# with partial mediation, but the indirect effect is NOT robust - a mild
# unmeasured confounder of the mediator-outcome relationship could account for
# it." We would avoid strong statements that the treatment 'works through' the
# mediator, call for measuring/adjusting more M-Y confounders, and treat the
# proportion mediated as fragile rather than established.

Python

# =============================================================================
# Chapter 17c, Exercise 4: Sensitivity Analysis - the E-value (Conceptual + code)
# How strong must unmeasured M-Y confounding be to explain away the indirect effect?
# =============================================================================
# The chapter obtains this "free" via R's CMAverse::cmsens(object, sens = "uc").
# Here we compute the E-value with the closed-form formula (VanderWeele & Ding
# 2017), which is exactly what cmsens() reports for an effect on the RR scale.

import numpy as np

np.random.seed(42)


# -----------------------------------------------------------------------------
# E-value formula (risk-ratio scale)
# -----------------------------------------------------------------------------
# For a point estimate RR (with RR >= 1):
#     E-value = RR + sqrt(RR * (RR - 1))
# It is the minimum strength of association (on the RR scale) that an unmeasured
# confounder would need with BOTH the mediator and the outcome, above and beyond
# measured covariates, to fully explain away the observed indirect effect.
# For a protective effect (RR < 1), first transform: RR <- 1 / RR.
def evalue(rr):
    if rr < 1:
        rr = 1 / rr                 # put protective effects on the >=1 scale
    return rr + np.sqrt(rr * (rr - 1))


# -----------------------------------------------------------------------------
# (a) Reason about the strength needed, with a worked calculation
# -----------------------------------------------------------------------------
# Suppose the indirect (mediated) effect from Exercise 1, re-expressed on a
# risk-ratio scale for a binary version of the outcome, corresponds to an
# indirect-effect RR of about 1.50. We compute its E-value:
rr_indirect = 1.50
ev = evalue(rr_indirect)

print("=== Exercise 4: E-value for the indirect effect ===\n")
print(f"Indirect-effect risk ratio (example): RR = {rr_indirect:.2f}")
print(f"E-value = RR + sqrt(RR*(RR-1)) = {ev:.3f}\n")
print("Interpretation: an unmeasured mediator-outcome confounder would need to be")
print("associated with BOTH the mediator and the outcome by a risk ratio of at")
print(f"least {ev:.2f} each (beyond measured covariates) to reduce the indirect")
print("effect to the null. Weaker confounding could shift but not erase it.\n")

# For reference, the E-value at a range of indirect-effect RRs:
print("E-value as the indirect effect grows:")
for rr in (1.1, 1.3, 1.5, 2.0, 3.0):
    print(f"  RR = {rr:.1f}  ->  E-value = {evalue(rr):.2f}")

# -----------------------------------------------------------------------------
# (b) How to temper the conclusion if the E-value were 1.3
# -----------------------------------------------------------------------------
# An E-value of 1.3 is SMALL. It says a fairly modest unmeasured mediator-outcome
# confounder - one associated with both the mediator and the outcome by only
# about a 1.3-fold risk ratio each - would be enough to explain away the entire
# indirect effect. Confounders of that magnitude are common and plausible in
# observational health data (e.g. an unmeasured lifestyle or comorbidity factor).
# So we would report the mechanistic claim cautiously: "the data are CONSISTENT
# with partial mediation, but the indirect effect is NOT robust - a mild
# unmeasured confounder of the mediator-outcome relationship could account for
# it." We would avoid strong statements that the treatment 'works through' the
# mediator, call for measuring/adjusting more M-Y confounders, and treat the
# proportion mediated as fragile rather than established.

23.7 Summary

Concept	Key point
Mediator	A variable on the causal path from exposure to outcome
Total = direct + indirect	A treatment’s overall effect splits into a direct part and a mediated part
Natural direct effect (NDE)	Exposure’s effect with the mediator held at its no-exposure value
Natural indirect effect (NIE)	Effect of just the exposure-driven change in the mediator
Controlled direct effect (CDE)	Exposure’s effect with the mediator fixed at a chosen value
Proportion mediated	Indirect / total; how much of the benefit runs through the mediator
Exposure-mediator interaction	Must be modelled; the modern estimators handle it, Baron-Kenny does not
Four confounding assumptions	No unmeasured confounding of A-Y, A-M, M-Y, and no exposure-affected M-Y confounder
Sensitivity analysis	E-value for mediation quantifies robustness to unmeasured confounding

TipKey Takeaways

• Mediation analysis opens the causal black box, decomposing a total effect into a direct part and an indirect (mediated) part to reveal mechanism.
• The modern causal framework (NDE, NIE, CDE) superseded the Baron and Kenny recipe because it stays valid with exposure-mediator interaction and non-linear outcomes.
• Identifying these effects requires four no-unmeasured-confounding assumptions — notably about the mediator-outcome relationship — so they are stronger than for a total effect.
• In R, regmedint is the practical first choice: it is on CRAN, handles exposure-mediator interaction, and covers linear, logistic, Poisson, Cox, and AFT outcomes. CMAverse adds a built-in sensitivity analysis but installs from GitHub, not CRAN (remotes::install_github("BS1125/CMAverse") — install.packages() will fail). The classic mediation package remains useful.
• Python’s mediation tooling is limited: the linear product-of-coefficients works, but for binary, survival, or interaction settings, prefer R — we recommend this honestly rather than overstate the ecosystem.
• Always pair a mediation estimate with a sensitivity analysis (E-value); the proportion mediated is a claim about mechanism that a sceptic must be able to weigh.

23.8 References and Further Reading

• For the definitive treatment of causal mediation, see VanderWeele (2015), with VanderWeele (2016) as the accessible overview.
• For the regression-based estimator that allows exposure–mediator interaction, see Valeri and VanderWeele (2013); for the simulation-based approach behind the mediation package, see Imai et al. (2010).
• For the historical approach and why it can mislead, see Baron and Kenny (1986) and Richiardi et al. (2013).
• For the potential-outcomes foundations underlying nested counterfactuals, see Hernán and Robins (2024).
• For software, see Shi et al. (2021) and Yoshida et al. (2024).

Baron, Reuben M, and David A Kenny. 1986. “The Moderator–Mediator Variable Distinction in Social Psychological Research: Conceptual, Strategic, and Statistical Considerations.” Journal of Personality and Social Psychology 51 (6): 1173–82. https://doi.org/10.1037/0022-3514.51.6.1173. The classic mediation procedure, still widely taught. Superseded for causal work by the counterfactual framework, since it can mislead when exposure and mediator interact or the outcome model is non-linear.

Hernán, Miguel A, and James M Robins. 2024. Causal Inference: What If. Chapman; Hall/CRC.

Imai, Kosuke, Luke Keele, and Dustin Tingley. 2010. “A General Approach to Causal Mediation Analysis.” Psychological Methods 15 (4): 309–34. https://doi.org/10.1037/a0020761. The basis of the mediation R package.

Richiardi, Lorenzo, Rino Bellocco, and Daniela Zugna. 2013. “Mediation Analysis in Epidemiology: Methods, Interpretation and Bias.” International Journal of Epidemiology 42 (5): 1511–19. https://doi.org/10.1093/ije/dyt127.

Shi, Baoyi, Christine Choirat, Brent A Coull, Tyler J VanderWeele, and Linda Valeri. 2021. “CMAverse: A Suite of Functions for Reproducible Causal Mediation Analyses.” Epidemiology 32 (5): e20–22. https://doi.org/10.1097/EDE.0000000000001378.

Valeri, Linda, and Tyler J VanderWeele. 2013. “Mediation Analysis Allowing for Exposure–Mediator Interactions and Causal Interpretation: Theoretical Assumptions and Implementation with SAS and SPSS Macros.” Psychological Methods 18 (2): 137–50. https://doi.org/10.1037/a0031034. The regression-based estimator behind regmedint and CMAverse.

VanderWeele, Tyler J. 2015. Explanation in Causal Inference: Methods for Mediation and Interaction. Oxford University Press. The definitive modern text on causal mediation.

VanderWeele, Tyler J. 2016. “Mediation Analysis: A Practitioner’s Guide.” Annual Review of Public Health 37: 17–32. https://doi.org/10.1146/annurev-publhealth-032315-021402.

Yoshida, Kazuki et al. 2024. regmedint: Regression-Based Causal Mediation Analysis with Interaction. R package. https://github.com/kaz-yos/regmedint.