# Weighted 4PL standard curve and interpolation of samples. R 4.3.3, base R (weighted least squares with optim).
d <- read.csv("plate.csv")
st <- d[d$type == "standard", ]
f4 <- function(x, bottom, top, ec50, slope) top + (bottom - top) / (1 + (x / ec50)^slope)
inv <- function(yv, p) p[["ec50"]] * ((p[["bottom"]] - p[["top"]]) / (yv - p[["top"]]) - 1)^(1 / p[["slope"]])
wss <- function(q) sum(((st$od - f4(st$conc_pg_ml, q[1], q[2], q[3], q[4])) / st$od)^2)   # weights 1/OD^2
o <- optim(c(0.05, 3, 100, 1), wss, method = "Nelder-Mead", control = list(maxit = 20000, reltol = 1e-14))
o <- optim(o$par, wss, method = "BFGS", control = list(maxit = 20000, reltol = 1e-15))
p <- setNames(o$par, c("bottom", "top", "ec50", "slope"))
cat(sprintf("Weighted 4PL: bottom %.4f, top %.4f, EC50 %.2f pg/mL, slope %.4f\n", p[["bottom"]], p[["top"]], p[["ec50"]], p[["slope"]]))
lev <- tapply(st$od, st$conc_pg_ml, mean); lev <- lev[order(-as.numeric(names(lev)))]
conc <- as.numeric(names(lev)); lloq <- min(conc); uloq <- max(conc)
cvs <- tapply(st$od, st$conc_pg_ml, function(v) 100 * sd(v) / mean(v))
cat("Back-calculated standards (mean of duplicate wells), % of nominal; CV of the duplicate ODs:\n")
for (i in seq_along(lev)) { cc <- conc[i]; acc <- 100 * inv(lev[[i]], p) / cc; lim <- if (cc %in% c(lloq, uloq)) 25 else 20
  cat(sprintf("  %8.2f pg/mL: %6.1f%%  %s (limit +/-%d%%); CV %.1f%%\n", cc, acc, if (abs(acc - 100) <= lim) "ok" else "OUTSIDE", lim, cvs[[as.character(cc)]])) }
lf <- coef(lm(od ~ conc_pg_ml, data = st))
cat("Straight line through all standards, back-calculated % of nominal:\n")
cat("  ", paste(sprintf("%s: %.0f%%", format(conc, drop0trailing = TRUE, trim = TRUE), 100 * ((lev - lf[[1]]) / lf[[2]]) / conc), collapse = ", "), "\n", sep = "")
lo_od <- f4(lloq, p[["bottom"]], p[["top"]], p[["ec50"]], p[["slope"]]); hi_od <- f4(uloq, p[["bottom"]], p[["top"]], p[["ec50"]], p[["slope"]])
flat <- p[["bottom"]] + 0.75 * (p[["top"]] - p[["bottom"]])   # example rule set in advance: upper quarter = plateau
cat(sprintf("Samples (mean OD of duplicate wells), reportable range %s to %s pg/mL:\n", format(lloq), format(uloq)))
sm <- tapply(d$od[d$type == "sample"], d$id[d$type == "sample"], mean)
for (k in names(sm)) { od <- sm[[k]]
  if (od < lo_od) cat(sprintf("  %s: OD %.4f -> below the LLOQ (report as < %s pg/mL)\n", k, od, format(lloq)))
  else if (od > hi_od) cat(sprintf("  %s: OD %.4f -> above the ULOQ (dilute and re-assay)\n", k, od))
  else if (od > flat) cat(sprintf("  %s: OD %.4f -> %.1f pg/mL, but on the upper plateau (dilute and re-assay)\n", k, od, inv(od, p)))
  else cat(sprintf("  %s: OD %.4f -> %.1f pg/mL\n", k, od, inv(od, p))) }