I² Statistic

# I² Statistic Calculation with 95% Confidence Interval
# Demonstrating heterogeneity quantification in psychotherapy meta-analysis

library(metafor)      # rma() for meta-analysis and confint() for I² CI
library(meta)         # forest() plotting
library(dplyr)
library(ggplot2)

# === STEP 1: Simulate Meta-Analytic Dataset ===
# In practice: data <- read.csv("meta_analysis_data.csv")
# Required: effect sizes (Hedges' g) and variances

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

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

# Sample sizes vary
n_treat <- sample(50:100, k, replace=TRUE)
n_control <- sample(50: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
observed_g <- rnorm(k, mean=true_effects, sd=sampling_se)
variance_g <- sampling_se^2

# Create dataset
meta_data <- data.frame(
  study_id = paste0("Study_", 1:k),
  author_year = paste0(LETTERS[1:k], " et al.(20", 12:23, ")"),
  therapy_type = sample(c("CBT", "Exposure", "Other"), k, replace=TRUE),
  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: Calculate Cochran's Q Statistic ===
# Q is foundation for I² calculation

# Fixed-effects weights
w_fixed <- 1 / meta_data$variance

# Fixed-effects pooled estimate
pooled_fixed <- sum(w_fixed * meta_data$hedges_g) / sum(w_fixed)

# Cochran's Q = Σw_i(y_i - ȳ)²
Q <- sum(w_fixed * (meta_data$hedges_g - pooled_fixed)^2)
df_Q <- k - 1  # Degrees of freedom

cat("\n=== Cochran's Q Statistic ===")
cat("\nQ(", df_Q, ") =", round(Q, 3))

# Q test p-value (chi-square distribution)
Q_pval <- 1 - pchisq(Q, df_Q)
cat("\np-value =", format.pval(Q_pval, digits=3))

if (Q_pval < 0.10) {
  cat("\n→ Significant heterogeneity detected(p < .10)")
} else {
  cat("\n→ No significant heterogeneity(p ≥ .10)")
}

# === STEP 3: Calculate I² Statistic ===
# I² = max(0, 100% × (Q - df) / Q)
# This is the percentage of total variance due to heterogeneity

I2_point <- max(0, 100 * (Q - df_Q) / Q)

cat("\n\n=== I² Statistic ===")
cat("\nI² =", round(I2_point, 2), "%")

# Interpretation category
if (I2_point < 25) {
  I2_category <- "low"
  I2_meaning <- "might not be important"
} else if (I2_point < 50) {
  I2_category <- "moderate"
  I2_meaning <- "may represent moderate heterogeneity"
} else if (I2_point < 75) {
  I2_category <- "substantial"
  I2_meaning <- "may represent substantial heterogeneity"
} else {
  I2_category <- "considerable"
  I2_meaning <- "represents considerable heterogeneity"
}

cat("\nCategory:", I2_category)
cat("\nInterpretation:", I2_meaning)

# === STEP 4: CRITICAL - Calculate 95% Confidence Interval for I² ===
# Uses non-central chi-square distribution method
# CI reflects uncertainty in I² estimate (especially important with small k)

# Function to calculate I² CI using non-central chi-square method
calculate_I2_CI <- function(Q, df, alpha=0.05) {
  # Lower bound: Find non-centrality parameter λ_L such that P(χ²_df(λ_L) > Q) = alpha/2
  # Upper bound: Find λ_U such that P(χ²_df(λ_U) < Q) = alpha/2
  
  # Lower CI for I²
  if (Q > df) {
    # Solve for λ_L: P(χ²_df(λ_L) ≥ Q) = 1 - alpha/2
    # This requires finding λ such that Q is at (1-alpha/2) quantile of χ²_df(λ)
    # Approximation: Use iterative search
    find_ncp_lower <- function(ncp) {
      pchisq(Q, df=df, ncp=ncp, lower.tail=FALSE) - (1 - alpha/2)
    }
    
    if (find_ncp_lower(0) < 0) {
      # Q is extreme; use root finding
      lambda_L <- tryCatch({
        uniroot(find_ncp_lower, c(0, Q))$root
      }, error = function(e) 0)
    } else {
      lambda_L <- 0
    }
    
    I2_lower <- max(0, 100 * lambda_L / Q)
  } else {
    I2_lower <- 0
  }
  
  # Upper CI for I²
  find_ncp_upper <- function(ncp) {
    pchisq(Q, df=df, ncp=ncp, lower.tail=TRUE) - alpha/2
  }
  
  lambda_U <- tryCatch({
    # Search up to 4×Q for upper bound
    uniroot(find_ncp_upper, c(0, max(Q*4, 100)))$root
  }, error = function(e) Q)  # If fails, use Q as conservative upper
  
  I2_upper <- min(100, 100 * lambda_U / Q)
  
  return(c(lower = I2_lower, upper = I2_upper))
}

I2_CI <- calculate_I2_CI(Q, df_Q)

cat("\n\n=== I² with 95% Confidence Interval ===")
cat("\nI² =", round(I2_point, 1), "%, 95% CI [", 
    round(I2_CI[1], 1), "%, ", round(I2_CI[2], 1), "%]")

CI_width <- I2_CI[2] - I2_CI[1]
cat("\nCI width =", round(CI_width, 1), "percentage points")

if (CI_width > 40) {
  cat("\n→ Wide CI indicates imprecise I² estimate(common with k < 10)")
} else if (CI_width > 25) {
  cat("\n→ Moderate CI width; I² estimate has some uncertainty")
} else {
  cat("\n→ Narrow CI indicates precise I² estimate")
}

cat("\n\nInterpretation:", I2_category, "heterogeneity(", I2_meaning, ")")
cat("\nTrue I² likely falls between", round(I2_CI[1], 1), "% and", 
    round(I2_CI[2], 1), "%")

# === STEP 5: Run Random-Effects Meta-Analysis for Complete Output ===
# metafor provides I² automatically with other heterogeneity statistics

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

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

# Extract statistics
pooled_g <- as.numeric(re_model$beta)
ci_lower <- re_model$ci.lb
ci_upper <- re_model$ci.ub
tau2 <- re_model$tau2
tau <- sqrt(tau2)
I2_metafor <- re_model$I2
H2 <- re_model$H2

cat("\n=== Comprehensive Heterogeneity Statistics ===")
cat("\nI² =", round(I2_metafor, 1), "% (from metafor)")
cat("\nτ² (tau-squared) =", round(tau2, 4))
cat("\nτ (tau, between-study SD) =", round(tau, 3))
cat("\nH² =", round(H2, 2))
cat("\nQ(", df_Q, ") =", round(Q, 2), ", p =", format.pval(Q_pval, digits=3))

# Verify I² formula: I² = (Q - df) / Q × 100%
I2_manual <- max(0, 100 * (Q - df_Q) / Q)
cat("\n\nVerification: I² = (Q - df) / Q × 100%")
cat("\n            = (", round(Q, 2), " - ", df_Q, ") / ", round(Q, 2), " × 100%")
cat("\n            =", round(I2_manual, 2), "%")

# Alternative formula: I² = (H² - 1) / H² × 100%
I2_from_H2 <- 100 * (H2 - 1) / H2
cat("\n\nAlternative formula: I² = (H² - 1) / H² × 100%")
cat("\n                   = (", round(H2, 2), " - 1) / ", round(H2, 2), " × 100%")
cat("\n                   =", round(I2_from_H2, 2), "%")

# === STEP 6: Get 95% CI for I² using metafor's confint() ===
# More rigorous than manual calculation

I2_confint <- confint(re_model, digits=2)

cat("\n\n=== metafor confint() for I² ===")
cat("\nI² =", round(I2_metafor, 1), "%")
cat("\n95% CI: [", round(I2_confint$random["I^2(%)", "ci.lb"], 1), "%, ",
    round(I2_confint$random["I^2(%)", "ci.ub"], 1), "%]")

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

title(main = paste0("Forest Plot: Psychotherapy for Anxiety(k=", k, " studies)\n",
                    "I² = ", round(I2_metafor, 1), "% (", I2_category, 
                    " heterogeneity), 95% CI [",
                    round(I2_confint$random["I^2(%)", "ci.lb"], 1), "%, ",
                    round(I2_confint$random["I^2(%)", "ci.ub"], 1), "%]"),
      cex.main = 0.9, line = 0.5)

# Add visual heterogeneity annotation
mtext(paste0("Q(", df_Q, ") = ", round(Q, 2), ", p ", 
             ifelse(Q_pval < 0.001, "< .001", paste0("= ", round(Q_pval, 3)))),
      side = 3, line = -1.5, cex = 0.8)

# === STEP 8: I² by Subgroup (Investigating Heterogeneity Source) ===
# When I² > 50%, investigate moderators

if (I2_metafor > 50) {
  cat("\n\n=== I² > 50%: Investigating Heterogeneity via Subgroup Analysis ===")
  
  # Subgroup analysis by therapy type
  subgroup_model <- rma(yi = hedges_g, vi = variance, data = meta_data,
                        mods = ~ therapy_type - 1, method = "REML")
  
  cat("\nSubgroup Analysis by Therapy Type:")
  print(subgroup_model)
  
  # Calculate I² within each subgroup
  for (therapy in unique(meta_data$therapy_type)) {
    subdata <- meta_data[meta_data$therapy_type == therapy, ]
    if (nrow(subdata) >= 3) {
      sub_model <- rma(yi = hedges_g, vi = variance, data = subdata, method = "REML")
      cat("\n", therapy, ":", nrow(subdata), "studies")
      cat("\n  I² =", round(sub_model$I2, 1), "%")
      cat("\n  Pooled g =", round(as.numeric(sub_model$beta), 2), 
          ", 95% CI [", round(sub_model$ci.lb, 2), ",", 
          round(sub_model$ci.ub, 2), "]\n")
    }
  }
  
  # Test between-group heterogeneity
  cat("\nBetween-group heterogeneity test:")
  cat("\nQ_between =", round(subgroup_model$QM, 2))
  cat("\np-value =", format.pval(subgroup_model$QMp, digits=3))
  
  if (subgroup_model$QMp < 0.05) {
    cat("\n→ Therapy type significantly moderates effect size")
    cat("\n(explains some heterogeneity)")
  } else {
    cat("\n→ No significant moderation by therapy type")
  }
}

# === STEP 9: I² Interpretation Guide ===
cat("\n\n=== I² INTERPRETATION GUIDE ===")
cat("\n\nObserved I² =", round(I2_metafor, 1), "%, 95% CI [",
    round(I2_confint$random["I^2(%)", "ci.lb"], 1), "%, ",
    round(I2_confint$random["I^2(%)", "ci.ub"], 1), "%]")
cat("\n\nMeaning:")
cat("\n- Approximately", round(I2_metafor, 0), 
    "% of total variability in effect sizes is due to")
cat("\n  true heterogeneity(between-study differences), not sampling error")
cat("\n-", round(100 - I2_metafor, 0), 
    "% of variability is attributable to sampling error(within-study variance)")

cat("\n\nClinical Implications:")
if (I2_metafor < 25) {
  cat("\n- Effects are relatively homogeneous across studies")
  cat("\n- Fixed-effect and random-effects models yield similar results")
  cat("\n- Pooled effect likely generalizes consistently")
} else if (I2_metafor < 50) {
  cat("\n- Moderate heterogeneity suggests some true variation in effects")
  cat("\n- Random-effects model preferred for pooling")
  cat("\n- Consider exploratory moderator analysis")
} else if (I2_metafor < 75) {
  cat("\n- Substantial heterogeneity indicates considerable true variation")
  cat("\n- INVESTIGATE sources via subgroup analysis or meta-regression")
  cat("\n- Report prediction interval prominently(not just CI)")
  cat("\n- Pooled estimate may mask important subgroup differences")
} else {
  cat("\n- Considerable heterogeneity suggests extreme variation")
  cat("\n- Question whether pooling is appropriate(may be comparing apples to oranges)")
  cat("\n- Focus on understanding WHY effects vary, not just pooled estimate")
  cat("\n- Narrative synthesis may be more informative than single pooled effect")
}

cat("\n\nStatistical Decisions:")
cat("\n- Model choice: Random-effects model STRONGLY preferred")
cat("\n(I² >", ifelse(I2_metafor > 25, "25% indicates meaningful heterogeneity)", "0%)"))
cat("\n- Moderator analysis:", 
    ifelse(I2_metafor > 50, "WARRANTED(I² > 50%)", 
           "Consider if theory-driven(I² < 50%)"))
cat("\n- Prediction interval: ESSENTIAL for interpreting generalizability")

# === STEP 10: Comparison with τ² ===
cat("\n\n=== I² vs. τ² (Complementary Heterogeneity Measures) ===")
cat("\n\nI² (relative measure):")
cat("\n  I² =", round(I2_metafor, 1), "% of total variance due to heterogeneity")
cat("\n  → Scale-free; comparable across meta-analyses")
cat("\n  → Interpreted via guidelines: <25% low, 25-50% moderate, etc.")

cat("\n\nτ² (absolute measure):")
cat("\n  τ² =", round(tau2, 4), "(between-study variance in g² units)")
cat("\n  τ =", round(tau, 3), "(between-study SD in g units)")
cat("\n  → Metric-dependent; magnitude depends on effect size scale")
cat("\n  → Interpret relative to pooled effect: g =", round(pooled_g, 2))
cat("\n     Effects vary approximately", round(pooled_g - tau, 2), "to", 
    round(pooled_g + tau, 2), "(θ ± τ)")

cat("\n\nRelationship:")
cat("\n  I² increases with τ² BUT also depends on study precision")
cat("\n  - Same τ² yields higher I² when studies are large(precise)")
cat("\n  - Same τ² yields lower I² when studies are small(imprecise)")
cat("\n  → Report BOTH: I² for interpretation category, τ² for absolute magnitude")

# === STEP 11: APA-Style Reporting ===
cat("\n\n=== APA-STYLE REPORT ===")
report <- paste0(
  "A random-effects meta-analysis of ", k, " RCTs(N = ", sum(meta_data$total_n),
  " participants)\nexamined psychotherapy interventions vs. control for anxiety disorders. ",
  "The pooled effect\nwas Hedges' g = ", round(pooled_g, 2), ", 95% CI [", 
  round(ci_lower, 2), ", ", round(ci_upper, 2), "],\n",
  "indicating a ", ifelse(abs(pooled_g) >= 0.8, "large", 
                          ifelse(abs(pooled_g) >= 0.5, "medium-to-large", "medium")),
  " effect favoring psychotherapy.\n\n",
  "Heterogeneity was ", I2_category, " (I² = ", round(I2_metafor, 1), 
  "%, 95% CI [", round(I2_confint$random["I^2(%)", "ci.lb"], 1), "%, ",
  round(I2_confint$random["I^2(%)", "ci.ub"], 1), "%];\n",
  "τ² = ", round(tau2, 3), ", τ = ", round(tau, 2), "; ",
  "Q(", df_Q, ") = ", round(Q, 2),
  ifelse(Q_pval < 0.001, ", p < .001", paste0(", p = ", round(Q_pval, 3))), ").\n\n",
  "The I² statistic indicates that approximately ", round(I2_metafor, 0),
  "% of the total variability\nin observed effect sizes is attributable to true heterogeneity(between-study differences)\n",
  "rather than sampling error. "
)

if (I2_metafor > 50) {
  report <- paste0(report, "Given the substantial heterogeneity, we conducted\n",
    "subgroup analysis by therapy type to investigate potential moderators. ")
  if (exists("subgroup_model") && subgroup_model$QMp < 0.05) {
    report <- paste0(report, "Therapy type significantly\nmoderated effect sizes ",
      "(Q_between = ", round(subgroup_model$QM, 2), ", p ", 
      ifelse(subgroup_model$QMp < 0.001, "< .001", 
             paste0("= ", round(subgroup_model$QMp, 3))),
      "),\nexplaining a portion of the observed heterogeneity. ")
  }
}

report <- paste0(report, "\n\nConclusion: Psychotherapy produces ",
  ifelse(abs(pooled_g) >= 0.8, "large", "medium-to-large"),
  " benefits for anxiety disorders on average,\nbut the ", I2_category, 
  " heterogeneity(I² = ", round(I2_metafor, 1), "%) indicates that effect\n",
  "magnitude varies considerably across studies. ")

if (I2_metafor > 50) {
  report <- paste0(report, 
    "The wide prediction interval and\nsubstantial I² suggest that treatment effects are context-dependent, warranting\n",
    "careful consideration of patient characteristics and intervention specifics when\n",
    "applying these findings to clinical practice. Future research should focus on\n",
    "identifying moderators to predict which patients benefit most from which therapies.")
} else {
  report <- paste0(report, 
    "The relatively\nconsistent effects across studies support broad applicability of psychotherapy\n",
    "for anxiety, though individual patient factors should always be considered.")
}

cat(report)

# === STEP 12: Visual Summary of I² with CI ===
par(mar=c(5,5,4,2))
barplot(I2_metafor, ylim=c(0, 100), 
        col="steelblue", border="darkblue", width=0.5,
        xlab="", ylab="I² (%)", cex.lab=1.2,
        main=paste0("I² Statistic with 95% Confidence Interval\n",
                    "Psychotherapy for Anxiety Meta-Analysis"),
        cex.main=1.1)

# Add CI error bars
arrows(x0=0.5, y0=I2_confint$random["I^2(%)", "ci.lb"],
       x1=0.5, y1=I2_confint$random["I^2(%)", "ci.ub"],
       angle=90, code=3, length=0.15, lwd=3, col="darkred")

# Add interpretation zones
abline(h=25, lty=2, col="gray50", lwd=1.5)
abline(h=50, lty=2, col="gray50", lwd=1.5)
abline(h=75, lty=2, col="gray50", lwd=1.5)

text(0.5, 12.5, "Low", cex=0.9, col="gray30")
text(0.5, 37.5, "Moderate", cex=0.9, col="gray30")
text(0.5, 62.5, "Substantial", cex=0.9, col="gray30")
text(0.5, 87.5, "Considerable", cex=0.9, col="gray30")

# Add I² value
text(0.5, I2_metafor + 8, 
     paste0("I² = ", round(I2_metafor, 1), "%\n",
            "95% CI [", round(I2_confint$random["I^2(%)", "ci.lb"], 1), "%, ",
            round(I2_confint$random["I^2(%)", "ci.ub"], 1), "%]"),
     cex=1.1, font=2, col="darkblue")