Three-Way Mixed ANOVA

# Three-way Mixed ANOVA: Treatment (between) × Time (within) × Gender (between)
# Clinical depression study: CBT effectiveness over time by gender

library(ez)          # ezANOVA for mixed designs
library(afex)        # Alternative mixed ANOVA
library(emmeans)     # Marginal means and contrasts
library(ggplot2)
library(dplyr)
library(tidyr)
library(car)         # Levene's test

# Simulate realistic clinical trial data
set.seed(2025)
n_per_cell <- 30  # 30 per Treatment×Gender combination
time_points <- c("Baseline", "Week8", "Week16")

# Create between-subjects design (Treatment × Gender)
data_between <- expand.grid(
  treatment = c("CBT", "Control"),
  gender = c("Male", "Female")
) %>%
  slice(rep(1:n(), each=n_per_cell)) %>%
  mutate(
    subject_id = 1:n(),
    # Baseline random intercepts (individual differences)
    baseline_bdi = rnorm(n(), mean=28, sd=8),
    baseline_bdi = pmax(10, pmin(50, baseline_bdi))  # Constrain
  )

# Generate longitudinal data (repeated measures on Time)
data_long <- data_between %>%
  crossing(time = factor(time_points, levels=time_points)) %>%
  mutate(
    time_numeric = as.numeric(time) - 1,  # 0, 1, 2
    
    # True model effects:
    # Treatment effect: CBT reduces by -6 pts per timepoint vs. -2 for Control
    # Gender effect: Females start 3 pts higher but no interaction
    # Individual variation: ICC ~0.65
    
    treatment_effect = ifelse(treatment == "CBT", -6, -2),
    gender_effect = ifelse(gender == "Female", 3, 0),
    
    # BDI = Baseline + Gender + Treatment×Time + random_error
    bdi_score = baseline_bdi + gender_effect + 
                treatment_effect * time_numeric +
                rnorm(n(), mean=0, sd=5),  # Within-subject error
    bdi_score = pmax(0, pmin(63, bdi_score))
  ) %>%
  select(subject_id, treatment, gender, time, bdi_score)

# Convert to factors
data_long$subject_id <- factor(data_long$subject_id)
data_long$treatment <- factor(data_long$treatment)
data_long$gender <- factor(data_long$gender)
data_long$time <- factor(data_long$time, levels=time_points)

cat("=== Data Structure ===")
cat("\nTotal participants:", length(unique(data_long$subject_id)))
cat("\nTotal observations:", nrow(data_long))
cat("\nDesign: 2 (Treatment) × 3 (Time) × 2 (Gender)\n")

# Descriptive statistics
cat("\n=== Descriptive Statistics ===")
desc_stats <- data_long %>%
  group_by(treatment, gender, time) %>%
  summarise(
    n = n(),
    M = mean(bdi_score),
    SD = sd(bdi_score),
    .groups = 'drop'
  )
print(desc_stats, n=12)

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

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

# 1. Normality of residuals
# Fit preliminary model to get residuals
model_prelim <- aov(bdi_score ~ treatment * gender * time + 
                    Error(subject_id/time), data=data_long)
residuals_long <- residuals(lm(bdi_score ~ treatment * gender * time, data=data_long))

cat("\n1. Normality of Residuals")
shapiro_test <- shapiro.test(sample(residuals_long, min(5000, length(residuals_long))))
cat("\n   Shapiro-Wilk: W=", round(shapiro_test$statistic, 4), 
    ", p=", round(shapiro_test$p.value, 3))
cat(ifelse(shapiro_test$p.value > 0.05, " ✓", " ⚠"), "\n")

par(mfrow=c(1,2))
qqnorm(residuals_long, main="Q-Q Plot: Residuals")
qqline(residuals_long, col="red", lwd=2)
hist(residuals_long, breaks=30, main="Histogram: Residuals", 
     xlab="Residuals", col="lightblue")

# 2. Homogeneity of variance (between-subjects factors)
cat("\n2. Homogeneity of Variance(Between-Subjects)")
levene_test <- leveneTest(bdi_score ~ treatment * gender, 
                          data=filter(data_long, time=="Baseline"))
cat("\n   Levene's test: F=", round(levene_test$`F value`[1], 3),
    ", p=", round(levene_test$`Pr(>F)`[1], 3))
cat(ifelse(levene_test$`Pr(>F)`[1] > 0.05, " ✓", " ⚠"), "\n")

par(mfrow=c(1,1))
boxplot(bdi_score ~ treatment * gender, 
        data=filter(data_long, time=="Baseline"),
        main="BDI by Treatment × Gender(Baseline)",
        xlab="Group", ylab="BDI-II Score", col=c("lightblue", "lightpink"))

# 3. Sphericity (within-subjects factor: Time)
cat("\n3. Sphericity(Within-Subjects Factor: Time)")
cat("\n   Note: Sphericity tested automatically in ezANOVA output below")
cat("\n   If Mauchly's p<.05 → use Greenhouse-Geisser or Huynh-Feldt correction\n")

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

cat("\n\n=== THREE-WAY MIXED ANOVA ===")

# Using ezANOVA (automatically tests sphericity and provides corrections)
anova_result <- ezANOVA(
  data = data_long,
  dv = bdi_score,
  wid = subject_id,
  within = time,
  between = .(treatment, gender),
  type = 3,
  detailed = TRUE,
  return_aov = TRUE
)

print(anova_result$ANOVA)

cat("\n=== Sphericity Tests(Mauchly's) ===")
if(!is.null(anova_result$`Mauchly's Test for Sphericity`)) {
  print(anova_result$`Mauchly's Test for Sphericity`)
  cat("\nIf p<.05, sphericity violated → use corrections below\n")
}

cat("\n=== Sphericity Corrections ===")
if(!is.null(anova_result$`Sphericity Corrections`)) {
  print(anova_result$`Sphericity Corrections`)
  cat("\nGreenhouse-Geisser(GGe): Conservative correction(use if ε<.75)")
  cat("\nHuynh-Feldt(HFe): Less conservative(use if ε>.75)\n")
}

# === STEP 3: Interpret Main Effects and Interactions ===

cat("\n\n=== INTERPRETATION OF RESULTS ===")

cat("\n1. MAIN EFFECT OF TREATMENT(between-subjects):")
cat("\n   Tests if CBT differs from Control averaged across Time and Gender")
cat("\n   If p<.05 → CBT and Control differ overall\n")

cat("\n2. MAIN EFFECT OF TIME(within-subjects):")
cat("\n   Tests if depression changes over time(averaged across Treatment and Gender)")
cat("\n   If p<.05 → Significant change from Baseline to Week 8 to Week 16")
cat("\n   Check sphericity correction! Use GG or HF p-value if Mauchly's p<.05\n")

cat("\n3. MAIN EFFECT OF GENDER(between-subjects):")
cat("\n   Tests if Males differ from Females averaged across Time and Treatment")
cat("\n   If p<.05 → Gender difference in overall depression levels\n")

cat("\n4. TREATMENT × TIME INTERACTION(critical for efficacy):")
cat("\n   Tests if rate of change over time differs between CBT and Control")
cat("\n   If p<.05 → CBT shows different trajectory than Control(expected!)")
cat("\n   This is the PRIMARY hypothesis in treatment trials\n")

cat("\n5. TREATMENT × GENDER INTERACTION:")
cat("\n   Tests if treatment effectiveness differs by gender(averaged over time)")
cat("\n   If p<.05 → CBT works better for one gender than the other\n")

cat("\n6. TIME × GENDER INTERACTION:")
cat("\n   Tests if rate of change over time differs by gender")
cat("\n   If p<.05 → Males and Females show different trajectories\n")

cat("\n7. TREATMENT × TIME × GENDER(three-way interaction):")
cat("\n   Tests if Treatment×Time interaction is moderated by Gender")
cat("\n   If p<.05 → CBT's time trajectory differs by gender(e.g., CBT works faster for females)")
cat("\n   Requires simple effects analysis to interpret\n")

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

cat("\n\n=== VISUALIZATION ===")

# Interaction plot: Treatment × Time by Gender
interaction_means <- data_long %>%
  group_by(treatment, gender, time) %>%
  summarise(M = mean(bdi_score), 
            SE = sd(bdi_score)/sqrt(n()),
            .groups='drop')

ggplot(interaction_means, aes(x=time, y=M, color=treatment, group=treatment)) +
  geom_line(size=1.2) +
  geom_point(size=3) +
  geom_errorbar(aes(ymin=M-1.96*SE, ymax=M+1.96*SE), width=0.1) +
  facet_wrap(~gender) +
  labs(title="Three-way Interaction: Treatment × Time × Gender",
       subtitle="Mean BDI-II Score ± 95% CI",
       x="Time Point", y="Depression(BDI-II)",
       color="Treatment") +
  theme_classic() +
  theme(legend.position="bottom") +
  scale_color_manual(values=c("CBT"="#2E86AB", "Control"="#A23B72"))

# Profile plot: Emphasize Treatment × Time interaction
ggplot(interaction_means, aes(x=time, y=M, color=treatment, linetype=gender, group=interaction(treatment, gender))) +
  geom_line(size=1.2) +
  geom_point(size=3) +
  labs(title="Depression Trajectory by Treatment, Time, and Gender",
       x="Time Point", y="Mean BDI-II Score",
       color="Treatment", linetype="Gender") +
  theme_classic()

# === STEP 5: Post-hoc Tests (if interactions significant) ===

cat("\n\n=== POST-HOC TESTS: Simple Effects Analysis ===")

# If Treatment×Time interaction is significant, test simple effects:
# Effect of Treatment at each Time point

cat("\nSimple effects: Effect of Treatment at each Time point\n")

emm <- emmeans(anova_result$aov, ~ treatment | time)
pairs_time <- pairs(emm, adjust="bonferroni")
print(pairs_time)

cat("\nInterpretation: Compare CBT vs. Control at Baseline, Week 8, Week 16")
cat("\nExpected: No difference at Baseline(randomization); CBT<Control at Weeks 8 & 16\n")

# Effect of Time within each Treatment
cat("\n\nSimple effects: Effect of Time within each Treatment\n")

emm_time <- emmeans(anova_result$aov, ~ time | treatment)
pairs_treatment <- pairs(emm_time, adjust="bonferroni")
print(pairs_treatment)

cat("\nInterpretation: Test if depression changes over time within CBT and Control separately")
cat("\nExpected: Significant decline in CBT; smaller/no decline in Control\n")

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

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

# Partial eta-squared from ezANOVA output
cat("\nPartial η² for each effect(from ANOVA table above)")
cat("\nInterpretation: .01=small, .06=medium, .14=large(Cohen, 1988)\n")

# Calculate effect size for key comparison: CBT vs. Control at Week 16
week16_data <- filter(data_long, time=="Week16")
cbt_16 <- filter(week16_data, treatment=="CBT")$bdi_score
control_16 <- filter(week16_data, treatment=="Control")$bdi_score

cohens_d <- (mean(cbt_16) - mean(control_16)) / 
            sqrt((var(cbt_16) + var(control_16)) / 2)

cat("\nCohen's d(CBT vs. Control at Week 16):", round(cohens_d, 2))
cat("\nInterpretation: .2=small, .5=medium, .8=large")
cat("\nExpected d ~ -0.6 to -0.8 (negative = CBT lower depression)\n")

# === APA-Style Reporting ===

cat("\n\n=== APA-STYLE REPORT ===")
cat("
A three-way mixed ANOVA was conducted to examine the effects of treatment 
(CBT vs. Control, between-subjects), time(Baseline, 8-week, 16-week, within-subjects), 
and gender(Male vs. Female, between-subjects) on depression severity(BDI-II scores). 
The sample included 120 participants(30 per Treatment×Gender combination) measured 
at three time points(N=360 observations).

Assumptions were evaluated: Normality of residuals was satisfactory(Shapiro-Wilk p>.05). 
Levene's test indicated homogeneity of variance for between-subjects factors(p=.XX). 
Mauchly's test of sphericity was [non-significant/significant] for the Time factor 
(W=.XX, p=.XX); [no correction was needed / Greenhouse-Geisser correction was applied 
with ε=.XX].

Results revealed a significant three-way interaction [F(2, 232)=X.XX, p=.XXX, 
partial η²=.XX] OR [The three-way interaction was not significant, F(2, 232)=X.XX, p=.XX]. 
The critical Treatment×Time interaction was significant [F(2, 232)=XX.XX, p<.001, 
partial η²=.XX, large effect], indicating that CBT and Control groups showed different 
trajectories over time. 

Simple effects analysis revealed that at baseline, CBT and Control did not differ 
(p=.XX), confirming successful randomization. At Week 8, CBT showed significantly 
lower depression than Control(M_CBT=XX.X, M_Control=XX.X, p<.001, d=X.XX). This 
difference increased by Week 16 (M_CBT=XX.X, M_Control=XX.X, p<.001, d=X.XX), 
demonstrating sustained and increasing treatment efficacy.

[If gender moderation:] The Treatment×Time×Gender interaction was significant(p=.XX), 
indicating that treatment trajectories differed by gender. Post-hoc analysis showed that 
[describe pattern].

Main effects: Time showed a significant main effect(F(2, 232)=XX.XX, p<.001, 
partial η²=.XX), indicating overall decrease in depression across time points. 
Treatment showed a significant main effect(F(1, 116)=XX.XX, p<.001, partial η²=.XX), 
with CBT participants showing lower overall depression than Controls. [Gender main 
effect results].

These findings support CBT as an effective treatment for depression, with therapeutic 
benefits emerging by 8 weeks and strengthening through 16 weeks(Cohen's d=-0.XX, 
large effect). Results are consistent with meta-analytic estimates of CBT efficacy 
(Cuijpers et al., 2016).
")

cat("\n\n=== KEY REPORTING ELEMENTS ===")
cat("
✓ Design clearly specified(2×3×2 mixed)
✓ Sample size and cell sizes reported
✓ Assumption checks reported(normality, homogeneity, sphericity)
✓ Sphericity corrections applied if needed(GG or HF)
✓ All main effects and interactions reported with F, df, p, partial η²
✓ Simple effects analysis for significant interactions
✓ Effect sizes for key comparisons(Cohen's d)
✓ Means, SDs, and confidence intervals for key cells
✓ Interpretation linked to research question and literature
")