# Clonogenic survival, sensitiser vs control: linear-quadratic fit per experiment and condition, SF2 compared across
# experiments (paired), versus one Poisson or quasi-Poisson model of all dishes pooled. R 4.3.3, base R.
d <- read.csv("colonies.csv")                        # one row per dish: experiment, condition, dose_gy, dish, seeded, colonies
f2 <- function(p) formatC(signif(p, 2), format = "g")
lq <- function(x) {
  m <- tapply(x$colonies / x$seeded, x$dose_gy, mean); pe <- m[["0"]]; sf <- m[names(m) != "0"] / pe; D <- as.numeric(names(sf))
  b <- coef(lm(-log(sf) ~ 0 + D + I(D^2))); c(pe = pe, a = b[[1]], b = b[[2]], sf2 = exp(-(2 * b[[1]] + 4 * b[[2]])))
}
ex <- sort(unique(d$experiment)); res <- list()
for (e in ex) for (cn in c("control", "sensitiser")) res[[paste(e, cn)]] <- lq(d[d$experiment == e & d$condition == cn, ])
for (e in ex) { a <- res[[paste(e, "control")]]; s <- res[[paste(e, "sensitiser")]]
  cat(sprintf("Experiment %d: control PE %.2f, alpha %.3f, beta %.3f, SF2 %.3f; sensitiser PE %.2f, SF2 %.3f\n", e, a[["pe"]], a[["a"]], a[["b"]], a[["sf2"]], s[["pe"]], s[["sf2"]])) }
r <- sapply(ex, function(e) log(res[[paste(e, "sensitiser")]][["sf2"]] / res[[paste(e, "control")]][["sf2"]])); tt <- t.test(r)
cat(sprintf("SF2, sensitiser / control, one ratio per experiment (paired, n = %d): %.2f (95%% CI %.2f to %.2f), p = %s\n", length(r), exp(mean(r)), exp(tt$conf.int[1]), exp(tt$conf.int[2]), f2(tt$p.value)))
# Dose enhancement ratio at SF = 0.1: dose for control / dose for sensitiser, from each experiment's LQ fits (>1 = sensitised)
d10 <- function(v) (-v[["a"]] + sqrt(v[["a"]]^2 - 4 * v[["b"]] * log(0.1))) / (2 * v[["b"]])
q <- sapply(ex, function(e) log(d10(res[[paste(e, "control")]]) / d10(res[[paste(e, "sensitiser")]]))); tq <- t.test(q)
cat(sprintf("Dose enhancement ratio at SF = 0.1, one per experiment (paired, n = %d): %.2f (95%% CI %.2f to %.2f), p = %s\n", length(q), exp(mean(q)), exp(tq$conf.int[1]), exp(tq$conf.int[2]), f2(tq$p.value)))
# Pooled Poisson model of all dishes: log mean = log seeded + a plating efficiency for each experiment and condition
# + dose and dose^2, plus condition x (dose, dose^2) in the larger model. Dishes are treated as independent counts.
d$s <- as.numeric(d$condition == "sensitiser"); d$pe <- interaction(d$experiment, d$condition)
m0 <- glm(colonies ~ 0 + pe + dose_gy + I(dose_gy^2) + offset(log(seeded)), family = poisson, data = d)
m1 <- update(m0, . ~ . + s:dose_gy + s:I(dose_gy^2))
lr <- deviance(m0) - deviance(m1); dfr <- df.residual(m1)
phi <- sum(residuals(m1, type = "pearson")^2) / dfr; Fq <- (lr / 2) / phi     # Pearson dispersion, as glm(family = quasipoisson)
cat(sprintf("All %d dishes pooled, Poisson model (plating efficiency per experiment and condition): likelihood-ratio chi2(2) = %.1f, p = %s; Pearson dispersion %.1f\n",
            nrow(d), lr, f2(pchisq(lr, 2, lower.tail = FALSE)), phi))
cat(sprintf("Same model with quasi-Poisson scaling: F(2, %d) = %.2f, p = %s\n", dfr, Fq, f2(pf(Fq, 2, dfr, lower.tail = FALSE))))