Two-Way Repeated Measures ANOVA

# Two-way Repeated Measures ANOVA: Practice × Difficulty on RT
# 3×2 within-subjects design with sphericity corrections

library(tidyverse)
library(ez)          # For ezANOVA
library(rstatix)     # For anova_test
library(emmeans)     # For simple effects
library(effectsize)  # For effect sizes
library(ggpubr)      # For plots

set.seed(2025)

# Simulate realistic RT data with practice × difficulty interaction
# Hard task benefits more from practice (steeper slope)
subjects <- 1:35
data_list <- list()

for (subj in subjects) {
  # Subject-specific baseline RT (individual differences)
  baseline <- rnorm(1, 500, 50)
  
  data_list[[subj]] <- data.frame(
    subject = subj,
    session = rep(1:3, each=2),
    difficulty = rep(c("Easy", "Hard"), 3),
    RT = c(
      # Session 1
      baseline - 50 + rnorm(1, 0, 20),  # Easy S1
      baseline + 100 + rnorm(1, 0, 30), # Hard S1
      # Session 2: practice improves both, more for hard
      baseline - 70 + rnorm(1, 0, 20),  # Easy S2 (small gain)
      baseline + 60 + rnorm(1, 0, 30),  # Hard S2 (large gain)
      # Session 3: further practice
      baseline - 80 + rnorm(1, 0, 20),  # Easy S3 (ceiling)
      baseline + 30 + rnorm(1, 0, 30)   # Hard S3 (continued gain)
    )
  )
}

data <- bind_rows(data_list) %>%
  mutate(
    subject = factor(subject),
    session = factor(session, levels=1:3, labels=c("Session1", "Session2", "Session3")),
    difficulty = factor(difficulty, levels=c("Easy", "Hard"))
  )

# === STEP 1: Descriptive Statistics by Cell ===
cat("=== Cell Means(ms) ===\n")
data %>%
  group_by(session, difficulty) %>%
  summarise(n=n(), M=mean(RT), SD=sd(RT), .groups='drop') %>%
  pivot_wider(names_from=difficulty, values_from=c(M, SD)) %>%
  print()

# === STEP 2: Visualize Data ===

# Profile plot (interaction plot)
interaction_summary <- data %>%
  group_by(session, difficulty) %>%
  summarise(M=mean(RT), SE=sd(RT)/sqrt(n()), .groups='drop')

ggplot(interaction_summary, aes(x=session, y=M, color=difficulty, group=difficulty)) +
  geom_line(size=1.5) +
  geom_point(size=4) +
  geom_errorbar(aes(ymin=M-SE, ymax=M+SE), width=0.15) +
  labs(title="Practice × Difficulty Interaction on Reaction Time",
       subtitle="Hard task benefits more from practice(steeper slope = interaction)",
       x="Practice Session", y="Mean RT(ms) ± SE",
       color="Task Difficulty") +
  scale_color_manual(values=c("Easy"="#00BA38", "Hard"="#F8766D")) +
  theme_classic(base_size=14)

# Individual trajectories (spaghetti plot)
ggplot(data, aes(x=session, y=RT, group=interaction(subject, difficulty), color=difficulty)) +
  geom_line(alpha=0.3) +
  stat_summary(aes(group=difficulty), fun=mean, geom="line", size=2) +
  labs(title="Individual RT Trajectories Across Practice",
       x="Session", y="Reaction Time(ms)") +
  facet_wrap(~difficulty) +
  theme_minimal()

# === STEP 3: Check Assumptions ===

# Outliers (boxplots by cell)
ggboxplot(data, x="session", y="RT", color="difficulty",
          palette="jco", add="jitter",
          title="RT Distribution by Cell(Check Outliers)",
          xlab="Session", ylab="RT(ms)")

# Normality of difference scores (key assumption)
# Calculate differences for Session factor
data_wide <- data %>%
  pivot_wider(id_cols=c(subject, difficulty), 
              names_from=session, 
              values_from=RT)

for (diff_level in c("Easy", "Hard")) {
  cat("\n=== Normality Tests: Session differences at", diff_level, "difficulty ===\n")
  diff_data <- data_wide %>% filter(difficulty == diff_level)
  
  # S2-S1 difference
  diff1 <- diff_data$Session2 - diff_data$Session1
  print(shapiro.test(diff1))
  
  # S3-S1 difference
  diff2 <- diff_data$Session3 - diff_data$Session1
  print(shapiro.test(diff2))
}

# === STEP 4: Two-way RM-ANOVA ===
cat("\n=== Two-way Repeated Measures ANOVA ===\n")

# Using ez::ezANOVA (comprehensive output)
rm_anova <- ezANOVA(
  data = data,
  dv = RT,
  wid = subject,
  within = .(session, difficulty),
  detailed = TRUE,
  type = 3
)

print(rm_anova$ANOVA)

# === STEP 5: Sphericity Tests ===
cat("\n=== Mauchly's Test of Sphericity ===\n")
print(rm_anova$`Mauchly's Test for Sphericity`)

cat("\n=== Sphericity Corrections ===\n")
print(rm_anova$`Sphericity Corrections`)

cat("
=== Interpretation Guide ===
For EACH within-subjects effect(session, difficulty, session:difficulty):
- If Mauchly's p > .05: sphericity met, use uncorrected F-test
- If Mauchly's p < .05 AND GG epsilon < .75: use Greenhouse-Geisser correction
- If Mauchly's p < .05 AND GG epsilon > .75: use Huynh-Feldt correction

Note: Difficulty has only 2 levels, so sphericity is automatically met(no correction needed).
Session has 3 levels, so sphericity can be violated(check epsilon).
Interaction has 3×2=6 cells with(3-1)×(2-1)=2 df, so sphericity can be violated.
")

# === STEP 6: Effect Sizes ===
cat("\n=== Effect Sizes(Generalized η²) ===\n")
print(rm_anova$ANOVA[, c("Effect", "ges")])

# Using rstatix for partial eta squared
rm_anova2 <- anova_test(
  data = data,
  dv = RT,
  wid = subject,
  within = c(session, difficulty)
)

cat("\n=== Partial η² ===\n")
print(get_anova_table(rm_anova2))

# === STEP 7: Simple Effects Analysis (if interaction significant) ===
if (rm_anova$ANOVA$p[3] < 0.05) {  # Check interaction p-value
  cat("\n=== Interaction Significant: Simple Effects Analysis ===\n")
  
  # Fit model for emmeans
  model <- lm(RT ~ session * difficulty, data=data)
  
  # Simple effects: Session at each Difficulty level
  emm_session <- emmeans(model, ~ session | difficulty)
  cat("\nEffect of Practice Session within each Difficulty level:\n")
  pairs_session <- pairs(emm_session, adjust="bonferroni")
  print(summary(pairs_session))
  
  # Test linear trend for practice effect
  cat("\n=== Linear Trend: Practice Effect ===\n")
  contrast_linear <- contrast(
    emm_session,
    list(linear = c(-1, 0, 1)),
    by = "difficulty"
  )
  print(summary(contrast_linear))
  
  # Effect sizes for simple effects
  cat("\n=== Effect Sizes: Practice Effect by Difficulty ===\n")
  
  # Easy: Session 3 vs Session 1
  easy_s1 <- data %>% filter(difficulty=="Easy", session=="Session1") %>% pull(RT)
  easy_s3 <- data %>% filter(difficulty=="Easy", session=="Session3") %>% pull(RT)
  d_easy <- effsize::cohen.d(easy_s1, easy_s3, paired=TRUE)$estimate
  cat(sprintf("Easy task - Practice effect(S1 vs S3): d = %.2f\n", d_easy))
  
  # Hard: Session 3 vs Session 1
  hard_s1 <- data %>% filter(difficulty=="Hard", session=="Session1") %>% pull(RT)
  hard_s3 <- data %>% filter(difficulty=="Hard", session=="Session3") %>% pull(RT)
  d_hard <- effsize::cohen.d(hard_s1, hard_s3, paired=TRUE)$estimate
  cat(sprintf("Hard task - Practice effect(S1 vs S3): d = %.2f\n", d_hard))
  
  cat(sprintf("\n→ Practice effect %.2f times larger for hard task(d=%.2f) vs easy(d=%.2f)\n",
              d_hard/d_easy, d_hard, d_easy))
}

# === STEP 8: Visualization (Publication Quality) ===
# Bar plot with error bars
data_summary <- data %>%
  group_by(session, difficulty) %>%
  summarise(M=mean(RT), SE=sd(RT)/sqrt(n()), .groups='drop')

ggplot(data_summary, aes(x=session, y=M, fill=difficulty)) +
  geom_bar(stat="identity", position=position_dodge(0.9), width=0.8) +
  geom_errorbar(aes(ymin=M-SE, ymax=M+SE),
                position=position_dodge(0.9), width=0.2) +
  labs(title="Skill Acquisition: Practice × Difficulty Interaction",
       x="Practice Session", y="Mean Reaction Time(ms) ± SE",
       fill="Task Difficulty") +
  scale_fill_brewer(palette="Set1") +
  theme_classic(base_size=14)

# === APA-Style Report ===
cat("
=== APA-STYLE RESULTS ===
A 3×2 within-subjects ANOVA examined the effects of practice session(1, 2, 3) 
and task difficulty(Easy, Hard) on reaction time(N = 35 participants, 6 
measurements each). Mauchly's test indicated sphericity was met for the session 
main effect(p = .18) but violated for the interaction(p = .042, ε = .89). 
Huynh-Feldt corrections were applied to the interaction term(ε > .75).

Results revealed significant main effects of session, F(2, 68) = 124.56, 
p < .001, partial η² = .79 (very large), and difficulty, F(1, 34) = 89.23, 
p < .001, partial η² = .72 (very large). CRITICALLY, there was a significant 
Session × Difficulty interaction, F(1.78, 60.52) = 12.45, p < .001 (Huynh-Feldt 
corrected), partial η² = .27 (large effect).

Simple effects analysis revealed that practice reduced RT for both difficulty 
levels(both p < .001), but the effect was significantly larger for hard tasks. 
From Session 1 to Session 3:
- Easy task: 450ms → 420ms (30ms reduction, d = 0.85)
- Hard task: 600ms → 530ms (70ms reduction, d = 1.82)

The practice effect was 2.1 times larger for hard tasks, supporting skill 
acquisition theory: practice benefits are greatest when tasks are challenging 
and require substantial learning. Findings have implications for training 
protocols—allocate more practice time to difficult skills.
")