# Planned n = 8 per group; 4 more per group added after seeing p = 0.08. R 4.3.3, base R.
d <- read.csv("rotarod.csv")
d$group <- factor(d$group, levels = c("control", "drug")); d$cohort <- factor(d$cohort)
for (lab in c("cohort 1", "both cohorts")) {
  x <- if (lab == "cohort 1") d[d$cohort == 1, ] else d
  m <- tapply(x$rotarod_s, x$group, mean); s <- tapply(x$rotarod_s, x$group, sd)
  cat(sprintf("Mean (SD), %s: control %.2f (%.2f), drug %.2f (%.2f)\n", lab, m[["control"]], s[["control"]], m[["drug"]], s[["drug"]]))
}
test <- function(x, welch = FALSE) {                 # Student's t-test (equal variances) unless welch = TRUE
  a <- x$rotarod_s[x$group == "control"]; b <- x$rotarod_s[x$group == "drug"]
  t <- t.test(b, a, var.equal = !welch)
  list(n = length(a), diff = mean(b) - mean(a), t = unname(t$statistic), p = t$p.value)
}
r <- test(d[d$cohort == 1, ])
cat(sprintf("Planned analysis, cohort 1 only: n = %d per group, difference %.2f s, t = %.2f, p = %.4f\n", r$n, r$diff, r$t, r$p))
r <- test(d)
cat(sprintf("After adding cohort 2 (unplanned): n = %d per group, difference %.2f s, t = %.2f, p = %.4f\n", r$n, r$diff, r$t, r$p))
cat(sprintf("Welch's t-test instead: cohort 1 p = %.4f; both cohorts p = %.4f\n", test(d[d$cohort == 1, ], TRUE)$p, test(d, TRUE)$p))
fit <- lm(rotarod_s ~ group + cohort, data = d)
ci <- confint(fit)["groupdrug", ]
cat(sprintf("Both cohorts, cohort as a block: difference %.2f s (95%% CI %.2f to %.2f), p = %.4f\n",
            coef(fit)[["groupdrug"]], ci[1], ci[2], summary(fit)$coefficients["groupdrug", 4]))
r <- test(d[d$cohort == 2, ])
cat(sprintf("Cohort 2 alone: n = %d per group, difference %.2f s, p = %.4f\n", r$n, r$diff, r$p))