Wald-Wolfowitz Runs Test (Two-Sample)

# ==============================================================================
# Wald-Wolfowitz Runs Test: Reaction Time Comparison
# Research Question: Do two training methods produce different RT distributions?
# ==============================================================================

library(tidyverse)      # Data manipulation and visualization
library(randtests)      # runs.test() for Wald-Wolfowitz
library(DescTools)      # WaldWolfTest() - alternative implementation
library(effectsize)     # Effect size calculations
library(ggpubr)         # Publication-ready plots

# Generate simulated data (n = 50, 25 per group)
set.seed(2025)
data <- data.frame(
  method = rep(c("Method_A", "Method_B"), each=25),
  reaction_time = c(
    rgamma(25, shape=4, scale=50),    # Method A: faster, less variable
    rgamma(25, shape=3, scale=70)     # Method B: slower, more variable
  )
)

data$method <- factor(data$method, levels = c("Method_A", "Method_B"))

# ==============================================================================
# STEP 1: DESCRIPTIVE STATISTICS
# ==============================================================================

cat("\n=== DESCRIPTIVE STATISTICS ===\n")

descriptives <- data %>%
  group_by(method) %>%
  summarise(
    n = n(),
    mean = mean(reaction_time),
    sd = sd(reaction_time),
    median = median(reaction_time),
    IQR = IQR(reaction_time),
    min = min(reaction_time),
    max = max(reaction_time)
  )
print(descriptives)

# ==============================================================================
# STEP 2: ASSUMPTION CHECKING
# ==============================================================================

cat("\n=== ASSUMPTION CHECKS ===\n")

# Assumption 1: Independence
cat("1. Independence: Design review confirms between-subjects(independent groups)\n")

# Assumption 2: Sample size (n₁, n₂ ≥ 10 for normal approximation)
cat("\n2. Sample Size Check:\n")
cat(sprintf("   Method A: n = %d(≥10 ✓)\n", sum(data$method == "Method_A")))
cat(sprintf("   Method B: n = %d(≥10 ✓)\n", sum(data$method == "Method_B")))
cat("   Normal approximation is valid.\n")

# Assumption 3: Check for ties
cat("\n3. Tied Values Check:\n")
n_ties <- sum(duplicated(data$reaction_time))
pct_ties <- (n_ties / nrow(data)) * 100
cat(sprintf("   Ties: %d(%.1f%% of data)\n", n_ties, pct_ties))
if(pct_ties < 10) {
  cat("   Ties are minimal(<10%), test is appropriate.\n")
} else {
  cat("   Warning: >10% ties may reduce power. Consider Mann-Whitney U.\n")
}

# ==============================================================================
# STEP 3: VISUALIZATIONS
# ==============================================================================

cat("\n=== GENERATING VISUALIZATIONS ===\n")

# Visualization 1: Boxplots
p1 <- ggplot(data, aes(x=method, y=reaction_time, fill=method)) +
  geom_boxplot(alpha=0.7) +
  geom_jitter(width=0.2, alpha=0.4, size=2) +
  labs(title="Reaction Time by Training Method",
       x="Training Method", y="Reaction Time(ms)") +
  theme_minimal() +
  theme(legend.position="none")
print(p1)

# Visualization 2: Density plots (check location AND scale differences)
p2 <- ggplot(data, aes(x=reaction_time, fill=method)) +
  geom_density(alpha=0.5) +
  labs(title="Distribution Comparison: Density Plots",
       subtitle="Runs test detects differences in location, scale, and shape",
       x="Reaction Time(ms)", y="Density") +
  theme_minimal()
print(p2)

# Visualization 3: ECDF plot (empirical cumulative distribution)
p3 <- ggplot(data, aes(x=reaction_time, color=method)) +
  stat_ecdf(geom="step", size=1.2) +
  labs(title="Empirical Cumulative Distribution Functions",
       subtitle="Separation indicates distribution differences",
       x="Reaction Time(ms)", y="Cumulative Probability") +
  theme_minimal()
print(p3)

# Visualization 4: Run sequence plot
cat("\n4. Creating Run Sequence Plot...\n")
# Combine and rank data
data_combined <- data %>%
  arrange(reaction_time) %>%
  mutate(
    rank = row_number(),
    group_code = ifelse(method == "Method_A", "A", "B")
  )

p4 <- ggplot(data_combined, aes(x=rank, y=reaction_time, color=method, shape=method)) +
  geom_point(size=3, alpha=0.7) +
  geom_line(aes(group=1), color="gray80", size=0.5) +
  labs(title="Run Sequence: Ordered Data by Rank",
       subtitle="Clustering of colors indicates distribution differences",
       x="Rank(ordered by RT)", y="Reaction Time(ms)") +
  theme_minimal()
print(p4)

# ==============================================================================
# STEP 4: WALD-WOLFOWITZ RUNS TEST
# ==============================================================================

cat("\n=== WALD-WOLFOWITZ RUNS TEST ===\n")

# Extract group data
group_a <- data %>% filter(method == "Method_A") %>% pull(reaction_time)
group_b <- data %>% filter(method == "Method_B") %>% pull(reaction_time)

# Method 1: randtests package
cat("\nMethod 1: randtests::runs.test\n")
runs_result <- runs.test(x=group_a, y=group_b, alternative="two.sided")
print(runs_result)

# Method 2: DescTools package (alternative)
cat("\nMethod 2: DescTools::WaldWolfTest\n")
ww_result <- WaldWolfTest(x=group_a, y=group_b)
print(ww_result)

# Extract key statistics
n1 <- length(group_a)
n2 <- length(group_b)
runs_observed <- runs_result$statistic
runs_expected <- 1 + (2*n1*n2)/(n1+n2)
runs_sd <- sqrt((2*n1*n2*(2*n1*n2 - n1 - n2)) / ((n1+n2)^2 * (n1+n2-1)))
z_stat <- (runs_observed - runs_expected) / runs_sd
p_value <- runs_result$p.value

cat("\n--- Test Summary ---\n")
cat(sprintf("n₁ (Method A): %d\n", n1))
cat(sprintf("n₂ (Method B): %d\n", n2))
cat(sprintf("Runs observed: %.0f\n", runs_observed))
cat(sprintf("Runs expected(H₀): %.2f\n", runs_expected))
cat(sprintf("Standard deviation: %.2f\n", runs_sd))
cat(sprintf("Z-statistic: %.3f\n", z_stat))
cat(sprintf("p-value: %.4f\n", p_value))
cat(sprintf("Result: %s at α = .05\n", 
    ifelse(p_value < 0.05, "SIGNIFICANT", "Not significant")))

# ==============================================================================
# STEP 5: EFFECT SIZE
# ==============================================================================

cat("\n=== EFFECT SIZE ===\n")

# Standardized runs effect size
standardized_runs <- (runs_observed - runs_expected) / runs_sd
cat(sprintf("Standardized runs statistic: %.3f\n", standardized_runs))
cat("Interpretation: |z| < 1.96 (small), 1.96-2.58 (medium), > 2.58 (large)\n")

# Cohen's d for comparison (location effect)
cohen_d <- (mean(group_a) - mean(group_b)) / 
           sqrt(((n1-1)*var(group_a) + (n2-1)*var(group_b)) / (n1+n2-2))
cat(sprintf("\nCohen's d(for comparison): %.3f\n", cohen_d))
cat("Note: Runs test is sensitive to ALL distribution differences, not just location.\n")

# Overlap coefficient
min_max_a <- max(group_a)
max_min_b <- min(group_b)
overlap <- sum(group_a > max_min_b & group_a < min_max_a) / (n1 + n2)
cat(sprintf("\nDistribution overlap: %.1f%%\n", overlap*100))

# ==============================================================================
# STEP 6: COMPARE WITH MANN-WHITNEY U TEST
# ==============================================================================

cat("\n=== COMPARISON WITH MANN-WHITNEY U TEST ===\n")
cat("Mann-Whitney tests primarily location shift; Runs test detects ALL differences.\n\n")

mw_result <- wilcox.test(group_a, group_b, alternative="two.sided")
cat("Mann-Whitney U Test:\n")
cat(sprintf("W = %.1f, p = %.4f\n", mw_result$statistic, mw_result$p.value))

cat("\nInterpretation:\n")
if(p_value < 0.05 && mw_result$p.value < 0.05) {
  cat("Both tests significant: distributions differ(likely in location AND scale/shape).\n")
} else if(p_value < 0.05 && mw_result$p.value >= 0.05) {
  cat("Runs significant, M-W not: distributions differ in scale/shape, not location.\n")
} else if(p_value >= 0.05 && mw_result$p.value < 0.05) {
  cat("M-W significant, Runs not: location shift present, but distributions similarly shaped.\n")
} else {
  cat("Both tests non-significant: distributions are similar.\n")
}

# ==============================================================================
# STEP 7: INTERPRETATION
# ==============================================================================

cat("\n=== INTERPRETATION ===\n")
cat(sprintf(
  "The Wald-Wolfowitz runs test showed %s difference between the two training\n",
  ifelse(p_value < 0.05, "a significant", "no significant")))
cat(sprintf(
  "methods, Z = %.3f, p = %.4f. The observed number of runs(%.0f) was %s\n",
  z_stat, p_value, runs_observed,
  ifelse(runs_observed < runs_expected, "fewer than expected", "similar to expected")))
cat(sprintf(
  "under the null hypothesis(expected: %.2f), indicating %s.\n",
  runs_expected,
  ifelse(runs_observed < runs_expected, "clustering(distribution differences)", "random intermixing")))

if(p_value < 0.05) {
  cat("\nThe distributions differ in location(Method A faster: M=%.1f vs M=%.1f)\n",
      mean(group_a), mean(group_b))
  cat("and/or variability(Method A SD=%.1f vs Method B SD=%.1f).\n",
      sd(group_a), sd(group_b))
  cat("This comprehensive distribution difference(detected by runs test) suggests\n")
  cat("Method A produces more consistent, faster reaction times than Method B.\n")
}

# ==============================================================================
# APA REPORTING TEMPLATE
# ==============================================================================

cat("\n=== APA-STYLE REPORTING ===\n")
cat("A Wald-Wolfowitz runs test was conducted to compare reaction time distributions\n")
cat("between two training methods(Method A: n=25, Method B: n=25). The test evaluates\n")
cat("whether the two samples are randomly intermixed(null hypothesis) or systematically\n")
cat(sprintf("clustered(alternative). Results showed %s, Z = %.3f,\n",
    ifelse(p_value < 0.05, "a significant difference", "no significant difference"),
    z_stat))
cat(sprintf("p = %.4f, with %.0f runs observed versus %.2f expected under random mixing.\n",
    p_value, runs_observed, runs_expected))
if(p_value < 0.05) {
  cat(sprintf("Method A(M=%.1f ms, SD=%.1f) produced faster and less variable reaction times\n",
      mean(group_a), sd(group_a)))
  cat(sprintf("than Method B(M=%.1f ms, SD=%.1f), d = %.2f (large effect). The runs test detected\n",
      mean(group_b), sd(group_b), abs(cohen_d)))
  cat("comprehensive distribution differences beyond location shift alone, suggesting\n")
  cat("Method A produces superior and more consistent performance.\n")
}