# Litter effects: a linear mixed model of every pup with litter as a random effect, sex and litter size as covariates,
# and Satterthwaite degrees of freedom. R 4.3.3 with lme4 and lmerTest. R only.
suppressPackageStartupMessages(library(lmerTest))
d <- read.csv("pups.csv"); d$litter <- paste(d$treatment, d$litter)
fit <- lmer(weight_g ~ treatment + sex + litter_size + (1 | litter), data = d)
s <- summary(fit)$coef["treatmenttreated", ]; q <- qt(0.975, s[["df"]])
vc <- as.data.frame(VarCorr(fit)); icc <- vc$vcov[1] / sum(vc$vcov)
cat(sprintf("Mixed model (litter random effect, sex and litter size): difference %.2f g (95%% CI %.2f to %.2f), Satterthwaite df %.1f, p = %s\n",
            s[["Estimate"]], s[["Estimate"]] - q * s[["Std. Error"]], s[["Estimate"]] + q * s[["Std. Error"]], s[["df"]], formatC(signif(s[["Pr(>|t|)"]], 2), format = "g")))
cat(sprintf("Share of the remaining variation that lies between litters (intraclass correlation): %.2f\n", icc))
