Eta-Squared (η²)

# Eta-squared effect sizes from one-way ANOVA
library(effectsize)
library(tidyverse)

# Simulate data
set.seed(2025)
data <- data.frame(
  method = rep(c("Lecture", "Active", "Flipped"), each=30),
  score = c(
    rnorm(30, 72, 12),  # Lecture: M=72, SD=12
    rnorm(30, 78, 11),  # Active: M=78, SD=11
    rnorm(30, 82, 10)   # Flipped: M=82, SD=10
  )
)

# ANOVA
model <- aov(score ~ method, data=data)
summary(model)

# Extract SS values
anova_summary <- summary(model)[[1]]
SS_effect <- anova_summary["method", "Sum Sq"]
SS_error <- anova_summary["Residuals", "Sum Sq"]
SS_total <- SS_effect + SS_error
df_effect <- anova_summary["method", "Df"]
MS_error <- anova_summary["Residuals", "Mean Sq"]

# Manual computation
eta_sq <- SS_effect / SS_total
partial_eta_sq <- SS_effect / (SS_effect + SS_error)  # Same as η² in one-way
omega_sq <- (SS_effect - df_effect * MS_error) / (SS_total + MS_error)

cat("=== Manual Calculation ===", "\n")
cat(sprintf("η² = %.3f (%.1f%% variance explained)\n", eta_sq, eta_sq*100))
cat(sprintf("ηp² = %.3f\n", partial_eta_sq))
cat(sprintf("ω² = %.3f (less biased estimate)\n", omega_sq))

# Using effectsize package (recommended)
cat("\n=== effectsize Package ===", "\n")
eta_squared(model, partial=FALSE)  # η²
eta_squared(model, partial=TRUE)   # ηp²
omega_squared(model)               # ω²

# Interpretation
cat("\n=== Cohen(1988) Benchmarks ===", "\n")
cat("Small: .01, Medium: .06, Large: .14\n")
cat(sprintf("Result: ω² = %.3f → %s effect\n", 
            omega_sq,
            ifelse(omega_sq >= .14, "LARGE",
                   ifelse(omega_sq >= .06, "Medium", "Small"))))

# Bootstrap CI for η²
library(boot)
boot_eta_sq <- function(data, indices) {
  d <- data[indices, ]
  model <- aov(score ~ method, data=d)
  SS <- summary(model)[[1]][, "Sum Sq"]
  SS[1] / sum(SS)
}

boot_results <- boot(data, boot_eta_sq, R=1000)
boot_ci <- boot.ci(boot_results, type="perc")
cat("\n95% CI for η²:", boot_ci$percent[4:5], "\n")

# APA report
cat("\n=== APA Report ===", "\n")
cat(sprintf("F(2, 87) = %.2f, p < .001, ω² = %.2f (large effect, 95%% CI [%.2f, %.2f])\n",
            anova_summary["method", "F value"],
            omega_sq,
            boot_ci$percent[4], boot_ci$percent[5]))