brycewang-stanford/mediation-meta-analyst

Mediation Meta-Analyst

Methods for synthesizing mediation effects across multiple studies

Use this skill when working on: meta-analysis of indirect effects, cross-study effect aggregation, heterogeneity assessment in mediation, individual participant data (IPD) meta-analysis, or systematic reviews of mediation studies.

---

Meta-Analysis Fundamentals

Why Meta-Analysis of Mediation is Challenging

| Challenge | Description | Solution Approach | |-----------|-------------|-------------------| | Non-normal effects | Product $ab$ is not normally distributed | Use appropriate pooling methods | | Correlated paths | $a$ and $b$ may be correlated | Model correlation structure | | Heterogeneity | Effects vary across studies | Random effects models | | Missing information | Studies report different statistics | Imputation or subset analysis | | Publication bias | Small studies with null effects unpublished | Sensitivity analysis |

Effect Size Measures

| Measure | Formula | Use Case | |---------|---------|----------| | Unstandardized $ab$ | $a \times b$ | Same scales across studies | | Partially standardized | $a \times b / SD_Y$ | Standardize by outcome only | | Fully standardized | $a' \times b'$ (standardized coefficients) | Different scales | | $R^2_{med}$ | Proportion of effect mediated | Bounded measure |

---

Aggregate Data Meta-Analysis

Fixed Effects Model

When assuming homogeneous true effects, the pooled estimate is:

$$\hat{\theta}_{FE} = \frac{\sum_i w_i \hat{\theta}_i}{\sum_i w_i}, \quad w_i = \frac{1}{\text{SE}_i^2}$$

This pooled estimate uses inverse-variance weights to optimally combine study-specific estimates.

Random Effects Model

When true effects vary across studies, the pooled estimate incorporates between-study variance:

$$\hat{\theta}_{RE} = \frac{\sum_i w_i^* \hat{\theta}_i}{\sum_i w_i^}, \quad w_i^ = \frac{1}{\text{SE}_i^2 + \hat{\tau}^2}$$

where $\hat{\tau}^2$ is the between-study variance (tau-squared). The pooled estimate under random effects provides a more generalizable result when heterogeneity is present.

Heterogeneity Assessment

Key statistics for assessing heterogeneity:

| Statistic | Formula | Interpretation | |-----------|---------|----------------| | Q statistic | $Q = \sum_i w_i(\hat{\theta}_i - \hat{\theta})^2$ | Test for heterogeneity | | I-squared ($I^2$) | $I^2 = \max(0, \frac{Q - (k-1)}{Q})$ | % variance due to heterogeneity | | tau-squared ($\tau^2$) | Between-study variance | Absolute heterogeneity | | H-squared | $H^2 = Q/(k-1)$ | Relative excess variance |

The I-squared statistic ranges from 0-100%: <25% indicates low heterogeneity, 25-75% moderate, and >75% high heterogeneity. The tau-squared provides the absolute magnitude of between-study variance.

R Implementation

#' Meta-Analysis of Indirect Effects
#'
#' @param effects Vector of indirect effect estimates
#' @param se Vector of standard errors
#' @param study_names Optional study identifiers
#' @param method "FE" for fixed effects, "RE" for random effects
#' @return Meta-analysis results
meta_indirect <- function(effects, se, study_names = NULL,
                           method = c("RE", "FE")) {
  method <- match.arg(method)
  k <- length(effects)

  if (is.null(study_names)) {
    study_names <- paste0("Study ", 1:k)
  }

  # Fixed effects weights
  w_fe <- 1 / se^2

  # Q statistic
  theta_fe <- sum(w_fe * effects) / sum(w_fe)
  Q <- sum(w_fe * (effects - theta_fe)^2)
  df <- k - 1

  # Heterogeneity
  I2 <- max(0, (Q - df) / Q)

  # DerSimonian-Laird tau^2 estimate
  c <- sum(w_fe) - sum(w_fe^2) / sum(w_fe)
  tau2 <- max(0, (Q - df) / c)

  if (method == "FE") {
    weights <- w_fe
    pooled <- theta_fe
    se_pooled <- sqrt(1 / sum(w_fe))
  } else {
    # Random effects weights
    weights <- 1 / (se^2 + tau2)
    pooled <- sum(weights * effects) / sum(weights)
    se_pooled <- sqrt(1 / sum(weights))
  }

  # Confidence interval
  ci <- pooled + c(-1.96, 1.96) * se_pooled

  # Test for heterogeneity
  p_het <- 1 - pchisq(Q, df)

  list(
    pooled_effect = pooled,
    se = se_pooled,
    ci = ci,
    z = pooled / se_pooled,
    p_value = 2 * pnorm(-abs(pooled / se_pooled)),
    heterogeneity = list(
      Q = Q,
      df = df,
      p = p_het,
      I2 = I2,
      tau2 = tau2
    ),
    study_data = data.frame(
      study = study_names,
      effect = effects,
      se = se,
      weight = weights / sum(weights)
    ),
    method = method
  )
}

---

Multivariate Meta-Analysis

Pooling Correlated Effects

When studies report both $a$ and $b$ paths:

$$\begin{pmatrix} \hat{a} \ \hat{b} \end{pmatrix} \sim N\left(\begin{pmatrix} a \ b \end{pmatrix}, \Sigma\right)$$

Two-Stage Approach

Stage 1: Extract path coefficients from each study Stage 2: Pool using multivariate random effects

#' Multivariate Meta-Analysis of Mediation Paths
#'
#' @param a_effects Vector of a path estimates
#' @param b_effects Vector of b path estimates
#' @param a_se Standard errors for a
#' @param b_se Standard errors for b
#' @param ab_cor Correlation between a and b estimates (often assumed 0)
#' @return Multivariate meta-analysis results
multivariate_meta_mediation <- function(a_effects, b_effects,
                                         a_se, b_se, ab_cor = 0) {
  library(metafor)

  k <- length(a_effects)

  # Construct variance-covariance matrices for each study
  V_list <- lapply(1:k, function(i) {
    cov_ab <- ab_cor * a_se[i] * b_se[i]
    matrix(c(a_se[i]^2, cov_ab, cov_ab, b_se[i]^2), 2, 2)
  })

  # Stack effects
  yi <- c(rbind(a_effects, b_effects))
  vi <- unlist(lapply(V_list, as.vector))

  # Create V matrix (block diagonal)
  V <- bldiag(V_list)

  # Fit multivariate model
  # Effect type indicator
  effect_type <- rep(c("a", "b"), k)
  study_id <- rep(1:k, each = 2)

  fit <- rma.mv(yi = yi, V = V,
                mods = ~ effect_type - 1,
                random = ~ effect_type | study_id,
                struct = "UN",
                data = data.frame(yi, effect_type, study_id))

  # Extract pooled estimates
  pooled_a <- coef(fit)["effect_typea"]
  pooled_b <- coef(fit)["effect_typeb"]

  # Compute indirect effect and CI via delta method
  vcov_pooled <- vcov(fit)
  indirect <- pooled_a * pooled_b

  # Delta method SE
  grad <- c(pooled_b, pooled_a)
  se_indirect <- sqrt(t(grad) %*% vcov_pooled %*% grad)

  list(
    pooled_a = pooled_a,
    pooled_b = pooled_b,
    pooled_indirect = indirect,
    se_indirect = as.numeric(se_indirect),
    ci_indirect = indirect + c(-1.96, 1.96) * as.numeric(se_indirect),
    model_fit = fit
  )
}

---

Individual Participant Data (IPD) Meta-Analysis

One-Stage Approach

Pool all data and fit single model with study-level random effects:

#' IPD Meta-Analysis for Mediation
#'
#' @param data Combined dataset with study indicator
#' @param study_var Name of study variable
#' @param treatment Name of treatment variable
#' @param mediator Name of mediator variable
#' @param outcome Name of outcome variable
#' @return IPD meta-analysis results
ipd_meta_mediation <- function(data, study_var, treatment, mediator, outcome) {
  library(lme4)

  # Mediator model with random slopes
  m_formula <- as.formula(paste(
    mediator, "~", treatment, "+ (1 +", treatment, "|", study_var, ")"
  ))
  m_model <- lmer(m_formula, data = data)

  # Outcome model with random slopes
  y_formula <- as.formula(paste(
    outcome, "~", treatment, "+", mediator,
    "+ (1 +", treatment, "+", mediator, "|", study_var, ")"
  ))
  y_model <- lmer(y_formula, data = data)

  # Extract fixed effects (pooled estimates)
  a <- fixef(m_model)[treatment]
  b <- fixef(y_model)[mediator]
  c_prime <- fixef(y_model)[treatment]

  # Indirect effect
  indirect <- a * b

  # Bootstrap for CI
  boot_indirect <- replicate(1000, {
    boot_idx <- sample(nrow(data), replace = TRUE)
    boot_data <- data[boot_idx, ]

    m_boot <- tryCatch(
      lmer(m_formula, data = boot_data),
      error = function(e) NULL
    )
    y_boot <- tryCatch(
      lmer(y_formula, data = boot_data),
      error = function(e) NULL
    )

    if (is.null(m_boot) || is.null(y_boot)) return(NA)

    fixef(m_boot)[treatment] * fixef(y_boot)[mediator]
  })

  boot_indirect <- boot_indirect[!is.na(boot_indirect)]

  list(
    pooled_a = a,
    pooled_b = b,
    pooled_c_prime = c_prime,
    pooled_indirect = indirect,
    pooled_total = indirect + c_prime,
    se_indirect = sd(boot_indirect),
    ci_indirect = quantile(boot_indirect, c(0.025, 0.975)),
    n_studies = length(unique(data[[study_var]])),
    n_total = nrow(data),
    m_model = m_model,
    y_model = y_model
  )
}

Two-Stage Approach

Estimate effects within each study, then pool:

#' Two-Stage IPD Meta-Analysis
#'
#' @param data Combined dataset
#' @param study_var Study identifier
#' @return Two-stage meta-analysis results
two_stage_ipd <- function(data, study_var, treatment, mediator, outcome) {

  studies <- unique(data[[study_var]])
  k <- length(studies)

  # Stage 1: Study-specific estimates
  study_results <- lapply(studies, function(s) {
    study_data <- data[data[[study_var]] == s, ]

    # Fit models
    m_model <- lm(as.formula(paste(mediator, "~", treatment)), data = study_data)
    y_model <- lm(as.formula(paste(outcome, "~", treatment, "+", mediator)),
                  data = study_data)

    a <- coef(m_model)[treatment]
    b <- coef(y_model)[mediator]

    # Delta method SE for indirect
    se_a <- sqrt(vcov(m_model)[treatment, treatment])
    se_b <- sqrt(vcov(y_model)[mediator, mediator])
    se_indirect <- sqrt(a^2 * se_b^2 + b^2 * se_a^2)

    data.frame(
      study = s,
      n = nrow(study_data),
      a = a,
      b = b,
      indirect = a * b,
      se_a = se_a,
      se_b = se_b,
      se_indirect = se_indirect
    )
  })

  study_df <- do.call(rbind, study_results)

  # Stage 2: Pool indirect effects
  meta_result <- meta_indirect(
    effects = study_df$indirect,
    se = study_df$se_indirect,
    study_names = study_df$study,
    method = "RE"
  )

  list(
    stage1 = study_df,
    stage2 = meta_result,
    pooled_indirect = meta_result$pooled_effect,
    ci = meta_result$ci,
    I2 = meta_result$heterogeneity$I2
  )
}

---

Publication Bias

Detection Methods

| Method | Description | Limitation | |--------|-------------|------------| | Funnel plot | SE vs effect plot | Visual, subjective | | Egger's test | Regression of effect on SE | Low power | | Trim-and-fill | Impute missing studies | Assumes specific mechanism | | PET-PEESE | Conditional regression | Requires assumptions | | Selection models | Model publication process | Complex, sensitive |

R Implementation

#' Publication Bias Assessment for Mediation Meta-Analysis
#'
#' @param effects Indirect effect estimates
#' @param se Standard errors
#' @return Publication bias diagnostics
publication_bias_mediation <- function(effects, se) {
  library(metafor)

  # Fit random effects model
  res <- rma(yi = effects, sei = se, method = "REML")

  # Funnel plot data
  funnel_data <- data.frame(
    effect = effects,
    se = se,
    precision = 1/se
  )

  # Egger's test
  egger <- regtest(res, model = "lm")

  # Trim and fill
  tf <- trimfill(res)

  # PET-PEESE (if significant, use PEESE; otherwise PET)
  pet <- lm(effects ~ se, weights = 1/se^2)
  peese <- lm(effects ~ I(se^2), weights = 1/se^2)

  pet_est <- coef(pet)[1]
  peese_est <- coef(peese)[1]

  # Use PEESE if PET significant, else PET
  if (coef(summary(pet))[2, 4] < 0.10) {
    adjusted_estimate <- peese_est
    method_used <- "PEESE"
  } else {
    adjusted_estimate <- pet_est
    method_used <- "PET"
  }

  list(
    original_estimate = coef(res),
    egger_test = list(
      z = egger$zval,
      p = egger$pval,
      interpretation = ifelse(egger$pval < 0.10,
                              "Evidence of funnel asymmetry",
                              "No strong evidence of asymmetry")
    ),
    trim_fill = list(
      original_k = res$k,
      imputed_k = tf$k0,
      adjusted_estimate = coef(tf),
      adjusted_ci = c(tf$ci.lb, tf$ci.ub)
    ),
    pet_peese = list(
      pet_estimate = pet_est,
      peese_estimate = peese_est,
      method_used = method_used,
      adjusted_estimate = adjusted_estimate
    ),
    funnel_data = funnel_data
  )
}

---

Moderator Analysis

Meta-Regression

Test whether study-level characteristics explain heterogeneity:

#' Meta-Regression for Mediation Effects
#'
#' @param effects Indirect effect estimates
#' @param se Standard errors
#' @param moderators Data frame of moderator variables
#' @return Meta-regression results
meta_regression_mediation <- function(effects, se, moderators) {
  library(metafor)

  # Build formula from moderator names
  mod_formula <- as.formula(paste("~", paste(names(moderators), collapse = " + ")))

  # Fit mixed-effects meta-regression
  res <- rma(yi = effects, sei = se,
             mods = mod_formula,
             data = moderators,
             method = "REML")

  # R^2 analog
  res_null <- rma(yi = effects, sei = se, method = "REML")
  R2 <- max(0, (res_null$tau2 - res$tau2) / res_null$tau2)

  # Test for residual heterogeneity
  QE_test <- list(
    QE = res$QE,
    df = res$k - res$p,
    p = res$QEp
  )

  list(
    coefficients = coef(summary(res)),
    tau2_residual = res$tau2,
    I2_residual = res$I2,
    R2 = R2,
    residual_heterogeneity = QE_test,
    model = res
  )
}

Subgroup Analysis

#' Subgroup Analysis for Mediation Meta-Analysis
#'
#' @param effects Indirect effect estimates
#' @param se Standard errors
#' @param subgroup Factor variable defining subgroups
#' @return Subgroup analysis results
subgroup_analysis <- function(effects, se, subgroup) {

  groups <- unique(subgroup)

  # Within-group estimates
  group_results <- lapply(groups, function(g) {
    idx <- subgroup == g
    meta_indirect(effects[idx], se[idx], method = "RE")
  })
  names(group_results) <- groups

  # Extract pooled estimates
  group_effects <- sapply(group_results, function(x) x$pooled_effect)
  group_se <- sapply(group_results, function(x) x$se)
  group_k <- sapply(group_results, function(x) length(x$study_data$effect))

  # Test for subgroup differences
  # Q_between = Q_total - sum(Q_within)
  overall <- meta_indirect(effects, se, method = "RE")
  Q_total <- overall$heterogeneity$Q

  Q_within <- sapply(group_results, function(x) x$heterogeneity$Q)
  Q_between <- Q_total - sum(Q_within)
  df_between <- length(groups) - 1
  p_between <- 1 - pchisq(Q_between, df_between)

  list(
    subgroup_estimates = data.frame(
      subgroup = groups,
      k = group_k,
      effect = group_effects,
      se = group_se,
      ci_lower = group_effects - 1.96 * group_se,
      ci_upper = group_effects + 1.96 * group_se
    ),
    test_for_differences = list(
      Q_between = Q_between,
      df = df_between,
      p = p_between,
      interpretation = ifelse(p_between < 0.05,
                              "Significant subgroup differences",
                              "No significant subgroup differences")
    ),
    group_results = group_results
  )
}

---

Reporting Checklist

PRISMA for Mediation Meta-Analysis

• [ ] Search strategy documented
• [ ] Inclusion/exclusion criteria specified
• [ ] Effect measure defined (standardized vs. unstandardized)
• [ ] Method for extracting/computing indirect effects stated
• [ ] Heterogeneity statistics reported ($I^2$, $\tau^2$, $Q$)
• [ ] Forest plot included
• [ ] Publication bias assessed
• [ ] Sensitivity analyses conducted
• [ ] GRADE assessment for certainty

Forest Plot Template

#' Create Forest Plot for Mediation Meta-Analysis
#'
#' @param meta_result Result from meta_indirect()
#' @return ggplot2 forest plot
forest_plot_mediation <- function(meta_result) {
  library(ggplot2)

  df <- meta_result$study_data
  df$ci_lower <- df$effect - 1.96 * df$se
  df$ci_upper <- df$effect + 1.96 * df$se

  # Add pooled estimate
  pooled <- data.frame(
    study = "Pooled",
    effect = meta_result$pooled_effect,
    se = meta_result$se,
    weight = NA,
    ci_lower = meta_result$ci[1],
    ci_upper = meta_result$ci[2]
  )
  df <- rbind(df, pooled)
  df$study <- factor(df$study, levels = rev(df$study))

  ggplot(df, aes(x = effect, y = study)) +
    geom_vline(xintercept = 0, linetype = "dashed", color = "gray50") +
    geom_point(aes(size = weight)) +
    geom_errorbarh(aes(xmin = ci_lower, xmax = ci_upper), height = 0.2) +
    geom_point(data = df[df$study == "Pooled", ],
               shape = 18, size = 5, color = "darkred") +
    labs(
      x = "Indirect Effect",
      y = "",
      title = "Forest Plot: Meta-Analysis of Indirect Effects",
      subtitle = sprintf("I² = %.1f%%, τ² = %.4f",
                         meta_result$heterogeneity$I2 * 100,
                         meta_result$heterogeneity$tau2)
    ) +
    theme_minimal() +
    theme(legend.position = "none")
}

---

References

Meta-Analysis Methods

• Borenstein, M., et al. (2009). Introduction to Meta-Analysis
• Higgins, J. P., & Green, S. (2011). Cochrane Handbook for Systematic Reviews

Mediation Meta-Analysis

• Cheung, M. W. L. (2015). Meta-Analysis: A Structural Equation Modeling Approach
• MacKinnon, D. P. (2008). Introduction to Statistical Mediation Analysis

Publication Bias

• Rothstein, H. R., et al. (2005). Publication Bias in Meta-Analysis
• Stanley, T. D., & Doucouliagos, H. (2014). Meta-regression approximations

Software

• Viechtbauer, W. (2010). Conducting meta-analyses in R with the metafor package
• Cheung, M. W. L. (2015). metaSEM: Meta-analysis using structural equation modeling

---

Version: 1.0.0 Created: 2025-12-09 Domain: Meta-analysis of mediation effects Applications: Systematic reviews, research synthesis, evidence aggregation