Omega-Squared (ω²)

# Omega-squared effect size from one-way ANOVA
library(effectsize)
library(tidyverse)

# Simulate data
set.seed(2025)
data <- data.frame(
  medication = rep(c("Med_A", "Med_B", "Med_C"), each=25),
  bp_reduction = c(
    rnorm(25, 12, 6),  # Med A: M=12, SD=6
    rnorm(25, 18, 7),  # Med B: M=18, SD=7
    rnorm(25, 16, 6.5) # Med C: M=16, SD=6.5
  )
)

# ANOVA
model <- aov(bp_reduction ~ medication, data=data)
summary(model)

# Extract components for manual calculation
anova_summary <- summary(model)[[1]]
SS_effect <- anova_summary["medication", "Sum Sq"]
SS_error <- anova_summary["Residuals", "Sum Sq"]
SS_total <- SS_effect + SS_error
df_effect <- anova_summary["medication", "Df"]
df_error <- anova_summary["Residuals", "Df"]
MS_error <- anova_summary["Residuals", "Mean Sq"]

# Manual computation of ω²
eta_sq <- SS_effect / SS_total
omega_sq <- (SS_effect - df_effect * MS_error) / (SS_total + MS_error)
bias <- eta_sq - omega_sq

cat("=== Manual Calculation ===", "\n")
cat(sprintf("SS_effect = %.2f\n", SS_effect))
cat(sprintf("SS_error = %.2f\n", SS_error))
cat(sprintf("SS_total = %.2f\n", SS_total))
cat(sprintf("df_effect = %d\n", df_effect))
cat(sprintf("MS_error = %.2f\n\n", MS_error))

cat(sprintf("η² = %.3f (%.1f%% sample variance)\n", eta_sq, eta_sq*100))
cat(sprintf("ω² = %.3f (%.1f%% population variance - LESS BIASED)\n", omega_sq, omega_sq*100))
cat(sprintf("Bias = %.3f (η² overestimates by %.1f percentage points)\n\n", bias, bias*100))

# Using effectsize package (recommended)
cat("=== effectsize Package ===", "\n")
eta_result <- eta_squared(model, partial=FALSE)
omega_result <- omega_squared(model)
print(eta_result)
print(omega_result)

# Cohen's f for power analysis
cohen_f <- sqrt(omega_sq / (1 - omega_sq))
cat(sprintf("\nCohen's f = %.3f\n", cohen_f))

# Interpretation using Cohen (1988) benchmarks
cat("\n=== Interpretation ===", "\n")
cat("Cohen(1988) benchmarks: ω² small=.01, medium=.06, large=.14\n")
if (omega_sq >= .14) {
  magnitude <- "LARGE"
} else if (omega_sq >= .06) {
  magnitude <- "Medium"
} else if (omega_sq >= .01) {
  magnitude <- "Small"
} else {
  magnitude <- "Negligible"
}
cat(sprintf("Result: ω² = %.3f → %s effect\n", omega_sq, magnitude))

# Bootstrap 95% CI for ω²
library(boot)
boot_omega_sq <- function(data, indices) {
  d <- data[indices, ]
  model <- aov(bp_reduction ~ medication, data=d)
  aov_summary <- summary(model)[[1]]
  SS_eff <- aov_summary[1, "Sum Sq"]
  SS_err <- aov_summary[2, "Sum Sq"]
  SS_tot <- SS_eff + SS_err
  df_eff <- aov_summary[1, "Df"]
  MS_err <- aov_summary[2, "Mean Sq"]
  (SS_eff - df_eff * MS_err) / (SS_tot + MS_err)
}

set.seed(2025)
boot_results <- boot(data, boot_omega_sq, R=1000)
boot_ci <- boot.ci(boot_results, type="perc")

cat("\n=== Bootstrap 95% CI ===", "\n")
cat(sprintf("ω² = %.3f, 95%% CI [%.3f, %.3f]\n", 
            omega_sq, boot_ci$percent[4], boot_ci$percent[5]))

# Diagnostic plots
par(mfrow=c(2,2))
plot(model)
par(mfrow=c(1,1))

# Sensitivity analysis: effect size with/without outliers
std_resid <- rstandard(model)
outliers <- abs(std_resid) > 3
if (sum(outliers) > 0) {
  data_no_outliers <- data[!outliers, ]
  model_no_outliers <- aov(bp_reduction ~ medication, data=data_no_outliers)
  omega_no_outliers <- omega_squared(model_no_outliers)
  cat("\n=== Sensitivity Analysis ===", "\n")
  cat(sprintf("ω² with outliers: %.3f\n", omega_sq))
  cat(sprintf("ω² without %d outlier(s): %.3f\n", sum(outliers), omega_no_outliers$Omega2))
} else {
  cat("\nNo outliers detected(|standardized residual| > 3)\n")
}

# APA-style report
F_val <- anova_summary["medication", "F value"]
p_val <- anova_summary["medication", "Pr(>F)"]

cat("\n=== APA Report ===", "\n")
cat(sprintf("Medication type significantly affected blood pressure reduction, "))
cat(sprintf("F(%d, %d) = %.2f, p = %.3f, ω² = %.2f, 95%% CI [%.2f, %.2f].\n",
            df_effect, df_error, F_val, p_val, omega_sq, 
            boot_ci$percent[4], boot_ci$percent[5]))
cat(sprintf("This represents a %s effect, with medication explaining approximately %.0f%% ",
            tolower(magnitude), omega_sq*100))
cat("of the population variance in BP reduction.\n")