# Organoid size, drug vs control: the same organoids analysed with organoids, wells or batches as n. R 4.3.3, base R.
d <- read.csv("organoids.csv")                       # one row per organoid: batch, condition, well, area_um2
f2 <- function(p) formatC(signif(p, 2), format = "g")
d$log_area <- log(d$area_um2)
w <- aggregate(log_area ~ batch + condition + well, data = d, FUN = mean); w$size <- aggregate(log_area ~ batch + condition + well, data = d, FUN = length)$log_area
b <- tapply(w$log_area, list(w$batch, w$condition), mean)
for (k in rownames(b)) cat(sprintf("Batch %s: geometric mean area control %s um2, drug %s um2 (ratio %.2f)\n", k,
                                   formatC(exp(b[k, "control"]), format = "f", digits = 0, big.mark = ","), formatC(exp(b[k, "drug"]), format = "f", digits = 0, big.mark = ","), exp(b[k, "drug"] - b[k, "control"])))
cat(sprintf("Organoids measured: %d (%d to %d per well)\n", nrow(d), min(w$size), max(w$size)))
ratio_test <- function(x, y, label) { t <- t.test(x, y)
  cat(sprintf("%s (%d vs %d), Welch t-test on log area: drug / control %.2f, p = %s\n", label, length(x), length(y), exp(mean(x) - mean(y)), f2(t$p.value))) }
ratio_test(d$log_area[d$condition == "drug"], d$log_area[d$condition == "control"], "A. Organoids as n")
ratio_test(w$log_area[w$condition == "drug"], w$log_area[w$condition == "control"], "B. Wells as n")
r <- b[, "drug"] - b[, "control"]; t <- t.test(r)
cat(sprintf("C. Batches as n, log ratio per batch (paired, n = %d): drug / control %.2f (95%% CI %.2f to %.2f), p = %s\n", length(r), exp(mean(r)), exp(t$conf.int[1]), exp(t$conf.int[2]), f2(t$p.value)))