# Tumour growth: vehicle vs drug, 8 mice each, measured every 3 days; mice culled at 1500 mm3. R 4.3.3, base R.
d <- read.csv("tumours.csv")
d$logv <- log(d$volume_mm3)
groups <- c("vehicle", "drug")
cat("Mice still measured on each day:\n")
for (g in groups) {
  cn <- table(d$day[d$group == g])
  cat(sprintf("  %-7s %s\n", g, paste(sprintf("d%s:%d", names(cn), as.integer(cn)), collapse = " ")))
}
welch <- function(a, b) { t <- t.test(b, a); c(unname(t$estimate[1] - t$estimate[2]), t$conf.int, t$p.value) }
last <- d[d$day == 21, ]
a <- last$logv[last$group == "vehicle"]; b <- last$logv[last$group == "drug"]
w <- welch(a, b)
cat(sprintf("Day 21, surviving mice only (biased): n = %d vehicle, %d drug; geometric means %.0f and %.0f mm3; T/C = %.2f, p = %.4f\n",
            length(a), length(b), exp(mean(a)), exp(mean(b)), exp(w[1]), w[4]))
cat("Separate test on each day (log volume, Welch):\n")
days <- sort(unique(d$day))
ps <- sapply(days, function(dd) { x <- d[d$day == dd, ]; welch(x$logv[x$group == "vehicle"], x$logv[x$group == "drug"])[4] })
cat("  ", paste(sprintf("d%d:%.3f", days, ps), collapse = " "), "\n", sep = "")
slope <- function(m) unname(coef(lm(logv ~ day, data = m))[2])
s <- sapply(split(d, d$mouse), slope); grp <- sapply(split(d$group, d$mouse), `[`, 1)
sv <- s[grp == "vehicle"]; sd_ <- s[grp == "drug"]
w <- welch(sv, sd_)
cat("Growth rate per mouse (slope of log volume per day, all measurements up to culling):\n")
cat(sprintf("  mean rate: vehicle %.4f, drug %.4f per day; doubling time %.1f and %.1f days\n",
            mean(sv), mean(sd_), log(2) / mean(sv), log(2) / mean(sd_)))
cat(sprintf("  difference (drug - vehicle) %.4f per day (95%% CI %.4f to %.4f), p = %s; rate ratio %.2f\n",
            w[1], w[2], w[3], formatC(signif(w[4], 2), format = "g"), mean(sd_) / mean(sv)))