Glass's Delta (Δ)

# Glass's Delta: Educational Intervention with Unequal Variances
library(ggplot2)
library(dplyr)
library(car)  # Levene's test

# Simulate data reflecting unequal variances
set.seed(2025)
data <- data.frame(
  group = c(rep("New_Method", 45), rep("Standard", 40)),
  math_score = c(
    rnorm(45, mean=72.5, sd=18.2),   # New: M=72.5, SD=18.2 (higher variance)
    rnorm(40, mean=63.8, sd=12.5)    # Standard: M=63.8, SD=12.5 (baseline)
  )
)

# === STEP 1: Descriptive Statistics ===
cat("=== Descriptive Statistics ===\n")
desc_stats <- data %>%
  group_by(group) %>%
  summarise(
    n = n(),
    M = mean(math_score),
    SD = sd(math_score),
    Min = min(math_score),
    Max = max(math_score),
    CV = SD/M  # Coefficient of variation
  )
print(desc_stats)

# === STEP 2: Test Variance Heterogeneity ===
cat("\n=== Variance Heterogeneity Tests ===\n")

# Levene's test
levene_result <- leveneTest(math_score ~ group, data=data)
print(levene_result)

# Variance ratio
SD_new <- sd(data$math_score[data$group == "New_Method"])
SD_standard <- sd(data$math_score[data$group == "Standard"])
var_ratio <- SD_new^2 / SD_standard^2

cat("\nVariance ratio(New/Standard):", round(var_ratio, 2))
cat("\nSD ratio(New/Standard):", round(SD_new/SD_standard, 2))

if (var_ratio > 2 | var_ratio < 0.5) {
  cat("\n→ Substantial variance heterogeneity detected.")
  cat("\n→ Glass's Delta recommended(use control SD).\n")
} else {
  cat("\n→ Variances approximately equal.")
  cat("\n→ Cohen's d(pooled SD) may be preferred.\n")
}

# === STEP 3: Visual Comparison ===
ggplot(data, aes(x=math_score, fill=group)) +
  geom_density(alpha=0.5) +
  geom_vline(data = desc_stats, aes(xintercept=M, color=group),
             linetype="dashed", size=1.2) +
  labs(title="Distribution of Math Scores by Teaching Method",
       subtitle=paste0("Note: New Method shows higher variance(SD=", 
                      round(SD_new, 1), ") vs Standard(SD=", 
                      round(SD_standard, 1), ")"),
       x="Math Test Score(0-100)", y="Density") +
  scale_fill_brewer(palette="Set1") +
  scale_color_brewer(palette="Set1") +
  theme_classic() +
  theme(legend.title=element_blank())

# Boxplots to check outliers
ggplot(data, aes(x=group, y=math_score, fill=group)) +
  geom_boxplot(alpha=0.6, outlier.color="red", outlier.size=3) +
  geom_jitter(width=0.1, alpha=0.3) +
  stat_summary(fun=mean, geom="point", size=4, color="blue", shape=18) +
  labs(title="Math Scores by Group(Blue diamond = mean)",
       x="Teaching Method", y="Math Score(0-100)") +
  scale_fill_brewer(palette="Set2") +
  theme_classic() +
  theme(legend.position="none")

# === STEP 4: Calculate Glass's Delta ===
cat("\n=== Glass's Delta Calculation ===\n")

M_new <- mean(data$math_score[data$group == "New_Method"])
M_standard <- mean(data$math_score[data$group == "Standard"])
n_new <- sum(data$group == "New_Method")
n_standard <- sum(data$group == "Standard")

# Glass's Delta: Standardize by CONTROL (Standard) SD only
glass_delta <- (M_new - M_standard) / SD_standard

cat("Mean New Method:", round(M_new, 2))
cat("\nMean Standard:", round(M_standard, 2))
cat("\nMean Difference:", round(M_new - M_standard, 2))
cat("\nControl SD(Standard):", round(SD_standard, 2))
cat("\nGlass's Δ:", round(glass_delta, 2), "\n")

# Interpretation
if (abs(glass_delta) < 0.2) {
  interpretation <- "negligible"
} else if (abs(glass_delta) < 0.5) {
  interpretation <- "small"
} else if (abs(glass_delta) < 0.8) {
  interpretation <- "medium"
} else {
  interpretation <- "large"
}
cat("\nEffect size:", interpretation, "(Cohen's benchmarks)\n")

# === STEP 5: Confidence Interval for Glass's Delta ===
# Approximate CI using SE formula
SE_delta <- sqrt((n_new + n_standard)/(n_new * n_standard) + 
                 glass_delta^2 / (2*n_standard))
CI_lower <- glass_delta - 1.96 * SE_delta
CI_upper <- glass_delta + 1.96 * SE_delta

cat("\n95% CI for Glass's Δ: [", round(CI_lower, 2), ",", 
    round(CI_upper, 2), "]\n")

# === STEP 6: Compare with Cohen's d ===
cat("\n=== Comparison with Cohen's d ===\n")

# Pooled SD for Cohen's d
SD_pooled <- sqrt(((n_new-1)*SD_new^2 + (n_standard-1)*SD_standard^2) / 
                    (n_new + n_standard - 2))
cohens_d <- (M_new - M_standard) / SD_pooled

cat("Cohen's d(pooled SD):", round(cohens_d, 2))
cat("\nGlass's Δ (control SD):", round(glass_delta, 2))
cat("\nDifference:", round(abs(cohens_d - glass_delta), 3))

cat("\n\nInterpretation:")
cat("\n- Glass's Δ =", round(glass_delta, 2), "uses only control SD =", 
    round(SD_standard, 2))
cat("\n- Cohen's d =", round(cohens_d, 2), "uses pooled SD =", 
    round(SD_pooled, 2))
cat("\n- With unequal variances, Glass's Δ preferred(stable control baseline)\n")

# === STEP 7: Multiple Treatment Groups Example ===
cat("\n=== Extension: Multiple Treatments vs. Control ===\n")

# Add a third treatment group
set.seed(2026)
data_multi <- rbind(
  data,
  data.frame(
    group = rep("Alternative_Method", 38),
    math_score = rnorm(38, mean=68.2, sd=15.7)
  )
)

# Calculate Glass's Delta for each treatment vs. control
treatments <- c("New_Method", "Alternative_Method")
control_mean <- mean(data_multi$math_score[data_multi$group == "Standard"])
control_sd <- sd(data_multi$math_score[data_multi$group == "Standard"])

results <- data.frame(
  Treatment = treatments,
  Delta = numeric(2),
  Interpretation = character(2),
  stringsAsFactors = FALSE
)

for (i in 1:length(treatments)) {
  treat_mean <- mean(data_multi$math_score[data_multi$group == treatments[i]])
  delta <- (treat_mean - control_mean) / control_sd
  results$Delta[i] <- delta
  
  if (abs(delta) < 0.2) interp <- "negligible"
  else if (abs(delta) < 0.5) interp <- "small"
  else if (abs(delta) < 0.8) interp <- "medium"
  else interp <- "large"
  
  results$Interpretation[i] <- interp
}

cat("\nGlass's Delta for each treatment(vs. Standard control):\n")
print(results)

cat("\nNote: All treatments standardized by same control SD =", 
    round(control_sd, 2))
cat("\nThis allows direct comparison across treatments.\n")

# === STEP 8: Sensitivity Analysis ===
cat("\n=== Sensitivity Analysis: Outlier Impact ===\n")

# Check for outliers in control group
control_data <- data$math_score[data$group == "Standard"]
control_z <- scale(control_data)
outliers <- which(abs(control_z) > 2.5)

if (length(outliers) > 0) {
  cat("\nOutliers detected in control group(|z| > 2.5):", length(outliers))
  
  # Recalculate without outliers
  control_no_outliers <- control_data[-outliers]
  SD_standard_robust <- sd(control_no_outliers)
  glass_delta_robust <- (M_new - M_standard) / SD_standard_robust
  
  cat("\nOriginal Glass's Δ:", round(glass_delta, 2), 
      "(SD =", round(SD_standard, 2), ")")
  cat("\nRobust Glass's Δ (outliers removed):", round(glass_delta_robust, 2),
      "(SD =", round(SD_standard_robust, 2), ")")
  cat("\nDifference:", round(abs(glass_delta - glass_delta_robust), 3), "\n")
} else {
  cat("\nNo extreme outliers in control group. Effect size stable.\n")
}

# === APA-Style Reporting ===
cat("\n=== APA-Style Report ===\n")
cat("Glass's Δ was calculated to quantify the magnitude of the new teaching\n")
cat("method's effect on mathematics achievement, standardized by the control\n")
cat("group's variability. Levene's test indicated significant variance\n")
cat("heterogeneity(F =", round(levene_result$`F value`[1], 2), 
    ", p =", round(levene_result$`Pr(>F)`[1], 3), "),\n")
cat("justifying Glass's Δ over Cohen's d. The new method group\n")
cat("(M =", round(M_new, 1), ", SD =", round(SD_new, 1), ", n =", n_new, ")\n")
cat("scored", round(M_new - M_standard, 1), "points higher than the standard\n")
cat("instruction group(M =", round(M_standard, 1), ", SD =", 
    round(SD_standard, 1), ", n =", n_standard, "),\n")
cat("Δ =", round(glass_delta, 2), ", 95% CI [", round(CI_lower, 2), ",", 
    round(CI_upper, 2), "].\n")
cat("This represents a", interpretation, "effect when standardized by the\n")
cat("control group's baseline variability, indicating the new method produced\n")
cat("meaningful improvement in mathematics achievement.\n")