Cochran's Q Test for Heterogeneity

# Cochran's Q Test for Heterogeneity Assessment
# Demonstrating Q statistic computation, interpretation, and integration with I² and τ²

library(metafor)      # rma() for meta-analysis, Q test built-in
library(meta)         # metagen() alternative
library(dplyr)
library(ggplot2)

# === STEP 1: Simulate Meta-Analytic Dataset ===
# In practice: data <- read.csv("meta_analysis_data.csv")
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 (I² ≈ 50%)
true_effects <- rnorm(k, mean = 0.70, sd = 0.20)

# 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

# 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 (true + sampling error)
observed_g <- rnorm(k, mean = true_effects, sd = sampling_se)
variance_g <- sampling_se^2

# Add therapy format moderator (individual vs. group)
therapy_format <- sample(c("Individual", "Group"), k, replace = TRUE, prob = c(0.6, 0.4))

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,
  therapy_format = therapy_format
)

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

# === STEP 2: Compute Cochran's Q Test (via metafor) ===
# Random-effects model automatically computes Q statistic
re_model <- rma(yi = hedges_g, vi = variance, data = meta_data, 
                method = "REML", slab = author_year)

print("\n\n=== COCHRAN'S Q TEST FOR HETEROGENEITY ===")
cat("\nTest of Heterogeneity:")
cat("\nQ statistic  =", round(re_model$QE, 2))
cat("\nDegrees of freedom(df) =", re_model$k - 1)
cat("\np-value =", format.pval(re_model$QEp, digits = 3))

# Critical value from χ² distribution
alpha <- 0.10  # Use .10, not .05, for heterogeneity tests
critical_value <- qchisq(1 - alpha, df = re_model$k - 1)
cat("\n\nCritical value(α = .10): χ²(", re_model$k - 1, ") =", round(critical_value, 2))

if (re_model$QEp < alpha) {
  q_conclusion <- "SIGNIFICANT heterogeneity detected"
  q_interpretation <- paste0(
    "Reject H₀ (p < .10). Effect sizes vary more than expected from sampling error alone.\n",
    "Random-effects meta-analysis is STRONGLY RECOMMENDED to account for between-study variance."
  )
} else {
  q_conclusion <- "No significant heterogeneity detected"
  q_interpretation <- paste0(
    "Fail to reject H₀ (p ≥ .10). Insufficient evidence that effects vary beyond sampling error.\n",
    "However, with k = ", k, ", Q test has LIMITED POWER. Random-effects model may still be preferred for generalization."
  )
}

cat("\n\nConclusion:", q_conclusion)
cat("\n\nInterpretation:", q_interpretation)

# === STEP 3: Extract and Interpret Heterogeneity Statistics ===
tau2 <- re_model$tau2
tau <- sqrt(tau2)
I2 <- re_model$I2
H2 <- re_model$H2

cat("\n\n=== HETEROGENEITY STATISTICS ===")
cat("\nτ² (tau-squared)  =", round(tau2, 4), "(between-study variance)")
cat("\nτ (tau)           =", round(tau, 3), "(between-study SD)")
cat("\nI²                =", round(I2, 1), "% (percent variance due to heterogeneity)")
cat("\nH²                =", round(H2, 2), "(variance inflation factor)")

if (I2 < 25) {
  I2_interp <- "low"
} else if (I2 < 50) {
  I2_interp <- "moderate"
} else if (I2 < 75) {
  I2_interp <- "substantial"
} else {
  I2_interp <- "considerable"
}

cat("\n\nI² Interpretation: Heterogeneity is", I2_interp)
cat("\n\nCRITICAL NOTE:")
cat("\n- Q statistic tests IF heterogeneity exists(hypothesis test, p-value)")
cat("\n- I² quantifies HOW MUCH heterogeneity(effect size, percentage)")
cat("\n- τ² quantifies heterogeneity in original metric units(variance)")
cat("\n- ALWAYS report all three: Q provides statistical test, I² and τ² provide magnitude\n")

# === STEP 4: Manual Calculation of Q (for pedagogical understanding) ===
cat("\n\n=== MANUAL CALCULATION OF Q STATISTIC(Step-by-Step) ===")

# Step 1: Fixed-effect weights (inverse variance)
w_fixed <- 1 / meta_data$variance
cat("\nStep 1: Calculate fixed-effect weights w_i = 1 / SE_i²")
cat("\nWeights(w_i):", paste(round(w_fixed, 2), collapse = ", "))

# Step 2: Pooled effect under fixed-effect model
pooled_fixed <- sum(w_fixed * meta_data$hedges_g) / sum(w_fixed)
cat("\n\nStep 2: Calculate pooled effect(fixed-effect)")
cat("\nθ̂ = Σ(w_i × θ_i) / Σw_i =", round(pooled_fixed, 3))

# Step 3: Squared deviations
deviations <- meta_data$hedges_g - pooled_fixed
squared_deviations <- deviations^2
weighted_squared_dev <- w_fixed * squared_deviations

cat("\n\nStep 3: Calculate squared deviations(θ_i - θ̂)²")
dev_table <- data.frame(
  Study = 1:k,
  Effect = round(meta_data$hedges_g, 3),
  Deviation = round(deviations, 3),
  Squared_Dev = round(squared_deviations, 4),
  Weight = round(w_fixed, 2),
  Weighted = round(weighted_squared_dev, 3)
)
print(head(dev_table, 5))
cat("\n...(showing first 5 of", k, "studies)")

# Step 4: Sum to get Q
Q_manual <- sum(weighted_squared_dev)
cat("\n\nStep 4: Sum weighted squared deviations")
cat("\nQ = Σ[w_i × (θ_i - θ̂)²] =", round(Q_manual, 2))

# Verify matches metafor output
cat("\n\nVerification: Manual Q =", round(Q_manual, 2), 
    "| metafor Q =", round(re_model$QE, 2))
cat("\nMatch:", ifelse(abs(Q_manual - re_model$QE) < 0.01, "YES ✓", "Check calculation"))

# === STEP 5: Power Analysis for Q Test ===
cat("\n\n=== POWER ANALYSIS FOR Q TEST ===")
cat("\nNumber of studies(k) =", k)
cat("\nDegrees of freedom(df) =", k - 1)

# Simulate power: probability of detecting heterogeneity given true I²
# Approximate using non-central χ² distribution
if (I2 > 0) {
  # Non-centrality parameter (NCP) approximation
  # NCP ≈ Q under alternative (depends on true τ² and study weights)
  # Rough approximation: NCP ≈ Q_observed if heterogeneity truly exists
  ncp_approx <- re_model$QE
  
  # Power: probability Q > critical value under alternative
  power_estimate <- 1 - pchisq(critical_value, df = k - 1, ncp = ncp_approx)
  
  cat("\n\nApproximate power to detect observed heterogeneity(I² =", round(I2, 1), "%):", 
      round(power_estimate * 100, 1), "%")
} else {
  cat("\n\nPower calculation not applicable(I² ≈ 0)")
}

cat("\n\nGeneral power guidelines for Q test:")
cat("\n- k < 10:  LOW power(~10-30% for moderate I² = 30-50%)")
cat("\n- k = 10-15: MODERATE power(~50-70% for moderate I²)")
cat("\n- k ≥ 20:  HIGH power(>80% for moderate I², >90% for substantial I²)")
cat("\n\nImplication with k =", k, ":")
if (k < 10) {
  cat(" LOW power. Non-significant Q doesn't prove homogeneity.")
  cat("\n  → DEFAULT to random-effects model given uncertainty.")
} else if (k < 20) {
  cat(" MODERATE power. Q test reasonably reliable.")
  cat("\n  → Use Q alongside I² and τ² for model selection.")
} else {
  cat(" HIGH power. Q test likely detects even small heterogeneity.")
  cat("\n  → Significant Q may reflect trivial heterogeneity; check I² for magnitude.")
}

# === STEP 6: Visual Assessment (Forest Plot) ===
par(mar = c(5, 4, 3, 2))
forest(re_model,
       xlab = "Hedges' g(CBT - Control)",
       header = c("Study", "g [95% CI]"),
       cex = 0.8,
       col = "steelblue",
       border = "steelblue")

mtext(paste0("Cochran's Q(", k - 1, ") = ", round(re_model$QE, 2), 
             ", p = ", format.pval(re_model$QEp, digits = 3),
             " | I² = ", round(I2, 1), "% (", I2_interp, " heterogeneity)"),
      side = 3, line = 0, cex = 0.9, font = 2)

# === STEP 7: Subgroup Analysis (Q Decomposition) ===
cat("\n\n=== SUBGROUP ANALYSIS: Q DECOMPOSITION ===")
cat("\nModerator: Therapy Format(Individual vs. Group CBT)\n")

# Subgroup meta-analysis
subgroup_model <- rma(yi = hedges_g, vi = variance, 
                      mods = ~ therapy_format - 1,  # Separate estimates per group
                      data = meta_data, method = "REML")

# Calculate Q statistics for subgroups manually
individual_data <- meta_data[meta_data$therapy_format == "Individual", ]
group_data <- meta_data[meta_data$therapy_format == "Group", ]

if (nrow(individual_data) >= 3 && nrow(group_data) >= 3) {
  # Individual CBT subgroup
  re_individual <- rma(yi = hedges_g, vi = variance, 
                       data = individual_data, method = "REML")
  Q_individual <- re_individual$QE
  df_individual <- re_individual$k - 1
  
  # Group CBT subgroup  
  re_group <- rma(yi = hedges_g, vi = variance, 
                  data = group_data, method = "REML")
  Q_group <- re_group$QE
  df_group <- re_group$k - 1
  
  # Q decomposition
  Q_within <- Q_individual + Q_group
  df_within <- df_individual + df_group
  Q_between <- re_model$QE - Q_within
  df_between <- (re_model$k - 1) - df_within
  p_between <- 1 - pchisq(Q_between, df = df_between)
  
  cat("\nQ Decomposition:")
  cat("\n  Q_total(", re_model$k - 1, ") =", round(re_model$QE, 2))
  cat("\n\n  Q_within(", df_within, ") =", round(Q_within, 2), 
      "(heterogeneity within subgroups)")
  cat("\n    - Q_individual(", df_individual, ") =", round(Q_individual, 2))
  cat("\n    - Q_group(", df_group, ") =", round(Q_group, 2))
  cat("\n\n  Q_between(", df_between, ") =", round(Q_between, 2), 
      "(heterogeneity between subgroups)")
  cat("\n    p-value =", format.pval(p_between, digits = 3))
  
  if (p_between < 0.10) {
    cat("\n\n  → SIGNIFICANT difference between subgroups(p < .10)")
    cat("\n    Therapy format is a significant moderator of CBT effects.")
  } else {
    cat("\n\n  → No significant difference between subgroups(p ≥ .10)")
    cat("\n    Therapy format does not significantly moderate CBT effects.")
  }
  
  # Subgroup estimates
  cat("\n\nSubgroup Estimates:")
  cat("\n  Individual CBT: g =", round(re_individual$beta[1], 3), 
      ", 95% CI [", round(re_individual$ci.lb, 3), ",", 
      round(re_individual$ci.ub, 3), "], I² =", round(re_individual$I2, 1), "%")
  cat("\n  Group CBT:      g =", round(re_group$beta[1], 3),
      ", 95% CI [", round(re_group$ci.lb, 3), ",",
      round(re_group$ci.ub, 3), "], I² =", round(re_group$I2, 1), "%")
} else {
  cat("\nInsufficient studies per subgroup for Q decomposition(need k ≥ 3 per group)")
}

# === STEP 8: Comparison with Fixed-Effect Model ===
cat("\n\n=== COMPARISON: FIXED vs. RANDOM EFFECTS ===")

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

cat("\nFixed-Effect Model(assumes τ² = 0):")
cat("\n  Pooled g =", round(as.numeric(fe_model$beta), 3))
cat("\n  95% CI: [", round(fe_model$ci.lb, 3), ",", round(fe_model$ci.ub, 3), "]")
cat("\n  SE =", round(fe_model$se, 4))

cat("\n\nRandom-Effects Model(estimates τ² from data):")
cat("\n  Pooled g =", round(as.numeric(re_model$beta), 3))
cat("\n  95% CI: [", round(re_model$ci.lb, 3), ",", round(re_model$ci.ub, 3), "]")
cat("\n  SE =", round(re_model$se, 4))
cat("\n  τ² =", round(tau2, 4))

ci_width_fe <- fe_model$ci.ub - fe_model$ci.lb
ci_width_re <- re_model$ci.ub - re_model$ci.lb
inflation <- (ci_width_re / ci_width_fe - 1) * 100

cat("\n\nDifference:")
cat("\n  Random-effects CI is", round(inflation, 1), "% wider than fixed-effect")
cat("\n(reflects additional uncertainty from between-study variance τ²)")

if (re_model$QEp < 0.10) {
  cat("\n\n→ Given significant Q test(p < .10) and I² =", round(I2, 1), "%,")
  cat("\n  RANDOM-EFFECTS model is STRONGLY RECOMMENDED.")
} else if (k < 10) {
  cat("\n\n→ Given k < 10 (low power for Q test),")
  cat("\n  RANDOM-EFFECTS model is RECOMMENDED for conservative generalization.")
} else {
  cat("\n\n→ Given non-significant Q test, FIXED-EFFECT model may be appropriate,")
  cat("\n  but RANDOM-EFFECTS often preferred for generalization beyond observed studies.")
}

# === STEP 9: Sensitivity Analysis (Influential Studies) ===
cat("\n\n=== SENSITIVITY ANALYSIS: INFLUENCE ON Q STATISTIC ===")

influence_Q <- numeric(k)
influence_I2 <- numeric(k)

for (i in 1:k) {
  # Remove study i
  mask <- (1:k) != i
  loo_model <- rma(yi = hedges_g[mask], vi = variance[mask], 
                   data = meta_data[mask, ], method = "REML")
  influence_Q[i] <- loo_model$QE
  influence_I2[i] <- loo_model$I2
}

influence_results <- data.frame(
  Study = meta_data$author_year,
  Q_without = round(influence_Q, 2),
  I2_without = round(influence_I2, 1),
  Q_change = round(influence_Q - re_model$QE, 2),
  I2_change = round(influence_I2 - I2, 1)
)

cat("\nLeave-One-Out Results(showing first 5):")
print(head(influence_results, 5))

max_Q_change <- max(abs(influence_results$Q_change))
most_influential <- which.max(abs(influence_results$Q_change))

cat("\n\nMost influential study:", meta_data$author_year[most_influential])
cat("\n  Removing this study changes Q by", round(influence_results$Q_change[most_influential], 2))
cat("\n  and changes I² by", round(influence_results$I2_change[most_influential], 1), "%")

if (max_Q_change > re_model$QE * 0.2) {
  cat("\n\n→ Q statistic SENSITIVE to individual studies(>20% change possible)")
  cat("\n  Report results with/without influential studies for transparency.")
} else {
  cat("\n\n→ Q statistic ROBUST to individual studies(<20% change across all leave-one-out)")
}

# === STEP 10: APA-Style Reporting ===
cat("\n\n=== APA-STYLE REPORTING TEMPLATE ===")
cat("\n\nCochran's Q test assessed heterogeneity in effect sizes across", k, "RCTs")
cat("\n(N =", sum(meta_data$total_n), "participants) examining CBT for depression.")
cat("\nThe test revealed", ifelse(re_model$QEp < 0.10, "significant", "non-significant"),
    "heterogeneity,")
cat("\nQ(", re_model$k - 1, ") =", round(re_model$QE, 2), ",",
    ifelse(re_model$QEp < 0.001, " p < .001", 
           paste0(" p = ", format.pval(re_model$QEp, digits = 3))),
    ", indicating that effect sizes")
cat("\nvaried", ifelse(re_model$QEp < 0.10, "significantly beyond", "within"),
    "sampling error expectations.")
cat("\n\nHeterogeneity statistics further quantified this variation: I² =", 
    round(I2, 1), "%")
cat("\n(", I2_interp, "heterogeneity) and τ² =", round(tau2, 3), 
    "(τ =", round(tau, 3), "),")
cat("\nindicating that approximately", round(I2, 0), "% of total variance was due to")
cat("\nbetween-study differences rather than sampling error.")

if (re_model$QEp < 0.10) {
  cat("\n\nGiven significant heterogeneity, a random-effects meta-analysis model was employed")
  cat("\nto account for both within- and between-study variance, yielding a pooled effect")
  cat("\nof g =", round(as.numeric(re_model$beta), 2), ", 95% CI [", 
      round(re_model$ci.lb, 2), ",", round(re_model$ci.ub, 2), "].")
} else {
  cat("\n\nWhile the Q test was non-significant, a random-effects model was retained")
  cat("\nfor conservative estimation and generalization beyond the observed studies,")
  cat("\nespecially given the limited power of the Q test with k =", k, "studies.")
}

cat("\n\nSensitivity analysis(leave-one-out) showed that no single study")
cat("\ndisproportionately influenced the heterogeneity assessment(Q range:")
cat("\n", round(min(influence_Q), 2), "to", round(max(influence_Q), 2), "),")
cat("\nsupporting the robustness of the heterogeneity conclusions.")

cat("\n\n[End of Analysis]\n")