Complete Mediation Analysis Workflow

# Load foundation package
library(mediationverse)  # Loads medfit automatically

# Load packages for this workflow
library(medsim)       # For data simulation
library(probmed)      # For P_med effect size
library(RMediation)   # For confidence intervals
library(medrobust)    # For sensitivity analysis

# Set seed for reproducibility
set.seed(2025)

# Population parameters
a_true <- 0.5 # X -> M effect
b_true <- 0.4 # M -> Y effect
c_prime_true <- 0.2 # X -> Y direct effect

# True indirect effect
indirect_true <- a_true * b_true
cat("True indirect effect:", indirect_true, "\n")

# Simulate data using medsim
sim_data <- simulate_mediation(
  n = 200,
  a = a_true,
  b = b_true,
  c_prime = c_prime_true,
  sigma_m = 1,
  sigma_y = 1
)

head(sim_data)

# Mediator model: M ~ X
fit_m <- lm(M ~ X, data = sim_data)

# Outcome model: Y ~ X + M
fit_y <- lm(Y ~ X + M, data = sim_data)

# View summaries
summary(fit_m)
summary(fit_y)

med_data <- extract_mediation(
  fit_m,
  model_y = fit_y,
  treatment = "X",
  mediator = "M"
)

# View the extracted structure
print(med_data)

# Access individual paths
cat("Estimated a path:", med_data@a_path, "\n")
cat("Estimated b path:", med_data@b_path, "\n")
cat("Estimated c' path:", med_data@c_prime, "\n")
cat("Estimated indirect effect:", med_data@a_path * med_data@b_path, "\n")

ci_dop <- ci(med_data, type = "dop", level = 0.95)
print(ci_dop)

ci_mc <- ci(med_data, type = "MC", n.mc = 10000, level = 0.95)
print(ci_mc)

# Create comparison table
ci_comparison <- data.frame(
  Method = c("DOP", "Monte Carlo"),
  Lower = c(ci_dop$CI[1], ci_mc$CI[1]),
  Upper = c(ci_dop$CI[2], ci_mc$CI[2]),
  Width = c(diff(ci_dop$CI), diff(ci_mc$CI))
)

print(ci_comparison)

pmed_result <- compute_pmed(med_data)
print(pmed_result)

# P_med interpretation:
# - 0.5 = no evidence of effect (coin flip)
# - > 0.8 = strong evidence of positive effect
# - < 0.2 = strong evidence of negative effect

# Compute sensitivity bounds
sens_bounds <- sensitivity_bounds(
  med_data,
  rho_range = c(-0.5, 0.5) # Range of confounding
)

print(sens_bounds)

# Falsification test
falsification <- falsify_mediation(med_data, threshold = 0.1)
print(falsification)

# Nonparametric bootstrap
boot_result <- bootstrap_mediation(
  med_data,
  n_boot = 2000,
  method = "nonparametric",
  ci_level = 0.95
)

print(boot_result)

# Access bootstrap distribution
hist(boot_result@boot_estimates,
  main = "Bootstrap Distribution of Indirect Effect",
  xlab = "Indirect Effect (a * b)",
  col = "steelblue",
  border = "white"
)
abline(v = indirect_true, col = "red", lwd = 2, lty = 2)
abline(v = boot_result@estimate, col = "blue", lwd = 2)
legend("topright",
  legend = c("True value", "Estimate"),
  col = c("red", "blue"),
  lwd = 2,
  lty = c(2, 1)
)

# Compile all results
results <- list(
  # Point estimates
  estimates = list(
    a = med_data@a_path,
    b = med_data@b_path,
    indirect = med_data@a_path * med_data@b_path,
    direct = med_data@c_prime
  ),

  # Confidence intervals
  ci_dop = ci_dop$CI,
  ci_mc = ci_mc$CI,
  ci_boot = c(boot_result@ci_lower, boot_result@ci_upper),

  # Effect size
  pmed = pmed_result$pmed,

  # Sensitivity
  sensitivity = sens_bounds
)

# Pretty print
cat("=== Mediation Analysis Results ===\n\n")
cat("Point Estimates:\n")
cat("  a path (X -> M):", round(results$estimates$a, 3), "\n")
cat("  b path (M -> Y):", round(results$estimates$b, 3), "\n")
cat("  c' path (X -> Y):", round(results$estimates$direct, 3), "\n")
cat("  Indirect effect:", round(results$estimates$indirect, 3), "\n\n")

cat("95% Confidence Intervals for Indirect Effect:\n")
cat("  DOP:", round(results$ci_dop, 3), "\n")
cat("  Monte Carlo:", round(results$ci_mc, 3), "\n")
cat("  Bootstrap:", round(results$ci_boot, 3), "\n\n")

cat("Effect Size:\n")
cat("  P_med:", round(results$pmed, 3), "\n\n")

library(lavaan)

# Specify mediation model
model <- "
  # Regressions
  M ~ a*X
  Y ~ b*M + c*X

  # Indirect effect
  indirect := a*b
"

# Fit model
fit_lavaan <- sem(model, data = sim_data)

# Extract mediation structure
med_data_lavaan <- extract_mediation(fit_lavaan)

# Use same analysis pipeline
ci_lavaan <- ci(med_data_lavaan, type = "dop")
print(ci_lavaan)

sessionInfo()