# Comparing a marker's fluorescence between WT and KO: events as n, per-mouse arithmetic mean, per-mouse median on the
# log scale, and the per-mouse median with staining day as a block. R 4.3.3, base R.
e <- read.csv("events.csv")                          # one row per gated event: day, mouse, genotype, fluorescence
f2 <- function(p) formatC(signif(p, 2), format = "g")
t0 <- t.test(e$fluorescence[e$genotype == "KO"], e$fluorescence[e$genotype == "WT"])
cat(sprintf("Events as n (%d vs %d events): Welch's t-test p = %s\n", sum(e$genotype == "KO"), sum(e$genotype == "WT"), f2(t0$p.value)))
m <- aggregate(fluorescence ~ day + mouse + genotype, data = e, FUN = function(v) c(mean = mean(v), median = median(v)))
m <- data.frame(m[, 1:3], mean = m$fluorescence[, "mean"], median = m$fluorescence[, "median"])
m <- m[order(m$day, m$mouse), ]; ko <- m$genotype == "KO"
t1 <- t.test(m$mean[ko], m$mean[!ko])
cat(sprintf("Per-mouse arithmetic mean (6 vs 6 mice): KO / WT %.2f, Welch's t-test p = %s\n", mean(m$mean[ko]) / mean(m$mean[!ko]), f2(t1$p.value)))
lm_ <- log(m$median); t2 <- t.test(lm_[ko], lm_[!ko])
cat(sprintf("Per-mouse median, log scale, ignoring day: KO / WT %.2f (95%% CI %.2f to %.2f), p = %s\n",
            exp(mean(lm_[ko]) - mean(lm_[!ko])), exp(t2$conf.int[1]), exp(t2$conf.int[2]), f2(t2$p.value)))
fit <- lm(log(median) ~ factor(genotype, levels = c("WT", "KO")) + factor(day), data = m)
cf <- summary(fit)$coefficients[2, ]; dfe <- df.residual(fit); tq <- qt(0.975, dfe)
cat(sprintf("Per-mouse median, log scale, staining day as a block: KO / WT %.2f (95%% CI %.2f to %.2f), t(%d) p = %s\n",
            exp(cf[1]), exp(cf[1] - tq * cf[2]), exp(cf[1] + tq * cf[2]), dfe, f2(cf[4])))
cat("Per-day WT median fluorescence (mean of the two WT mice): ",
    paste(sapply(1:3, function(dd) sprintf("day %d %.0f", dd, mean(m$median[!ko & m$day == dd]))), collapse = "; "), "\n", sep = "")