One-Way Repeated Measures ANOVA

# One-way Repeated Measures ANOVA: Cortisol across 4 daily timepoints
# Based on diurnal cortisol patterns in mindfulness research

library(tidyverse)   # Data manipulation
library(ez)          # For ezANOVA (comprehensive RM-ANOVA)
library(rstatix)     # For get_summary_stats, anova_test
library(effectsize)  # For effect sizes
library(ggpubr)      # For ggqqplot

# Simulate realistic cortisol data (or load: data <- read.csv("cortisol.csv"))
set.seed(2025)
subject_ids <- rep(1:40, each=4)
timepoint <- rep(c("Awakening", "30min", "Noon", "Bedtime"), 40)

# Realistic cortisol pattern: high at awakening, peak at +30min, decline by bedtime
data_wide <- data.frame(
  subject = 1:40,
  Awakening = rnorm(40, mean=0.52, sd=0.12),
  ThirtyMin = rnorm(40, mean=0.68, sd=0.15),  # +30% awakening response
  Noon = rnorm(40, mean=0.24, sd=0.08),
  Bedtime = rnorm(40, mean=0.08, sd=0.04)
)

# Convert to long format (required for RM-ANOVA)
data_long <- data_wide %>%
  pivot_longer(cols = c(Awakening, ThirtyMin, Noon, Bedtime),
               names_to = "timepoint",
               values_to = "cortisol") %>%
  mutate(timepoint = factor(timepoint, 
                            levels=c("Awakening", "ThirtyMin", "Noon", "Bedtime")))

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

# 1. Descriptive statistics per timepoint
data_long %>%
  group_by(timepoint) %>%
  get_summary_stats(cortisol, type = "mean_sd")

# 2. Visualize: Profile plot (means across time)
data_long %>%
  group_by(timepoint) %>%
  summarise(M = mean(cortisol), SE = sd(cortisol)/sqrt(n())) %>%
  ggplot(aes(x=timepoint, y=M, group=1)) +
  geom_line(color="steelblue", size=1.2) +
  geom_point(size=3) +
  geom_errorbar(aes(ymin=M-SE, ymax=M+SE), width=0.2) +
  labs(title="Diurnal Cortisol Pattern(Post-Intervention)",
       x="Time of Day", y="Mean Cortisol(μg/dL) ± SE") +
  theme_classic()

# 3. Check outliers (boxplots per timepoint)
ggboxplot(data_long, x="timepoint", y="cortisol", add="jitter",
          title="Cortisol by Timepoint", xlab="Time", ylab="Cortisol(μg/dL)")

# 4. Check normality of DIFFERENCE SCORES
diff_data <- data_wide %>%
  mutate(
    diff_Awk_30min = ThirtyMin - Awakening,
    diff_Awk_Noon = Noon - Awakening,
    diff_Awk_Bed = Bedtime - Awakening,
    diff_30min_Noon = Noon - ThirtyMin,
    diff_30min_Bed = Bedtime - ThirtyMin,
    diff_Noon_Bed = Bedtime - Noon
  )

# Q-Q plots for key differences
par(mfrow=c(2,3))
for (diff_var in c("diff_Awk_30min", "diff_Awk_Noon", "diff_Awk_Bed",
                   "diff_30min_Noon", "diff_30min_Bed", "diff_Noon_Bed")) {
  qqnorm(diff_data[[diff_var]], main=diff_var)
  qqline(diff_data[[diff_var]])
}

# Shapiro-Wilk on difference scores
shapiro.test(diff_data$diff_Awk_30min)  # p > .05 OK
shapiro.test(diff_data$diff_Awk_Noon)   # p > .05 OK

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

# Method 1: Using ez::ezANOVA (comprehensive output)
rm_anova <- ezANOVA(
  data = data_long,
  dv = cortisol,
  wid = subject,
  within = timepoint,
  detailed = TRUE,
  type = 3
)

print(rm_anova)

# Output includes:
# - ANOVA table with F, p-value, generalized eta squared (ges)
# - Mauchly's test of sphericity
# - Greenhouse-Geisser and Huynh-Feldt corrections (if sphericity violated)

# Method 2: Using rstatix (tidy output)
rm_anova2 <- anova_test(
  data = data_long,
  dv = cortisol,
  wid = subject,
  within = timepoint
)

get_anova_table(rm_anova2)

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

# If p < .05 (sphericity violated), use corrected results:
cat("\n=== Sphericity Corrections ===")
print(rm_anova$`Sphericity Corrections`)

# Interpretation:
# - If Mauchly's p > .05: use uncorrected ANOVA results
# - 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

# === STEP 4: Effect Size ===
eta_squared(rm_anova2$ANOVA, partial=TRUE)
# partial η² for within-subjects factor

# === STEP 5: Post-hoc Pairwise Comparisons (if p < .05) ===
posthoc <- data_long %>%
  pairwise_t_test(
    cortisol ~ timepoint,
    paired = TRUE,
    p.adjust.method = "bonferroni"  # or "holm"
  )

print(posthoc)

# Compact display
posthoc %>%
  select(group1, group2, p.adj, p.adj.signif) %>%
  arrange(p.adj)

# === STEP 6: Visualize Results ===
# Bar plot with error bars
data_summary <- data_long %>%
  group_by(timepoint) %>%
  summarise(
    M = mean(cortisol),
    SE = sd(cortisol)/sqrt(n())
  )

ggplot(data_summary, aes(x=timepoint, y=M, fill=timepoint)) +
  geom_bar(stat="identity", width=0.6, alpha=0.8) +
  geom_errorbar(aes(ymin=M-1.96*SE, ymax=M+1.96*SE), width=0.2) +
  labs(title="Diurnal Cortisol Pattern After Mindfulness Training",
       x="Time of Day", y="Mean Cortisol(μg/dL) ± 95% CI") +
  scale_fill_brewer(palette="Blues") +
  theme_classic() +
  theme(legend.position="none")

# === APA-Style Reporting ===
cat("\n=== APA Results ===")
cat("
A one-way repeated measures ANOVA was conducted to examine diurnal cortisol
patterns following 8 weeks of mindfulness training(N=40). Mauchly's test 
indicated that the assumption of sphericity was met, χ²(5) = 8.32, p = .14.
There was a significant effect of time of day on cortisol levels, 
F(3, 117) = 287.45, p < .001, partial η² = .88 (large effect).

Post-hoc pairwise comparisons with Bonferroni correction revealed:
- Cortisol significantly increased from awakening(M=0.52, SD=0.12) to 
  +30 minutes(M=0.68, SD=0.15), p < .001 (cortisol awakening response).
- Cortisol then declined significantly from +30min to noon(M=0.24, SD=0.08), 
  p < .001.
- Bedtime cortisol(M=0.08, SD=0.04) was significantly lower than all other 
  timepoints, all p < .001.

These results demonstrate a normal diurnal cortisol rhythm with expected
awakening response and evening nadir, consistent with healthy HPA axis
function following mindfulness training.
")