Mediation analysis for biostat in R and explain how it works and give a tutorial on example data set with a mix of numer

Mediation Analysis in R: A Complete Tutorial

What Mediation Analysis Actually Does

Mediation analysis answers the question: "How does X affect Y?" — specifically, whether the effect of a predictor (X) on an outcome (Y) runs through an intermediate variable (M, the mediator).

The classic decomposition splits the total effect of X on Y into:

• Direct effect (ADE): X → Y, holding the mediator fixed
• Indirect effect (ACME): X → M → Y, the part that flows through the mediator
• Total effect = Direct + Indirect

        M (mediator)
       ↗   ↘
      a     b
     ↗       ↘
X ─────c'──────→ Y
       (direct)

Indirect (ACME) = a × b
Total effect    = c' + (a × b)

Key idea: We fit two models — one for the mediator, one for the outcome — then use them to simulate what the outcome would be under different hypothetical exposure values.

The Two-Model Framework (Baron–Kenny → modern approach)

The old Baron & Kenny (1986) steps are still taught, but the modern standard is the potential-outcomes / counterfactual framework implemented in the mediation package (Imai, Keele, Tingley). It:

• Works with any model type (linear, logistic, etc.)
• Handles continuous and categorical exposures/mediators
• Uses bootstrapping or quasi-Bayesian simulation for confidence intervals
• Gives valid inference even with interactions and nonlinearities

Setup

install.packages("mediation")
library(mediation)
set.seed(2024)

Building an Example Dataset (mixed numeric + categorical)

Let's simulate a realistic biostatistics scenario:

Research question: Does a health program (categorical treatment) reduce blood pressure (numeric outcome), and is this effect mediated by BMI (numeric mediator)? We'll also include covariates: age (numeric) and sex (categorical).

n <- 500

# Covariates
age <- round(rnorm(n, 50, 10))
sex <- factor(sample(c("Female", "Male"), n, replace = TRUE))

# Treatment: categorical (0 = control, 1 = program), depends partly on age
program <- rbinom(n, 1, plogis(-2 + 0.03 * age))
program <- factor(program, labels = c("Control", "Program"))

# Mediator: BMI — reduced by the program
# (Program lowers BMI; age and sex also matter)
bmi <- 28 - 2.5 * (program == "Program") +
       0.05 * age +
       1.2 * (sex == "Male") +
       rnorm(n, 0, 2)

# Outcome: systolic blood pressure
# Program has a small DIRECT effect + a large effect THROUGH bmi
sbp <- 90 +
       1.5 * bmi +                       # mediator → outcome (path b)
       -3 * (program == "Program") +     # direct effect (path c')
       0.4 * age +
       2 * (sex == "Male") +
       rnorm(n, 0, 5)

dat <- data.frame(program, bmi, sbp, age, sex)
head(dat)

Here the true structure is:

• Program lowers BMI (path a)
• BMI raises SBP (path b) → so program indirectly lowers SBP
• Program also directly lowers SBP by ~3 units (path c')

Step 1: Fit the Mediator Model

Model the mediator (BMI) as a function of the treatment + covariates.

model.M <- lm(bmi ~ program + age + sex, data = dat)
summary(model.M)

Look at the coefficient for programProgram — this is path a (effect of treatment on the mediator). You should see it's near −2.5.

Step 2: Fit the Outcome Model

Model the outcome (SBP) as a function of the treatment and the mediator + covariates.

model.Y <- lm(sbp ~ program + bmi + age + sex, data = dat)
summary(model.Y)

• The bmi coefficient is path b (~1.5).
• The programProgram coefficient here is the direct effect c' (~−3).

Step 3: Run the Mediation Analysis

mediate() combines both models and simulates counterfactuals. You must tell it the treatment and mediator variable names.

med.out <- mediate(model.M, model.Y,
                   treat = "program",
                   mediator = "bmi",
                   boot = TRUE,         # bootstrap CIs
                   sims = 1000)         # number of simulations

summary(med.out)

Interpreting the Output

You'll get a table like this (values approximate):

                Estimate  95% CI Lower  95% CI Upper  p-value
ACME              -3.8        -4.6         -3.0        <2e-16   ← indirect (a*b)
ADE               -3.1        -4.1         -2.1        <2e-16   ← direct (c')
Total Effect      -6.9        -7.9         -5.9        <2e-16   ← total
Prop. Mediated     0.55        0.44         0.66       <2e-16   ← % via BMI

How to read it:

Both CIs exclude 0 → significant mediation. The mediator explains about half the total effect.

plot(med.out)   # visualize ACME, ADE, Total with CIs

Handling a Treatment × Mediator Interaction

If the mediator's effect differs by treatment group, include an interaction and let mediate() account for it:

model.Y2 <- lm(sbp ~ program * bmi + age + sex, data = dat)

med.out2 <- mediate(model.M, model.Y2,
                    treat = "program",
                    mediator = "bmi",
                    boot = TRUE, sims = 1000)
summary(med.out2)

With interactions, mediate() reports ACME/ADE for control and treated groups separately (they can differ).

Categorical Outcomes or Mediators

The framework generalizes cleanly:

# Binary outcome (e.g., hypertension yes/no) → use glm
model.Y.bin <- glm(hypertension ~ program + bmi + age + sex,
                   data = dat, family = binomial)

med.bin <- mediate(model.M, model.Y.bin,
                   treat = "program", mediator = "bmi",
                   sims = 1000)

mediate() automatically handles the nonlinear link — effects are reported on the appropriate scale.

Sensitivity Analysis (crucial!)

Mediation assumes no unmeasured confounding between mediator and outcome. This is untestable, so check robustness:

sens <- medsens(med.out, rho.by = 0.1)
summary(sens)
plot(sens)

The parameter ρ (rho) represents correlation between the error terms of the two models (i.e., hidden confounding). The output tells you how strong a confounder would need to be to nullify your ACME. If ACME only disappears at implausibly large ρ, your conclusion is robust.

Key Assumptions to Report

1. No unmeasured X→M confounding
2. No unmeasured M→Y confounding (check with medsens)
3. No unmeasured X→Y confounding
4. No M→Y confounder affected by X (the trickiest one)
5. Correct model specification (linearity, link functions)

Randomizing X (as in a trial) satisfies #1 and #3, but #2 and #4 still require careful covariate adjustment or sensitivity analysis.

Quick Reference Workflow

# 1. Mediator model
model.M <- lm(M ~ X + covariates, data = d)

# 2. Outcome model
model.Y <- lm(Y ~ X + M + covariates, data = d)

# 3. Mediate
out <- mediate(model.M, model.Y, treat = "X", mediator = "M",
               boot = TRUE, sims = 1000)
summary(out)

# 4. Sensitivity
summary(medsens(out))

A note on honesty: The numeric outputs above are approximate expectations from the simulation design, not values I actually ran. When you run the code with set.seed(2024), your estimates will be close but not identical to the tables shown — always report your actual output.

Would you like me to extend this to a categorical mediator (e.g., a 3-level mediator), a multiple-mediator setup, or show the older Baron–Kenny hand-calculation for comparison?