# Mixed models for the SuperPlots design (both conditions in every experiment), no true average effect.
# Scenario A of simulation.py: 4 experiments, 50 cells per condition per experiment, cell SD 2,
# day-to-day baseline SD 1, drug effect varying between experiments with SD 0.5 around 0.
# Compares false-positive rates at alpha = 0.05 for:
#   m1: speed ~ condition + (1 | experiment)            (random intercept only)
#   m2: speed ~ condition + (1 | experiment/condition)  (adds a condition-by-experiment term)
#   paired t-test on the experiment means
# R 4.3.3, lme4, lmerTest (Satterthwaite df). Run: Rscript mixed_model_simulation.R > mixed_model_simulation_output.txt
suppressPackageStartupMessages({ library(lme4); library(lmerTest) })
cat("R", as.character(getRversion()), "| lme4", as.character(packageVersion("lme4")),
    "| lmerTest", as.character(packageVersion("lmerTest")), "\n")
set.seed(20261110)
n_sim <- 1000; n_exp <- 4; n_cells <- 50
hits <- c(m1 = 0, m2 = 0, paired = 0); singular <- c(m1 = 0, m2 = 0)
for (s in seq_len(n_sim)) {
  day <- rnorm(n_exp, 0, 1); eff <- rnorm(n_exp, 0, 0.5)
  d <- expand.grid(cell = seq_len(n_cells), condition = c("control", "drug"), experiment = paste0("E", 1:n_exp))
  e <- as.integer(d$experiment)
  d$speed <- 10 + day[e] + ifelse(d$condition == "drug", eff[e], 0) + rnorm(nrow(d), 0, 2)
  m1 <- suppressMessages(lmer(speed ~ condition + (1 | experiment), data = d))
  m2 <- suppressMessages(lmer(speed ~ condition + (1 | experiment/condition), data = d))
  singular["m1"] <- singular["m1"] + isSingular(m1)
  singular["m2"] <- singular["m2"] + isSingular(m2)
  hits["m1"] <- hits["m1"] + (coef(summary(m1))["conditiondrug", "Pr(>|t|)"] < 0.05)
  hits["m2"] <- hits["m2"] + (coef(summary(m2))["conditiondrug", "Pr(>|t|)"] < 0.05)
  mm <- tapply(d$speed, list(d$experiment, d$condition), mean)
  hits["paired"] <- hits["paired"] + (t.test(mm[, "drug"], mm[, "control"], paired = TRUE)$p.value < 0.05)
}
cat("simulations:", n_sim, "\n")
cat("false-positive rate:\n"); print(round(hits / n_sim, 3))
cat("singular fits:\n"); print(round(singular / n_sim, 3))