# Calcium imaging responses, WT vs KO: neurons as n, mice as n, and the effect of selecting "responsive" neurons on the same
# trials that measure them. R 4.3.3, base R.
d <- read.csv("responses.csv")                       # one row per trial: mouse, genotype, neuron, trial, dff
f2 <- function(p) formatC(signif(p, 2), format = "g")
nm <- aggregate(dff ~ mouse + genotype + neuron, data = d, FUN = mean)       # one response per neuron
ko <- nm$genotype == "KO"; t1 <- t.test(nm$dff[ko], nm$dff[!ko])
cat(sprintf("Neurons as n (%d vs %d neurons): KO - WT %.3f, Welch's t-test p = %s\n", sum(ko), sum(!ko), mean(nm$dff[ko]) - mean(nm$dff[!ko]), f2(t1$p.value)))
mm <- aggregate(dff ~ mouse + genotype, data = nm, FUN = mean); km <- mm$genotype == "KO"; t2 <- t.test(mm$dff[km], mm$dff[!km])
cat(sprintf("Mice as n (5 vs 5 mice, mean of each mouse's neurons): KO - WT %.3f (95%% CI %.3f to %.3f), p = %s\n",
            mean(mm$dff[km]) - mean(mm$dff[!km]), t2$conf.int[1], t2$conf.int[2], f2(t2$p.value)))
respond <- function(x) {
  a <- aggregate(dff ~ mouse + genotype + neuron, data = x, FUN = function(v) c(m = mean(v), p = t.test(v, mu = 0)$p.value))
  data.frame(a[, 1:3], mean = a$dff[, "m"], p = a$dff[, "p"])
}
allt <- respond(d); sel <- allt$p < 0.05 & allt$mean > 0
cat(sprintf("Responsive on all 20 trials: %d of %d neurons (%.1f%%); their mean ΔF/F on the same trials %.3f; mean ΔF/F of all neurons %.3f (true average 0.05)\n",
            sum(sel), nrow(allt), 100 * mean(sel), mean(allt$mean[sel]), mean(allt$mean)))
so <- respond(d[d$trial %% 2 == 1, ]); se <- respond(d[d$trial %% 2 == 0, ]); s_odd <- so$p < 0.05 & so$mean > 0
cat(sprintf("Selected on odd trials (%d neurons): mean ΔF/F on odd trials %.3f, on held-out even trials %.3f\n", sum(s_odd), mean(so$mean[s_odd]), mean(se$mean[s_odd])))
allt$resp <- sel; pm <- aggregate(resp ~ mouse + genotype, data = allt, FUN = mean); kp <- pm$genotype == "KO"
t3 <- t.test(pm$resp[kp], pm$resp[!kp])
cat(sprintf("Share of responsive neurons per mouse: WT %.1f%%, KO %.1f%%; mice as n, Welch's t-test p = %s\n", 100 * mean(pm$resp[!kp]), 100 * mean(pm$resp[kp]), f2(t3$p.value)))