Hedges' g

# Hedges' g: Small Sample Clinical Trial - Mindfulness for Anxiety
# Demonstrates bias correction in small samples (n=28)

library(effsize)     # cohen.d() function
library(ggplot2)     # visualization
library(dplyr)       # data manipulation
library(car)         # Levene's test
library(moments)     # skewness/kurtosis

# ============================================================================
# SIMULATE DATA: Small RCT (Intervention n=15, Control n=13)
# ============================================================================
set.seed(2025)
data <- data.frame(
  group = c(rep("Intervention", 15), rep("Control", 13)),
  anxiety = c(
    rnorm(15, mean=42.3, sd=9.2),   # Intervention: M=42.3, SD=9.2
    rnorm(13, mean=51.7, sd=8.5)    # Control: M=51.7, SD=8.5
  )
)

cat("======================================================================\n")
cat("HEDGES' G: BIAS-CORRECTED EFFECT SIZE FOR SMALL SAMPLES\n")
cat("======================================================================\n\n")

# ============================================================================
# STEP 1: DESCRIPTIVE STATISTICS
# ============================================================================
cat("STEP 1: DESCRIPTIVE STATISTICS\n")
cat("------------------------------\n")
desc_stats <- data %>%
  group_by(group) %>%
  summarise(
    n = n(),
    M = mean(anxiety),
    SD = sd(anxiety),
    Min = min(anxiety),
    Max = max(anxiety),
    Skewness = skewness(anxiety),
    Kurtosis = kurtosis(anxiety) - 3  # Excess kurtosis
  )
print(desc_stats)

n_intervention <- sum(data$group == "Intervention")
n_control <- sum(data$group == "Control")
total_n <- nrow(data)

cat("\nTotal N =", total_n, "(SMALL SAMPLE → Hedges' g recommended)\n\n")

# ============================================================================
# STEP 2: DIAGNOSTIC CHECKS
# ============================================================================
cat("STEP 2: DIAGNOSTIC CHECKS\n")
cat("-------------------------\n\n")

# 2.1 Normality checks
cat("2.1 Normality Assessment:\n")
shapiro_intervention <- shapiro.test(data$anxiety[data$group == "Intervention"])
shapiro_control <- shapiro.test(data$anxiety[data$group == "Control"])

cat("Intervention: Shapiro-Wilk W =", round(shapiro_intervention$statistic, 3),
    ", p =", round(shapiro_intervention$p.value, 3))
if (shapiro_intervention$p.value > 0.05) {
  cat(" → Normal\n")
} else {
  cat(" → Non-normal(consider robust methods)\n")
}

cat("Control:      Shapiro-Wilk W =", round(shapiro_control$statistic, 3),
    ", p =", round(shapiro_control$p.value, 3))
if (shapiro_control$p.value > 0.05) {
  cat(" → Normal\n")
} else {
  cat(" → Non-normal(consider robust methods)\n")
}

cat("Interpretation: p > .05 suggests approximate normality(acceptable for g)\n\n")

# 2.2 Homogeneity of variance
cat("2.2 Homogeneity of Variance:\n")
levene_test <- leveneTest(anxiety ~ group, data=data)
cat("Levene's Test: F =", round(levene_test$`F value`[1], 3),
    ", p =", round(levene_test$`Pr(>F)`[1], 3))
if (levene_test$`Pr(>F)`[1] > 0.05) {
  cat(" → Variances equal(pooled SD appropriate)\n")
} else {
  cat(" → Variances unequal(consider Glass's delta or Welch correction)\n")
}

SD_intervention <- sd(data$anxiety[data$group == "Intervention"])
SD_control <- sd(data$anxiety[data$group == "Control"])
variance_ratio <- SD_intervention^2 / SD_control^2
cat("Variance ratio(Intervention/Control) =", round(variance_ratio, 2))
if (variance_ratio > 0.5 & variance_ratio < 2.0) {
  cat(" → Within acceptable range(0.5-2.0)\n\n")
} else {
  cat(" → Outside acceptable range(pooled SD may be inappropriate)\n\n")
}

# 2.3 Outlier detection
cat("2.3 Outlier Detection:\n")
outliers <- data %>%
  group_by(group) %>%
  mutate(
    z_score = (anxiety - mean(anxiety)) / sd(anxiety),
    is_outlier = abs(z_score) > 3
  ) %>%
  filter(is_outlier)

if (nrow(outliers) == 0) {
  cat("No extreme outliers detected(|z| > 3)\n\n")
} else {
  cat("WARNING:", nrow(outliers), "extreme outliers detected:\n")
  print(outliers[, c("group", "anxiety", "z_score")])
  cat("Consider sensitivity analysis or robust methods\n\n")
}

# ============================================================================
# STEP 3: VISUAL COMPARISON
# ============================================================================
cat("STEP 3: CREATING VISUALIZATIONS\n")
cat("--------------------------------\n\n")

# 3.1 Overlaid density plots
p1 <- ggplot(data, aes(x=anxiety, fill=group)) +
  geom_density(alpha=0.5) +
  geom_vline(data = data %>% group_by(group) %>% summarise(M=mean(anxiety)),
             aes(xintercept=M, color=group), linetype="dashed", linewidth=1) +
  labs(title="Distribution of Anxiety Scores by Group(Small Sample)",
       subtitle=paste0("N = ", total_n, " (Intervention n=", n_intervention, 
                       ", Control n=", n_control, ")"),
       x="STAI Anxiety Score(20-80)", y="Density") +
  scale_fill_brewer(palette="Set1") +
  scale_color_brewer(palette="Set1") +
  theme_classic() +
  theme(legend.position="top")
print(p1)

# 3.2 Boxplots with individual points (helpful for small n)
p2 <- ggplot(data, aes(x=group, y=anxiety, fill=group)) +
  geom_boxplot(alpha=0.6, width=0.5, outlier.shape=NA) +
  geom_jitter(width=0.15, alpha=0.6, size=2.5) +
  stat_summary(fun=mean, geom="point", size=4, color="red", shape=18) +
  labs(title="Anxiety Scores: Intervention vs Control",
       subtitle="Red diamonds = means; dots = individual observations",
       x="Group", y="STAI Anxiety Score(20-80)") +
  scale_fill_brewer(palette="Set1") +
  theme_classic() +
  theme(legend.position="none")
print(p2)

# 3.3 Q-Q plots for normality
par(mfrow=c(1,2))
qqnorm(data$anxiety[data$group == "Intervention"], main="Q-Q Plot: Intervention")
qqline(data$anxiety[data$group == "Intervention"], col="red")
qqnorm(data$anxiety[data$group == "Control"], main="Q-Q Plot: Control")
qqline(data$anxiety[data$group == "Control"], col="blue")
par(mfrow=c(1,1))

# ============================================================================
# STEP 4: CALCULATE COHEN'S D (UNCORRECTED)
# ============================================================================
cat("\nSTEP 4: COHEN'S D(UNCORRECTED)\n")
cat("-------------------------------\n")

M_intervention <- mean(data$anxiety[data$group == "Intervention"])
M_control <- mean(data$anxiety[data$group == "Control"])

# Pooled standard deviation
SD_pooled <- sqrt(((n_intervention-1)*SD_intervention^2 + 
                    (n_control-1)*SD_control^2) / 
                   (n_intervention + n_control - 2))

cohen_d <- (M_intervention - M_control) / SD_pooled

cat("Mean Intervention:", round(M_intervention, 2), "\n")
cat("Mean Control:     ", round(M_control, 2), "\n")
cat("Mean Difference:  ", round(M_intervention - M_control, 2), "\n")
cat("SD Intervention:  ", round(SD_intervention, 2), "\n")
cat("SD Control:       ", round(SD_control, 2), "\n")
cat("SD Pooled:        ", round(SD_pooled, 2), "\n")
cat("\nCohen's d(uncorrected):", round(cohen_d, 3), "\n")

# Interpretation
if (abs(cohen_d) < 0.2) {
  d_interpretation <- "negligible"
} else if (abs(cohen_d) < 0.5) {
  d_interpretation <- "small"
} else if (abs(cohen_d) < 0.8) {
  d_interpretation <- "medium"
} else {
  d_interpretation <- "large"
}
cat("Effect size magnitude:", d_interpretation, "(Cohen, 1988)\n\n")

# ============================================================================
# STEP 5: CALCULATE HEDGES' G (BIAS-CORRECTED)
# ============================================================================
cat("STEP 5: HEDGES' G(BIAS-CORRECTED)\n")
cat("-----------------------------------\n")

# Degrees of freedom
df <- n_intervention + n_control - 2

# Hedges' correction factor J
J <- 1 - (3 / (4*df - 1))

# Hedges' g
hedges_g <- cohen_d * J

# Bias magnitude
bias <- cohen_d - hedges_g
bias_percent <- (bias / cohen_d) * 100

cat("Degrees of freedom(df):", df, "\n")
cat("Correction factor J:     ", round(J, 5), "\n")
cat("\nCohen's d(uncorrected): ", round(cohen_d, 3), "\n")
cat("Hedges' g(corrected):   ", round(hedges_g, 3), "\n")
cat("\nBias reduction:          ", round(bias, 4), 
    " (", round(bias_percent, 2), "%)\n")

if (abs(cohen_d) < 0.2) {
  g_interpretation <- "negligible"
} else if (abs(cohen_d) < 0.5) {
  g_interpretation <- "small"
} else if (abs(cohen_d) < 0.8) {
  g_interpretation <- "medium"
} else {
  g_interpretation <- "large"
}
cat("Effect size magnitude:   ", g_interpretation, "(Cohen, 1988)\n\n")

cat("INTERPRETATION:\n")
cat("With small sample(N=", total_n, ", df=", df, "), the correction factor J=", 
    round(J, 3), "\n")
cat("reduces Cohen's d by", round(bias_percent, 1), 
    "%, yielding Hedges' g=", round(hedges_g, 2), ".\n")
cat("This correction is IMPORTANT to avoid overestimating the population effect.\n\n")

# ============================================================================
# STEP 6: CONFIDENCE INTERVALS
# ============================================================================
cat("STEP 6: CONFIDENCE INTERVALS\n")
cat("----------------------------\n")

# CI for Cohen's d (using effsize package)
cohen_result <- cohen.d(anxiety ~ group, data=data)
cat("\nCohen's d: ", round(cohen_result$estimate, 3), 
    ", 95% CI [", round(cohen_result$conf.int[1], 3), 
    ", ", round(cohen_result$conf.int[2], 3), "]\n")

# Approximate CI for Hedges' g (apply J correction to CI bounds)
# Note: This is an approximation; exact CI requires more complex calculations
hedges_ci_lower <- cohen_result$conf.int[1] * J
hedges_ci_upper <- cohen_result$conf.int[2] * J

cat("Hedges' g: ", round(hedges_g, 3), 
    ", 95% CI [", round(hedges_ci_lower, 3), 
    ", ", round(hedges_ci_upper, 3), "] (approximate)\n\n")

cat("INTERPRETATION:\n")
if (hedges_ci_lower * hedges_ci_upper > 0) {
  cat("CI does NOT include zero → Effect is statistically significant.\n")
} else {
  cat("CI includes zero → Effect is NOT statistically significant.\n")
  cat("Small sample may lack power to detect true effect.\n")
}

ci_width <- hedges_ci_upper - hedges_ci_lower
cat("CI width =", round(ci_width, 2), "→", 
    ifelse(ci_width < 0.5, "Precise estimate", 
           ifelse(ci_width < 1.0, "Moderate precision", 
                  "Wide CI(imprecise, need larger sample)")))
cat("\n\n")

# ============================================================================
# STEP 7: ALTERNATIVE EFFECT SIZES
# ============================================================================
cat("STEP 7: ALTERNATIVE EFFECT SIZES\n")
cat("---------------------------------\n")

# Glass's delta (standardize by control SD only)
glass_delta <- (M_intervention - M_control) / SD_control
glass_delta_corrected <- glass_delta * J  # Can also apply Hedges correction

cat("Glass's Δ (control SD):          ", round(glass_delta, 3), "\n")
cat("Glass's Δ (bias-corrected):      ", round(glass_delta_corrected, 3), "\n")

# Unstandardized difference with 95% CI
t_result <- t.test(anxiety ~ group, data=data)
cat("\nUnstandardized mean difference:  ", round(M_intervention - M_control, 2), 
    " STAI points\n")
cat("95% CI for mean difference:      [", round(t_result$conf.int[1], 2), 
    ", ", round(t_result$conf.int[2], 2), "]\n")
cat("t-test: t(", round(t_result$parameter, 1), ") = ", round(t_result$statistic, 2),
    ", p = ", round(t_result$p.value, 4), "\n\n")

# ============================================================================
# STEP 8: EFFECT SIZE COMPARISON VISUALIZATION
# ============================================================================
cat("STEP 8: EFFECT SIZE COMPARISON\n")
cat("------------------------------\n\n")

effect_sizes <- data.frame(
  Metric = c("Cohen's d\n(uncorrected)", 
             "Hedges' g\n(corrected)",
             "Glass's Δ\n(control SD)"),
  Value = c(abs(cohen_d), abs(hedges_g), abs(glass_delta)),
  CI_lower = c(abs(cohen_result$conf.int[1]), abs(hedges_ci_lower), NA),
  CI_upper = c(abs(cohen_result$conf.int[2]), abs(hedges_ci_upper), NA)
)

p3 <- ggplot(effect_sizes, aes(x=Metric, y=Value, fill=Metric)) +
  geom_bar(stat="identity", alpha=0.8, width=0.6) +
  geom_errorbar(aes(ymin=CI_lower, ymax=CI_upper), width=0.2, linewidth=1) +
  geom_hline(yintercept=c(0.2, 0.5, 0.8), linetype="dashed", 
             color="gray50", alpha=0.7) +
  annotate("text", x=3.3, y=0.2, label="Small(0.2)", hjust=0, size=3) +
  annotate("text", x=3.3, y=0.5, label="Medium(0.5)", hjust=0, size=3) +
  annotate("text", x=3.3, y=0.8, label="Large(0.8)", hjust=0, size=3) +
  labs(title="Effect Size Comparison: Cohen's d vs Hedges' g",
       subtitle=paste0("Small sample(N=", total_n, ") → ", 
                       round(bias_percent, 1), "% bias reduction with Hedges' g"),
       x="", y="Effect Size(absolute value)") +
  ylim(0, max(effect_sizes$Value, na.rm=TRUE) * 1.2) +
  scale_fill_brewer(palette="Set2") +
  theme_classic() +
  theme(legend.position="none",
        axis.text.x = element_text(size=10))
print(p3)

# ============================================================================
# STEP 9: COMPREHENSIVE SUMMARY
# ============================================================================
cat("\n======================================================================\n")
cat("COMPREHENSIVE SUMMARY\n")
cat("======================================================================\n\n")

cat("Sample: N =", total_n, "(Intervention n=", n_intervention, 
    ", Control n=", n_control, ")\n")
cat("Outcome: STAI Anxiety(20-80, higher = more anxiety)\n\n")

cat("DESCRIPTIVE STATISTICS:\n")
cat("Intervention: M =", round(M_intervention, 1), ", SD =", round(SD_intervention, 1), "\n")
cat("Control:      M =", round(M_control, 1), ", SD =", round(SD_control, 1), "\n")
cat("Difference:   M =", round(M_intervention - M_control, 1), " points\n\n")

cat("EFFECT SIZES:\n")
cat("Cohen's d(uncorrected):  ", round(cohen_d, 3), " (", d_interpretation, ")\n")
cat("Hedges' g(corrected):    ", round(hedges_g, 3), " (", g_interpretation, ")\n")
cat("Bias reduction:           ", round(bias, 4), " (", round(bias_percent, 1), "%)\n")
cat("95% CI for Hedges' g:     [", round(hedges_ci_lower, 2), ", ", 
    round(hedges_ci_upper, 2), "]\n\n")

cat("STATISTICAL SIGNIFICANCE:\n")
cat("t(", round(t_result$parameter, 0), ") = ", round(t_result$statistic, 2),
    ", p = ", round(t_result$p.value, 4))
if (t_result$p.value < 0.001) {
  cat(" ***\n")
} else if (t_result$p.value < 0.01) {
  cat(" **\n")
} else if (t_result$p.value < 0.05) {
  cat(" *\n")
} else {
  cat(" (ns)\n")
}

cat("\nINTERPRETATION:\n")
cat("The mindfulness intervention showed", 
    ifelse(abs(M_intervention - M_control) > 5, "substantially", "moderately"),
    "lower anxiety\n")
cat("compared to control(Hedges' g =", round(hedges_g, 2), 
    ", 95% CI [", round(hedges_ci_lower, 2), ", ", round(hedges_ci_upper, 2), "]).\n")
cat("This represents a", tolower(g_interpretation), "effect(Cohen, 1988).\n")
cat("The bias correction(J =", round(J, 3), ") reduced the uncorrected Cohen's d\n")
cat("by", round(bias_percent, 1), "%, which is IMPORTANT in this small sample(N=", 
    total_n, ").\n")
cat("Hedges' g is the appropriate effect size for meta-analysis and provides\n")
cat("an unbiased estimate of the population standardized mean difference.\n\n")

# ============================================================================
# STEP 10: APA-STYLE REPORTING
# ============================================================================
cat("======================================================================\n")
cat("APA-STYLE REPORT\n")
cat("======================================================================\n\n")

cat("Hedges' g was calculated to quantify the effect of mindfulness\n")
cat("intervention on anxiety, with bias correction applied due to small sample\n")
cat("size(N = ", total_n, "). The intervention group(M = ", round(M_intervention, 1),
    ", SD = ", round(SD_intervention, 1), ", n = ", n_intervention, ")\n")
cat("showed lower anxiety than the control group(M = ", round(M_control, 1),
    ", SD = ", round(SD_control, 1), ", n = ", n_control, "),\n")
cat("g = ", round(hedges_g, 2), ", 95% CI [", round(hedges_ci_lower, 2), ", ",
    round(hedges_ci_upper, 2), "]. This represents a ", tolower(g_interpretation),
    "\n")
cat("effect(Cohen, 1988). The bias correction factor(J = ", round(J, 3),
    ") reduced\n")
cat("Cohen's d by ", round(bias_percent, 1), "% (from d = ", round(cohen_d, 2),
    " to g = ", round(hedges_g, 2), "),\n")
cat("providing an unbiased estimate of the population effect size. The\n")
cat("intervention produced a clinically meaningful reduction of approximately\n")
cat(round(abs(M_intervention - M_control), 1), " points on the STAI(20-80 scale).\n\n")

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