# Litter effects: a treatment given to pregnant dams, pups weighed at weaning. Pups as n, litter means, and litter
# means adjusted for litter size. R 4.3.3, base R.
d <- read.csv("pups.csv")                            # one row per pup: treatment, litter, litter_size, pup, sex, weight_g
f2 <- function(p) formatC(signif(p, 2), format = "g")
m <- aggregate(cbind(weight_g, litter_size) ~ treatment + litter, data = d, FUN = mean)
for (g in c("control", "treated")) { x <- m[m$treatment == g, ]
  cat(sprintf("%s: %d litters, %d pups, mean litter size %.1f, mean of litter means %.2f g\n", g, nrow(x), sum(d$treatment == g), mean(x$litter_size), mean(x$weight_g))) }
pt <- d$weight_g[d$treatment == "treated"]; pc <- d$weight_g[d$treatment == "control"]; t <- t.test(pt, pc)
cat(sprintf("A. Pups as n (%d vs %d), Welch t-test: difference %.2f g, p = %s\n", length(pt), length(pc), mean(pt) - mean(pc), f2(t$p.value)))
tr <- m$weight_g[m$treatment == "treated"]; co <- m$weight_g[m$treatment == "control"]; t <- t.test(tr, co)
cat(sprintf("B. Litters as n (%d vs %d), Welch t-test on litter means: difference %.2f g (95%% CI %.2f to %.2f), p = %s\n", length(tr), length(co), mean(tr) - mean(co), t$conf.int[1], t$conf.int[2], f2(t$p.value)))
f <- lm(weight_g ~ treatment + litter_size, data = m); s <- summary(f)$coef["treatmenttreated", ]; ci <- confint(f)["treatmenttreated", ]
cat(sprintf("C. Litters as n, litter means adjusted for litter size (linear model): difference %.2f g (95%% CI %.2f to %.2f), p = %s\n", s[["Estimate"]], ci[1], ci[2], f2(s[["Pr(>|t|)"]])))