# Simulation for /guides/repeated-measures-anova-vs-mixed-model/ (R 4.3.3, lme4/lmerTest/afex).
# Design as in the worked example: 2 groups x 6 mice x 4 weeks, mouse SD 1.5 g, residual SD 0.6 g.
# Seed 20261110 fixed in advance. 2,000 datasets per scenario.
# Run from this folder: Rscript simulation.R > simulation_output.txt
suppressPackageStartupMessages({ library(lmerTest); library(afex) })
set.seed(20261110)
nsim <- 2000

sim_data <- function(slope_vehicle, slope_drug, n = 6) {
  rows <- list()
  for (g in c("vehicle", "drug")) for (i in seq_len(n)) {
    base <- 25 + rnorm(1, 0, 1.5); s <- if (g == "vehicle") slope_vehicle else slope_drug
    rows[[length(rows) + 1]] <- data.frame(mouse = paste0(substr(g, 1, 1), i), group = g, week = 0:3,
                                           weight = base + s * (0:3) + rnorm(4, 0, 0.6))
  }
  d <- do.call(rbind, rows)
  d$group <- factor(d$group, levels = c("vehicle", "drug")); d$week_f <- factor(d$week); d$mouse <- factor(d$mouse)
  d
}

run <- function(slope_vehicle, slope_drug, drop_one) {
  res <- matrix(NA, nsim, 5, dimnames = list(NULL, c("naive_group", "naive_int", "lmm_group", "lmm_int", "rm_int")))
  for (k in seq_len(nsim)) {
    d <- sim_data(slope_vehicle, slope_drug)
    if (drop_one) d <- d[-sample(which(d$week > 0), 1), ]   # one value missing completely at random
    a <- anova(lm(weight ~ group * week_f, data = d))
    res[k, "naive_group"] <- a["group", "Pr(>F)"]; res[k, "naive_int"] <- a["group:week_f", "Pr(>F)"]
    m <- suppressMessages(lmer(weight ~ group * week_f + (1 | mouse), data = d))
    b <- anova(m); res[k, "lmm_group"] <- b["group", "Pr(>F)"]; res[k, "lmm_int"] <- b["group:week_f", "Pr(>F)"]
    r <- suppressWarnings(suppressMessages(aov_ez("mouse", "weight", d, between = "group", within = "week_f")))
    res[k, "rm_int"] <- summary(r)$pval.adjustments["group:week_f", "Pr(>F[GG])"]
  }
  round(100 * colMeans(res < 0.05), 1)
}

cat("Share of datasets with p < 0.05 (%), 2,000 per row\n")
cat("\nNo true effect (both groups gain 1.0 g/week), complete data:\n");  print(run(1.0, 1.0, FALSE))
cat("\nNo true effect, one value missing at random:\n");                print(run(1.0, 1.0, TRUE))
cat("\nTrue interaction (vehicle 1.0, drug 0.4 g/week), complete data:\n"); print(run(1.0, 0.4, FALSE))
cat("\nTrue interaction, one value missing at random:\n");              print(run(1.0, 0.4, TRUE))
cat("\nColumns: naive = two-way ANOVA ignoring mice; lmm = lmer + Satterthwaite; rm_int = RM-ANOVA, GG-corrected, complete cases only.\n")