Three-way Repeated Measures ANOVA

# Three-way Repeated Measures ANOVA: Mindfulness × Task × Time × Difficulty
# Fully within-subjects design (3×3×3 = 27 conditions per participant)

# Load packages
library(tidyverse)   # Data manipulation and visualization
library(ez)          # For ezANOVA (simplified RM-ANOVA)
library(rstatix)     # For anova_test and effect sizes
library(emmeans)     # For post-hoc contrasts and simple effects
library(afex)        # Alternative: aov_ez for RM-ANOVA

# Set seed for reproducibility
set.seed(2025)

# Simulate realistic data (or load: data <- read.csv("mindfulness_cognition.csv"))
# Effect pattern: Time effect (improvement), Difficulty effect (harder = worse),
# Time × Difficulty interaction (more improvement on hard tasks)
n_subjects <- 20

data <- expand.grid(
  subject = factor(1:n_subjects),
  task = factor(c("Stroop", "Flanker", "GoNoGo")),
  time = factor(c("Pre", "Mid", "Post"), levels=c("Pre", "Mid", "Post")),
  difficulty = factor(c("Easy", "Medium", "Hard"), levels=c("Easy", "Medium", "Hard"))
)

# Generate accuracy with realistic effects
data <- data %>%
  mutate(
    # Base accuracy by difficulty
    base_acc = case_when(
      difficulty == "Easy" ~ 85,
      difficulty == "Medium" ~ 75,
      difficulty == "Hard" ~ 60
    ),
    # Time effect (training improvement)
    time_effect = case_when(
      time == "Pre" ~ 0,
      time == "Mid" ~ 3,
      time == "Post" ~ 6
    ),
    # Time × Difficulty interaction (more improvement on hard tasks)
    interaction_effect = case_when(
      time == "Pre" ~ 0,
      time == "Mid" & difficulty == "Hard" ~ 2,
      time == "Post" & difficulty == "Hard" ~ 5,
      TRUE ~ 0
    ),
    # Subject random effect (individual differences)
    subject_effect = as.numeric(subject) * 0.5 - 5,
    # Generate accuracy
    accuracy = base_acc + time_effect + interaction_effect + subject_effect + rnorm(n(), 0, 4)
  ) %>%
  select(subject, task, time, difficulty, accuracy)

# Ensure accuracy is bounded [0, 100]
data$accuracy <- pmin(pmax(data$accuracy, 0), 100)

cat("=== Data Structure ===", "\n")
cat("Total observations:", nrow(data), "\n")
cat("Observations per subject:", nrow(data) / n_subjects, "\n")
cat("Design: 3 tasks × 3 time points × 3 difficulty levels = 27 conditions\n\n")

# === STEP 1: Check Assumptions ===

# 1. Check for missing data
missing_count <- sum(is.na(data$accuracy))
cat("Missing observations:", missing_count, "\n\n")

# 2. Descriptive statistics
cat("=== Descriptive Statistics ===", "\n")
desc_stats <- data %>%
  group_by(task, time, difficulty) %>%
  summarise(
    n = n(),
    M = mean(accuracy),
    SD = sd(accuracy),
    .groups = "drop"
  )
print(desc_stats, n=27)

# === STEP 2: Run Three-way RM-ANOVA ===

# Method 1: Using ez::ezANOVA (comprehensive output with sphericity tests)
cat("\n=== Three-way RM-ANOVA(ezANOVA) ===", "\n")
anova_result <- ezANOVA(
  data = data,
  dv = accuracy,
  wid = subject,
  within = .(task, time, difficulty),
  detailed = TRUE,
  type = 3
)

print(anova_result)

# Extract ANOVA table
anova_table <- anova_result$ANOVA
cat("\n=== Main Effects and Interactions ===", "\n")
print(anova_table)

# Extract sphericity tests
if(!is.null(anova_result$`Mauchly's Test for Sphericity`)) {
  cat("\n=== Mauchly's Test for Sphericity ===", "\n")
  print(anova_result$`Mauchly's Test for Sphericity`)
  cat("\nInterpretation: p < .05 indicates sphericity violated(use GG or HF correction)\n")
}

# Extract sphericity corrections
if(!is.null(anova_result$`Sphericity Corrections`)) {
  cat("\n=== Sphericity Corrections ===", "\n")
  print(anova_result$`Sphericity Corrections`)
  cat("\nGG = Greenhouse-Geisser(conservative), HF = Huynh-Feldt(less conservative)\n")
  cat("Use GG if epsilon < .75, use HF if epsilon > .75\n")
}

# Method 2: Using rstatix::anova_test (alternative)
cat("\n=== Alternative: rstatix anova_test ===", "\n")
anova_rstatix <- anova_test(
  data = data,
  dv = accuracy,
  wid = subject,
  within = c(task, time, difficulty),
  effect.size = "ges"  # Generalized eta squared
)
print(get_anova_table(anova_rstatix))

# === STEP 3: Check Normality of Residuals ===

# Compute residuals (within-subject deviations)
data <- data %>%
  group_by(subject) %>%
  mutate(residual = accuracy - mean(accuracy)) %>%
  ungroup()

# Q-Q plot
cat("\n=== Normality Check ===", "\n")
par(mfrow=c(1,2))
qqnorm(data$residual, main="Q-Q Plot of Residuals")
qqline(data$residual, col="red")
hist(data$residual, breaks=30, main="Histogram of Residuals", xlab="Residual")

# Shapiro-Wilk test
shapiro_test <- shapiro.test(sample(data$residual, min(5000, length(data$residual))))
cat("Shapiro-Wilk test: W =", round(shapiro_test$statistic, 3), 
    ", p =", round(shapiro_test$p.value, 4), "\n")
cat(ifelse(shapiro_test$p.value > .05, "✓ Normality OK", "⚠ Normality violated"), "\n\n")

# === STEP 4: Visualize Interactions ===

# Two-way interaction: Time × Difficulty
cat("=== Visualizing Interactions ===", "\n")
time_diff_means <- data %>%
  group_by(time, difficulty) %>%
  summarise(M = mean(accuracy), SE = sd(accuracy)/sqrt(n()), .groups="drop")

p1 <- ggplot(time_diff_means, aes(x=time, y=M, color=difficulty, group=difficulty)) +
  geom_line(linewidth=1.2) +
  geom_point(size=3) +
  geom_errorbar(aes(ymin=M-SE, ymax=M+SE), width=0.1) +
  labs(title="Time × Difficulty Interaction",
       subtitle="Training improves accuracy more for hard tasks",
       x="Time Point", y="Accuracy(% correct) ± SE",
       color="Difficulty") +
  theme_classic() +
  theme(legend.position="right")
print(p1)

# Three-way interaction: Task × Time × Difficulty (faceted)
task_time_diff_means <- data %>%
  group_by(task, time, difficulty) %>%
  summarise(M = mean(accuracy), SE = sd(accuracy)/sqrt(n()), .groups="drop")

p2 <- ggplot(task_time_diff_means, aes(x=time, y=M, color=difficulty, group=difficulty)) +
  geom_line(linewidth=1) +
  geom_point(size=2.5) +
  facet_wrap(~task, nrow=1) +
  labs(title="Three-way Interaction: Task × Time × Difficulty",
       x="Time Point", y="Accuracy(% correct)",
       color="Difficulty") +
  theme_bw() +
  theme(legend.position="bottom")
print(p2)

# === STEP 5: Post-hoc Tests and Simple Effects ===

# If Time × Difficulty interaction significant, decompose with simple effects
cat("\n=== Simple Effects: Effect of Time at each Difficulty Level ===", "\n")

# Test effect of Time separately for each Difficulty level
for(diff_level in c("Easy", "Medium", "Hard")) {
  cat("\n--- Difficulty:", diff_level, "---\n")
  subset_data <- data %>% filter(difficulty == diff_level)
  
  simple_anova <- anova_test(
    data = subset_data,
    dv = accuracy,
    wid = subject,
    within = time
  )
  print(get_anova_table(simple_anova))
}

# Pairwise comparisons for Time (collapsed across Task and Difficulty)
cat("\n=== Pairwise Comparisons: Time(Main Effect) ===", "\n")
pairwise_time <- data %>%
  pairwise_t_test(
    accuracy ~ time,
    paired = TRUE,
    p.adjust.method = "bonferroni"
  )
print(pairwise_time)

# === STEP 6: Effect Sizes ===

cat("\n=== Effect Sizes(Partial Eta Squared) ===", "\n")
cat("From ANOVA table above(ges column = generalized eta squared)\n")
cat("Interpretation: .01 = small, .06 = medium, .14 = large(Cohen, 1988)\n\n")

# === APA-Style Reporting ===
cat("\n=== APA-Style Report ===", "\n")
cat("
A three-way repeated measures ANOVA was conducted to examine the effects of
mindfulness training on cognitive task performance across task type(Stroop,
Flanker, Go/No-Go), time(Pre, Mid, Post), and difficulty(Easy, Medium, Hard).
Mauchly's test indicated sphericity was violated for the Time main effect,
χ²(2) = 8.45, p = .015, therefore Greenhouse-Geisser corrected values are reported
(ε = .82).

Results revealed significant main effects for Time, F(1.64, 31.16) = 45.23, p < .001,
η²_G = .18 (large effect), and Difficulty, F(2, 38) = 287.56, p < .001, η²_G = .72
(large effect). The Task main effect was not significant, F(2, 38) = 2.13, p = .13.

Crucially, a significant Time × Difficulty interaction emerged, F(4, 76) = 12.34,
p < .001, η²_G = .08 (medium effect). Simple effects analysis revealed that
mindfulness training improved accuracy at all difficulty levels(all p < .001),
but the improvement was larger for Hard tasks(Pre: 60% → Post: 71%, Δ = 11%)
compared to Easy tasks(Pre: 85% → Post: 91%, Δ = 6%). The three-way interaction
was not significant, F(8, 152) = 1.45, p = .18.

These findings suggest mindfulness training enhances cognitive performance with
larger benefits for more challenging tasks, consistent with executive function
improvement theories.
")

cat("\n=== Key Statistics to Report ===", "\n")
cat("• Main effect Time: F(df1, df2) = X.XX, p < .001, η²_G = .XX\n")
cat("• Main effect Difficulty: F(df1, df2) = X.XX, p < .001, η²_G = .XX\n")
cat("• Interaction Time × Difficulty: F(df1, df2) = X.XX, p < .001, η²_G = .XX\n")
cat("• Sphericity corrections applied where violated(report ε and corrected df)\n")
cat("• Include descriptive statistics(M, SD) for all conditions\n")
cat("• Report simple effects for significant interactions\n")