Mediation analysis for biostat public health data in R different methods but most efficient code possible and tutorial a

Mediation analysis in public health examines whether an intermediate variable (mediator, $M$) explains the pathway between an exposure/intervention ($X$) and a health outcome ($Y$), while controlling for confounders ($C$).

Below is a complete, copy-paste-ready guide using simulated public health data across the two most established, efficient packages in R: mediation (Causal Mediation / Potential Outcomes framework) and lavaan (Structural Equation Modeling).

Step 0: Generate Simulated Public Health Data

We will simulate a cross-sectional study:

• Exposure ($X$): Weekly physical activity (hours/week).
• Mediator ($M$): Body Mass Index ($\text{kg/m}^2$).
• Outcome ($Y$): Systolic Blood Pressure (mmHg).
• Confounders ($C$): Age (years) and Sex ($0 = \text{Female}, 1 = \text{Male}$).

# Install required packages if needed
# install.packages(c("mediation", "lavaan", "tidyverse"))

library(tidyverse)
library(mediation)
library(lavaan)

set.seed(42)
n <- 500

# Confounders
age <- rnorm(n, mean = 50, sd = 10)
sex <- rbinom(n, size = 1, prob = 0.5)

# Exposure (Physical Activity in hours/week)
exercise <- pmax(0, rnorm(n, mean = 3.5, sd = 2) - 0.02 * age)

# Mediator (BMI): negatively impacted by exercise, positively by age
bmi <- 28 - (0.6 * exercise) + (0.08 * age) + (0.5 * sex) + rnorm(n, 0, 2)

# Outcome (Systolic BP): reduced by exercise directly and indirectly via lower BMI
sbp <- 110 + (1.2 * bmi) - (1.5 * exercise) + (0.3 * age) + (2 * sex) + rnorm(n, 0, 8)

ph_data <- data.frame(exercise, bmi, sbp, age, sex)
head(ph_data)

Method 1: The Causal Mediation Framework (mediation package)

Best for: Epidemiological studies, non-linear models (e.g., logistic regression for binary disease outcomes), and testing exposure-mediator interactions or sensitivity to unmeasured confounding.

Code:

# 1. Fit the Mediator Model (M ~ X + Confounders)
model_m <- lm(bmi ~ exercise + age + sex, data = ph_data)

# 2. Fit the Outcome Model (Y ~ X + M + Confounders)
model_y <- lm(sbp ~ exercise + bmi + age + sex, data = ph_data)

# 3. Estimate mediation effects with Quasi-Bayesian Monte Carlo or Nonparametric Bootstrap
set.seed(42)
causal_med <- mediate(
  model.m = model_m,
  model.y = model_y,
  treat = "exercise",
  mediator = "bmi",
  boot = TRUE,          # Set to TRUE for nonparametric bootstrap CIs (standard in biostat)
  sims = 1000           # 1000-5000 iterations recommended for publications
)

# 4. View results
summary(causal_med)

How to Read the Output:

• ACME (Average Causal Mediation Effect / Indirect Effect): The change in SBP attributed solely to changes in BMI caused by exercise.
• ADE (Average Direct Effect): The change in SBP caused by exercise through all pathways other than BMI.
• Total Effect: Total change in SBP per unit increase in exercise ($\text{ACME} + \text{ADE}$).
• Prop. Mediated: Percentage of the total effect explained by BMI ($\frac{\text{ACME}}{\text{Total Effect}}$).

Essential Public Health Step: Sensitivity Analysis

Epidemiological mediation relies on the assumption of no unmeasured confounding. The medsens() function evaluates how robust your findings are to violations of this assumption:

sens <- medsens(causal_med, rho.by = 0.1, effect.type = "indirect")
summary(sens)
plot(sens) # Shows the correlation (rho) needed between residuals to nullify the effect

Method 2: Structural Equation Modeling (lavaan package)

Best for: Computation speed, multiple mediators simultaneously, latent variables, and continuous outcomes.

Code:

# 1. Define the model using path labels (a, b, c_prime)
sem_model <- '
  # Mediator path (a)
  bmi ~ a*exercise + age + sex

  # Outcome path (b and c_prime)
  sbp ~ b*bmi + c_prime*exercise + age + sex

  # Derived parameters
  indirect   := a * b
  direct     := c_prime
  total      := c_prime + (a * b)
  prop_med   := indirect / total
'

# 2. Fit SEM with robust standard errors or bootstrapping
fit_sem <- sem(
  model = sem_model,
  data = ph_data,
  se = "bootstrap",     # Bootstrap handles non-normal indirect effect distributions
  bootstrap = 1000
)

# 3. View clean parameter estimates and 95% Confidence Intervals
parameterEstimates(fit_sem, ci = TRUE, boot.ci.type = "perc") %>%
  filter(op == ":=") %>%
  select(label = lhs, est, se, pvalue, ci.lower, ci.upper)

Method Comparison Summary

Key Epidemiological Reporting Checklist

When reporting these results in a manuscript:

1. Report Effect Sizes & 95% Bootstrap CIs: Present the indirect (ACME), direct (ADE), total effect, and proportion mediated.
2. State Confounder Adjustments Clearly: Address baseline confounding of the $X \to Y$, $X \to M$, and $M \to Y$ relationships.
3. Verify Temporal Precedence: Ensure $X$ precedes $M$, and $M$ precedes $Y$ in study timing to avoid reverse causality.