Meta-Regression (Meta-Analytic Regression)

# Meta-Regression: Moderators of MBI Effects on Anxiety
# Demonstrates testing continuous and categorical moderators

library(metafor)      # rma() for meta-regression
library(dplyr)
library(ggplot2)

# === STEP 1: Simulate Meta-Regression Dataset ===
# In practice: data <- read.csv("meta_regression_data.csv")
# Required: study_id, effect_size, variance, moderators (duration, age, sample_type)

set.seed(2025)
k <- 28  # Number of studies

# Moderators
duration <- sample(4:12, k, replace=TRUE)  # Intervention duration (4-12 weeks)
mean_age <- runif(k, min=24, max=65)       # Mean age (24-65 years)
sample_type <- sample(c("Clinical", "Non-clinical"), k, replace=TRUE, prob=c(0.6, 0.4))

# Generate effect sizes with moderator relationships
# True model: θ = 0.30 + 0.04*duration + 0.005*age + 0.25*clinical + u + e
# β₀ = 0.30 (intercept, effect at duration=0, age=0, non-clinical)
# β_duration = 0.04 (longer interventions → larger effects)
# β_age = 0.005 (older participants → slightly larger effects)
# β_clinical = 0.25 (clinical samples → larger effects than non-clinical)

clinical_dummy <- ifelse(sample_type == "Clinical", 1, 0)
true_effects <- 0.30 + 0.04*duration + 0.005*mean_age + 0.25*clinical_dummy

# Add residual between-study heterogeneity (τ²_residual ≈ 0.04)
true_effects <- true_effects + rnorm(k, 0, sd=0.20)

# Sample sizes
n_treat <- sample(30:80, k, replace=TRUE)
n_control <- sample(30:80, k, replace=TRUE)
total_n <- n_treat + n_control

# Sampling standard errors
sampling_se <- sqrt((n_treat + n_control)/(n_treat * n_control) + 
                     true_effects^2 / (2*(n_treat + n_control)))

# Observed effect sizes
observed_g <- rnorm(k, mean=true_effects, sd=sampling_se)
variance_g <- sampling_se^2

meta_data <- data.frame(
  study_id = paste0("Study_", 1:k),
  author_year = paste0(LETTERS[1:k], " (", 2010:2024, ")")[1:k],
  hedges_g = observed_g,
  variance = variance_g,
  se = sqrt(variance_g),
  duration_weeks = duration,
  mean_age = mean_age,
  sample_type = sample_type,
  clinical_dummy = clinical_dummy,
  n_treatment = n_treat,
  n_control = n_control,
  total_n = total_n
)

print("=== Meta-Regression Dataset ===")
print(head(meta_data, 10))
cat("\nTotal k =", k, "studies, N =", sum(meta_data$total_n), "participants")
cat("\nModerators: Duration(", range(duration)[1], "-", range(duration)[2], 
    "weeks), Age(", round(min(mean_age)), "-", round(max(mean_age)), "years)")
cat("\nSample type: Clinical =", sum(clinical_dummy), ", Non-clinical =", 
    sum(1-clinical_dummy))

# === STEP 2: Baseline Random-Effects Meta-Analysis (No Moderators) ===
# Assess total heterogeneity before moderator analysis

baseline_model <- rma(yi = hedges_g, vi = variance, data = meta_data, 
                      method = "REML", slab = author_year)

cat("\n\n=== BASELINE: Random-Effects Meta-Analysis(No Moderators) ===")
print(baseline_model)

tau2_total <- baseline_model$tau2
I2_total <- baseline_model$I2
Q_total <- baseline_model$QE
Q_pval_total <- baseline_model$QEp

cat("\n=== Total Heterogeneity(Before Moderators) ===")
cat("\nτ²_total =", round(tau2_total, 4))
cat("\nI²_total =", round(I2_total, 1), "%")
cat("\nQ(", baseline_model$k-1, ") =", round(Q_total, 2), ", p =", 
    format.pval(Q_pval_total, digits=3))

if (I2_total < 50) {
  cat("\n\nWARNING: Moderate heterogeneity(I²<50%). Meta-regression may be premature.")
  cat("\nConsider whether moderator analysis is justified.")
} else {
  cat("\n\nSubstantial heterogeneity detected(I²≥50%). Meta-regression warranted.")
}

# === STEP 3: Meta-Regression with Multiple Moderators ===
# Model: θ = β₀ + β₁*duration + β₂*age + β₃*clinical + u + e

meta_reg_model <- rma(yi = hedges_g, vi = variance, 
                      mods = ~ duration_weeks + mean_age + clinical_dummy,
                      data = meta_data, 
                      method = "REML",
                      slab = author_year)

cat("\n\n=== META-REGRESSION: Multivariable Model ===")
print(meta_reg_model)

# Extract results
beta_intercept <- meta_reg_model$beta[1]
beta_duration <- meta_reg_model$beta[2]
beta_age <- meta_reg_model$beta[3]
beta_clinical <- meta_reg_model$beta[4]

ci_duration <- c(meta_reg_model$ci.lb[2], meta_reg_model$ci.ub[2])
ci_age <- c(meta_reg_model$ci.lb[3], meta_reg_model$ci.ub[3])
ci_clinical <- c(meta_reg_model$ci.lb[4], meta_reg_model$ci.ub[4])

p_duration <- meta_reg_model$pval[2]
p_age <- meta_reg_model$pval[3]
p_clinical <- meta_reg_model$pval[4]

# QM test (omnibus test for all moderators)
QM <- meta_reg_model$QM
QM_df <- meta_reg_model$m  # Number of moderators
QM_pval <- meta_reg_model$QMp

# Residual heterogeneity
tau2_residual <- meta_reg_model$tau2
I2_residual <- meta_reg_model$I2
Q_residual <- meta_reg_model$QE
Q_residual_pval <- meta_reg_model$QEp

cat("\n=== Moderator Coefficients ===")
cat("\nIntercept(β₀) =", round(beta_intercept, 3))
cat("\nDuration(β₁) =", round(beta_duration, 3), 
    ", 95% CI [", round(ci_duration[1], 3), ",", round(ci_duration[2], 3), "],",
    "p =", format.pval(p_duration, digits=3))
cat("\nAge(β₂) =", round(beta_age, 4), 
    ", 95% CI [", round(ci_age[1], 4), ",", round(ci_age[2], 4), "],",
    "p =", format.pval(p_age, digits=3))
cat("\nClinical vs. Non-clinical(β₃) =", round(beta_clinical, 3), 
    ", 95% CI [", round(ci_clinical[1], 3), ",", round(ci_clinical[2], 3), "],",
    "p =", format.pval(p_clinical, digits=3))

cat("\n\n=== Omnibus Test(QM) ===")
cat("\nQM(", QM_df, ") =", round(QM, 2), ", p =", format.pval(QM_pval, digits=3))
if (QM_pval < 0.05) {
  cat("\n→ At least one moderator significantly explains heterogeneity")
} else {
  cat("\n→ Moderators collectively do not explain significant heterogeneity")
}

# === STEP 4: Calculate R² (% Variance Explained) ===
# R² = (τ²_total - τ²_residual) / τ²_total
R2 <- max(0, 100 * (tau2_total - tau2_residual) / tau2_total)

cat("\n\n=== Variance Explained(R²) ===")
cat("\nτ²_total(baseline) =", round(tau2_total, 4))
cat("\nτ²_residual(after moderators) =", round(tau2_residual, 4))
cat("\nR² =", round(R2, 1), "% of between-study variance explained")

if (R2 < 25) {
  R2_interp <- "minimal"
} else if (R2 < 50) {
  R2_interp <- "moderate"
} else if (R2 < 75) {
  R2_interp <- "substantial"
} else {
  R2_interp <- "extensive"
}
cat("\nInterpretation: Moderators explain", R2_interp, "heterogeneity")

cat("\n\n=== Residual Heterogeneity ===")
cat("\nI²_residual =", round(I2_residual, 1), "%")
cat("\nQ_residual(", meta_reg_model$k - QM_df - 1, ") =", round(Q_residual, 2),
    ", p =", format.pval(Q_residual_pval, digits=3))

if (I2_residual > 50) {
  cat("\n→ Substantial residual heterogeneity remains(I²>50%); additional moderators needed")
} else if (I2_residual > 25) {
  cat("\n→ Moderate residual heterogeneity remains")
} else {
  cat("\n→ Minimal residual heterogeneity; moderators explain most variance")
}

# === STEP 5: Univariate Meta-Regressions (Compare to Multivariable) ===
cat("\n\n=== UNIVARIATE META-REGRESSIONS(for comparison) ===")

# Duration only
uni_duration <- rma(yi = hedges_g, vi = variance, mods = ~ duration_weeks,
                    data = meta_data, method = "REML")
cat("\nDuration(univariate): β =", round(coef(uni_duration)[2], 3),
    ", p =", format.pval(uni_duration$pval[2], digits=3),
    ", R² =", round(max(0, 100*(tau2_total - uni_duration$tau2)/tau2_total), 1), "%")

# Age only
uni_age <- rma(yi = hedges_g, vi = variance, mods = ~ mean_age,
               data = meta_data, method = "REML")
cat("\nAge(univariate): β =", round(coef(uni_age)[2], 4),
    ", p =", format.pval(uni_age$pval[2], digits=3),
    ", R² =", round(max(0, 100*(tau2_total - uni_age$tau2)/tau2_total), 1), "%")

# Clinical vs. non-clinical only
uni_clinical <- rma(yi = hedges_g, vi = variance, mods = ~ clinical_dummy,
                    data = meta_data, method = "REML")
cat("\nClinical vs. Non-clinical(univariate): β =", round(coef(uni_clinical)[2], 3),
    ", p =", format.pval(uni_clinical$pval[2], digits=3),
    ", R² =", round(max(0, 100*(tau2_total - uni_clinical$tau2)/tau2_total), 1), "%")

cat("\n\nNote: Compare univariate vs. multivariable estimates to assess confounding.")
cat("\nLarge changes suggest moderators are correlated(confounding/multicollinearity).")

# === STEP 6: Multicollinearity Check ===
cat("\n\n=== MULTICOLLINEARITY DIAGNOSTICS ===")

# Correlation matrix
cor_matrix <- cor(meta_data[, c("duration_weeks", "mean_age", "clinical_dummy")])
cat("\nCorrelation Matrix:")
print(round(cor_matrix, 3))

if (any(abs(cor_matrix[upper.tri(cor_matrix)]) > 0.70)) {
  cat("\nWARNING: High correlation detected(r > .70). Consider removing one moderator.")
} else {
  cat("\nNo severe multicollinearity detected(all r < .70).")
}

# Variance Inflation Factors (VIF) - approximate
# VIF = 1 / (1 - R²_j), where R²_j is R² from regressing moderator j on other moderators
vif_duration <- 1 / (1 - summary(lm(duration_weeks ~ mean_age + clinical_dummy, 
                                     data=meta_data))$r.squared)
vif_age <- 1 / (1 - summary(lm(mean_age ~ duration_weeks + clinical_dummy, 
                                data=meta_data))$r.squared)
vif_clinical <- 1 / (1 - summary(lm(clinical_dummy ~ duration_weeks + mean_age, 
                                     data=meta_data))$r.squared)

cat("\nVariance Inflation Factors(VIF):")
cat("\n  Duration: VIF =", round(vif_duration, 2))
cat("\n  Age: VIF =", round(vif_age, 2))
cat("\n  Clinical: VIF =", round(vif_clinical, 2))

if (max(vif_duration, vif_age, vif_clinical) > 5) {
  cat("\nWARNING: VIF > 5 detected. Problematic multicollinearity. Consider univariate models.")
} else if (max(vif_duration, vif_age, vif_clinical) > 3) {
  cat("\nModerate multicollinearity(VIF 3-5). Interpret coefficients cautiously.")
} else {
  cat("\nAcceptable multicollinearity(VIF < 3).")
}

# === STEP 7: Bubble Plots (Moderator Relationships) ===
par(mfrow=c(2,2), mar=c(4,4,3,2))

# Duration vs. Effect Size
weights_meta <- 1 / (meta_data$variance + tau2_residual)
weights_scaled <- (weights_meta / max(weights_meta)) * 500

plot(meta_data$duration_weeks, meta_data$hedges_g,
     pch=19, cex=sqrt(weights_scaled/100), col=rgb(0,0,1,0.5),
     xlab="Intervention Duration(weeks)", ylab="Hedges' g",
     main="Duration vs. Effect Size")
abline(a = beta_intercept + beta_age*mean(meta_data$mean_age) + 
           beta_clinical*mean(meta_data$clinical_dummy),
       b = beta_duration, col="red", lwd=2)
text(max(meta_data$duration_weeks)*0.7, max(meta_data$hedges_g)*0.9,
     paste0("β = ", round(beta_duration, 3), ", p = ", 
            format.pval(p_duration, digits=2)), cex=0.9)

# Age vs. Effect Size
plot(meta_data$mean_age, meta_data$hedges_g,
     pch=19, cex=sqrt(weights_scaled/100), col=rgb(0,0,1,0.5),
     xlab="Mean Age(years)", ylab="Hedges' g",
     main="Age vs. Effect Size")
abline(a = beta_intercept + beta_duration*mean(meta_data$duration_weeks) + 
           beta_clinical*mean(meta_data$clinical_dummy),
       b = beta_age, col="red", lwd=2)
text(mean(meta_data$mean_age), max(meta_data$hedges_g)*0.9,
     paste0("β = ", round(beta_age, 4), ", p = ", 
            format.pval(p_age, digits=2)), cex=0.9)

# Clinical vs. Non-clinical (Boxplot)
boxplot(hedges_g ~ sample_type, data=meta_data,
        xlab="Sample Type", ylab="Hedges' g",
        main="Clinical vs. Non-clinical",
        col=c("lightblue", "lightgreen"))
points(jitter(as.numeric(factor(meta_data$sample_type)), 0.2), 
       meta_data$hedges_g,
       pch=19, cex=sqrt(weights_scaled/100), col=rgb(0,0,0,0.4))
text(1.5, max(meta_data$hedges_g)*0.9,
     paste0("β = ", round(beta_clinical, 3), ", p = ", 
            format.pval(p_clinical, digits=2)), cex=0.9)

# Residual Plot
residuals_meta <- residuals(meta_reg_model)
fitted_meta <- fitted(meta_reg_model)
plot(fitted_meta, residuals_meta,
     pch=19, cex=sqrt(weights_scaled/100), col=rgb(0,0,1,0.5),
     xlab="Fitted Values", ylab="Residuals",
     main="Residual Plot")
abline(h=0, col="red", lty=2, lwd=2)

par(mfrow=c(1,1))

# === STEP 8: Influence Analysis (Leave-One-Out) ===
cat("\n\n=== INFLUENCE ANALYSIS(Leave-One-Out) ===")

influence_results <- list()
for (i in 1:k) {
  # Remove study i
  data_loo <- meta_data[-i, ]
  
  # Re-run meta-regression
  model_loo <- rma(yi = hedges_g, vi = variance, 
                   mods = ~ duration_weeks + mean_age + clinical_dummy,
                   data = data_loo, method = "REML")
  
  influence_results[[i]] <- data.frame(
    study_removed = meta_data$author_year[i],
    beta_duration = coef(model_loo)[2],
    beta_age = coef(model_loo)[3],
    beta_clinical = coef(model_loo)[4],
    R2 = max(0, 100*(tau2_total - model_loo$tau2)/tau2_total)
  )
}

influence_df <- do.call(rbind, influence_results)

cat("\nEffect size range across leave-one-out:")
cat("\n  Duration: β ∈ [", round(min(influence_df$beta_duration), 3), ",", 
    round(max(influence_df$beta_duration), 3), "]")
cat("\n  Age: β ∈ [", round(min(influence_df$beta_age), 4), ",", 
    round(max(influence_df$beta_age), 4), "]")
cat("\n  Clinical: β ∈ [", round(min(influence_df$beta_clinical), 3), ",", 
    round(max(influence_df$beta_clinical), 3), "]")
cat("\n  R²: ", round(min(influence_df$R2), 1), "% to", round(max(influence_df$R2), 1), "%")

# Check if any study substantially changes results
if (max(abs(influence_df$beta_duration - beta_duration)) > abs(beta_duration)*0.3) {
  cat("\nWARNING: Removing ≥1 study changes Duration coefficient by >30%. Investigate influential studies.")
} else {
  cat("\n→ Moderator estimates robust to individual study removal")
}

# === STEP 9: APA-Style Reporting ===
cat("\n\n=== APA-STYLE REPORT ===")
cat("\nA random-effects meta-regression of", k, "RCTs(N =", sum(meta_data$total_n),
    "participants) examined\nwhether intervention duration, participant age, and sample type(clinical vs. non-clinical)\nmoderate mindfulness-based intervention effects on anxiety.\n\nBaseline random-effects meta-analysis(no moderators) revealed substantial heterogeneity\n(I² =", round(I2_total, 1), "%, τ² =", round(tau2_total, 3), 
    ", Q(", baseline_model$k-1, ") =", round(Q_total, 2),
    ", p", ifelse(Q_pval_total < 0.001, " < .001", paste0(" = ", round(Q_pval_total, 3))),
    "),\nwarranting moderator analysis.\n\nMeta-regression results(k/p =", round(k/QM_df, 1), 
    ") showed that moderators collectively explained\nsignificant heterogeneity(QM(", QM_df, ") =", round(QM, 2), ", p", 
    ifelse(QM_pval < 0.001, " < .001", paste0(" = ", round(QM_pval, 3))), 
    "), accounting for\nR² =", round(R2, 1), "% of between-study variance.\n\nIntervention duration significantly predicted effect sizes(β =", round(beta_duration, 3),
    ", 95% CI [", round(ci_duration[1], 3), ",", round(ci_duration[2], 3), "], p",
    ifelse(p_duration < 0.001, " < .001", paste0(" = ", round(p_duration, 3))),
    "):\neach additional week of intervention was associated with a", 
    round(beta_duration, 3), "increase in effect size.\n")

if (p_age < 0.05) {
  cat("Participant age also significantly predicted effects(β =", round(beta_age, 4),
      ", 95% CI [", round(ci_age[1], 4), ",", round(ci_age[2], 4), "], p =",
      round(p_age, 3), "),\nsuggesting", ifelse(beta_age > 0, "larger", "smaller"), 
      "effects in studies with older participants.\n")
} else {
  cat("Participant age did not significantly predict effects(β =", round(beta_age, 4),
      ", p =", round(p_age, 3), ").\n")
}

if (p_clinical < 0.05) {
  cat("Sample type moderated effects(β =", round(beta_clinical, 3),
      ", 95% CI [", round(ci_clinical[1], 3), ",", round(ci_clinical[2], 3), "], p",
      ifelse(p_clinical < 0.001, " < .001", paste0(" = ", round(p_clinical, 3))),
      "):\nclinical samples showed", ifelse(beta_clinical > 0, "larger", "smaller"),
      "effects than non-clinical samples.\n")
} else {
  cat("Sample type did not significantly moderate effects(β =", round(beta_clinical, 3),
      ", p =", round(p_clinical, 3), ").\n")
}

cat("\nResidual heterogeneity after accounting for moderators was", 
    ifelse(I2_residual > 50, "substantial", ifelse(I2_residual > 25, "moderate", "minimal")),
    "\n(I²_residual =", round(I2_residual, 1), "%, τ²_residual =", round(tau2_residual, 4),
    "), suggesting")
if (I2_residual > 50) {
  cat("\nadditional unmeasured moderators contribute to between-study variance.")
} else {
  cat("\nmoderators explain most heterogeneity.")
}

cat("\n\nInfluence analysis confirmed robustness: moderator estimates varied minimally across\nleave-one-out analyses, indicating no single study disproportionately influenced findings.\n\nConclusion: Longer interventions and", ifelse(p_clinical < 0.05, "clinical samples", ""),
    "were associated with larger MBI\neffects on anxiety, explaining", round(R2, 1), 
    "% of heterogeneity. These findings are ASSOCIATIONAL\n(not causal) and require replication.")

if (R2 < 50) {
  cat(" Additional moderators(e.g., intervention format,\ntherapist training, outcome measure) should be explored in future meta-analyses.")
}
cat("\n")