suppressPackageStartupMessages(library(drc))
d <- read.csv("viability.csv")
m <- aggregate(viability ~ experiment + compound + conc_M, data = d, FUN = mean)   # average the triplicate wells
fits <- do.call(rbind, lapply(split(m, list(m$experiment, m$compound)), function(g) {
  f <- drm(viability ~ conc_M, data = g, fct = LL.4())         # four-parameter log-logistic
  data.frame(experiment = g$experiment[1], compound = g$compound[1], pIC50 = -log10(coef(f)[["e:(Intercept)"]]))
}))
w <- reshape(fits, idvar = "experiment", timevar = "compound", direction = "wide")
w <- w[order(w$experiment), ]
cat("pIC50 per experiment (-log10 IC50 in mol/L):\n")
for (i in seq_len(nrow(w))) cat(sprintf("  experiment %d: A %.2f, B %.2f, A - B %.2f\n", w$experiment[i], w$pIC50.A[i], w$pIC50.B[i], w$pIC50.A[i] - w$pIC50.B[i]))
cat(sprintf("Geometric mean IC50: A %.0f nM, B %.0f nM\n", 10^(-mean(w$pIC50.A)) * 1e9, 10^(-mean(w$pIC50.B)) * 1e9))
diff <- w$pIC50.A - w$pIC50.B
t <- t.test(diff)
cat(sprintf("Paired (same plates): IC50 ratio B/A = %.2f (95%% CI %.1f to %.1f), t = %.2f, df = %.0f, p = %.3f\n",
            10^mean(diff), 10^t$conf.int[1], 10^t$conf.int[2], t$statistic, t$parameter, t$p.value))
wt <- t.test(w$pIC50.A, w$pIC50.B)
cat(sprintf("Welch, pairing ignored: IC50 ratio B/A = %.2f (95%% CI %.1f to %.1f), t = %.2f, df = %.1f, p = %.3f\n",
            10^(mean(w$pIC50.A) - mean(w$pIC50.B)), 10^wt$conf.int[1], 10^wt$conf.int[2], wt$statistic, wt$parameter, wt$p.value))
# A common shortcut, shown for comparison: every well pooled, F test of one shared IC50
full <- drm(viability ~ conc_M, compound, data = d, fct = LL.4())
shared <- drm(viability ~ conc_M, compound, data = d, fct = LL.4(), pmodels = list(~compound - 1, ~compound - 1, ~compound - 1, ~1))
a <- anova(shared, full, details = FALSE)
pv <- a[2, "p value"]
cat(sprintf("Shortcut, all %d wells pooled, F test of one shared IC50: F(%d, %d) = %.1f, %s\n", nrow(d), a$Df[2], a$ModelDf[2], a[2, "F value"],
            if (pv < 1e-15) "p < 1e-15" else sprintf("p = %.1e", pv)))