# Histopathology scores (0-4), treated vs control: fields as n, a t-test on mouse scores, a rank test on mouse scores,
# a proportional-odds (ordinal) model, and the mean of the field scores per mouse. R 4.3.3 with MASS (shipped with R).
suppressPackageStartupMessages(library(MASS))
d <- read.csv("scores.csv")                          # one row per field: group, mouse, field, score (0-4)
f2 <- function(p) formatC(signif(p, 2), format = "g")
m <- aggregate(score ~ group + mouse, data = d, FUN = median)          # one score per mouse: the median of its 9 fields
for (g in sort(unique(m$group))) { x <- m[m$group == g, ]
  cat(sprintf("%s: mouse median scores %s; median of mice %g\n", g, paste(sapply(0:4, function(k) sprintf("%d: %d", k, sum(x$score == k))), collapse = ", "), median(x$score))) }
a <- d$score[d$group == "treated"]; b <- d$score[d$group == "control"]
cat(sprintf("A. Fields as n (%d vs %d), Mann-Whitney test: p = %s\n", length(a), length(b), f2(suppressWarnings(wilcox.test(a, b, exact = FALSE, correct = TRUE))$p.value)))
mt <- m$score[m$group == "treated"]; mc <- m$score[m$group == "control"]; t <- t.test(mt, mc)
cat(sprintf("B. Mouse scores treated as numbers, Welch t-test: mean %.2f vs %.2f, p = %s\n", mean(mt), mean(mc), f2(t$p.value)))
cat(sprintf("C. Mouse scores (%d vs %d), Mann-Whitney test: p = %s\n", length(mt), length(mc), f2(suppressWarnings(wilcox.test(mt, mc, exact = FALSE, correct = TRUE))$p.value)))
mf <- aggregate(score ~ group + mouse, data = d, FUN = mean); t <- t.test(mf$score[mf$group == "treated"], mf$score[mf$group == "control"])
cat(sprintf("D. Mean of the 9 field scores per mouse, Welch t-test: %.2f vs %.2f, p = %s\n", mean(mf$score[mf$group == "treated"]), mean(mf$score[mf$group == "control"]), f2(t$p.value)))
m$y <- factor(m$score, ordered = TRUE); m$trt <- as.numeric(m$group == "treated")
f1 <- polr(y ~ trt, data = m, method = "logistic"); f0 <- polr(y ~ 1, data = m, method = "logistic")
lr <- deviance(f0) - deviance(f1)
cat(sprintf("E. Proportional-odds model of mouse scores: odds ratio for a higher score, treated vs control %.2f, likelihood-ratio chi2(1) = %.2f, p = %s\n",
            exp(coef(f1)[["trt"]]), lr, f2(pchisq(lr, 1, lower.tail = FALSE))))