# Bliss independence and highest single agent (HSA) for a single-dose combination, per plate and across experiments.
# R 4.3.3, base R.
d <- read.csv("viability.csv", check.names = FALSE)  # columns: experiment, treatment (control, A, B, A+B), well, signal
ctl <- tapply(d$signal[d$treatment == "control"], d$experiment[d$treatment == "control"], mean)
d$sf <- d$signal / ctl[as.character(d$experiment)]
plate <- function(e) { x <- d[d$experiment == e, ]; list(sf = sapply(c("A", "B", "A+B"), function(g) mean(x$sf[x$treatment == g])), x = x) }
f2 <- function(p) formatC(signif(p, 2), format = "g")

r <- plate(1); sf <- r$sf; x <- r$x
bliss <- sf[["A"]] * sf[["B"]]; hsa <- min(sf[["A"]], sf[["B"]])
cat(sprintf("Experiment 1: surviving fraction A %.3f, B %.3f, A+B %.3f; Bliss expectation %.3f; HSA (best single drug) %.3f\n",
            sf[["A"]], sf[["B"]], sf[["A+B"]], bliss, hsa))
best <- if (sf[["A"]] < sf[["B"]]) "A" else "B"
t1 <- t.test(x$sf[x$treatment == "A+B"], x$sf[x$treatment == best])
t2 <- t.test(x$sf[x$treatment == "A+B"], mu = bliss)
cat(sprintf("  Wells as n, HSA: A+B vs %s alone, Welch's t-test p = %s\n", best, f2(t1$p.value)))
cat(sprintf("  Wells as n, Bliss: A+B wells vs the Bliss expectation, one-sample t-test p = %s\n", f2(t2$p.value)))

cat("Each experiment: observed / Bliss-expected surviving fraction [observed / HSA]:\n")
lb <- c(); lh <- c()
for (e in sort(unique(d$experiment))) {
  sf <- plate(e)$sf
  lb <- c(lb, log(sf[["A+B"]]) - log(sf[["A"]]) - log(sf[["B"]])); lh <- c(lh, log(sf[["A+B"]]) - log(min(sf[["A"]], sf[["B"]])))
  cat(sprintf("  Experiment %d: %.2f [%.2f]\n", e, exp(tail(lb, 1)), exp(tail(lh, 1))))
}
for (lab in c("Bliss", "HSA")) {
  v <- if (lab == "Bliss") lb else lh; tt <- t.test(v)
  cat(sprintf("Across %d experiments, observed / %s expected: %.2f (95%% CI %.2f to %.2f), one-sample t-test on the log ratios p = %s\n",
              length(v), lab, exp(mean(v)), exp(tt$conf.int[1]), exp(tt$conf.int[2]), f2(tt$p.value)))
}