brycewang-stanford 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