# Two treatments x two sexes, 5 mice per cell. R 4.3.3, base R only.
d <- read.csv("grip.csv")
d$sex <- factor(d$sex, levels = c("female", "male")); d$treatment <- factor(d$treatment, levels = c("vehicle", "drug"))
g <- function(s, k) d$grip_g[d$sex == s & d$treatment == k]
cat("Cell means (g):\n")
for (s in c("female", "male")) for (k in c("vehicle", "drug"))
  cat(sprintf("  %-6s %-7s %.2f (n = %d)\n", s, k, mean(g(s, k)), length(g(s, k))))
cat("Separate t-tests in each sex (the reading error):\n")
for (s in c("male", "female")) {
  t <- t.test(g(s, "drug"), g(s, "vehicle"), var.equal = TRUE)
  cat(sprintf("  %s: difference %.2f g, t = %.2f, p = %.4f\n", s, mean(g(s, "drug")) - mean(g(s, "vehicle")), t$statistic, t$p.value))
}
t <- t.test(d$grip_g[d$treatment == "drug"], d$grip_g[d$treatment == "vehicle"], var.equal = TRUE)
cat(sprintf("Pooled t-test ignoring sex: t = %.2f, p = %.4f\n", t$statistic, t$p.value))
a <- car::Anova(lm(grip_g ~ treatment * sex, data = d), type = 2)   # Type II (all types agree here)
cat("Two-way ANOVA (treatment x sex):\n")
for (i in 1:3) cat(sprintf("  %-15s F(1, %d) = %.2f, p = %.4f\n", c("treatment", "sex", "treatment x sex")[i], a$Df[4], a$`F value`[i], a$`Pr(>F)`[i]))
n <- 5; dfe <- nrow(d) - 4; mse <- a$`Sum Sq`[4] / a$Df[4]; tq <- qt(0.975, dfe)
eff <- c(female = mean(g("female", "drug")) - mean(g("female", "vehicle")), male = mean(g("male", "drug")) - mean(g("male", "vehicle")))
rows <- list(list("treatment effect, averaged over sexes", mean(eff), sqrt(mse / n)),
             list("treatment effect in females", eff[["female"]], sqrt(2 * mse / n)),
             list("treatment effect in males", eff[["male"]], sqrt(2 * mse / n)),
             list("difference between sexes in the effect (male - female)", eff[["male"]] - eff[["female"]], sqrt(4 * mse / n)))
cat("Effects with 95% CI (g):\n")
for (r in rows) cat(sprintf("  %s: %.2f (%.2f to %.2f)\n", r[[1]], r[[2]], r[[2]] - tq * r[[3]], r[[2]] + tq * r[[3]]))