Begg's Rank Correlation Test for Funnel Plot Asymmetry

# Begg's Rank Correlation Test for Publication Bias
# Exercise Intervention for Depression Meta-Analysis Example

library(metafor)      # For meta-analysis and ranktest()
library(meta)         # For alternative implementations
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 <- 22  # Number of studies

# Simulate publication bias scenario:
# True mean effect θ = 0.50 (moderate exercise benefit)
# Small studies with null/negative results underrepresented

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

# Sample sizes: heterogeneous (realistic for exercise trials)
# Mix of small pilot studies and larger trials
n_treat <- c(sample(15:35, 12, replace=TRUE),   # Small studies
              sample(35:70, 8, replace=TRUE),    # Medium studies
              sample(70:150, 2, replace=TRUE))   # Large studies
n_control <- c(sample(15:35, 12, replace=TRUE),
               sample(35:70, 8, replace=TRUE),
               sample(70:150, 2, 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:
# Smaller studies with weaker effects less likely published
# Publication probability increases with effect size, decreases with SE
pub_prob <- plogis(1.8 * observed_g - 2.5 * sampling_se + 0.3)
published <- rbinom(k, 1, prob=pub_prob) == 1

# Ensure minimum k ≥ 15 for Begg's test
if (sum(published) < 15) {
  # Add some null studies to reach minimum
  unpublished_indices <- which(!published)
  additional <- sample(unpublished_indices, 15 - sum(published))
  published[additional] <- TRUE
}

# 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", 5:23)[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[, c("author_year", "hedges_g", "se", "total_n")])
cat("\nPublished studies: k =", k_published, "(out of", k, "conducted)\n")
cat("Total N =", sum(meta_data_published$total_n), "participants\n")

# Check precision variability (important for Begg's test)
se_range <- range(meta_data_published$se)
se_cv <- sd(meta_data_published$se) / mean(meta_data_published$se)
cat("\nPrecision variability check:")
cat("\nSE range:", round(se_range[1], 3), "to", round(se_range[2], 3))
cat("\nSE coefficient of variation:", round(se_cv, 3))
if (se_cv < 0.30) {
  cat(" [WARNING: Limited precision heterogeneity may reduce Begg's test power]\n")
} else {
  cat(" [Adequate precision variability for Begg's test]\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))

if (I2 > 75) {
  cat("\n[NOTE: Substantial heterogeneity may confound bias detection]")
}

# === 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(Exercise Effect)",
       ylab = "Standard Error",
       main = "Funnel Plot",
       back = "white",
       shade = "white")
abline(v = pooled_g, col="red", lwd=2, lty=2)

# Funnel plot with effect vs. variance (Begg's test uses variance)
plot(meta_data_published$variance, meta_data_published$hedges_g,
     xlab = "Variance",
     ylab = "Hedges' g",
     main = "Effect Size vs. Variance\n(Begg's Test Relationship)",
     pch = 19, col = "steelblue", cex = 1.5)
abline(h = pooled_g, col="red", lwd=2, lty=2)

# Add lowess smoother to visualize trend
lines(lowess(meta_data_published$variance, meta_data_published$hedges_g),
      col = "darkgreen", lwd = 2)
legend("topright", 
       c("Studies", "Pooled effect", "Lowess smooth"),
       col = c("steelblue", "red", "darkgreen"),
       pch = c(19, NA, NA),
       lty = c(NA, 2, 1),
       lwd = c(NA, 2, 2))

par(mfrow=c(1,1))

cat("\n\n=== Funnel Plot Visual Assessment ===")
cat("\nStandard funnel: Look for asymmetry(missing studies in lower corners)")
cat("\nEffect vs. Variance: Begg's test assesses if larger variances(small studies)")
cat("\nassociate with larger effects. Positive trend suggests publication bias.\n")

# === STEP 4: Begg's Rank Correlation Test ===
# Tests Kendall's tau correlation between effect sizes and variances
# H₀: τ = 0 (no correlation)
# Hₐ: τ ≠ 0 (correlation exists)

begg_test <- ranktest(re_model)
# ranktest() in metafor implements Begg & Mazumdar (1994) test

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

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

# Extract test statistic details
# ranktest returns Kendall's tau and associated p-value

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

# Interpret tau magnitude
if (abs(begg_tau) < 0.2) {
  tau_interp <- "small(weak correlation)"
} else if (abs(begg_tau) < 0.4) {
  tau_interp <- "moderate"
} else {
  tau_interp <- "large(strong correlation)"
}

cat("\nTau magnitude:", tau_interp)

# Interpret direction
if (begg_tau > 0) {
  cat("\nDirection: POSITIVE correlation")
  cat("\n→ Small studies(high variance) tend to show LARGER effects")
  cat("\n(Consistent with typical publication bias pattern)")
} else if (begg_tau < 0) {
  cat("\nDirection: NEGATIVE correlation")
  cat("\n→ Small studies(high variance) tend to show SMALLER effects")
  cat("\n(Unusual pattern; investigate further)")
} else {
  cat("\nDirection: No correlation(τ ≈ 0)")
}

# Interpret using α = 0.10 threshold (recommended)
cat("\n\n=== INTERPRETATION(α = 0.10 threshold) ===")
if (begg_p < 0.10) {
  cat("\n✓ SIGNIFICANT rank correlation detected(p < .10)")
  cat("\n→ Statistically significant association between effect sizes and variances")
  cat("\n→ Possible publication bias or small-study effects")
  
  if (begg_tau > 0) {
    cat("\n→ Positive τ: Small studies show larger effects(typical bias pattern)")
    cat("\n  Small null studies may be missing from published literature")
  }
  
  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 rank correlation detected(p ≥ .10)")
  cat("\n→ Limited statistical evidence of association")
  cat("\n→ However, this does NOT prove absence of publication bias")
  cat("\n(Begg's test has modest power, especially with k =", k_published, ")")
  
  cat("\n\nConclusion: No significant evidence of rank correlation.")
  cat("\nAbsence of evidence is not evidence of absence.")
  cat("\nBias may exist but test lacks power to detect.")
}

# Power consideration
cat("\n\n=== POWER CONSIDERATION ===")
if (k_published < 10) {
  cat("\nWARNING: k < 10 studies. Begg's test has VERY LOW POWER(<15%).")
  cat("\nTest is unreliable. Do not rely on this result.")
} else if (k_published < 15) {
  cat("\nCAUTION: k < 15 studies. Begg's test has LOW TO MODEST POWER(~30-50%).")
  cat("\nBegg's test requires larger k than Egger's for equivalent power.")
  cat("\nInterpret with caution. Non-significant may reflect inadequate power.")
} else if (k_published < 20) {
  cat("\nAdequate but not ideal sample size(k = 15-19).")
  cat("\nBegg's test has moderate power(~60-75%) for moderate bias.")
  cat("\nResults more reliable but still consider power limitations.")
} else {
  cat("\nGood sample size(k ≥ 20) for Begg's test.")
  cat("\nTest has reasonable power(>80%) to detect moderate bias.")
}

# === STEP 5: Egger's Regression Test (Comparison) ===
# More powerful parametric alternative

egger_test <- regtest(re_model, model="lm", predictor="sei")

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

egger_intercept <- egger_test$est
egger_p <- egger_test$pval

cat("\n=== Egger's Test Results ===")
cat("\nIntercept(β₀) =", round(egger_intercept, 3))
cat("\np-value =", round(egger_p, 4))

# Compare Begg's and Egger's
cat("\n\n=== COMPARISON: BEGG'S vs. EGGER'S TEST ===")
cat("\nBegg's test:  τ =", round(begg_tau, 3), ", p =", round(begg_p, 3))
cat("\nEgger's test: β₀ =", round(egger_intercept, 3), ", p =", round(egger_p, 3))

# Convergence assessment
begg_sig <- begg_p < 0.10
egger_sig <- egger_p < 0.10

if (begg_sig & egger_sig) {
  cat("\n\n→ CONVERGENCE: Both tests significant(p < .10)")
  cat("\n  Strong evidence of funnel plot asymmetry")
  cat("\n  Consistent signal across parametric and non-parametric methods")
  cat("\n  HIGH CONFIDENCE in asymmetry detection")
  
} else if (!begg_sig & !egger_sig) {
  cat("\n\n→ CONVERGENCE: Both tests non-significant(p ≥ .10)")
  cat("\n  Limited statistical evidence of asymmetry from either method")
  cat("\n  However, both tests may lack power(especially Begg's)")
  cat("\n  Cannot rule out bias; absence of evidence ≠ evidence of absence")
  
} else if (egger_sig & !begg_sig) {
  cat("\n\n→ DIVERGENCE: Egger's significant, Begg's not significant")
  cat("\n  COMMON pattern: Egger's has higher power than Begg's")
  cat("\n  Egger's detected asymmetry; Begg's may be underpowered")
  cat("\n  INTERPRETATION: Moderate evidence of asymmetry")
  cat("\n  Prioritize Egger's result IF no extreme outliers")
  
} else {  # begg_sig & !egger_sig (RARE)
  cat("\n\n→ DIVERGENCE: Begg's significant, Egger's not significant")
  cat("\n  UNUSUAL pattern: Suggests potential outliers affecting Egger's")
  cat("\n  Begg's(non-parametric) may be more robust here")
  cat("\n  INVESTIGATE: Check for influential outliers in Egger's regression")
}

# Methodological comparison
cat("\n\n=== METHODOLOGICAL COMPARISON ===")
cat("\n┌─────────────────┬───────────────┬───────────────┐")
cat("\n│ Feature         │ Begg's Test   │ Egger's Test  │")
cat("\n├─────────────────┼───────────────┼───────────────┤")
cat("\n│ Method          │ Rank corr.    │ Regression    │")
cat("\n│ Power           │ Lower         │ Higher        │")
cat("\n│ Robustness      │ More robust   │ Less robust   │")
cat("\n│ Outliers        │ Handles well  │ Sensitive     │")
cat("\n│ Min. k          │ 15-20         │ 10-15         │")
cat("\n│ Binary outcomes │ Better        │ Poor          │")
cat("\n│ Result(p-val)  │", sprintf("%-13s", round(begg_p, 3)), "│", sprintf("%-13s", round(egger_p, 3)), "│")
cat("\n└─────────────────┴───────────────┴───────────────┘\n")

# === STEP 6: Outlier Check (Influence on Tests) ===
cat("\n=== OUTLIER ASSESSMENT ===")

# Identify potential outliers
residuals <- meta_data_published$hedges_g - pooled_g
std_residuals <- residuals / meta_data_published$se
outliers <- abs(std_residuals) > 3

cat("\nOutlier detection(|standardized residual| > 3):")
if (any(outliers)) {
  cat("\n", sum(outliers), "potential outlier(s) detected:\n")
  for (i in which(outliers)) {
    cat("  -", meta_data_published$author_year[i], 
        ": g =", round(meta_data_published$hedges_g[i], 3),
        ", std. resid =", round(std_residuals[i], 2), "\n")
  }
  
  # Sensitivity analysis: Begg's test without outliers
  if (sum(!outliers) >= 10) {
    cat("\nSensitivity analysis(excluding outliers):\n")
    
    model_no_outliers <- rma(yi = hedges_g, vi = variance,
                              data = meta_data_published[!outliers,],
                              method = "REML")
    begg_no_outliers <- ranktest(model_no_outliers)
    egger_no_outliers <- regtest(model_no_outliers, model="lm", predictor="sei")
    
    cat("Begg's test(no outliers):  τ =", round(begg_no_outliers$tau, 3),
        ", p =", round(begg_no_outliers$pval, 3), "\n")
    cat("Egger's test(no outliers): β₀ =", round(egger_no_outliers$est, 3),
        ", p =", round(egger_no_outliers$pval, 3), "\n")
    
    # Compare with full sample
    cat("\nComparison:")
    cat("\nBegg's p-value:  ", round(begg_p, 3), "(full) vs.",
        round(begg_no_outliers$pval, 3), "(no outliers)")
    cat("\nEgger's p-value: ", round(egger_p, 3), "(full) vs.",
        round(egger_no_outliers$pval, 3), "(no outliers)")
    
    # Assess sensitivity
    begg_changes <- (begg_p < 0.10) != (begg_no_outliers$pval < 0.10)
    egger_changes <- (egger_p < 0.10) != (egger_no_outliers$pval < 0.10)
    
    if (egger_changes & !begg_changes) {
      cat("\n\n→ Egger's test result CHANGES with outlier removal")
      cat("\n  Begg's test result STABLE(more robust to outliers)")
      cat("\n  Supports using Begg's result when outliers present")
    } else if (begg_changes & !egger_changes) {
      cat("\n\n→ Begg's test result changes(unusual)")
      cat("\n  Extreme outliers may affect rank ordering")
    } else if (begg_changes & egger_changes) {
      cat("\n\n→ Both tests change with outlier removal")
      cat("\n  Results FRAGILE; conclusions depend on outliers")
      cat("\n  Report both analyses(with/without outliers)")
    } else {
      cat("\n\n→ Both tests ROBUST to outliers")
      cat("\n  Conclusions do not depend on extreme studies")
    }
  }
  
} else {
  cat("\nNo extreme outliers detected(all |std. residuals| < 3)")
  cat("\nBegg's and Egger's test results not confounded by outliers.")
}

# === STEP 7: Trim-and-Fill Analysis ===
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=== Bias Impact Assessment ===")
cat("\nImputed missing studies(k₀) =", k_imputed)
cat("\nUnadjusted estimate: g =", round(pooled_g, 3), 
    "[95% CI:", round(ci_lower, 3), ",", round(ci_upper, 3), "]")
cat("\nAdjusted estimate:   g =", round(g_adjusted, 3),
    "[95% CI:", round(ci_adj_lower, 3), ",", round(ci_adj_upper, 3), "]")

if (k_imputed > 0) {
  diff <- pooled_g - g_adjusted
  pct_change <- (diff / pooled_g) * 100
  
  cat("\nDifference:          Δg =", round(diff, 3))
  cat("\nPercent change:      ", round(abs(pct_change), 1), "%")
  
  if (abs(pct_change) < 10) {
    cat("\n\nInterpretation: Modest bias impact(<10% change)")
    cat("\n→ Conclusions relatively robust despite potential bias")
  } else if (abs(pct_change) < 25) {
    cat("\n\nInterpretation: Moderate bias impact(", round(abs(pct_change), 1), "% change)")
    cat("\n→ Non-trivial overestimation; interpret with caution")
  } else {
    cat("\n\nInterpretation: Substantial bias impact(", round(abs(pct_change), 1), "% change)")
    cat("\n→ Major overestimation; conclusions may be fragile")
  }
  
  # Check if still significant
  if (ci_adj_lower > 0) {
    cat("\n→ Adjusted CI excludes zero: Effect remains significant")
  } else {
    cat("\n→ WARNING: Adjusted CI includes zero")
    cat("\n  Bias-correction eliminates statistical significance")
  }
  
} else {
  cat("\n\nInterpretation: No missing studies imputed")
  cat("\nTrim-and-fill suggests minimal bias(or bias on opposite side)")
}

# === 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 complementary methods. ",
  "Visual inspection of the funnel plot suggested ",
  ifelse(begg_p < 0.10 | egger_p < 0.10, 
         "potential asymmetry, with possible missing studies in regions of non-significance. ",
         "approximate symmetry, though formal statistical testing yielded mixed results. "),
  "\n\nBegg's rank correlation test(non-parametric) ",
  ifelse(begg_p < 0.10, "detected significant", "did not detect significant"),
  " correlation between effect sizes and variances(Kendall's τ = ", 
  round(begg_tau, 3), ", p = ", round(begg_p, 3),
  " at the liberal α = .10 threshold recommended for bias detection). ",
  ifelse(begg_tau > 0 & begg_p < 0.10,
         "The positive correlation indicates small studies tended to show larger treatment effects, consistent with possible publication bias where small null studies remain unpublished.",
         ifelse(begg_p >= 0.10,
                paste0("However, Begg's test has modest power, particularly with k = ", k_published, 
                       ", so this result does not definitively rule out publication bias."),
                "The correlation pattern warrants further investigation.")),
  "\n\nEgger's regression test(parametric, higher power) ",
  ifelse(egger_p < 0.10, "also detected", "did not detect"),
  " significant asymmetry(intercept = ", round(egger_intercept, 3),
  ", p = ", round(egger_p, 3), "). ",
  ifelse(begg_sig == egger_sig,
         "The convergence between Begg's and Egger's tests strengthens confidence in the bias assessment conclusion. ",
         ifelse(egger_sig & !begg_sig,
                "The discrepancy(Egger's significant, Begg's not) reflects Egger's higher statistical power. Begg's non-significant result likely reflects inadequate power rather than true absence of bias. ",
                ifelse(!egger_sig & begg_sig,
                       "The discrepancy(Begg's significant, Egger's not) suggests potential outliers may affect Egger's parametric regression. Begg's robust non-parametric result may be more reliable here. ",
                       "Both tests suggest limited evidence of asymmetry. "))),
  "\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 ", ifelse(k_imputed == 1, "this study", "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 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 and raising concerns about effect robustness.",
                "")),
  ifelse(any(outliers),
         paste0("\n\nOutlier analysis identified ", sum(outliers), 
                " potential outlier stud", ifelse(sum(outliers)==1, "y", "ies"),
                " with extreme effect sizes. Sensitivity analysis excluding outliers showed ",
                ifelse(begg_changes | egger_changes,
                       "that bias test results were sensitive to these studies, indicating some fragility. ",
                       "that both Begg's and Egger's tests were robust to outlier removal. ")),
         ""),
  "\n\nConclusion: ",
  ifelse((begg_p < 0.10 | egger_p < 0.10) & k_imputed > 0 & abs(pct_change) >= 20,
         "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 appropriate caution. Prioritizing evidence from large, high-quality studies is recommended.",
         ifelse((begg_p < 0.10 | egger_p < 0.10) & k_imputed > 0 & abs(pct_change) < 20,
                "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 in interpreting results.",
                ifelse(begg_p >= 0.10 & 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 < 20, "modest power with k < 20 studies(particularly for Begg's test). ", "available power. "),
                              "Comprehensive search strategies including gray literature and trial registries strengthen confidence in findings."),
                       "Publication bias assessment yielded mixed results requiring careful interpretation and triangulation across methods.")))
)

cat(report)

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