# AR(1)-correlated data for /guides/repeated-measures-anova-vs-mixed-model/ (R 4.3.3).
# Same design as simulation.R (2 groups x 6 mice x 4 weeks, complete data, no true group x week effect),
# but the within-mouse correlation falls with distance: corr = 0.9^|lag|, total SD 1.6 g, so the
# random-intercept (compound symmetry) assumption of the mixed model in the worked example is wrong.
# Seed 20261112 fixed in advance. 2,000 datasets.
# Run from this folder: Rscript simulation_ar1.R > simulation_ar1_output.txt
suppressPackageStartupMessages({ library(lmerTest); library(afex); library(nlme); library(MASS) })
set.seed(20261112)
nsim <- 2000
Sig <- outer(0:3, 0:3, function(i, j) 0.9^abs(i - j)) * 1.6^2
p <- matrix(NA, nsim, 4, dimnames = list(NULL, c("lmm_cs_satterthwaite", "rm_anova_gg", "gls_ar1", "slope_welch")))
for (k in seq_len(nsim)) {
  y <- mvrnorm(12, rep(0, 4), Sig) + matrix(25 + (0:3), 12, 4, byrow = TRUE)
  d <- data.frame(mouse = factor(rep(1:12, each = 4)), group = factor(rep(c("vehicle", "drug"), each = 24)),
                  week = rep(0:3, 12), weight = as.vector(t(y)))
  d$week_f <- factor(d$week)
  m <- suppressMessages(lmer(weight ~ group * week_f + (1 | mouse), data = d))
  p[k, 1] <- anova(m)["group:week_f", "Pr(>F)"]
  r <- suppressWarnings(suppressMessages(aov_ez("mouse", "weight", d, between = "group", within = "week_f")))
  p[k, 2] <- summary(r)$pval.adjustments["group:week_f", "Pr(>F[GG])"]
  g <- gls(weight ~ group * week_f, data = d, correlation = corAR1(form = ~ week | mouse), method = "REML")
  p[k, 3] <- anova(g, type = "marginal")["group:week_f", "p-value"]
  s <- sapply(split(d, d$mouse), function(x) coef(lm(weight ~ week, data = x))[["week"]])
  grp <- sapply(split(d, d$mouse), function(x) as.character(x$group[1]))
  p[k, 4] <- t.test(s[grp == "drug"], s[grp == "vehicle"])$p.value
}
cat("No true effect, AR(1) correlation 0.9^lag, complete data: share with p < 0.05 (%), 2,000 datasets\n")
print(round(100 * colMeans(p < 0.05), 1))
cat("Monte Carlo SE at 5%: about 0.5 percentage points\n")