Egger's Regression Test

# Egger's Regression Test for Publication Bias
# Antidepressant vs. Placebo Meta-Analysis Example

library(metafor)      # For meta-analysis and regtest()
library(meta)         # For funnel plot enhancements
library(dplyr)
library(ggplot2)

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

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

# Simulate publication bias scenario:
# True mean effect θ = 0.40 (moderate antidepressant effect)
# Small studies with null/negative results are missing (publication bias)

# Generate true effects with moderate heterogeneity
true_mean <- 0.40
tau <- 0.18  # Between-study SD
true_effects <- rnorm(k, mean=true_mean, sd=tau)

# Sample sizes: vary considerably (realistic for pharma trials)
n_treat <- c(sample(30:60, 10, replace=TRUE),   # Small studies
              sample(60:120, 10, replace=TRUE),  # Medium studies
              sample(120:250, 5, replace=TRUE))  # Large studies
n_control <- c(sample(30:60, 10, replace=TRUE),
               sample(60:120, 10, replace=TRUE),
               sample(120:250, 5, replace=TRUE))

# Sampling standard errors (larger for small studies)
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

# SIMULATE PUBLICATION BIAS:
# Remove small studies with null/negative results (g < 0.2)
# Probability of publication decreases with smaller effects and larger SE
pub_prob <- plogis(2 * observed_g - 3 * sampling_se + 0.5)
published <- rbinom(k, 1, prob=pub_prob) == 1

# Create published sample (biased)
meta_data_published <- data.frame(
  study_id = paste0("Study_", which(published)),
  author_year = paste0(LETTERS[which(published)], " et al.(20", 
                       sprintf("%02d", 10:24)[which(published)], ")"),
  hedges_g = observed_g[published],
  variance = variance_g[published],
  se = sqrt(variance_g[published]),
  n_treatment = n_treat[published],
  n_control = n_control[published],
  total_n = (n_treat + n_control)[published]
)

k_published <- nrow(meta_data_published)

print("=== Published Studies Dataset(After Publication Bias) ===")
print(meta_data_published)
cat("\nPublished studies: k =", k_published, "(out of", k, "conducted)\n")
cat("Total N =", sum(meta_data_published$total_n), "participants\n")

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

print("\n=== Random-Effects Meta-Analysis Results ===")
print(re_model)

pooled_g <- as.numeric(re_model$beta)
ci_lower <- re_model$ci.lb
ci_upper <- re_model$ci.ub
p_value <- re_model$pval

I2 <- re_model$I2
tau2 <- re_model$tau2
Q <- re_model$QE
Q_pval <- re_model$QEp

cat("\n=== Pooled Effect(Potentially Biased) ===")
cat("\nHedges' g =", round(pooled_g, 3))
cat("\n95% CI: [", round(ci_lower, 3), ",", round(ci_upper, 3), "]")
cat("\np-value:", format.pval(p_value, digits=3))
cat("\n\nHeterogeneity: I² =", round(I2, 1), "%, τ² =", round(tau2, 4))

# === STEP 3: Funnel Plot (Visual Inspection) ===
par(mfrow=c(1,2), mar=c(5,4,3,2))

# Standard funnel plot
funnel(re_model, 
       xlab = "Hedges' g",
       ylab = "Standard Error",
       main = "Funnel Plot",
       back = "white",
       shade = "white")

# Add reference line
abline(v = pooled_g, col="red", lwd=2, lty=2)

# Contour-enhanced funnel plot (distinguishes bias from heterogeneity)
funnel(re_model,
       xlab = "Hedges' g",
       ylab = "Standard Error", 
       main = "Contour-Enhanced Funnel Plot",
       back = "white",
       shade = c("white", "lightgray", "darkgray"),
       level = c(0.10, 0.05, 0.01))
legend("topright", c("p > .10", ".05 < p < .10", ".01 < p < .05", "p < .01"),
       fill = c("white", "lightgray", "darkgray", "darkgray"),
       cex = 0.7)

par(mfrow=c(1,1))

cat("\n\n=== Funnel Plot Visual Assessment ===")
cat("\nVisual inspection: Look for asymmetry(missing studies in bottom-right/left)")
cat("\nContour-enhanced plot: If missing studies cluster in non-significant")
cat("\nregions(white area, p>.10), suggests publication bias.")
cat("\nIf missing studies spread across significance contours, suggests heterogeneity.\n")

# === STEP 4: Egger's Regression Test ===
# Regress standardized effect (effect/SE) on precision (1/SE)
# H₀: Intercept = 0 (no asymmetry)
# Hₐ: Intercept ≠ 0 (asymmetry present)

egger_test <- regtest(re_model, model="lm", predictor="sei")
# Note: predictor="sei" uses standard error as predictor (equivalent to precision in regression)
# This is the standard Egger's test

print("\n=== EGGER'S REGRESSION TEST ===")
print(egger_test)

egger_intercept <- egger_test$est
egger_se <- egger_test$se
egger_z <- egger_test$zval  # Actually t-statistic for intercept
egger_p <- egger_test$pval
ci_egger_lower <- egger_intercept - 1.96 * egger_se
ci_egger_upper <- egger_intercept + 1.96 * egger_se

cat("\n=== Egger's Test Interpretation ===")
cat("\nRegression Intercept(β₀) =", round(egger_intercept, 3))
cat("\nStandard Error =", round(egger_se, 3))
cat("\n95% CI for intercept: [", round(ci_egger_lower, 3), ",", 
    round(ci_egger_upper, 3), "]")
cat("\nt-statistic =", round(egger_z, 3))
cat("\np-value =", round(egger_p, 4))

# Interpret using α = 0.10 threshold (recommended)
cat("\n\n=== INTERPRETATION(α = 0.10 threshold) ===")
if (egger_p < 0.10) {
  cat("\n✓ SIGNIFICANT asymmetry detected(p < .10)")
  cat("\n→ Funnel plot shows significant asymmetry")
  cat("\n→ Possible publication bias or small-study effects")
  
  if (egger_intercept > 0) {
    cat("\n→ Positive intercept: Small studies show LARGER effects")
    cat("\n(Typical publication bias pattern: small null studies missing)")
  } else {
    cat("\n→ Negative intercept: Small studies show SMALLER effects")
    cat("\n(Unusual pattern; investigate further)")
  }
  
  cat("\n\nConclusion: Evidence suggests potential publication bias.")
  cat("\nPooled effect estimate may be overestimated.")
  cat("\nConduct bias-correction analyses(trim-and-fill, PET-PEESE).")
  
} else {
  cat("\n✗ No significant asymmetry detected(p ≥ .10)")
  cat("\n→ Limited statistical evidence of funnel plot asymmetry")
  cat("\n→ However, this does NOT prove absence of publication bias")
  cat("\n(Test may be underpowered with k =", k_published, ")")
  
  cat("\n\nConclusion: No significant evidence of asymmetry, but absence")
  cat("\nof evidence is not evidence of absence. Bias may still exist.")
}

# Power consideration
cat("\n\n=== POWER CONSIDERATION ===")
if (k_published < 10) {
  cat("\nWARNING: k < 10 studies. Egger's test has VERY LOW POWER(<20%).")
  cat("\nTest is unreliable with this sample size. Do not rely on this result.")
} else if (k_published < 15) {
  cat("\nCAUTION: k < 15 studies. Egger's test has MODEST POWER(~40-60%).")
  cat("\nInterpret with caution. Non-significant result may reflect low power.")
} else {
  cat("\nAdequate sample size(k ≥ 15) for Egger's test.")
  cat("\nTest has reasonable power(~70-80%) to detect moderate bias.")
}

# === STEP 5: Begg's Rank Correlation Test (Non-Parametric Sensitivity) ===
# Less powerful than Egger's but more robust to outliers

begg_test <- ranktest(re_model)

print("\n\n=== BEGG'S RANK CORRELATION TEST(Sensitivity Check) ===")
print(begg_test)

begg_tau <- begg_test$tau
begg_p <- begg_test$pval

cat("\n=== Begg's Test Interpretation ===")
cat("\nKendall's tau =", round(begg_tau, 3))
cat("\np-value =", round(begg_p, 4))

if (begg_p < 0.10) {
  cat("\n✓ Significant rank correlation(p < .10)")
  cat("\n→ Corroborates Egger's test: asymmetry detected")
} else {
  cat("\n✗ No significant rank correlation(p ≥ .10)")
  cat("\n→ Note: Begg's test has lower power than Egger's test")
}

cat("\n\nComparison: Egger's p =", round(egger_p, 3), 
    ", Begg's p =", round(begg_p, 3))
if (egger_p < 0.10 & begg_p >= 0.10) {
  cat("\n→ Egger's significant but Begg's not: Egger's more powerful for continuous outcomes")
} else if (egger_p >= 0.10 & begg_p < 0.10) {
  cat("\n→ Begg's significant but Egger's not: Unusual; check for outliers affecting Egger's")
} else if (egger_p < 0.10 & begg_p < 0.10) {
  cat("\n→ Both significant: Strong convergent evidence of asymmetry")
} else {
  cat("\n→ Both non-significant: Limited evidence of asymmetry from either test")
}

# === STEP 6: Trim-and-Fill Analysis (Estimate Missing Studies) ===
# Imputes potentially missing studies and adjusts pooled effect

taf <- trimfill(re_model)

print("\n\n=== TRIM-AND-FILL ANALYSIS ===")
print(taf)

k_imputed <- taf$k0
g_adjusted <- as.numeric(taf$beta)
ci_adj_lower <- taf$ci.lb
ci_adj_upper <- taf$ci.ub

cat("\n=== Trim-and-Fill Results ===")
cat("\nImputed missing studies(k₀) =", k_imputed)
cat("\nAdjusted Hedges' g =", round(g_adjusted, 3))
cat("\nAdjusted 95% CI: [", round(ci_adj_lower, 3), ",", 
    round(ci_adj_upper, 3), "]")

cat("\n\n=== Bias Impact Assessment ===")
cat("\nUnadjusted estimate: g =", round(pooled_g, 3))
cat("\nAdjusted estimate:   g =", round(g_adjusted, 3))
cat("\nDifference:          Δg =", round(pooled_g - g_adjusted, 3))
pct_change <- ((pooled_g - g_adjusted) / pooled_g) * 100
cat("\nPercent change:      ", round(abs(pct_change), 1), "%")

if (k_imputed == 0) {
  cat("\n\nInterpretation: No missing studies imputed.")
  cat("\nTrim-and-fill suggests minimal bias(or bias on left side of funnel).")
} else if (abs(pct_change) < 10) {
  cat("\n\nInterpretation:", k_imputed, "missing studies imputed.")
  cat("\nAdjusted estimate differs by <10% from unadjusted.")
  cat("\n→ Modest impact of potential bias; conclusions relatively robust.")
} else if (abs(pct_change) < 25) {
  cat("\n\nInterpretation:", k_imputed, "missing studies imputed.")
  cat("\nAdjusted estimate differs by", round(abs(pct_change), 1), "% from unadjusted.")
  cat("\n→ Moderate impact of potential bias; interpret with caution.")
} else {
  cat("\n\nInterpretation:", k_imputed, "missing studies imputed.")
  cat("\nAdjusted estimate differs by", round(abs(pct_change), 1), "% from unadjusted.")
  cat("\n→ Substantial impact of potential bias; conclusions may be fragile.")
  cat("\n→ Pooled effect may be considerably overestimated.")
}

# Check if adjusted estimate still significant
if (ci_adj_lower > 0) {
  cat("\n→ Adjusted CI still excludes zero: Effect remains statistically significant.")
} else if (ci_adj_upper < 0) {
  cat("\n→ Adjusted CI still excludes zero(negative): Harmful effect remains significant.")
} else {
  cat("\n→ Adjusted CI includes zero: Effect no longer statistically significant.")
  cat("\n  WARNING: Bias-correction eliminates significance; findings may be spurious.")
}

# Funnel plot with trim-and-fill imputed studies
par(mar=c(5,4,3,2))
funnel(taf,
       xlab = "Hedges' g",
       ylab = "Standard Error",
       main = paste0("Trim-and-Fill: ", k_imputed, " Imputed Studies"),
       back = "white",
       col = c("blue", "red"),
       pch = c(19, 17))
legend("topright", c("Observed studies", "Imputed studies"),
       col = c("blue", "red"), pch = c(19, 17), cex=0.9)

# === STEP 7: Influence Diagnostics for Egger's Regression ===
# Check if outliers distort Egger's test

cat("\n\n=== INFLUENCE DIAGNOSTICS FOR EGGER'S REGRESSION ===")

# Create regression data
precision <- 1 / meta_data_published$se
standardized_effect <- meta_data_published$hedges_g / meta_data_published$se

# Fit Egger's regression manually to get diagnostics
egger_lm <- lm(standardized_effect ~ precision)

# Cook's distance
cooks_d <- cooks.distance(egger_lm)
influential <- cooks_d > 4/k_published

cat("\nCook's Distance(identifies influential studies):")
for (i in 1:k_published) {
  cat("\n ", meta_data_published$author_year[i], ": D =", 
      round(cooks_d[i], 3),
      ifelse(influential[i], " [INFLUENTIAL]", ""))
}

if (any(influential)) {
  cat("\n\nWARNING:", sum(influential), "influential study(ies) detected(Cook's D > 4/k).")
  cat("\nThese studies may disproportionately affect Egger's test result.")
  cat("\nConsider leave-one-out sensitivity analysis.")
} else {
  cat("\n\nNo highly influential outliers detected.")
  cat("\nEgger's test result appears robust to individual studies.")
}

# Leave-one-out sensitivity for Egger's test
cat("\n\n=== LEAVE-ONE-OUT SENSITIVITY FOR EGGER'S TEST ===")

loo_egger_p <- numeric(k_published)
loo_egger_intercept <- numeric(k_published)

for (i in 1:k_published) {
  # Remove study i
  loo_model <- rma(yi = hedges_g, vi = variance, 
                    data = meta_data_published[-i,],
                    method = "REML")
  loo_test <- regtest(loo_model, model="lm", predictor="sei")
  loo_egger_p[i] <- loo_test$pval
  loo_egger_intercept[i] <- loo_test$est
}

cat("\nEgger's p-value range(leave-one-out):", 
    round(min(loo_egger_p), 4), "to", round(max(loo_egger_p), 4))
cat("\nFull model p-value:", round(egger_p, 4))

# Check if any single study changes conclusion
if (egger_p < 0.10) {
  # Originally significant
  n_nonsig <- sum(loo_egger_p >= 0.10)
  if (n_nonsig > 0) {
    cat("\n\nWARNING: Removing", n_nonsig, "study(ies) makes Egger's test non-significant.")
    cat("\nEgger's test result is FRAGILE; depends on specific studies included.")
    cat("\nStudies causing shift to non-significance:")
    for (i in which(loo_egger_p >= 0.10)) {
      cat("\n -", meta_data_published$author_year[i], 
          "(p changes from", round(egger_p, 3), "to", round(loo_egger_p[i], 3), ")")
    }
  } else {
    cat("\n\nEgger's test remains significant(p < .10) across all leave-one-out analyses.")
    cat("\nResult is ROBUST to removal of any single study.")
  }
} else {
  # Originally non-significant
  n_sig <- sum(loo_egger_p < 0.10)
  if (n_sig > 0) {
    cat("\n\nNote: Removing", n_sig, "study(ies) makes Egger's test significant.")
    cat("\nThese studies may be suppressing detection of asymmetry.")
  } else {
    cat("\n\nEgger's test remains non-significant across all leave-one-out analyses.")
    cat("\nConsistently no evidence of asymmetry.")
  }
}

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

report <- paste0(
  "Publication bias was assessed using multiple methods. ",
  "Visual inspection of the funnel plot suggested ",
  ifelse(egger_p < 0.10, "asymmetry, with potential missing studies in regions of non-significance. ",
         "approximate symmetry, though formal statistical testing is necessary. "),
  "\n\nEgger's regression test ",
  ifelse(egger_p < 0.10, "detected significant", "did not detect significant"),
  " funnel plot asymmetry(intercept = ", round(egger_intercept, 3),
  ", 95% CI [", round(ci_egger_lower, 3), ", ", round(ci_egger_upper, 3),
  "], p = ", round(egger_p, 3),
  ifelse(egger_p < 0.10, 
         " at the liberal α = .10 threshold recommended for bias detection). The positive intercept indicates small studies showed larger treatment effects than large studies, consistent with possible publication bias.",
         "). Using the recommended liberal α = .10 threshold for bias detection, this result suggests limited statistical evidence of funnel plot asymmetry."),
  "\n\nBegg's rank correlation test(non-parametric sensitivity check) ",
  ifelse(begg_p < 0.10, "also detected", "did not detect"),
  " significant asymmetry(Kendall's tau = ", round(begg_tau, 3),
  ", p = ", round(begg_p, 3), "). ",
  ifelse((egger_p < 0.10 & begg_p < 0.10),
         "The convergence of Egger's and Begg's tests provides stronger evidence of asymmetry. ",
         ifelse((egger_p < 0.10 & begg_p >= 0.10),
                "The discrepancy between tests(Egger's significant, Begg's not) reflects Egger's greater power for continuous outcomes, though Begg's test is more robust to outliers. ",
                "")),
  "\n\nTrim-and-fill analysis estimated ", k_imputed,
  ifelse(k_imputed == 0, " missing studies",
         ifelse(k_imputed == 1, " missing study", " missing studies")),
  ifelse(k_imputed > 0,
         paste0(". Imputing these studies yielded an adjusted pooled effect of g = ",
                round(g_adjusted, 3), " (95% CI [", round(ci_adj_lower, 3), ", ",
                round(ci_adj_upper, 3), "]), compared to the unadjusted estimate of g = ",
                round(pooled_g, 3), " (95% CI [", round(ci_lower, 3), ", ",
                round(ci_upper, 3), "]), representing a ",
                round(abs(pct_change), 1), "% ",
                ifelse(pct_change > 0, "reduction", "increase"), "."),
         paste0(", suggesting that any bias, if present, may favor the null rather than the alternative hypothesis, or that bias is minimal.")),
  ifelse(k_imputed > 0 & ci_adj_lower > 0,
         " Importantly, the adjusted estimate remained statistically significant, suggesting conclusions are relatively robust despite potential bias.",
         ifelse(k_imputed > 0 & ci_adj_upper > 0 & ci_adj_lower <= 0,
                " However, the adjusted confidence interval included zero, indicating that bias-correction eliminated statistical significance. This raises concerns about the robustness of the treatment effect.",
                "")),
  "\n\nInfluence diagnostics revealed ",
  ifelse(any(influential),
         paste0(sum(influential), " influential study(ies) (Cook's D > 4/k) that may disproportionately affect Egger's test. "),
         "no highly influential outliers in Egger's regression. "),
  "Leave-one-out sensitivity analysis showed Egger's p-value ranged from ",
  round(min(loo_egger_p), 3), " to ", round(max(loo_egger_p), 3),
  " when removing each study sequentially, indicating the result is ",
  ifelse((egger_p < 0.10 & all(loo_egger_p < 0.10)) | 
           (egger_p >= 0.10 & all(loo_egger_p >= 0.10)),
         "robust", "somewhat fragile"),
  " to the inclusion of individual studies.",
  "\n\nConclusion: ",
  ifelse(egger_p < 0.10 & k_imputed > 0 & abs(pct_change) >= 25,
         "Evidence suggests possible publication bias, with substantial impact on the pooled effect estimate. The adjusted estimate should be considered alongside the unadjusted estimate, and conclusions should be interpreted with caution. Prioritizing evidence from large, high-quality studies is recommended.",
         ifelse(egger_p < 0.10 & k_imputed > 0 & abs(pct_change) < 25,
                "Evidence suggests possible publication bias, though bias-correction methods indicate modest impact on conclusions. The pooled effect estimate appears relatively robust, but potential bias should be acknowledged.",
                ifelse(egger_p >= 0.10,
                       paste0("Limited statistical evidence of publication bias was detected, though this does not prove absence of bias given ",
                              ifelse(k_published < 15, "modest power with k < 15 studies. ", "available power. "),
                              "Comprehensive search strategies including gray literature and trial registries strengthen confidence in findings."),
                       "Publication bias assessment yielded mixed results requiring careful interpretation.")))
)

cat(report)

cat("\n\n========================================\n")
cat("=== END OF ANALYSIS ===")
cat("\n========================================\n")