# Anxiety tests, knockout vs wild-type: eight measures from the elevated plus maze (EPM) and open field (OF), and four
# ways to draw a conclusion about anxiety. R 4.3.3, base R.
d <- read.csv("anxiety.csv")                         # one row per mouse
f2 <- function(p) formatC(signif(p, 2), format = "g")
ANX <- c("pct_open_time", "pct_open_entries", "open_entries", "pct_centre_time", "centre_entries")   # higher = less anxious
ACT <- c("closed_entries", "distance_m", "rearing")                                                  # activity
ko <- d[d$genotype == "knockout", ]; wt <- d[d$genotype == "wild-type", ]
p <- sapply(c(ANX, ACT), function(cl) t.test(ko[[cl]], wt[[cl]])$p.value)
for (cl in c(ANX, ACT)) cat(sprintf("%s: wild-type %.1f, knockout %.1f, Welch p = %s\n", cl, mean(wt[[cl]]), mean(ko[[cl]]), f2(p[[cl]])))
pa <- p[ANX]
cat(sprintf("A. Any of the %d anxiety measures with p < 0.05: %s (%d of %d)\n", length(ANX), if (any(pa < 0.05)) "yes" else "no", sum(pa < 0.05), length(ANX)))
cat(sprintf("B. Holm-adjusted over the %d anxiety measures: smallest adjusted p = %s\n", length(ANX), f2(min(p.adjust(pa, "holm")))))
t <- t.test(ko$pct_open_time, wt$pct_open_time)
cat(sprintf("C. Primary outcome chosen in advance, %% time in open arms: difference %.1f points (95%% CI %.1f to %.1f), p = %s\n",
            mean(ko$pct_open_time) - mean(wt$pct_open_time), t$conf.int[1], t$conf.int[2], f2(t$p.value)))
Z <- c("pct_open_time", "pct_open_entries", "pct_centre_time")               # composite of the three percentage measures
z <- Reduce(`+`, lapply(Z, function(cl) (d[[cl]] - mean(wt[[cl]])) / sd(wt[[cl]]))) / length(Z)   # z-scores against wild-type
t <- t.test(z[d$genotype == "knockout"], z[d$genotype == "wild-type"])
cat(sprintf("D. Composite anxiety score (mean z of %s): difference %.2f SD (95%% CI %.2f to %.2f), p = %s\n",
            paste(Z, collapse = ", "), mean(z[d$genotype == "knockout"]) - mean(z[d$genotype == "wild-type"]), t$conf.int[1], t$conf.int[2], f2(t$p.value)))