# Bacterial burden (CFU per spleen), treated vs control: raw CFU, log10 CFU with values below the limit of detection
# set to LOD/2, a rank test, and a censored (Tobit) model. R 4.3.3 with the survival package (shipped with R).
suppressPackageStartupMessages(library(survival))
d <- read.csv("cfu.csv")                             # one row per mouse: group, mouse, dilution (10^-d plated), colonies
f2 <- function(p) formatC(signif(p, 2), format = "g")
VOL_ML <- 0.1; HOMOG_ML <- 1.0                       # 100 uL plated; organ homogenised in 1 mL
d$cfu <- d$colonies * 10^d$dilution / VOL_ML * HOMOG_ML
d$lod <- 10^d$dilution / VOL_ML * HOMOG_ML          # one colony on the plate read: 10 CFU per organ on the undiluted plate
d$below <- d$colonies == 0; LOD <- if (any(d$below)) max(d$lod[d$below]) else min(d$lod)
d$log_half <- log10(ifelse(d$below, d$lod / 2, d$cfu))
for (g in sort(unique(d$group))) { x <- d[d$group == g, ]
  det <- x$cfu[!x$below]
  cat(sprintf("%s: %d mice, %d below the limit of detection (%g CFU); detected mice range %g to %g CFU\n", g, nrow(x), sum(x$below), LOD, min(det), max(det))) }
c <- d$group == "control"; tr <- d$group == "treated"
t <- t.test(d$cfu[tr], d$cfu[c])
cat(sprintf("A. Raw CFU, Welch t-test: mean treated %.0f, control %.0f, p = %s\n", mean(d$cfu[tr]), mean(d$cfu[c]), f2(t$p.value)))
t <- t.test(d$log_half[tr], d$log_half[c])
cat(sprintf("B. log10 CFU, below-LOD values set to LOD/2, Welch t-test: difference %.2f log10 (95%% CI %.2f to %.2f), p = %s\n",
            mean(d$log_half[tr]) - mean(d$log_half[c]), t$conf.int[1], t$conf.int[2], f2(t$p.value)))
u <- suppressWarnings(wilcox.test(d$cfu[tr], d$cfu[c], exact = FALSE, correct = TRUE))
cat(sprintf("C. Mann-Whitney test on CFU (below-LOD values tied at the bottom): p = %s\n", f2(u$p.value)))
d$y <- log10(ifelse(d$below, d$lod, d$cfu)); d$trt <- as.numeric(tr)
f1 <- survreg(Surv(y, !below, type = "left") ~ trt, data = d, dist = "gaussian")
f0 <- survreg(Surv(y, !below, type = "left") ~ 1, data = d, dist = "gaussian")
lr <- 2 * (f1$loglik[2] - f0$loglik[2])
b <- coef(f1)[["trt"]]; se <- sqrt(vcov(f1)["trt", "trt"])
cat(sprintf("D. Censored (Tobit) model of log10 CFU, LOD treated as censoring: difference %.2f log10 (95%% Wald CI %.2f to %.2f), likelihood-ratio chi2(1) = %.2f, p = %s\n",
            b, b - 1.96 * se, b + 1.96 * se, lr, f2(pchisq(lr, 1, lower.tail = FALSE))))