# Dendritic spine density, knockout vs wild-type: segments, neurons or mice as n, and a Poisson model of all segments
# with log(length) as an offset. R 4.3.3, base R.
d <- read.csv("spines.csv")                          # one row per segment: genotype, mouse, neuron, segment, length_um, spines
f2 <- function(p) formatC(signif(p, 2), format = "g")
d$density10 <- 10 * d$spines / d$length_um           # spines per 10 um
nr <- aggregate(cbind(spines, length_um) ~ genotype + mouse + neuron, data = d, FUN = sum); nr$density10 <- 10 * nr$spines / nr$length_um
ms <- aggregate(cbind(spines, length_um) ~ genotype + mouse, data = d, FUN = sum); ms$density10 <- 10 * ms$spines / ms$length_um
for (g in sort(unique(ms$genotype))) { x <- ms[ms$genotype == g, ]
  cat(sprintf("%s: %d mice, %d spines on %.0f um; mouse densities %.1f to %.1f per 10 um\n", g, nrow(x), sum(x$spines), sum(x$length_um), min(x$density10), max(x$density10))) }
welch <- function(df, label) { k <- df$density10[df$genotype == "knockout"]; w <- df$density10[df$genotype == "wild-type"]; t <- t.test(k, w)
  cat(sprintf("%s (%d vs %d), Welch t-test: knockout / wild-type %.2f, p = %s\n", label, length(k), length(w), mean(k) / mean(w), f2(t$p.value))) }
welch(d, "A. Segments as n"); welch(nr, "B. Neurons as n")
k <- ms$density10[ms$genotype == "knockout"]; w <- ms$density10[ms$genotype == "wild-type"]; t <- t.test(k, w)
cat(sprintf("C. Mice as n (%d vs %d), Welch t-test on spines per 10 um: difference %.2f (95%% CI %.2f to %.2f), p = %s\n", length(k), length(w), mean(k) - mean(w), t$conf.int[1], t$conf.int[2], f2(t$p.value)))
d$ko <- as.numeric(d$genotype == "knockout")
f1 <- glm(spines ~ ko + offset(log(length_um)), family = poisson, data = d); f0 <- glm(spines ~ 1 + offset(log(length_um)), family = poisson, data = d)
lr <- deviance(f0) - deviance(f1)
cat(sprintf("D. Poisson model of all segments, log(length) offset, mice ignored: rate ratio %.2f, likelihood-ratio chi2(1) = %.2f, p = %s\n", exp(coef(f1)[["ko"]]), lr, f2(pchisq(lr, 1, lower.tail = FALSE))))