# Mixed model for cells nested in mice: everything a report needs. R 4.3.3, lme4 1.1-35.1, lmerTest 3.1-3.
suppressMessages(library(lmerTest))
d <- read.csv("cells.csv"); d$group <- factor(d$group, levels = c("vehicle", "drug"))
fit <- lmer(y ~ group + (1 | mouse), data = d, REML = TRUE)
cat(sprintf("Data: %d cells from %d mice (%d per group)\n", nrow(d), length(unique(d$mouse)), length(unique(d$mouse[d$group == "drug"]))))
cat(sprintf("Singular fit: %s\n", isSingular(fit)))
co <- summary(fit, ddf = "Satterthwaite")$coefficients["groupdrug", ]
ci <- confint(fit, parm = "groupdrug", method = "Wald")
vc <- as.data.frame(VarCorr(fit)); sm <- vc$sdcor[vc$grp == "mouse"]; se <- vc$sdcor[vc$grp == "Residual"]
cat(sprintf("Fixed effect, drug - vehicle: %.3f (SE %.3f; normal-approximation (Wald z) 95%% CI %.3f to %.3f)\n", co["Estimate"], co["Std. Error"], ci[1], ci[2]))
cat(sprintf("Satterthwaite: t(%.1f) = %.2f, p = %.4f; 95%% CI with these df %.3f to %.3f\n", co["df"], co["t value"], co["Pr(>|t|)"],
            co["Estimate"] - qt(0.975, co["df"]) * co["Std. Error"], co["Estimate"] + qt(0.975, co["df"]) * co["Std. Error"]))
co2 <- summary(fit, ddf = "Kenward-Roger")$coefficients["groupdrug", ]
cat(sprintf("Kenward-Roger: t(%.1f) = %.2f, p = %.4f\n", co2["df"], co2["t value"], co2["Pr(>|t|)"]))
cat(sprintf("Random effects: between-mouse SD %.3f; residual (between-cell) SD %.3f; ICC %.3f\n", sm, se, sm^2 / (sm^2 + se^2)))
vf <- var(as.vector(model.matrix(fit) %*% fixef(fit)))
cat(sprintf("R-squared (Nakagawa & Schielzeth): marginal %.3f, conditional %.3f\n", vf / (vf + sm^2 + se^2), (vf + sm^2) / (vf + sm^2 + se^2)))
m <- aggregate(y ~ mouse + group, data = d, mean)
tt <- t.test(y ~ factor(group, levels = c("drug", "vehicle")), data = m, var.equal = TRUE)
cat(sprintf("Check: t-test on the 16 mouse means: difference %.3f, t(%d) = %.2f, p = %.4f\n", mean(m$y[m$group == "drug"]) - mean(m$y[m$group == "vehicle"]), tt$parameter, tt$statistic, tt$p.value))