Within-Within Subjects ANOVA

# Two-way RM-ANOVA: Task × Time on Cognitive Performance
# Fully within-subjects design (2×3 = 6 conditions per participant)

# Load packages
library(tidyverse)   # Data manipulation and ggplot2
library(ez)          # For ezANOVA
library(rstatix)     # For anova_test
library(emmeans)     # For post-hoc contrasts

# Set seed
set.seed(2025)

# Simulate realistic data
n_subjects <- 30

data <- expand.grid(
  subject = factor(1:n_subjects),
  task = factor(c("Memory", "Attention")),
  time = factor(c("30min", "60min", "90min"), levels = c("30min", "60min", "90min"))
)

# Generate accuracy with realistic effects
data <- data %>%
  mutate(
    # Base accuracy by task
    base_acc = ifelse(task == "Memory", 72, 75),
    
    # Time effect (caffeine peaks at 60min)
    time_effect = case_when(
      time == "30min" ~ 0,
      time == "60min" ~ 5,
      time == "90min" ~ 3
    ),
    
    # Task × Time interaction (caffeine more beneficial for Attention)
    interaction_effect = case_when(
      task == "Attention" & time == "60min" ~ 4,
      task == "Attention" & time == "90min" ~ 2,
      TRUE ~ 0
    ),
    
    # Subject random effect
    subject_effect = as.numeric(subject) * 0.4 - 6,
    
    # Generate accuracy
    accuracy = base_acc + time_effect + interaction_effect + subject_effect + rnorm(n(), 0, 3.5)
  ) %>%
  mutate(accuracy = pmin(pmax(accuracy, 0), 100)) %>%
  select(subject, task, time, accuracy)

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

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

# 1. Missing data
cat("Missing observations:", sum(is.na(data$accuracy)), "\n\n")

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

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

cat("\n=== Two-way RM-ANOVA(ezANOVA) ===", "\n")
anova_result <- ezANOVA(
  data = data,
  dv = accuracy,
  wid = subject,
  within = .(task, time),
  detailed = TRUE,
  type = 3
)

print(anova_result)

# Check sphericity
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\n")
}

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

# === STEP 3: Check Normality ===

cat("=== Normality Check ===", "\n")
# Compute within-subject residuals
data <- data %>%
  group_by(subject) %>%
  mutate(residual = accuracy - mean(accuracy)) %>%
  ungroup()

# Q-Q plot and histogram
par(mfrow=c(1,2))
qqnorm(data$residual, main="Q-Q Plot of Residuals")
qqline(data$residual, col="red")
hist(data$residual, breaks=25, main="Histogram of Residuals", xlab="Residual")

# Shapiro-Wilk test
shapiro_test <- shapiro.test(data$residual)
cat("\nShapiro-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 Interaction ===

cat("=== Interaction Plot ===", "\n")
task_time_means <- data %>%
  group_by(task, time) %>%
  summarise(M = mean(accuracy), SE = sd(accuracy)/sqrt(n()), .groups="drop")

p1 <- ggplot(task_time_means, aes(x=time, y=M, color=task, group=task)) +
  geom_line(linewidth=1.2) +
  geom_point(size=3) +
  geom_errorbar(aes(ymin=M-SE, ymax=M+SE), width=0.1) +
  labs(title="Task × Time Interaction: Caffeine Effects on Performance",
       subtitle="Caffeine more beneficial for Attention at peak(60min)",
       x="Time Since Caffeine Ingestion", 
       y="Accuracy(% correct) ± SE",
       color="Task Type") +
  theme_classic() +
  theme(legend.position="right")
print(p1)

# === STEP 5: Simple Effects Analysis (if interaction significant) ===

if(anova_result$ANOVA$`Pr(>F)`[3] < .05) {  # If interaction p < .05
  cat("\n=== Simple Effects: Effect of Time for each Task ===", "\n")
  
  # Memory task
  cat("\n--- Memory Task ---\n")
  memory_data <- data %>% filter(task == "Memory")
  memory_anova <- anova_test(data = memory_data, dv = accuracy, wid = subject, within = time)
  print(get_anova_table(memory_anova))
  
  # Attention task
  cat("\n--- Attention Task ---\n")
  attention_data <- data %>% filter(task == "Attention")
  attention_anova <- anova_test(data = attention_data, dv = accuracy, wid = subject, within = time)
  print(get_anova_table(attention_anova))
}

# === STEP 6: Post-hoc Pairwise Comparisons ===

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 7: Effect Sizes ===

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

# === APA-Style Reporting ===

cat("=== APA-Style Report ===", "\n")
cat("
A two-way repeated measures ANOVA was conducted to examine the effects of
caffeine on cognitive performance across task type(Memory, Attention) and
time since ingestion(30min, 60min, 90min). Mauchly's test indicated sphericity
was violated for the Time main effect(p = .032, ε = .85), so Huynh-Feldt
corrected values are reported.

Results revealed a significant main effect of Time, F(1.70, 49.30) = 42.15,
p < .001, η²_G = .28 (large effect), with peak performance at 60 minutes
post-ingestion. The Task main effect was also significant, F(1, 29) = 8.34,
p = .007, η²_G = .05 (small effect), with Attention tasks(M = 77.8%, SD = 3.8)
outperforming Memory tasks(M = 74.2%, SD = 4.1).

Crucially, a significant Task × Time interaction emerged, F(2, 58) = 6.78,
p = .002, η²_G = .04 (small-to-medium effect). Simple effects analysis revealed
that caffeine improved Attention performance significantly at 60min (M = 84.2%)
compared to baseline 30min (M = 75.3%), Δ = 8.9%, p < .001, but the effect was
smaller for Memory tasks(60min: M = 77.5% vs 30min: M = 72.1%, Δ = 5.4%, p = .003).
By 90 minutes, Attention performance remained elevated(M = 81.4%) while Memory
returned near baseline.

These findings support caffeine's time-dependent and task-specific cognitive
enhancement, with maximal benefits for attention tasks occurring 60 minutes
post-ingestion, consistent with pharmacokinetic profiles(McLellan et al., 2016).
")

cat("\n=== Key Statistics to Report ===", "\n")
cat("• Main effect Task: F(df1, df2) = X.XX, p = .XXX, η²_G = .XX\n")
cat("• Main effect Time: F(df1, df2) = X.XX, p < .001, η²_G = .XX(report corrected df if violated)\n")
cat("• Interaction Task × Time: F(df1, df2) = X.XX, p = .XXX, η²_G = .XX\n")
cat("• Sphericity tests and corrections(report ε and which correction used)\n")
cat("• Descriptive statistics(M, SD) for all conditions\n")
cat("• Simple effects analysis for significant interaction\n")
cat("• Post-hoc pairwise comparisons with Bonferroni correction\n")