Leave-One-Out Meta-Analysis (Influence Diagnostics)

# Leave-One-Out Meta-Analysis: Influence Diagnostics for CBT Depression Meta-Analysis
# Systematically assess impact of each study on pooled estimate

library(metafor)      # leave1out() for systematic sensitivity
library(dplyr)
library(ggplot2)
library(tidyr)

# === STEP 1: Simulate Meta-Analytic Dataset ===
# In practice: data <- read.csv("meta_analysis_data.csv")
# Required: study_id, effect_size (Hedges' g), variance, sample sizes

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

# Simulate effect sizes with moderate heterogeneity
# True effects vary: mean θ=0.70, between-study SD τ=0.20
true_effects <- rnorm(k, mean=0.70, sd=0.20)

# One influential study with larger effect (Study 5)
true_effects[5] <- 1.10  # Outlier-ish but not extreme

# Sample sizes vary
n_treat <- sample(40:100, k, replace=TRUE)
n_control <- sample(40:100, k, replace=TRUE)
total_n <- n_treat + n_control

# One study with large sample (Study 3) - influential due to precision
n_treat[3] <- 150
n_control[3] <- 150
total_n[3] <- 300

# Observed effect sizes (true effect + sampling error)
sampling_se <- sqrt((n_treat + n_control)/(n_treat * n_control) + 
                     true_effects^2 / (2*(n_treat + n_control)))
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], " et al.(20", 10:24, ")"),
  hedges_g = observed_g,
  variance = variance_g,
  se = sqrt(variance_g),
  n_treatment = n_treat,
  n_control = n_control,
  total_n = total_n
)

print("=== Meta-Analytic Dataset ===")
print(meta_data)
cat("\nTotal N =", sum(meta_data$total_n), "participants across", k, "studies\n")

# === STEP 2: Primary Random-Effects Meta-Analysis ===
re_model <- rma(yi = hedges_g, vi = variance, data = meta_data, 
                method = "REML", slab = author_year)

print("\n=== PRIMARY RANDOM-EFFECTS META-ANALYSIS ===")
print(re_model)

# Extract primary results
pooled_g <- as.numeric(re_model$beta)
ci_lower <- re_model$ci.lb
ci_upper <- re_model$ci.ub
ci_width <- ci_upper - ci_lower
se_pooled <- re_model$se
p_value <- re_model$pval

# Heterogeneity
tau2 <- re_model$tau2
tau <- sqrt(tau2)
I2 <- re_model$I2
Q <- re_model$QE
Q_pval <- re_model$QEp

cat("\n=== PRIMARY MODEL RESULTS ===")
cat("\nPooled Hedges' g =", round(pooled_g, 3))
cat("\n95% CI: [", round(ci_lower, 3), ",", round(ci_upper, 3), "]")
cat("\nSE =", round(se_pooled, 4))
cat("\np-value =", format.pval(p_value, digits=3))
cat("\n\nHeterogeneity: I² =", round(I2, 1), "%, τ² =", round(tau2, 4))
cat("\nCochran's Q(", re_model$k-1, ") =", round(Q, 2), ", p =",
    format.pval(Q_pval, digits=3))

if (I2 < 25) {
  heterogeneity_interp <- "low"
} else if (I2 < 50) {
  heterogeneity_interp <- "moderate"
} else if (I2 < 75) {
  heterogeneity_interp <- "substantial"
} else {
  heterogeneity_interp <- "considerable"
}
cat("\nInterpretation:", heterogeneity_interp, "heterogeneity\n")

# === STEP 3: LEAVE-ONE-OUT SENSITIVITY ANALYSIS ===
cat("\n\n=== LEAVE-ONE-OUT SENSITIVITY ANALYSIS ===")
cat("\nSystematically removing each of", k, "studies and re-estimating pooled effect...\n\n")

# metafor's leave1out() function conducts k meta-analyses, each excluding one study
loo_results <- leave1out(re_model, digits=3)

print(loo_results)

# Extract leave-one-out results for detailed analysis
loo_df <- data.frame(
  study_removed = meta_data$author_year,
  study_id = meta_data$study_id,
  original_effect = meta_data$hedges_g,
  original_weight = weights(re_model),
  pooled_g_loo = loo_results$estimate,
  se_loo = loo_results$se,
  ci_lower_loo = loo_results$ci.lb,
  ci_upper_loo = loo_results$ci.ub,
  ci_width_loo = loo_results$ci.ub - loo_results$ci.lb,
  pval_loo = loo_results$pval,
  Q_loo = loo_results$Q,
  Qp_loo = loo_results$Qp,
  tau2_loo = loo_results$tau2,
  I2_loo = loo_results$I2,
  H2_loo = loo_results$H2
)

# Calculate change metrics
loo_df <- loo_df %>%
  mutate(
    delta_g = pooled_g - pooled_g_loo,  # Change in pooled estimate
    delta_g_abs = abs(delta_g),
    delta_g_standardized = delta_g / se_pooled,  # Standardized change (in SE units)
    delta_ci_width = ci_width - ci_width_loo,
    delta_I2 = I2 - I2_loo,
    delta_tau2 = tau2 - tau2_loo,
    sig_changed = (ci_lower > 0 & ci_lower_loo < 0) | (ci_upper < 0 & ci_upper_loo > 0)
  )

cat("\n=== LEAVE-ONE-OUT SUMMARY STATISTICS ===")
cat("\nRange of pooled estimates: [", round(min(loo_df$pooled_g_loo), 3), ",",
    round(max(loo_df$pooled_g_loo), 3), "]")
cat("\nFull model pooled estimate:", round(pooled_g, 3))
cat("\nMaximum change in estimate(Δμ̂): ", round(max(loo_df$delta_g_abs), 3))
cat("\nMaximum standardized change: ", round(max(abs(loo_df$delta_g_standardized)), 2), "SE")
cat("\n\nRange of I²: [", round(min(loo_df$I2_loo), 1), "%-",
    round(max(loo_df$I2_loo), 1), "%]")
cat("\nFull model I²:", round(I2, 1), "%")
cat("\nMaximum change in I² (ΔI²): ", round(max(abs(loo_df$delta_I2)), 1), "%")

# Identify most influential study
most_influential_idx <- which.max(loo_df$delta_g_abs)
most_influential <- loo_df[most_influential_idx, ]

cat("\n\n=== MOST INFLUENTIAL STUDY ===")
cat("\nStudy:", as.character(most_influential$study_removed))
cat("\nOriginal effect size: g =", round(most_influential$original_effect, 3))
cat("\nWeight in full model:", round(most_influential$original_weight, 1), "%")
cat("\nPooled estimate with study: g =", round(pooled_g, 3))
cat("\nPooled estimate without study: g =", round(most_influential$pooled_g_loo, 3))
cat("\nChange in estimate(Δμ̂):", round(most_influential$delta_g, 3))
cat("\nStandardized change:", round(most_influential$delta_g_standardized, 2), "SE")
cat("\nChange in I²:", round(most_influential$delta_I2, 1), "%")

# Interpret influence magnitude
if (max(loo_df$delta_g_abs) < 0.10 * se_pooled) {
  influence_interp <- "negligible(Δμ̂ < 10% SE) - ROBUST"
} else if (max(loo_df$delta_g_abs) < 0.20 * se_pooled) {
  influence_interp <- "small(Δμ̂ < 20% SE) - Generally robust"
} else if (max(loo_df$delta_g_abs) < se_pooled) {
  influence_interp <- "moderate(Δμ̂ < 1 SE) - Some influence, investigate"
} else {
  influence_interp <- "large(Δμ̂ ≥ 1 SE) - Highly influential, warrants investigation"
}

cat("\n\nOverall influence interpretation:", influence_interp)

# Check for significance instability
if (any(loo_df$sig_changed)) {
  cat("\n\nWARNING: Statistical significance changed for some iterations!")
  sig_changed_studies <- loo_df$study_removed[loo_df$sig_changed]
  cat("\nFragile significance when removing:", paste(sig_changed_studies, collapse=", "))
} else {
  cat("\n\nStatistical significance stable across all leave-one-out iterations.")
}

# === STEP 4: Visualize Leave-One-Out Results ===

# Plot 1: Pooled estimates with CIs
par(mfrow=c(2,2), mar=c(4,4,3,2))

# Effect size influence
plot(1:k, loo_df$pooled_g_loo,
     ylim = range(c(loo_df$ci_lower_loo, loo_df$ci_upper_loo, pooled_g)),
     xlab = "Study Removed(Number)", 
     ylab = "Pooled Hedges' g",
     main = "Leave-One-Out: Pooled Effect Estimates",
     pch = 19, col = "steelblue", cex=1.2)

# Add CIs
segments(1:k, loo_df$ci_lower_loo, 1:k, loo_df$ci_upper_loo, col="steelblue")

# Add full model estimate
abline(h = pooled_g, col="red", lwd=2, lty=2)
abline(h = ci_lower, col="red", lwd=1, lty=3)
abline(h = ci_upper, col="red", lwd=1, lty=3)
abline(h = 0, col="gray", lty=2)

legend("topright", c("Full model", "95% CI"), 
       col=c("red", "red"), lty=c(2,3), lwd=c(2,1), cex=0.7)

# Highlight most influential
points(most_influential_idx, most_influential$pooled_g_loo, 
       pch=8, col="darkred", cex=2)

# Plot 2: Change in estimate
plot(1:k, loo_df$delta_g,
     xlab = "Study Removed(Number)",
     ylab = "Change in Pooled g(Δμ̂)",
     main = "Influence: Change in Estimate",
     pch = 19, col = "darkgreen", cex=1.2)
abline(h = 0, col="gray", lwd=2, lty=2)

# Reference lines at ±10% SE and ±20% SE
abline(h = c(-0.20*se_pooled, -0.10*se_pooled, 0.10*se_pooled, 0.20*se_pooled),
       col="orange", lty=3)
text(k*0.8, 0.15*se_pooled, "±10% SE", col="orange", cex=0.7)
text(k*0.8, 0.25*se_pooled, "±20% SE", col="orange", cex=0.7)

# Highlight most influential
points(most_influential_idx, most_influential$delta_g, 
       pch=8, col="darkred", cex=2)

# Plot 3: Heterogeneity (I²) changes
plot(1:k, loo_df$I2_loo,
     xlab = "Study Removed(Number)",
     ylab = "I² (%)",
     main = "Leave-One-Out: Heterogeneity(I²)",
     pch = 19, col = "purple", cex=1.2,
     ylim = c(0, max(loo_df$I2_loo, I2)*1.1))

# Add full model I²
abline(h = I2, col="red", lwd=2, lty=2)
legend("topright", "Full model I²", col="red", lty=2, lwd=2, cex=0.7)

# Highlight most influential
points(most_influential_idx, most_influential$I2_loo, 
       pch=8, col="darkred", cex=2)

# Plot 4: p-value stability
plot(1:k, loo_df$pval_loo,
     xlab = "Study Removed(Number)",
     ylab = "p-value",
     main = "Leave-One-Out: Significance Stability",
     pch = 19, col = "navy", cex=1.2,
     ylim = c(0, max(loo_df$pval_loo, p_value)*1.2))

# Add significance threshold
abline(h = 0.05, col="red", lwd=2, lty=2)
abline(h = p_value, col="blue", lwd=1, lty=3)
legend("topright", c("α = .05", "Full model p"), 
       col=c("red", "blue"), lty=c(2,3), lwd=c(2,1), cex=0.7)

par(mfrow=c(1,1))

# === STEP 5: Forest Plot with Leave-One-Out Results ===
cat("\n\n=== GENERATING FOREST PLOT WITH LEAVE-ONE-OUT ===")

# Create forest plot showing primary analysis + leave-one-out summaries
par(mar=c(5,4,3,2))

# Forest plot of primary meta-analysis
forest(re_model,
       xlab = "Hedges' g(CBT - Control)",
       header = c("Study", "g [95% CI]"),
       cex = 0.75,
       col = "blue")

mtext(paste0("Primary Random-Effects Model: g = ", round(pooled_g, 2),
             ", 95% CI [", round(ci_lower, 2), ", ", round(ci_upper, 2), "]"),
      side=3, line=1, cex=0.85, font=2)
mtext(paste0("Leave-One-Out Range: [", round(min(loo_df$pooled_g_loo), 2),
             ", ", round(max(loo_df$pooled_g_loo), 2), "] | ",
             "Max Δμ̂ = ", round(max(loo_df$delta_g_abs), 3)),
      side=3, line=0.2, cex=0.75)

# === STEP 6: Influence Diagnostics (Cook's Distance, DFBETAS) ===
cat("\n=== ADVANCED INFLUENCE DIAGNOSTICS ===")

# metafor::influence() provides comprehensive diagnostics
inf <- influence(re_model)

cat("\nCook's Distance(influence metric):")
print(round(inf$inf$cook.d, 4))

# Identify studies with high Cook's distance (threshold: 4/k)
cooks_threshold <- 4 / k
high_cooks <- which(inf$inf$cook.d > cooks_threshold)

if (length(high_cooks) > 0) {
  cat("\nStudies exceeding Cook's distance threshold(4/k =", 
      round(cooks_threshold, 3), "):")
  cat("\n", paste(meta_data$author_year[high_cooks], collapse=", "))
} else {
  cat("\nNo studies exceed Cook's distance threshold(4/k =", 
      round(cooks_threshold, 3), ")")
}

# DFBETAS (standardized change in coefficient)
cat("\n\nDFBETAS(standardized influence):")
print(round(inf$inf$dfbs, 3))

# Studentized residuals (outlier detection)
cat("\n\nStudentized Residuals(outlier detection):")
rstudent_vals <- rstudent(re_model)
print(round(rstudent_vals$z, 3))

outliers <- which(abs(rstudent_vals$z) > 2.5)
if (length(outliers) > 0) {
  cat("\nPotential outliers(|z| > 2.5):")
  cat("\n", paste(meta_data$author_year[outliers], collapse=", "))
} else {
  cat("\nNo outliers detected(all |z| ≤ 2.5)")
}

# === STEP 7: Baujat Plot (Heterogeneity Contribution vs. Influence) ===
cat("\n\n=== BAUJAT PLOT: Heterogeneity Contribution vs. Influence ===")

par(mar=c(5,5,3,2))
baujat(re_model, 
       symbol = "slab",
       xlab = "Contribution to Q(Heterogeneity)",
       ylab = "Influence on Pooled Estimate",
       main = "Baujat Plot: Identifying Influential Studies")

# === STEP 8: Detailed Leave-One-Out Table ===
cat("\n\n=== DETAILED LEAVE-ONE-OUT TABLE ===")

# Create publication-ready table
loo_table <- loo_df %>%
  select(study_removed, pooled_g_loo, ci_lower_loo, ci_upper_loo,
         pval_loo, I2_loo, tau2_loo, delta_g, delta_I2) %>%
  mutate(
    CI_95 = paste0("[", round(ci_lower_loo, 2), ", ", round(ci_upper_loo, 2), "]"),
    pooled_g_loo = round(pooled_g_loo, 3),
    pval_loo = format.pval(pval_loo, digits=3, eps=0.001),
    I2_loo = paste0(round(I2_loo, 1), "%"),
    tau2_loo = round(tau2_loo, 4),
    delta_g = round(delta_g, 3),
    delta_I2 = paste0(round(delta_I2, 1), "%")
  ) %>%
  select(study_removed, pooled_g_loo, CI_95, pval_loo, 
         I2_loo, delta_g, delta_I2)

colnames(loo_table) <- c("Study Removed", "Pooled g", "95% CI", 
                          "p-value", "I²", "Δμ̂", "ΔI²")

print(loo_table, row.names=FALSE)

# === STEP 9: Interpretation and Recommendations ===
cat("\n\n=== INTERPRETATION AND RECOMMENDATIONS ===")

cat("\n\n1. ROBUSTNESS ASSESSMENT:")
if (max(loo_df$delta_g_abs) < 0.10 * se_pooled) {
  cat("\n   ✓ Pooled estimate is ROBUST. No individual study substantially affects results.")
  cat("\n   ✓ Maximum change(Δμ̂ =", round(max(loo_df$delta_g_abs), 3), 
      ") is < 10% of SE(0.10 × ", round(se_pooled, 3), " =", round(0.10*se_pooled, 3), ").")
  cat("\n   → Conclusions are trustworthy and not dependent on single studies.")
} else if (max(loo_df$delta_g_abs) < 0.20 * se_pooled) {
  cat("\n   ✓ Pooled estimate is generally ROBUST with minor influence.")
  cat("\n   • Maximum change(Δμ̂ =", round(max(loo_df$delta_g_abs), 3),
      ") is 10-20% of SE.")
  cat("\n   • Study '", as.character(most_influential$study_removed), 
      "' shows moderate influence but does not change conclusions.")
  cat("\n   → Investigate characteristics of influential study(sample size, effect magnitude).")
} else {
  cat("\n   ⚠ Pooled estimate shows SUBSTANTIAL INFLUENCE from individual studies.")
  cat("\n   • Maximum change(Δμ̂ =", round(max(loo_df$delta_g_abs), 3),
      ") exceeds 20% of SE.")
  cat("\n   • Study '", as.character(most_influential$study_removed),
      "' is highly influential(Δμ̂ =", round(most_influential$delta_g, 3), ").")
  cat("\n   → INVESTIGATE: Is study an outlier? High quality? Distinct population?")
  cat("\n   → Report results both WITH and WITHOUT influential study for transparency.")
}

cat("\n\n2. HETEROGENEITY STABILITY:")
if (max(abs(loo_df$delta_I2)) < 10) {
  cat("\n   ✓ Heterogeneity estimates stable(max ΔI² =", 
      round(max(abs(loo_df$delta_I2)), 1), "%).")
  cat("\n   → No single study disproportionately contributes to heterogeneity.")
} else if (max(abs(loo_df$delta_I2)) < 20) {
  cat("\n   • Heterogeneity moderately affected by some studies(max ΔI² =",
      round(max(abs(loo_df$delta_I2)), 1), "%).")
  cat("\n   → Some studies contribute more to heterogeneity; consider subgroup analysis.")
} else {
  cat("\n   ⚠ Heterogeneity substantially affected(max ΔI² =",
      round(max(abs(loo_df$delta_I2)), 1), "%).")
  most_hetero_idx <- which.max(abs(loo_df$delta_I2))
  cat("\n   • Removing '", as.character(loo_df$study_removed[most_hetero_idx]),
      "' changes I² by", round(loo_df$delta_I2[most_hetero_idx], 1), "%.")
  cat("\n   → Study may represent distinct population or be outlier; investigate.")
}

cat("\n\n3. SIGNIFICANCE STABILITY:")
if (!any(loo_df$sig_changed)) {
  cat("\n   ✓ Statistical significance is STABLE across all leave-one-out iterations.")
  if (ci_lower > 0) {
    cat("\n   ✓ Pooled effect remains significant even when removing any single study.")
  } else {
    cat("\n   • Pooled effect remains non-significant even when removing any single study.")
  }
  cat("\n   → Conclusion does not depend critically on individual studies(robust).")
} else {
  cat("\n   ⚠ Statistical significance is FRAGILE.")
  cat("\n   • Removing", paste(loo_df$study_removed[loo_df$sig_changed], collapse=" or "),
      "changes significance.")
  cat("\n   → Evidence is not robust; conclusions should be tempered.")
  cat("\n   → Additional studies needed to establish effect more definitively.")
}

cat("\n\n4. RECOMMENDATIONS:")
cat("\n   a) Investigate influential studies:")
cat("\n      - Assess study quality(risk of bias) independently")
cat("\n      - Examine why influential(large n? extreme effect? unique population?)")
cat("\n      - DO NOT automatically exclude—influential ≠ biased")

if (length(outliers) > 0) {
  cat("\n   b) Outliers detected:")
  cat("\n      - Studies:", paste(meta_data$author_year[outliers], collapse=", "))
  cat("\n      - Assess whether outliers represent true population variation or data errors")
  cat("\n      - Consider subgroup analysis or meta-regression to explore differences")
}

if (max(abs(loo_df$delta_I2)) > 15) {
  cat("\n   c) Heterogeneity investigation:")
  cat("\n      - Conduct subgroup analysis or meta-regression")
  cat("\n      - Identify moderators(population, intervention characteristics, study quality)")
}

cat("\n   d) Transparent reporting:")
cat("\n      - Report leave-one-out range of pooled estimates")
cat("\n      - Identify most influential study and investigate characteristics")
cat("\n      - If highly influential study found, report with/without analyses")

# === STEP 10: APA-Style Reporting ===
cat("\n\n=== APA-STYLE REPORT ===")

report <- paste0(
  "A leave-one-out sensitivity analysis was conducted to assess robustness of the pooled ",
  "effect estimate. Systematically removing each of the ", k, " studies one at a time and ",
  "re-estimating the pooled effect revealed that estimates ranged from g = ",
  round(min(loo_df$pooled_g_loo), 2), " to g = ", round(max(loo_df$pooled_g_loo), 2),
  " (full model: g = ", round(pooled_g, 2), "). The maximum change in the pooled estimate ",
  "was Δμ̂ = ", round(max(loo_df$delta_g_abs), 3), " (",
  round(max(loo_df$delta_g_abs)/se_pooled * 100, 0), "% of SE), observed when removing ",
  as.character(most_influential$study_removed), ".\n\n"
)

if (max(loo_df$delta_g_abs) < 0.10 * se_pooled) {
  report <- paste0(report,
    "This small change(< 10% of SE) indicates that no single study disproportionately ",
    "influenced the pooled effect, supporting robustness of the meta-analytic conclusion. "
  )
} else if (max(loo_df$delta_g_abs) < 0.20 * se_pooled) {
  report <- paste0(report,
    "This moderate change(10-20% of SE) suggests some influence from ",
    as.character(most_influential$study_removed), ", but the pooled effect remained ",
    ifelse(most_influential$ci_lower_loo > 0, "significant and ", ""),
    "substantively similar(g = ", round(most_influential$pooled_g_loo, 2), "). ",
    "Investigation revealed this study had ",
    ifelse(most_influential$original_weight > median(loo_df$original_weight),
           "high precision(large sample)",
           "an extreme effect size"),
    ", explaining its influence. "
  )
} else {
  report <- paste0(report,
    "This substantial change(> 20% of SE) indicates that ",
    as.character(most_influential$study_removed), " is highly influential. ",
    "Removing this study yielded g = ", round(most_influential$pooled_g_loo, 2),
    ", 95% CI [", round(most_influential$ci_lower_loo, 2), ", ",
    round(most_influential$ci_upper_loo, 2), "]. Investigation of study characteristics ",
    "is warranted to understand whether influence reflects legitimate precision(large sample), ",
    "outlier status(extreme effect), or distinct population. Results are reported both ",
    "with and without this study for transparency. "
  )
}

if (!any(loo_df$sig_changed)) {
  report <- paste0(report,
    "Statistical significance remained stable across all leave-one-out iterations, with ",
    "confidence intervals consistently excluding zero. "
  )
} else {
  report <- paste0(report,
    "Notably, statistical significance was fragile: removing ",
    paste(loo_df$study_removed[loo_df$sig_changed], collapse=" or "),
    " rendered the pooled effect non-significant(CI crossing zero). This indicates that ",
    "conclusions are not robust to individual studies, suggesting need for additional ",
    "evidence before making strong claims. "
  )
}

report <- paste0(report,
  "Heterogeneity(I²) ranged from ", round(min(loo_df$I2_loo), 1), "% to ",
  round(max(loo_df$I2_loo), 1), "% across leave-one-out iterations(full model: ",
  round(I2, 1), "%). "
)

if (max(abs(loo_df$delta_I2)) > 20) {
  most_hetero_idx <- which.max(abs(loo_df$delta_I2))
  report <- paste0(report,
    "Removing ", as.character(loo_df$study_removed[most_hetero_idx]),
    " reduced I² by ", abs(round(loo_df$delta_I2[most_hetero_idx], 1)),
    "%, indicating this study contributed disproportionately to heterogeneity, ",
    "possibly representing a distinct population warranting subgroup analysis.\n\n"
  )
} else {
  report <- paste0(report,
    "No single study disproportionately contributed to heterogeneity(max ΔI² = ",
    round(max(abs(loo_df$delta_I2)), 1), "%).\n\n"
  )
}

report <- paste0(report,
  "Conclusion: Leave-one-out sensitivity analysis ",
  ifelse(max(loo_df$delta_g_abs) < 0.10 * se_pooled && !any(loo_df$sig_changed),
         "demonstrated robust pooled effect estimates, with no individual study ",
         "identified influential studies that warrant investigation. While the pooled effect "),
  ifelse(max(loo_df$delta_g_abs) < 0.10 * se_pooled && !any(loo_df$sig_changed),
         "critically affecting conclusions. The meta-analytic finding is trustworthy.",
         "remains meaningful, transparency about influence patterns is essential for interpretation.")
)

cat("\n", report, "\n")