# C. elegans lifespan, treated vs control: Kaplan-Meier medians, log-rank tests pooled and stratified by trial, and
# analyses with plates or trials as the units. R 4.3.3 with the survival package (shipped with R).
suppressPackageStartupMessages(library(survival))
d <- read.csv("lifespan.csv")                        # one row per worm: trial, condition, plate, day (of death or censoring), died (1) or censored (0)
f2 <- function(p) formatC(signif(p, 2), format = "g")
km_median <- function(x) unname(summary(survfit(Surv(day, died) ~ 1, data = x))$table["median"])
for (t in sort(unique(d$trial))) {
  g <- d[d$trial == t, ]
  cat(sprintf("Trial %d: median lifespan control %s days, treated %s days (%d worms censored)\n", t,
              format(km_median(g[g$condition == "control", ])), format(km_median(g[g$condition == "treated", ])), sum(g$died == 0)))
}
c1 <- survdiff(Surv(day, died) ~ condition, data = d)$chisq
cat(sprintf("A. All worms pooled (%d vs %d), log-rank test: chi2(1) = %.2f, p = %s\n", sum(d$condition == "treated"), sum(d$condition == "control"), c1, f2(pchisq(c1, 1, lower.tail = FALSE))))
c2 <- survdiff(Surv(day, died) ~ condition + strata(trial), data = d)$chisq
cat(sprintf("B. Log-rank test stratified by trial: chi2(1) = %.2f, p = %s\n", c2, f2(pchisq(c2, 1, lower.tail = FALSE))))
pl <- do.call(rbind, lapply(split(d, d$plate), function(x) data.frame(trial = x$trial[1], condition = x$condition[1], med = km_median(x))))
fit <- lm(log(med) ~ condition + factor(trial), data = pl); s <- coef(summary(fit))["conditiontreated", ]; ci <- confint(fit)["conditiontreated", ]
cat(sprintf("C. Plate medians as n (%d vs %d), log scale with trial as a block: treated / control %.2f (95%% CI %.2f to %.2f), p = %s\n",
            nrow(pl) / 2, nrow(pl) / 2, exp(s[["Estimate"]]), exp(ci[1]), exp(ci[2]), f2(s[["Pr(>|t|)"]])))
tr <- sapply(c("control", "treated"), function(cn) sapply(sort(unique(d$trial)), function(t) km_median(d[d$trial == t & d$condition == cn, ])))
r <- log(tr[, "treated"] / tr[, "control"]); t <- t.test(r)
cat(sprintf("D. Trial medians, log ratio per trial (paired, n = %d): treated / 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)))