# Myelin g-ratio, knockout vs wild-type: axons as n, a per-axon regression of g-ratio on axon diameter with axons as n,
# the mean g-ratio per mouse, and each mouse's g-ratio at a common axon diameter (1 um) from a line of fibre on axon
# diameter. R 4.3.3, base R.
d <- read.csv("axons.csv")                           # one row per myelinated axon: genotype, mouse, axon, axon_um, fibre_um
f2 <- function(p) formatC(signif(p, 2), format = "g")
d$g <- d$axon_um / d$fibre_um
D0 <- 1.0                                            # common axon diameter for the adjusted comparison
m <- do.call(rbind, lapply(split(d, paste(d$genotype, d$mouse)), function(x) { b <- coef(lm(fibre_um ~ axon_um, data = x))   # line of fibre on axon diameter
  data.frame(genotype = x$genotype[1], mouse = x$mouse[1], g_mean = mean(x$g), g_at_1um = D0 / (b[[1]] + b[[2]] * D0), diam = median(x$axon_um), n = nrow(x)) }))
for (g in sort(unique(m$genotype))) { x <- m[m$genotype == g, ]
  cat(sprintf("%s: %d mice, %d axons; mean g-ratio %.3f; median axon diameter %.2f um\n", g, nrow(x), sum(x$n), mean(x$g_mean), mean(x$diam))) }
k <- d$g[d$genotype == "knockout"]; w <- d$g[d$genotype == "wild-type"]; t <- t.test(k, w)
cat(sprintf("A. Axons as n (%d vs %d), Welch t-test on g-ratio: difference %.3f, p = %s\n", length(k), length(w), mean(k) - mean(w), f2(t$p.value)))
d$ko <- as.numeric(d$genotype == "knockout"); s <- coef(summary(lm(g ~ ko + axon_um, data = d)))["ko", ]
cat(sprintf("B. Regression of g-ratio on genotype and axon diameter, axons as n: genotype difference %.3f, p = %s\n", s[["Estimate"]], f2(s[["Pr(>|t|)"]])))
for (v in c("g_mean", "g_at_1um")) {
  k <- m[[v]][m$genotype == "knockout"]; w <- m[[v]][m$genotype == "wild-type"]; t <- t.test(k, w)
  lab <- if (v == "g_mean") "C. Mice as n, mean g-ratio per mouse" else "D. Mice as n, each mouse's g-ratio at 1 um axon diameter"
  cat(sprintf("%s (%d vs %d), Welch t-test: difference %.3f (95%% CI %.3f to %.3f), p = %s\n", lab, length(k), length(w), mean(k) - mean(w), t$conf.int[1], t$conf.int[2], f2(t$p.value)))
}