Mediation Analysis with Bootstrap Confidence Intervals: A Worked Example for Behavioral Research
You hypothesize that an intervention (X) improves an outcome (Y), and you suspect a specific mechanism (M) carries the effect. The reviewer comment will be predictable: “Test the indirect effect with bootstrapped confidence intervals.” The Sobel test and Baron & Kenny’s causal-steps method are not what they want. They want a×b with an empirical CI.
This walks through a complete mediation analysis with bootstrap CIs on a worked example: stress (X) → sleep quality (M) → cognitive performance (Y). Numbers are fabricated to keep the focus on method, but the workflow and reporting are exactly what you would submit. Code is in R using the mediation package; the structure transfers to Hayes’ PROCESS macro and Python’s pingouin.mediation_analysis.
The Model and What You’re Estimating
The simple mediation model decomposes the total effect of X on Y into a direct effect (c′) and an indirect effect (a×b), where M sits in the causal chain:
$Y = c'X + bM + e_Y, \quad M = aX + e_M$The interesting quantity is the indirect effect, a×b — the part of X’s influence on Y that flows through M. The total effect c equals c′ + a×b. The Baron & Kenny “causal steps” method (35,000+ citations, originally from 1986) required four conditions and never tested a×b directly. The modern standard, from Preacher & Hayes (2004, 2008), tests a×b with bootstrapped confidence intervals — works with smaller samples, has higher power, and matches the inferential question.
Total effect: c (X → Y, ignoring M)
Direct effect: c′ (X → Y, controlling for M)
Indirect effect: a×b (X → M → Y)
Identity: c = c′ + a×b
Setup: Variables and Sample
For the worked example, assume a sample of n = 200 graduate students measured at one time point (we’ll come back to the cross-sectional problem). Variables:
- X (Stress): Perceived Stress Scale, continuous, range 0–40
- M (Sleep quality): PSQI total score, continuous, range 0–21 (lower = better)
- Y (Cognitive performance): composite z-score from three tasks
Assumed regression results from fitting the two equations:
- a-path (Stress → Sleep): a = 0.18, SE = 0.04, p < .001
- b-path (Sleep → Cognition, controlling for Stress): b = −0.22, SE = 0.06, p < .001
- c′-path (Stress → Cognition, controlling for Sleep): c′ = −0.09, SE = 0.05, p = .062
- c-path (Stress → Cognition, total): c = −0.13, SE = 0.05, p = .010
Point estimate of the indirect effect: a×b = 0.18 × (−0.22) = −0.0396. Sign means higher stress → worse sleep → worse cognition. The point estimate by itself is not enough — you need an inference about whether it differs from zero.
Why Sobel Is Wrong (and Bootstrap Is Right)
The Sobel test computes a standard error for a×b assuming the sampling distribution of the product is normal. It isn’t. The product of two normal variables has a skewed distribution, especially when one path is small. Sobel-test CIs are symmetric around the point estimate, but the true sampling distribution often is not. Result: under-coverage of the true effect, and reviewers in 2026 will flag the use.
Reporting a Sobel z-test or a Baron & Kenny “full vs. partial mediation” label as the main inference. Reviewers in psychology, behavioral medicine, and biostatistics now expect bootstrapped CIs around a×b. Even when Sobel agrees with the bootstrap, the methods section should describe the bootstrap procedure.
Bootstrapping resamples your data with replacement many times (typically 5,000), recomputes a×b in each resample, and uses the empirical distribution to build a confidence interval. The bias-corrected and accelerated (BCa) variant adjusts for skew and median bias and is the default in current mediation and boot packages.
Worked Example: R Code with mediation Package
The minimum-viable workflow assumes a data frame d with columns stress, sleep, and cognition:
library(mediation)
# Fit the two component models
model.M <- lm(sleep ~ stress, data = d)
model.Y <- lm(cognition ~ stress + sleep, data = d)
# Bootstrap the indirect effect
set.seed(2026)
med <- mediate(
model.m = model.M,
model.y = model.Y,
treat = "stress",
mediator = "sleep",
boot = TRUE,
boot.ci.type = "bca",
sims = 5000
)
summary(med)
Key arguments: boot = TRUE activates bootstrapping (default is quasi-Bayesian Monte Carlo, which is also valid but less commonly requested by reviewers); boot.ci.type = "bca" requests the bias-corrected and accelerated CI; sims = 5000 is the conventional minimum for stable CIs (some methodologists recommend 10,000 for publication).
The output reports four quantities:
- ACME (Average Causal Mediation Effect): the indirect effect a×b with its bootstrap CI
- ADE (Average Direct Effect): c′ with its CI
- Total Effect: c with its CI
- Prop. Mediated: a×b divided by c
For the worked example, the bootstrap output for ACME would be approximately: estimate = −0.040, 95% CI [−0.072, −0.014], p = .003. Because the CI excludes zero, the indirect effect is significant. The direct effect (ADE) is c′ = −0.09, 95% CI [−0.18, 0.01] — CI includes zero, so the direct effect is not detectable. The total effect c is significant.
“Stress predicted lower cognitive performance, c = −0.13, 95% CI [−0.23, −0.03]. This effect was carried by sleep quality: the indirect effect was significant, a×b = −0.040, 95% BCa bootstrap CI [−0.072, −0.014], 5,000 resamples. The direct effect of stress controlling for sleep was not detectable, c′ = −0.09, 95% CI [−0.18, 0.01].”
What “Full vs. Partial Mediation” Misses
Older papers report “full mediation” when c′ is non-significant and “partial mediation” when c′ is significant. This labeling has fallen out of favor. The distinction depends entirely on the power to detect c′ — underpowered studies will look like “full mediation” even when a meaningful direct effect exists. The current recommendation: report the magnitude of a×b and c′ with their CIs, and let the reader judge proportion-mediated. Don’t binarize.
The prop. mediated quantity is reasonable to report when it’s well-estimated, but its CI is often wide (it’s a ratio of two effects, each with sampling error). When c is near zero, prop-mediated estimates are unstable.
The Cross-Sectional Problem
Mediation models are inherently causal. They claim X causes M causes Y, in that temporal order. Cross-sectional data — everything measured at the same time — cannot establish temporal precedence. Cross-sectional mediation estimates can be substantially biased relative to longitudinal mediation parameters, even under ideal conditions (Maxwell & Cole, 2007).
What this means in practice: if your data are cross-sectional, you must explicitly disclaim causal interpretation in the methods and limitations sections. “The mediation model was estimated cross-sectionally; the indirect effect is consistent with sleep mediating stress’s effect on cognition but does not establish temporal causal ordering. A longitudinal design measuring stress before sleep before cognition would be required for causal inference.” Reviewers in clinical psychology and behavioral medicine routinely insist on this language.
If your study collects X, M, and Y at three time points (T1 → T2 → T3), use longitudinal mediation models — cross-lagged panel models or latent change score models. The lavaan package handles these. Cross-sectional models with three-wave data are wasting design strength you already paid for.
Sanity Checks Before Submitting
Independent checks that catch common errors:
- Coefficient sign of a×b matches the sign of c. If a is positive, b is negative, and c is positive, your indirect-effect sign disagrees with the total-effect sign. This signals a suppression effect (real but unusual) or a model misspecification (more common). Investigate before reporting.
- Indirect-effect magnitude is plausible. If a×b is larger in magnitude than c, the direct and indirect effects are working in opposite directions. Again: legitimate (inconsistent mediation) but rare; usually a sign that a confounder is missing.
- Bootstrap CI does not straddle the point estimate symmetrically. Asymmetric CIs around a×b are expected. If the BCa CI is symmetric, double-check that
boot.ci.typewas set; the default in some routines is “perc” (percentile), which is less precise. - Number of bootstrap resamples reported. Methods sections should state “5,000 bootstrap resamples” (or 10,000) and the CI type (BCa). Missing either is grounds for revision.
Reporting Template for the Manuscript
What the results section should contain, in order:
- Brief description of the mediation hypothesis with a slot diagram or path equation
- Total effect c with CI
- a-path coefficient with CI; b-path coefficient with CI
- Indirect effect a×b with BCa bootstrap CI, number of resamples, and the package/version used
- Direct effect c′ with CI
- Causal disclaimer if data are cross-sectional
For posts on the upstream questions — how to power a study that includes a mediator and checking assumptions before regression-based mediation models — the workflow connects directly. Both linear models in the worked example carry standard regression assumptions (linearity, normality of residuals, homoscedasticity); diagnose those on each model separately before trusting the bootstrap CI.
When to Use Something Other Than Simple Mediation
The model above (Hayes’ PROCESS Model 4) is the simplest case — one X, one M, one Y. Common extensions:
- Multiple mediators in parallel (PROCESS Model 4 with multiple Ms): bootstrap each indirect effect and the total. Useful when several mechanisms operate jointly.
- Serial mediation (Model 6): X → M1 → M2 → Y. The two-link indirect path through M1 and M2 must be bootstrapped.
- Moderated mediation (Model 7, 8, 14): the indirect effect depends on a moderator. Conditional indirect effects are reported at specific moderator values.
- Structural equation models with latent variables: switch to
lavaan. Note thatlavaan’s default for indirect-effect SEs is the delta method, not bootstrap — you must explicitly requestse = "bootstrap"with at least 1,000 resamples.
For exploratory analyses with many candidate mediators, the indirect effects are subject to the same multiple-comparisons concerns as any other family of tests. Use FDR adjustment across the indirect-effect CIs, the same way you would adjust pairwise post-hoc tests.
If you want to walk through the analysis with your own data without writing R code, GraphHelix surfaces the bootstrap configuration (number of resamples, CI type, seed) as visible parameters — the same configuration choices that belong in a methods section. A reviewer asking “how many bootstrap resamples?” gets a one-line answer instead of a re-run.