Fixed-Effects Meta-Analysis

# Fixed-Effects Meta-Analysis: Homogeneous CBT Protocol
# Demonstrating when fixed-effects is appropriate and comparison with random-effects

library(metafor)      # rma() for meta-analysis
library(meta)         # forest(), funnel() plotting
library(dplyr)
library(ggplot2)

# === STEP 1: Simulate Meta-Analytic Dataset ===
# Scenario: Highly homogeneous studies (same CBT manual, similar populations)
# Low heterogeneity expected (τ ≈ 0.05, I² ≈ 18%)

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

# Simulate effect sizes with MINIMAL heterogeneity (fixed-effects scenario)
# True effects vary slightly: mean θ=0.72, small between-study SD τ=0.05
true_effects <- rnorm(k, mean=0.72, sd=0.05)  # Minimal heterogeneity

# Sample sizes vary across studies
n_treat <- sample(50:80, k, replace=TRUE)
n_control <- sample(50:80, k, replace=TRUE)
total_n <- n_treat + n_control

# 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:21, ")"),
  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(Homogeneous Studies) ===")
print(meta_data)
cat("\nTotal N =", sum(total_n), "participants across", k, "studies")
cat("\nScenario: Identical CBT-Beck 12-week protocol; similar populations(MDD, HRSD>18)\n")

# === STEP 2A: Fixed-Effects Meta-Analysis (Inverse Variance) ===
# Weights = 1/SE² (no τ² component)
fe_model <- rma(yi = hedges_g, vi = variance, data = meta_data, 
                method = "FE", slab = author_year)

print("\n=== FIXED-EFFECTS Meta-Analysis Results ===")
print(fe_model)

# Extract key statistics
pooled_g_fe <- as.numeric(fe_model$beta)
ci_lower_fe <- fe_model$ci.lb
ci_upper_fe <- fe_model$ci.ub
p_value_fe <- fe_model$pval

# Heterogeneity statistics (calculated but not used in fixed-effects model)
Q <- fe_model$QE             # Cochran's Q statistic
Q_pval <- fe_model$QEp       # Q test p-value
I2 <- fe_model$I2            # % variance due to heterogeneity
H2 <- fe_model$H2            # Ratio of total to sampling variance

cat("\n=== Pooled Effect(Fixed-Effects) ===")
cat("\nHedges' g =", round(pooled_g_fe, 3))
cat("\n95% CI: [", round(ci_lower_fe, 3), ",", round(ci_upper_fe, 3), "]" )
cat("\np-value:", format.pval(p_value_fe, digits=3))

cat("\n\n=== Heterogeneity Statistics(Test of Fixed-Effects Assumption) ===")
cat("\nCochran's Q(", fe_model$k-1, ") =", round(Q, 2), ", p =", 
    format.pval(Q_pval, digits=3))
cat("\nI² =", round(I2, 1), "%")
cat("\nH² =", round(H2, 2))

if (I2 < 25) {
  heterogeneity_interp <- "low"
  fe_appropriate <- TRUE
} else if (I2 < 50) {
  heterogeneity_interp <- "moderate"
  fe_appropriate <- FALSE
} else if (I2 < 75) {
  heterogeneity_interp <- "substantial"
  fe_appropriate <- FALSE
} else {
  heterogeneity_interp <- "considerable"
  fe_appropriate <- FALSE
}

cat("\nInterpretation: Heterogeneity is", heterogeneity_interp)

if (Q_pval > 0.10 && I2 < 25) {
  cat("\n→ Q test non-significant(p >", round(Q_pval, 3), ") AND I² < 25%")
  cat("\n→ FIXED-EFFECTS MODEL APPROPRIATE: Homogeneity assumption supported")
  cat("\n→ Studies appear to share a common true effect(τ² ≈ 0)\n")
} else if (Q_pval <= 0.10) {
  cat("\n→ Q test SIGNIFICANT(p ≤ .10): Heterogeneity detected")
  cat("\n→ FIXED-EFFECTS MODEL INAPPROPRIATE: Use random-effects instead")
  cat("\n→ Studies do NOT share a common true effect(τ² > 0)\n")
} else {
  cat("\n→ Q test non-significant but I² =", round(I2, 1), "% (borderline)")
  cat("\n→ Consider random-effects for robustness and generalizability\n")
}

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

pooled_g_re <- as.numeric(re_model$beta)
ci_lower_re <- re_model$ci.lb
ci_upper_re <- re_model$ci.ub
tau2 <- re_model$tau2
tau <- sqrt(tau2)

cat("\n=== RANDOM-EFFECTS Meta-Analysis Results(Comparison) ===")
cat("\nPooled Hedges' g =", round(pooled_g_re, 3))
cat("\n95% CI: [", round(ci_lower_re, 3), ",", round(ci_upper_re, 3), "]")
cat("\nτ² =", round(tau2, 4), ", τ =", round(tau, 3))

# Prediction interval (random-effects only)
pi <- predict(re_model, digits=3)
cat("\n95% Prediction Interval: [", round(pi$pi.lb, 3), ",", round(pi$pi.ub, 3), "]")

cat("\n\n=== COMPARISON: Fixed vs. Random Effects ===")
cat("\nFixed-Effect:  g =", round(pooled_g_fe, 3),
    ", 95% CI [", round(ci_lower_fe, 3), ",", round(ci_upper_fe, 3), "]")
cat("\nRandom-Effect: g =", round(pooled_g_re, 3),
    ", 95% CI [", round(ci_lower_re, 3), ",", round(ci_upper_re, 3), "]")

ci_width_fe <- ci_upper_fe - ci_lower_fe
ci_width_re <- ci_upper_re - ci_lower_re

cat("\n\nCI width: Fixed =", round(ci_width_fe, 3), 
    ", Random =", round(ci_width_re, 3))

if (abs(pooled_g_fe - pooled_g_re) < 0.05 && ci_width_fe < ci_width_re * 1.1) {
  cat("\n→ Estimates very similar; low heterogeneity confirms fixed-effects appropriate")
  cat("\n→ Fixed-effects provides more precise estimate(narrower CI) due to τ²≈0")
} else if (ci_width_re > ci_width_fe * 1.3) {
  cat("\n→ Random-effects CI substantially wider(accounts for τ²)")
  cat("\n→ Heterogeneity present; random-effects more appropriate for honest uncertainty")
} else {
  cat("\n→ Moderate difference; both models yield similar conclusions")
}

if (fe_appropriate) {
  cat("\n\nRECOMMENDATION: Use FIXED-EFFECTS model(I²<25%, Q p>.10)")
  cat("\n- Studies are homogeneous(same protocol, population)")
  cat("\n- Inference is conditional: estimate applies to THESE studies")
  cat("\n- For generalization, still prefer random-effects(accounts for uncertainty)\n")
} else {
  cat("\n\nRECOMMENDATION: Use RANDOM-EFFECTS model(I²≥25% or Q p≤.10)")
  cat("\n- Heterogeneity detected; fixed-effects assumption violated")
  cat("\n- Random-effects accounts for between-study variance(τ²)")
  cat("\n- Enables generalization via prediction interval\n")
}

# === STEP 3: Forest Plot (Fixed-Effects with Comparison) ===
par(mfrow=c(2,1), mar=c(4,4,3,2))

# Fixed-effects forest plot
forest(fe_model, 
       xlab = "Hedges' g(CBT - Control)",
       slab = meta_data$author_year,
       header = c("Study", "g [95% CI]"),
       cex = 0.75,
       col = "darkblue",
       border = "darkblue",
       lwd = 2)
mtext(paste0("Fixed-Effects Model: g = ", round(pooled_g_fe, 2), 
             ", 95% CI [", round(ci_lower_fe, 2), ", ", round(ci_upper_fe, 2), "]"),
      side=3, line=1.5, cex=0.85, font=2)
mtext(paste0("I² = ", round(I2, 1), "% (", heterogeneity_interp, 
             "); Q(", fe_model$k-1, ") p = ", round(Q_pval, 3)),
      side=3, line=0.3, cex=0.75)

# Random-effects forest plot (for comparison)
forest(re_model, 
       xlab = "Hedges' g(CBT - Control)",
       slab = meta_data$author_year,
       header = c("Study", "g [95% CI]"),
       cex = 0.75,
       addpred = TRUE,  # Add prediction interval
       col = "darkgreen",
       border = "darkgreen",
       lwd = 2)
mtext(paste0("Random-Effects Model: g = ", round(pooled_g_re, 2), 
             ", 95% CI [", round(ci_lower_re, 2), ", ", round(ci_upper_re, 2), "]"),
      side=3, line=1.5, cex=0.85, font=2)
mtext(paste0("τ² = ", round(tau2, 4), "; 95% PI [", 
             round(pi$pi.lb, 2), ", ", round(pi$pi.ub, 2), "]"),
      side=3, line=0.3, cex=0.75)

par(mfrow=c(1,1))

# === STEP 4: Study Weights Comparison ===
cat("\n=== Study Weights: Fixed vs. Random Effects ===")

weights_fe <- weights(fe_model)
weights_re <- weights(re_model)

weights_df <- data.frame(
  study = meta_data$author_year,
  n_total = meta_data$total_n,
  se = meta_data$se,
  weight_fe = weights_fe,
  weight_re = weights_re,
  diff = weights_fe - weights_re
)

print(weights_df)

cat("\nNOTE: Fixed-effects weights = 100% × (1/SE²) / Σ(1/SE²)")
cat("\n      Random-effects weights = 100% × (1/(SE²+τ²)) / Σ(1/(SE²+τ²))")
cat("\n      With low τ², weights are similar; large studies dominate in both models\n")

# Visualize weights
par(mfrow=c(1,2), mar=c(5,4,3,2))

# Fixed-effects weights
barplot(weights_fe, names.arg=1:k, 
        main="Fixed-Effects Weights",
        xlab="Study", ylab="Weight(%)",
        col="steelblue", border="black")
abline(h=mean(weights_fe), lty=2, col="red", lwd=2)

# Random-effects weights
barplot(weights_re, names.arg=1:k,
        main="Random-Effects Weights",
        xlab="Study", ylab="Weight(%)",
        col="darkgreen", border="black")
abline(h=mean(weights_re), lty=2, col="red", lwd=2)

par(mfrow=c(1,1))

# === STEP 5: Funnel Plot & Publication Bias ===
par(mfrow=c(1,2))

# Funnel plot (fixed-effects)
funnel(fe_model, 
       xlab = "Hedges' g",
       ylab = "Standard Error",
       main = "Funnel Plot(Fixed-Effects)",
       back = "white",
       shade = "white")

# Egger's regression test
egger_test <- regtest(fe_model, model="lm")
cat("\n\n=== Publication Bias Assessment ===")
cat("\nEgger's Regression Test(Fixed-Effects):")
cat("\n  Intercept =", round(egger_test$zval, 3))
cat("\n  p-value =", format.pval(egger_test$pval, digits=3))
if (egger_test$pval < 0.10) {
  cat("\n  → Significant asymmetry detected(p<.10); publication bias possible")
  cat("\n  → WARNING: Fixed-effects especially vulnerable(large studies dominate)\n")
} else {
  cat("\n  → No significant asymmetry(p≥.10); limited evidence of bias\n")
}

# Trim-and-fill analysis
taf <- trimfill(fe_model)
cat("\nTrim-and-Fill Analysis(Fixed-Effects):")
cat("\n  Imputed studies(k₀) =", taf$k0)
if (taf$k0 > 0) {
  cat("\n  Adjusted g =", round(as.numeric(taf$beta), 3))
  cat("\n  Adjusted 95% CI: [", round(taf$ci.lb, 3), ",", round(taf$ci.ub, 3), "]")
  cat("\n  → Bias correction attenuates effect by", 
      round((pooled_g_fe - as.numeric(taf$beta))/pooled_g_fe * 100, 1), "%\n")
} else {
  cat("\n  → No missing studies imputed; no evidence of bias\n")
}

# Funnel plot with trim-and-fill
funnel(taf, 
       xlab = "Hedges' g",
       ylab = "Standard Error",
       main = "Trim-and-Fill Adjusted",
       back = "white",
       shade = "white",
       col = c("blue", "red"))
if (taf$k0 > 0) {
  legend("topright", c("Observed", "Imputed"), 
         col=c("blue", "red"), pch=19, cex=0.8)
}

par(mfrow=c(1,1))

# === STEP 6: Influence Analysis (Leave-One-Out) ===
influence_results <- leave1out(fe_model, digits=3)

cat("\n=== Influence Analysis(Leave-One-Out, Fixed-Effects) ===")
print(influence_results)

cat("\nEffect size range(leave-one-out):", 
    round(min(influence_results$estimate), 3), "to", 
    round(max(influence_results$estimate), 3))

if (max(abs(influence_results$estimate - pooled_g_fe)) < 0.1) {
  cat("\n→ Pooled estimate is ROBUST(minimal change when removing any study)\n")
} else {
  cat("\n→ Pooled estimate is SENSITIVE to individual studies; caution warranted\n")
}

# Plot influence
par(mar=c(5,4,3,2))
plot(1:k, influence_results$estimate, 
     ylim = range(c(influence_results$ci.lb, influence_results$ci.ub)),
     xlab = "Study Removed", ylab = "Pooled g(Fixed-Effects)",
     main = "Influence Analysis: Leave-One-Out",
     pch = 19, col = "darkblue", cex=1.2)
abline(h = pooled_g_fe, lty=2, col="red", lwd=2)
segments(1:k, influence_results$ci.lb, 1:k, influence_results$ci.ub, 
         col="darkblue", lwd=1.5)
legend("topright", "Full model", lty=2, col="red", lwd=2, cex=0.9)

# === STEP 7: Cumulative Meta-Analysis ===
cat("\n\n=== Cumulative Meta-Analysis(Fixed-Effects) ===")
cumul_results <- cumul(fe_model, order=order(meta_data$author_year))
print(cumul_results)

par(mar=c(5,4,3,2))
forest(cumul_results,
       xlab="Cumulative Hedges' g",
       header=c("Study Added", "g [95% CI]"),
       cex=0.75)
mtext("Cumulative Meta-Analysis(Fixed-Effects)",
      side=3, line=1, cex=0.9, font=2)

# === STEP 8: APA-Style Reporting ===
cat("\n\n=== APA-STYLE REPORT ===")
cat("\nA fixed-effects meta-analysis of", k, "RCTs(N =", sum(total_n), 
    "participants) examined\na standardized CBT protocol(Beck 12-week manual) vs. waitlist control for major\ndepression. Heterogeneity testing supported the fixed-effects assumption: Cochran's\nQ(", fe_model$k-1, ") =", round(Q, 2), ", p =", round(Q_pval, 3), 
    "(non-significant), I² =", round(I2, 1), "%\n(low heterogeneity), indicating studies share a common true effect.\n")

cat("\nThe pooled effect(inverse variance weights) was Hedges' g =", 
    round(pooled_g_fe, 2), ",\n95% CI [", round(ci_lower_fe, 2), ",", round(ci_upper_fe, 2), "], p",
    ifelse(p_value_fe < 0.001, " < .001", paste0(" = ", round(p_value_fe, 3))),
    ", indicating a\n", ifelse(abs(pooled_g_fe) < 0.5, "small to medium", 
              ifelse(abs(pooled_g_fe) < 0.8, "medium to large", "large")),
    " effect favoring CBT. This estimate represents the common effect\nacross the", k, 
    "included studies using this specific protocol.\n")

cat("\nFor comparison, a random-effects model yielded g =", round(pooled_g_re, 2),
    ", 95% CI\n[", round(ci_lower_re, 2), ",", round(ci_upper_re, 2), 
    "], τ² =", round(tau2, 4), 
    ", nearly identical to fixed-effects,\nconfirming minimal heterogeneity. ",
    "The random-effects 95% prediction interval\n[", round(pi$pi.lb, 2), ",", 
    round(pi$pi.ub, 2), "] suggests that a new study using this protocol\nwould ",
    ifelse(pi$pi.lb > 0, "consistently show beneficial effects.",
           "show variable effects across contexts."))

if (egger_test$pval < 0.10) {
  cat("\n\nPublication bias assessment revealed funnel plot asymmetry(Egger's test p",
      ifelse(egger_test$pval < 0.001, " < .001", 
             paste0(" = ", round(egger_test$pval, 3))),
      "),\nsuggesting potential small-study effects.")
  if (taf$k0 > 0) {
    cat(" Trim-and-fill analysis estimated", taf$k0, 
        "missing\nstudies, with adjusted g =", round(as.numeric(taf$beta), 2),
        ", 95% CI [", round(taf$ci.lb, 2), ",", round(taf$ci.ub, 2), "],\n",
        ifelse(abs(as.numeric(taf$beta)) > abs(pooled_g_fe)*0.7,
               "which remains substantial.", 
               "indicating potential overestimation."))
  }
} else {
  cat("\n\nPublication bias assessment showed no significant funnel plot asymmetry\n(Egger's test p =", 
      round(egger_test$pval, 3), "), with limited evidence of bias.")
}

cat("\n\nInfluence analysis(leave-one-out) confirmed robustness: pooled estimate\nranged from", 
    round(min(influence_results$estimate), 2), "to", 
    round(max(influence_results$estimate), 2), 
    "across analyses, indicating\nno single study disproportionately influenced results.\n")

cat("\nConclusion: For the specific standardized CBT-Beck 12-week protocol in adults\nwith MDD, a", 
    ifelse(abs(pooled_g_fe) >= 0.8, "large", 
           ifelse(abs(pooled_g_fe) >= 0.5, "medium to large", "medium")),
    " pooled effect(g =", round(pooled_g_fe, 2), 
    ") was observed with\nlow heterogeneity. Fixed-effects model was appropriate given homogeneous studies.\nInference is conditional: this estimate applies to the included studies using\nthis specific protocol. For broader generalization to other CBT approaches or\npopulations, random-effects meta-analysis with diverse studies is preferred.\n")