# Values below the limit of detection (LOD = 5 pg/mL): detected = 0 means "< 5" (recorded as 5). R 4.3.3 with survival.
suppressMessages(library(survival))
LOD <- 5
d <- read.csv("cytokine.csv"); d$group <- factor(d$group, levels = c("control", "treated"))
for (k in levels(d$group)) { x <- d[d$group == k, ]
  cat(sprintf("%s: n = %d, below the LOD = %d, log-SD of detected values = %.3f\n", k, nrow(x), sum(x$detected == 0), sd(log(x$conc_pg_ml[x$detected == 1])))) }
sub <- ifelse(d$detected == 1, d$conc_pg_ml, LOD / 2)
lc <- log(sub[d$group == "control"]); lt <- log(sub[d$group == "treated"])
cat(sprintf("LOD/2 substitution, Welch's t-test on log values: ratio of geometric means %.2f, log-SD control %.3f, p = %.4f\n",
            exp(mean(lt) - mean(lc)), sd(lc), t.test(lt, lc)$p.value))
cat(sprintf("Mann-Whitney test (values below the LOD tied at the bottom): p = %.4f\n",
            suppressWarnings(wilcox.test(d$conc_pg_ml[d$group == "treated"], d$conc_pg_ml[d$group == "control"], exact = FALSE, correct = TRUE))$p.value))
fit <- survreg(Surv(conc_pg_ml, detected, type = "left") ~ group, data = d, dist = "lognormal")
b1 <- coef(fit)[["grouptreated"]]; s1 <- sqrt(vcov(fit)["grouptreated", "grouptreated"])
cat(sprintf("Left-censored log-normal model (maximum likelihood): ratio of geometric means %.2f (95%% CI %.2f to %.2f), log-SD %.3f, p = %.4f\n",
            exp(b1), exp(b1 - 1.96 * s1), exp(b1 + 1.96 * s1), fit$scale, 2 * pnorm(-abs(b1 / s1))))