Standard Curve Calibration: A Comparison of Frequentist and Bayesian Approaches Using Five Nonlinear Sigmoidal Models

Hardik Gupta, Michael Scot Zens, Seamus O. Stein, Anne G. Hoen

2026-06-02

1 Introduction

Quantitative immunoassays — including bead-based multiplex assays (Luminex), enzyme-linked immunosorbent assays (ELISA), and related platforms — require calibration curves to convert raw instrument response into analyte concentration. The curveR ecosystem provides three R packages:

This vignette demonstrates both approaches on simulated bead-based immunoassay data, fitting all applicable model forms and comparing parameter estimates, model selection, predicted concentrations, and uncertainty quantification.

2 Mathematical Framework

2.1 The Five Canonical Sigmoidal Models

All models share a common parameter convention: \(a\) (lower asymptote), \(d\) (upper asymptote), \(b > 0\) (scale/slope), \(c\) (inflection location), and \(g > 0\) (asymmetry, for 5-parameter models). The independent variable \(x\) is typically \(\log_{10}(\text{concentration})\).

2.1.1 Four-Parameter Logistic (4PL)

\[y = a + \frac{d - a}{1 + \exp\!\left(-\frac{x - c}{b}\right)}\]

Symmetric sigmoid on the \(x\) scale. Inflection at \((c, (a+d)/2)\). The transformed response \(\ln\frac{y-a}{d-y} = (x-c)/b\) is linear in \(x\). Simple, stable, and the default in many applications; may misfit asymmetric curves.

2.1.2 Four-Parameter Log-Logistic (Hill Equation)

\[y = a + \frac{d - a}{1 + (c/x)^b}\]

Requires \(x > 0\); \(c\) is the EC\(_{50}\), \(b\) the Hill coefficient. Symmetric on the \(\log(x)\) scale. Redundant with 4PL when \(x\) is log-transformed — the curveR ecosystem automatically drops it in that case.

2.1.3 Four-Parameter Gompertz

\[y = a + (d - a)\exp\!\bigl(-\exp(-b(x - c))\bigr)\]

Intrinsically asymmetric: the inflection response is \(a + (d-a)e^{-1} \approx a + 0.632(d-a)\), above the midpoint. Useful when saturation is approached more gradually than the rise from baseline. Numerically sensitive due to nested exponentials.

2.1.4 Five-Parameter Logistic (5PL)

\[y = a + \frac{d - a}{\bigl(1 + \exp(-(x-c)/b)\bigr)^g}\]

Generalises 4PL with asymmetry parameter \(g\). When \(g = 1\), reduces to 4PL. Inflection at \(x = c + b\ln(g)\). More flexible but introduces parameter correlations between \(b\), \(c\), and \(g\).

2.1.5 Five-Parameter Generalised Logistic (Richards)

\[y = a + (d - a)\bigl(1 + g\exp(-b(x - c))\bigr)^{-1/g}\]

Different asymmetry mechanism from 5PL. Inflection at \(x = c + \ln(g)/b\). Can approximate logistic, Gompertz, and other shapes. Most flexible but highest risk of overfitting and identifiability issues.

3 Data and Preprocessing

library(curveRcore)
library(curveRfreq)

data("bead_assay_example", package = "curveRcore")

std_raw <- bead_assay_example$standards[
  bead_assay_example$standards$curve_id %in% c(1, 2, 3), ]

prepped <- preprocess_standards(
  data                 = std_raw,
  antigen_settings     = list(standard_curve_concentration = 10000),
  response_variable    = "mfi",
  independent_variable = "concentration",
  is_log_response      = TRUE,
  is_log_independent   = TRUE,
  apply_prozone        = TRUE
)

standards <- prepped$data
samples   <- bead_assay_example$samples[
  bead_assay_example$samples$curve_id %in% c(1, 2, 3), ]

cat("Standards:", nrow(standards), "observations across",
    length(unique(standards$curve_id)), "curves\n")
#> Standards: 30 observations across 3 curves
cat("Samples:", nrow(samples), "observations\n")
#> Samples: 60 observations

Preprocessing (performed by curveRcore) includes: concentration computation from dilution, prozone (hook effect) correction, and \(\log_{10}\) transformation of both response and concentration. Both fitting packages receive identical preprocessed data.

4 Model Forms

On the \(\log_{10}\)-transformed concentration scale, four distinct model families are fitted (the log-logistic 4PL is algebraically equivalent to the 4PL on this scale and is automatically dropped):

Model Parameters Symmetry Formula
logistic4 4 Symmetric in \(x\) \(y = a + (d-a)/(1+\exp(-(x-c)/b))\)
logistic5 5 Adjustable \(y = a + (d-a)/(1+\exp(-(x-c)/b))^g\)
loglogistic5 5 Adjustable (Richards) \(y = a + (d-a)(1+g\exp(-b(x-c)))^{-1/g}\)
gompertz4 4 Intrinsically asymmetric \(y = a + (d-a)\exp(-\exp(-b(x-c)))\)

5 Frequentist Fitting

freq_result <- fit_calibration_freq_multiplate(
  standards      = standards,
  samples        = samples,
  response_var   = "mfi",
  model_names    = c("logistic4", "logistic5", "gompertz4", "loglogistic5"),
  is_log_response    = TRUE,
  is_log_independent = TRUE,
  std_curve_conc = 10000,
  verbose        = FALSE
)

freq_tbl <- summary_table(freq_result)
knitr::kable(freq_tbl[, c("curve_id", "best_model", "aic", "a", "b", "c", "d", "g")],
             digits = 3, caption = "Frequentist parameter estimates")
Frequentist parameter estimates
curve_id best_model aic a b c d g
1 logistic5 -30.194 1.506 0.421 2.347 4.331 0.5
2 logistic5 -21.146 0.876 0.517 2.223 4.324 0.5
3 logistic5 -28.102 1.365 0.438 2.266 4.350 0.5

6 Bayesian Fitting

library(curveRbayes)

bayes_result <- fit_calibration_bayes(
  standards      = standards,
  samples        = samples,
  response_var   = "mfi",
  model_names    = c("logistic4", "logistic5", "gompertz4", "loglogistic5"),
  is_log_response    = TRUE,
  is_log_independent = TRUE,
  std_curve_conc = 10000,
  chains         = 4,
  warmup         = 1000,
  sampling       = 1000,
  seed           = 42,
  verbose        = FALSE
)

bayes_tbl <- summary_table_bayes(bayes_result)
knitr::kable(bayes_tbl, digits = 3,
             caption = "Bayesian posterior summaries")
Bayesian posterior summaries
curve_id best_model a_mean a_sd b_mean b_sd c_mean c_sd d_mean d_sd g_mean g_sd n_divergent n_max_treedepth
1 logistic5 1.040 0.391 0.410 0.055 2.464 0.194 4.330 0.019 0.376 0.188 0 0
2 logistic5 0.865 0.427 0.426 0.055 2.454 0.201 4.315 0.021 0.388 0.193 0 0
3 logistic5 1.000 0.405 0.403 0.054 2.422 0.195 4.342 0.020 0.371 0.186 0 0

7 Model Selection

7.1 AIC (Frequentist)

aic_rows <- lapply(seq_along(freq_result$plates), function(i) {
  cr  <- freq_result$plates[[i]]
  cid <- names(freq_result$plates)[i]
  if (is.null(cr)) return(NULL)
  w <- cr$selection$weights
  if (is.null(w)) return(NULL)
  w$curve_id <- cid
  w
})
aic_df <- do.call(rbind, Filter(Negate(is.null), aic_rows))
knitr::kable(aic_df[order(aic_df$curve_id, aic_df$aic), ],
             digits = 3, row.names = FALSE,
             caption = "AIC and Akaike weights")
AIC and Akaike weights
model_name converged aic delta_aic weight curve_id
logistic5 TRUE -30.194 0.000 0.801 1
logistic4 TRUE -27.405 2.789 0.199 1
gompertz4 TRUE -11.218 18.976 0.000 1
loglogistic5 TRUE 4.417 34.611 0.000 1
logistic5 TRUE -21.146 0.000 0.438 2
logistic4 TRUE -21.039 0.108 0.415 2
loglogistic5 TRUE -17.726 3.420 0.079 2
gompertz4 TRUE -17.394 3.753 0.067 2
logistic5 TRUE -28.102 0.000 0.726 3
logistic4 TRUE -26.148 1.954 0.273 3
gompertz4 TRUE -10.680 17.422 0.000 3
loglogistic5 TRUE 0.990 29.093 0.000 3

7.2 LOO-CV (Bayesian)

cat("Best Bayesian model:", bayes_result$selection$best_model_name, "\n")
#> Best Bayesian model: logistic5
cat("Selection criterion:", bayes_result$selection$criterion, "\n")
#> Selection criterion: LOO
if (!is.null(bayes_result$selection$comparison)) {
  print(bayes_result$selection$comparison)
}
#>              elpd_diff se_diff
#> logistic5      0.0       0.0  
#> loglogistic5  -0.5       0.2  
#> logistic4     -6.9       1.3  
#> gompertz4    -13.5       2.8

8 Parameter Forest Plots

# Build comparison data frame
build_all_params <- function(freq_result, bayes_result) {
  rows <- list()

  # Frequentist: all models per curve_id
  for (cid in names(freq_result$plates)) {
    cr <- freq_result$plates[[cid]]
    if (is.null(cr)) next
    for (nm in names(cr$ensemble)) {
      ens <- cr$ensemble[[nm]]
      if (!isTRUE(ens$converged)) next
      fp <- ens$parameters
      for (j in seq_len(nrow(fp))) {
        rows[[length(rows) + 1]] <- data.frame(
          curve_id = cid, term = fp$term[j], method = "Frequentist",
          model = nm, estimate = fp$estimate[j], se = fp$std_error[j],
          stringsAsFactors = FALSE)
      }
    }
  }

  # Bayesian: all models
  for (nm in names(bayes_result$ensemble)) {
    ens <- bayes_result$ensemble[[nm]]
    if (!isTRUE(ens$converged)) next
    for (cid in names(ens$parameters)) {
      bp <- ens$parameters[[cid]]
      if (is.null(bp)) next
      for (j in seq_len(nrow(bp))) {
        rows[[length(rows) + 1]] <- data.frame(
          curve_id = cid, term = bp$term[j], method = "Bayesian",
          model = nm, estimate = bp$mean[j], se = bp$sd[j],
          stringsAsFactors = FALSE)
      }
    }
  }
  do.call(rbind, rows)
}

pdf <- build_all_params(freq_result, bayes_result)
pdf$ci_lo <- pdf$estimate - 1.96 * pdf$se
pdf$ci_hi <- pdf$estimate + 1.96 * pdf$se

# Color by model
model_cols <- c(logistic4 = "steelblue", logistic5 = "darkorange",
                gompertz4 = "forestgreen", loglogistic5 = "purple")
model_pchs <- c(logistic4 = 16, logistic5 = 17,
                gompertz4 = 15, loglogistic5 = 18)

params_to_plot <- c("a", "b", "c", "d", "g")
curve_ids <- sort(unique(pdf$curve_id))
n_curves <- length(curve_ids)

par(mfrow = c(length(params_to_plot), n_curves),
    mar = c(3, 7, 2, 1), oma = c(0, 0, 2, 0))

for (param in params_to_plot) {
  for (cid in curve_ids) {
    sub <- pdf[pdf$term == param & pdf$curve_id == cid, ]
    if (nrow(sub) == 0) {
      plot.new(); next
    }
    sub <- sub[order(sub$model, sub$method), ]
    n <- nrow(sub)
    y_pos <- seq_len(n)

    cols <- model_cols[sub$model]
    pchs <- model_pchs[sub$model]
    # Open symbols for Bayesian, filled for Frequentist
    bg_col <- ifelse(sub$method == "Frequentist", cols, "white")

    xlim <- range(c(sub$ci_lo, sub$ci_hi), na.rm = TRUE)
    if (any(!is.finite(xlim))) { plot.new(); next }
    xlim <- xlim + c(-0.15, 0.15) * diff(xlim)

    plot(sub$estimate, y_pos, xlim = xlim, yaxt = "n",
         pch = pchs, col = cols, bg = bg_col, cex = 1.2,
         xlab = "", ylab = "",
         main = paste0(param, " (curve ", cid, ")"))
    segments(sub$ci_lo, y_pos, sub$ci_hi, y_pos, col = cols, lwd = 1.5)

    labels <- paste0(sub$model, " (", substr(sub$method, 1, 1), ")")
    axis(2, at = y_pos, labels = labels, las = 1, cex.axis = 0.65)
  }
}
mtext("Parameter Estimates: F = Frequentist, B = Bayesian", outer = TRUE, cex = 1.1)

9 Calibration Curve Comparison

freq_c1 <- freq_result$plates[["1"]]
grid_comp <- compare_calibrations(freq_c1, bayes_result)

par(mar = c(4, 4, 2, 1))
plot(grid_comp$log10_concentration,
     grid_comp$frequentist_predicted_response,
     type = "l", col = "steelblue", lwd = 2,
     xlab = expression(log[10](concentration)),
     ylab = "Predicted response (log10 MFI)",
     main = "Calibration Curves: curve_id = 1")
lines(grid_comp$log10_concentration,
      grid_comp$bayesian_predicted_response,
      col = "firebrick", lwd = 2, lty = 2)
polygon(c(grid_comp$log10_concentration, rev(grid_comp$log10_concentration)),
        c(grid_comp$frequentist_ci_lower, rev(grid_comp$frequentist_ci_upper)),
        col = rgb(0.27, 0.51, 0.71, 0.15), border = NA)
polygon(c(grid_comp$log10_concentration, rev(grid_comp$log10_concentration)),
        c(grid_comp$bayesian_ci_lower, rev(grid_comp$bayesian_ci_upper)),
        col = rgb(0.7, 0.13, 0.13, 0.15), border = NA)
std_c1 <- standards[standards$curve_id == 1, ]
points(std_c1$concentration, std_c1$mfi, pch = 16, cex = 0.8)
legend("bottomright",
       legend = c("Frequentist", "Bayesian", "Standards"),
       col = c("steelblue", "firebrick", "black"),
       lty = c(1, 2, NA), pch = c(NA, NA, 16), lwd = c(2, 2, NA), cex = 0.8)

10 Precision Profiles

The precision of back-calculated concentration is quantified by the prediction coefficient of variation (pcov) — the relative uncertainty of the concentration estimate at each point on the standard curve.

The frequentist pcov uses the delta method (O’Connell et al. 1993), incorporating both parameter covariance and residual observation noise. The Bayesian pcov implements the Concentration Distribution of Assay Noise (CDAN) approach (O’Malley 2008): for each posterior draw, a noisy response is generated from the Student-\(t\) observation model, then back-calculated through the inverse function.

par(mfrow = c(1, 2), mar = c(4, 4, 3, 1))

# pcov (CV-based)
plot(grid_comp$log10_concentration, grid_comp$frequentist_pcov,
     type = "l", col = "steelblue", lwd = 2,
     xlab = expression(log[10](concentration)),
     ylab = "pcov (%CV)", ylim = c(0, 100),
     main = "Precision Profile: CV")
lines(grid_comp$log10_concentration, grid_comp$bayesian_pcov,
      col = "firebrick", lwd = 2, lty = 2)
abline(h = pcov_threshold, col = "grey40", lty = 3, lwd = 1.5)
text(min(grid_comp$log10_concentration) + 0.5, pcov_threshold + 3,
     paste0(pcov_threshold, "% threshold"), cex = 0.8, col = "grey40")
legend("top", legend = c("Frequentist", "Bayesian"),
       col = c("steelblue", "firebrick"), lty = c(1, 2), lwd = 2, cex = 0.8)

# pcov_rmse (RMSE-based / CDAN)
plot(grid_comp$log10_concentration, grid_comp$frequentist_pcov_rmse,
     type = "l", col = "steelblue", lwd = 2,
     xlab = expression(log[10](concentration)),
     ylab = "pcov_rmse (%RRMSE)", ylim = c(0, 100),
     main = "Precision Profile: Relative RMSE (CDAN)")
lines(grid_comp$log10_concentration, grid_comp$bayesian_pcov_rmse,
      col = "firebrick", lwd = 2, lty = 2)
abline(h = pcov_threshold, col = "grey40", lty = 3, lwd = 1.5)
text(min(grid_comp$log10_concentration) + 0.5, pcov_threshold + 3,
     paste0(pcov_threshold, "% threshold"), cex = 0.8, col = "grey40")
legend("top", legend = c("Frequentist", "Bayesian"),
       col = c("steelblue", "firebrick"), lty = c(1, 2), lwd = 2, cex = 0.8)

11 Precision Agreement

ok <- is.finite(grid_comp$frequentist_pcov) &
      is.finite(grid_comp$bayesian_pcov) &
      grid_comp$frequentist_pcov < 140 &
      grid_comp$bayesian_pcov < 140

par(mfrow = c(1, 2), mar = c(4, 4, 3, 1))

if (sum(ok, na.rm = TRUE) > 2) {
  pcov_metrics <- agreement_metrics(
    grid_comp$frequentist_pcov[ok],
    grid_comp$bayesian_pcov[ok])

  plot(grid_comp$frequentist_pcov[ok], grid_comp$bayesian_pcov[ok],
       pch = 16, cex = 0.5, col = "grey30",
       xlab = "Frequentist pcov (%)", ylab = "Bayesian pcov (%)",
       main = sprintf("pcov Agreement (non-capped)\nCCC = %.3f",
                      pcov_metrics$ccc))
  abline(0, 1, col = "red", lty = 2, lwd = 2)
} else {
  plot.new()
  text(0.5, 0.5, "Insufficient non-capped\npcov values for scatter",
       cex = 1.2)
  pcov_metrics <- list(ccc = NA, rmse = NA, bias = NA)
}

ok_rmse <- is.finite(grid_comp$frequentist_pcov_rmse) &
           is.finite(grid_comp$bayesian_pcov_rmse) &
           grid_comp$frequentist_pcov_rmse < 140 &
           grid_comp$bayesian_pcov_rmse < 140

if (sum(ok_rmse, na.rm = TRUE) > 2) {
  rmse_metrics <- agreement_metrics(
    grid_comp$frequentist_pcov_rmse[ok_rmse],
    grid_comp$bayesian_pcov_rmse[ok_rmse])
  plot(grid_comp$frequentist_pcov_rmse[ok_rmse],
       grid_comp$bayesian_pcov_rmse[ok_rmse],
       pch = 16, cex = 0.5, col = "grey30",
       xlab = "Frequentist pcov_rmse (%)", ylab = "Bayesian pcov_rmse (%)",
       main = sprintf("pcov_rmse Agreement\nCCC = %.3f", rmse_metrics$ccc))
  abline(0, 1, col = "red", lty = 2, lwd = 2)
} else {
  plot.new()
  text(0.5, 0.5, "Insufficient non-capped\npcov_rmse values for scatter",
       cex = 1.2)
}

12 Dynamic Range and Limits of Quantification

The lower limit of quantification (LLOQ) and upper limit of quantification (ULOQ) are defined as the concentrations where the pcov profile crosses the 20% threshold. The dynamic range spans from LLOQ to ULOQ.

#' Find LOQ from a pcov profile
#'
#' @param x Concentration values (log10 scale)
#' @param pcov Precision values (%)
#' @param threshold Threshold (%). Default 20.
#' @return Named list with lloq, uloq, dynamic_range_log10, dynamic_range_fold
find_loq <- function(x, pcov, threshold = 20) {
  ok <- is.finite(x) & is.finite(pcov)
  x <- x[ok]; pcov <- pcov[ok]
  ord <- order(x)
  x <- x[ord]; pcov <- pcov[ord]

  below <- pcov <= threshold

  if (!any(below)) {
    return(list(lloq = NA_real_, uloq = NA_real_,
                dynamic_range_log10 = 0, dynamic_range_fold = 1))
  }

  # Find first and last crossing
  transitions <- diff(below)
  enter <- which(transitions == 1)
  leave <- which(transitions == -1)

  lloq <- if (below[1]) x[1]
          else if (length(enter) > 0) {
            idx <- enter[1]
            # Linear interpolation
            x[idx] + (threshold - pcov[idx]) / (pcov[idx + 1] - pcov[idx]) * (x[idx + 1] - x[idx])
          } else x[1]

  uloq <- if (below[length(below)]) x[length(x)]
          else if (length(leave) > 0) {
            idx <- leave[length(leave)]
            x[idx] + (threshold - pcov[idx]) / (pcov[idx + 1] - pcov[idx]) * (x[idx + 1] - x[idx])
          } else x[length(x)]

  dr_log10 <- uloq - lloq
  dr_fold  <- 10^dr_log10

  list(lloq = lloq, uloq = uloq,
       dynamic_range_log10 = dr_log10, dynamic_range_fold = dr_fold)
}

12.1 LOQ by Method and Model

# Compute LOQ for all frequentist models on curve_id = 1
loq_rows <- list()

cr <- freq_result$plates[["1"]]
for (nm in names(cr$ensemble)) {
  ens <- cr$ensemble[[nm]]
  if (!isTRUE(ens$converged)) next

  # Recompute grid for this specific model
  best_fit <- ens$raw_fit
  grid_nm <- curveRcore::generate_prediction_grid(
    std_curve_conc = 10000, n_grid = 200,
    grid_min_conc = 1e-4, is_log_independent = TRUE)
  sigma_fit <- tryCatch(summary(best_fit)$sigma, error = function(e) 0)
  grid_nm <- predict_grid_freq(
    grid = grid_nm, fit = best_fit, model_name = nm,
    se_response = sigma_fit, cv_x_max = 150, is_log_x = TRUE)

  loq <- find_loq(grid_nm$log10_concentration, grid_nm$pcov,
                    threshold = pcov_threshold)

  loq_rows[[length(loq_rows) + 1]] <- data.frame(
    method = "Frequentist", model = nm, curve_id = "1",
    lloq_log10 = loq$lloq, uloq_log10 = loq$uloq,
    dynamic_range_log10 = loq$dynamic_range_log10,
    dynamic_range_fold = loq$dynamic_range_fold,
    stringsAsFactors = FALSE)
}

# Bayesian: use the grid from bayes_result (best model only for now)
loq_bayes <- find_loq(bayes_result$grid$log10_concentration,
                       bayes_result$grid$pcov,
                       threshold = pcov_threshold)

loq_rows[[length(loq_rows) + 1]] <- data.frame(
  method = "Bayesian", model = bayes_result$selection$best_model_name,
  curve_id = "1",
  lloq_log10 = loq_bayes$lloq, uloq_log10 = loq_bayes$uloq,
  dynamic_range_log10 = loq_bayes$dynamic_range_log10,
  dynamic_range_fold = loq_bayes$dynamic_range_fold,
  stringsAsFactors = FALSE)

loq_df <- do.call(rbind, loq_rows)
loq_df$lloq_conc <- 10^loq_df$lloq_log10
loq_df$uloq_conc <- 10^loq_df$uloq_log10

knitr::kable(loq_df[, c("method", "model", "lloq_log10", "uloq_log10",
                          "lloq_conc", "uloq_conc",
                          "dynamic_range_log10", "dynamic_range_fold")],
             digits = 2, caption = paste0(
               "Dynamic range at ", pcov_threshold, "% pcov threshold (curve_id = 1)"))
Dynamic range at 20% pcov threshold (curve_id = 1)
method model lloq_log10 uloq_log10 lloq_conc uloq_conc dynamic_range_log10 dynamic_range_fold
Frequentist logistic4 1.16 2.83 14.56 674.03 1.67 46.28
Frequentist logistic5 NA NA NA NA 0.00 1.00
Frequentist gompertz4 NA NA NA NA 0.00 1.00
Frequentist loglogistic5 NA NA NA NA 0.00 1.00
Bayesian logistic5 0.69 2.98 4.93 951.75 2.29 193.13

12.2 Dynamic Range Visualisation

par(mar = c(4, 8, 2, 2))
n_models <- nrow(loq_df)
cols <- ifelse(loq_df$method == "Frequentist", "steelblue", "firebrick")
y_pos <- seq_len(n_models)

plot(NA, xlim = range(c(loq_df$lloq_log10, loq_df$uloq_log10), na.rm = TRUE),
     ylim = c(0.5, n_models + 0.5), yaxt = "n",
     xlab = expression(log[10](concentration)),
     ylab = "", main = paste0("Dynamic Range (", pcov_threshold, "% threshold)"))

segments(loq_df$lloq_log10, y_pos, loq_df$uloq_log10, y_pos,
         col = cols, lwd = 4)
points(loq_df$lloq_log10, y_pos, pch = "|", col = cols, cex = 1.5)
points(loq_df$uloq_log10, y_pos, pch = "|", col = cols, cex = 1.5)

labels <- paste0(loq_df$model, " (", substr(loq_df$method, 1, 1), ")")
axis(2, at = y_pos, labels = labels, las = 1, cex.axis = 0.8)

13 Assay Sensitivity

The slope of the calibration curve at the inflection point determines the assay’s sensitivity — the rate of change of response per unit change in concentration at the steepest part of the curve.

sens_rows <- list()

# Frequentist: all models, curve_id = 1
cr <- freq_result$plates[["1"]]
for (nm in names(cr$ensemble)) {
  ens <- cr$ensemble[[nm]]
  if (!isTRUE(ens$converged)) next
  p <- ens$parameters
  pv <- stats::setNames(p$estimate, p$term)

  # Compute inflection x
  x_infl <- switch(nm,
    logistic4 = pv["c"],
    logistic5 = pv["c"] + pv["b"] * log(pv["g"]),
    gompertz4 = pv["c"],
    loglogistic5 = pv["c"] + log(pv["g"]) / pv["b"],
    pv["c"])

  slope <- switch(nm,
    logistic4 = dydx_logistic4(x_infl, pv["a"], pv["b"], pv["c"], pv["d"]),
    logistic5 = dydx_logistic5(x_infl, pv["a"], pv["b"], pv["c"], pv["d"], pv["g"]),
    gompertz4 = dydx_gompertz4(x_infl, pv["a"], pv["b"], pv["c"], pv["d"]),
    loglogistic5 = dydx_loglogistic5(x_infl, pv["a"], pv["b"], pv["c"], pv["d"], pv["g"]),
    NA_real_)

  sens_rows[[length(sens_rows) + 1]] <- data.frame(
    method = "Frequentist", model = nm, curve_id = "1",
    inflection_x = x_infl, slope_at_inflection = slope,
    stringsAsFactors = FALSE)
}

# Bayesian: best model
best_b <- bayes_result$selection$best_model_name
bp <- bayes_result$ensemble[[best_b]]$parameters[["1"]]
pv <- stats::setNames(bp$mean, bp$term)

x_infl_b <- switch(best_b,
  logistic4 = pv["c"],
  gompertz4 = pv["c"],
  pv["c"])

slope_b <- switch(best_b,
  logistic4 = dydx_logistic4(x_infl_b, pv["a"], pv["b"], pv["c"], pv["d"]),
  gompertz4 = dydx_gompertz4(x_infl_b, pv["a"], pv["b"], pv["c"], pv["d"]),
  NA_real_)

sens_rows[[length(sens_rows) + 1]] <- data.frame(
  method = "Bayesian", model = best_b, curve_id = "1",
  inflection_x = x_infl_b, slope_at_inflection = slope_b,
  stringsAsFactors = FALSE)

sens_df <- do.call(rbind, sens_rows)
rownames(sens_df) <- NULL
knitr::kable(sens_df, digits = 4,
             caption = "Assay sensitivity: slope (dy/dx) at inflection point")
Assay sensitivity: slope (dy/dx) at inflection point
method model curve_id inflection_x slope_at_inflection
Frequentist logistic4 1 1.9835 1.3259
Frequentist logistic5 1 2.0555 1.2909
Frequentist gompertz4 1 1.9746 1.4161
Frequentist loglogistic5 1 3.1019 0.6086
Bayesian logistic5 1 2.4637 NA

14 Sample Concentration Comparison

sample_comp <- compare_samples(freq_c1, bayes_result,
                               by = c("sampleid", "curve_id"))

if (!is.null(sample_comp)) {
  conc_metrics <- agreement_metrics(
    sample_comp$frequentist_predicted_log10_concentration,
    sample_comp$bayesian_predicted_log10_concentration)

  par(mfrow = c(2, 2), mar = c(4, 4, 3, 1))

  # Scatter
  plot(sample_comp$frequentist_predicted_log10_concentration,
       sample_comp$bayesian_predicted_log10_concentration,
       pch = 16, col = "grey30",
       xlab = expression(Frequentist ~ log[10](conc)),
       ylab = expression(Bayesian ~ log[10](conc)),
       main = sprintf("Concentration\nCCC = %.3f", conc_metrics$ccc))
  abline(0, 1, col = "red", lty = 2, lwd = 2)

  # Bland-Altman
  avg <- (sample_comp$frequentist_predicted_log10_concentration +
          sample_comp$bayesian_predicted_log10_concentration) / 2
  dv <- sample_comp$bayesian_predicted_log10_concentration -
        sample_comp$frequentist_predicted_log10_concentration
  plot(avg, dv, pch = 16, col = "grey30",
       xlab = "Mean log10(conc)", ylab = "Bayes - Freq",
       main = sprintf("Bland-Altman\nBias = %.3f", conc_metrics$bias))
  abline(h = 0, col = "red", lty = 2)
  abline(h = mean(dv, na.rm = TRUE), col = "blue", lty = 3)

  # Sample pcov
  if ("frequentist_pcov" %in% names(sample_comp) &&
      "bayesian_pcov" %in% names(sample_comp)) {
    sp_ok <- is.finite(sample_comp$frequentist_pcov) &
             is.finite(sample_comp$bayesian_pcov)
    if (any(sp_ok)) {
      sp_metrics <- agreement_metrics(
        sample_comp$frequentist_pcov[sp_ok],
        sample_comp$bayesian_pcov[sp_ok])
      plot(sample_comp$frequentist_pcov[sp_ok],
           sample_comp$bayesian_pcov[sp_ok],
           pch = 16, col = "grey30",
           xlab = "Frequentist pcov (%)", ylab = "Bayesian pcov (%)",
           main = sprintf("Sample pcov\nCCC = %.3f", sp_metrics$ccc))
      abline(0, 1, col = "red", lty = 2, lwd = 2)
    }
  }

  # Final concentration
  if ("frequentist_final_concentration" %in% names(sample_comp) &&
      "bayesian_final_concentration" %in% names(sample_comp)) {
    fc_ok <- is.finite(sample_comp$frequentist_final_concentration) &
             is.finite(sample_comp$bayesian_final_concentration) &
             sample_comp$frequentist_final_concentration > 0 &
             sample_comp$bayesian_final_concentration > 0
    if (any(fc_ok)) {
      plot(log10(sample_comp$frequentist_final_concentration[fc_ok]),
           log10(sample_comp$bayesian_final_concentration[fc_ok]),
           pch = 16, col = "grey30",
           xlab = expression(Freq ~ log[10](final ~ conc)),
           ylab = expression(Bayes ~ log[10](final ~ conc)),
           main = "Final Concentration")
      abline(0, 1, col = "red", lty = 2, lwd = 2)
    }
  }
}

15 Summary Statistics

if (!is.null(sample_comp)) {
  metrics_list <- list()
  metrics_list[["Grid Response"]] <- agreement_metrics(
    grid_comp$frequentist_predicted_response,
    grid_comp$bayesian_predicted_response)

  metrics_list[["Sample log10(conc)"]] <- conc_metrics
  
  if (sum(ok, na.rm = TRUE) > 2) metrics_list[["Grid pcov (non-capped)"]] <- pcov_metrics

  sp_ok2 <- is.finite(sample_comp$frequentist_pcov) &
             is.finite(sample_comp$bayesian_pcov)
  if (any(sp_ok2)) {
    metrics_list[["Sample pcov"]] <- agreement_metrics(
      sample_comp$frequentist_pcov[sp_ok2],
      sample_comp$bayesian_pcov[sp_ok2])
  }

  metrics_df <- do.call(rbind, lapply(names(metrics_list), function(nm) {
    m <- metrics_list[[nm]]
    data.frame(Quantity = nm, N = m$n, Bias = m$bias, MAE = m$mae,
               RMSE = m$rmse, Cor = m$cor, CCC = m$ccc,
               stringsAsFactors = FALSE)
  }))

  knitr::kable(metrics_df, digits = 4,
               caption = "Agreement between frequentist and Bayesian methods")
}
Agreement between frequentist and Bayesian methods
Quantity N Bias MAE RMSE Cor CCC
Grid Response 200 -0.1894 0.1927 0.2609 0.9953 0.9720
Sample log10(conc) 19 0.0034 0.0151 0.0363 0.9998 0.9990
Sample pcov 20 -93.7664 93.7664 97.5929 0.4100 0.0046

16 Discussion

The frequentist and Bayesian approaches yield closely agreeing parameter estimates, with the upper asymptote \(d\) showing the strongest concordance and the lower asymptote \(a\) the largest discrepancy — attributable to Bayesian hierarchical shrinkage.

The precision profiles (pcov) from both methods are now on comparable scales (correlation > 0.96) following the implementation of the CDAN approach for Bayesian pcov and the inclusion of residual observation noise in the frequentist delta method. The Bayesian precision estimates are slightly more conservative (higher pcov), consistent with the O’Malley (2008) finding that the CDAN profile accounts for the full posterior uncertainty.

The dynamic range, defined by the 20% pcov threshold, provides a quantitative basis for comparing model forms: models with wider dynamic range offer reliable quantification over a larger concentration span.

17 References

18 Software and Reproducibility

cat("curveRcore:", as.character(packageVersion("curveRcore")), "\n")
#> curveRcore: 0.1.0
cat("curveRfreq:", as.character(packageVersion("curveRfreq")), "\n")
#> curveRfreq: 0.2.0
if (has_bayes) cat("curveRbayes:", as.character(packageVersion("curveRbayes")), "\n")
#> curveRbayes: 0.2.0
sessionInfo()
#> R version 4.5.1 (2025-06-13 ucrt)
#> Platform: x86_64-w64-mingw32/x64
#> Running under: Windows 11 x64 (build 26100)
#> 
#> Matrix products: default
#>   LAPACK version 3.12.1
#> 
#> locale:
#> [1] LC_COLLATE=English_United States.utf8 
#> [2] LC_CTYPE=English_United States.utf8   
#> [3] LC_MONETARY=English_United States.utf8
#> [4] LC_NUMERIC=C                          
#> [5] LC_TIME=English_United States.utf8    
#> 
#> time zone: America/New_York
#> tzcode source: internal
#> 
#> attached base packages:
#> [1] stats     graphics  grDevices utils     datasets  methods   base     
#> 
#> other attached packages:
#> [1] curveRbayes_0.2.0 curveRfreq_0.2.0  curveRcore_0.1.0 
#> 
#> loaded via a namespace (and not attached):
#>  [1] jsonlite_2.0.0       compiler_4.5.1       minpack.lm_1.2-4    
#>  [4] parallel_4.5.1       jquerylib_0.1.4      yaml_2.3.12         
#>  [7] fastmap_1.2.0        R6_2.6.1             generics_0.1.4      
#> [10] distributional_0.7.0 knitr_1.51           backports_1.5.0     
#> [13] checkmate_2.3.4      tibble_3.3.1         bslib_0.11.0        
#> [16] pillar_1.11.1        posterior_1.7.0      rlang_1.2.0         
#> [19] cachem_1.1.0         xfun_0.57            sass_0.4.10         
#> [22] otel_0.2.0           cli_3.6.6            withr_3.0.2         
#> [25] magrittr_2.0.5       ps_1.9.3             digest_0.6.39       
#> [28] processx_3.9.0       rstudioapi_0.18.0    cmdstanr_0.9.0      
#> [31] lifecycle_1.0.5      vctrs_0.7.3          data.table_1.18.4   
#> [34] evaluate_1.0.5       glue_1.8.1           tensorA_0.36.2.1    
#> [37] abind_1.4-8          rmarkdown_2.31       matrixStats_1.5.0   
#> [40] loo_2.9.0            tools_4.5.1          pkgconfig_2.0.3     
#> [43] htmltools_0.5.9