curveRcore: Model Mathematics and Preprocessing Pipeline

Overview

curveRcore provides the shared mathematical foundation for the curveR ecosystem. Both curveRfreq (frequentist) and curveRbayes (Bayesian) depend on this package for model functions, inverse functions, gradient closures, data transformations, and the calibration_result output class.

This vignette demonstrates:

The five canonical model forms
Inverse (back-calculation) functions
The data preprocessing pipeline
Prediction grid generation

The Five Model Forms

All models use the convention: a = lower asymptote, d = upper asymptote, b > 0 = slope, c = inflection point, g > 0 = asymmetry. All are monotonically increasing when a < d and b > 0.

x <- seq(-2, 6, length.out = 300)
par(mfrow = c(3, 2), mar = c(4, 4, 2, 1))

# Shared parameters
a <- 100; b <- 1.0; c_val <- 2; d <- 20000; g <- 1.3

plot(x, logistic4(x, a, b, c_val, d), type = "l", lwd = 2,
     main = "logistic4 (4PL)", ylab = "Response", xlab = "log10(conc)")

plot(x, logistic5(x, a, b, c_val, d, g), type = "l", lwd = 2,
     main = "logistic5 (5PL)", ylab = "Response", xlab = "log10(conc)")

plot(x, loglogistic5(x, a, b, c_val, d, g), type = "l", lwd = 2,
     main = "loglogistic5 (Richards)", ylab = "Response", xlab = "log10(conc)")

plot(x, gompertz4(x, a, b, c_val, d), type = "l", lwd = 2,
     main = "gompertz4", ylab = "Response", xlab = "log10(conc)")

# loglogistic4 requires x > 0 (raw concentration)
x_pos <- 10^seq(-2, 4, length.out = 300)
plot(log10(x_pos), loglogistic4(x_pos, a, b = 2, c = 10, d), type = "l",
     lwd = 2, main = "loglogistic4 (Hill, raw conc)",
     ylab = "Response", xlab = "log10(conc)")

plot.new()
legend("center", bty = "n", cex = 1.2,
       legend = c("a = lower asymptote (baseline)",
                  "d = upper asymptote (saturation)",
                  "b > 0 = slope (all models)",
                  "c = inflection point",
                  "g > 0 = asymmetry (5-param only)"))

Round-Trip: Forward then Inverse

Every model has an analytical inverse. The round-trip property inv(forward(x)) == x is guaranteed within the dynamic range.

x_test <- seq(-1, 5, length.out = 20)
y <- logistic4(x_test, a = 100, b = 1.5, c = 2, d = 20000)
x_recovered <- inv_logistic4(y, a = 100, b = 1.5, c = 2, d = 20000)

all.equal(x_test, x_recovered)
#> [1] TRUE

First Derivatives (dy/dx)

All first derivatives are positive (monotonically increasing curves). The derivative peaks near the inflection point.

x <- seq(-2, 6, length.out = 300)
dydx <- dydx_logistic4(x, a = 100, b = 1.5, c = 2, d = 20000)

par(mfrow = c(1, 1))
plot(x, dydx, type = "l", lwd = 2,
     main = "dy/dx for logistic4 -- peaks at inflection",
     ylab = "dy/dx", xlab = "log10(conc)")
abline(v = 2, lty = 2, col = "gray50")
text(2.1, max(dydx) * 0.9, "c = 2 (inflection)", adj = 0)

Settings Constructors

curveRcore provides typed constructors for the three settings objects that both curveRfreq and curveRbayes use:

ac <- new_antigen_constraints(
  antigen        = "alpha",
  l_asy_min      = 0,
  l_asy_max      = 0,
  l_asy_method   = "default",
  std_curve_conc = 10000,
  pcov_threshold = 15
)

sp <- new_study_params(
  is_log_response    = TRUE,
  is_log_independent = TRUE,
  apply_prozone      = TRUE,
  blank_option       = "ignored"
)

fo <- new_fit_options(n_grid = 200, cv_x_max = 150)

str(ac)
#> List of 7
#>  $ antigen                     : chr "alpha"
#>  $ l_asy_min_constraint        : num 0
#>  $ l_asy_max_constraint        : num 0
#>  $ l_asy_constraint_method     : chr "default"
#>  $ standard_curve_concentration: num 10000
#>  $ pcov_threshold              : num 15
#>  $ std_error_blank             : NULL
#>  - attr(*, "class")= chr [1:2] "antigen_constraints" "list"
str(sp)
#> List of 4
#>  $ is_log_response   : logi TRUE
#>  $ is_log_independent: logi TRUE
#>  $ apply_prozone     : logi TRUE
#>  $ blank_option      : chr "ignored"
#>  - attr(*, "class")= chr [1:2] "study_params" "list"
str(fo)
#> List of 5
#>  $ model_names  : chr [1:5] "logistic4" "loglogistic4" "gompertz4" "logistic5" ...
#>  $ n_grid       : int 200
#>  $ cv_x_max     : num 150
#>  $ grid_min_conc: num 1e-04
#>  $ grid_max_conc: NULL
#>  - attr(*, "class")= chr [1:2] "fit_options" "list"

Prediction Grid

Both packages generate their prediction grid using the same function, ensuring the concentration points align for cross-method comparison:

grid <- generate_prediction_grid(ac, fo, is_log_independent = TRUE)

head(grid)
#>   log10_concentration concentration     x_fit
#> 1           -4.000000  0.0001000000 -4.000000
#> 2           -3.959799  0.0001096986 -3.959799
#> 3           -3.919598  0.0001203378 -3.919598
#> 4           -3.879397  0.0001320088 -3.879397
#> 5           -3.839196  0.0001448118 -3.839196
#> 6           -3.798995  0.0001588565 -3.798995
cat("Grid range: ", range(grid$log10_concentration), "\n")
#> Grid range:  -4 4
cat("Grid points:", nrow(grid), "\n")
#> Grid points: 200

params <- c(a = 100, b = 1.0, c = 2, d = 20000)
y_pred <- predict_grid_response(grid, "logistic4", params)

plot(grid$log10_concentration, y_pred, type = "l", lwd = 2,
     main = "Prediction grid: logistic4",
     xlab = "log10(concentration)", ylab = "Predicted response")

The calibration_result Output Class

Both curveRfreq and curveRbayes return calibration_result objects with identical structure, enabling direct comparison:

meta <- list(
  method = "frequentist",
  package = "curveRfreq",
  antigen = "alpha",
  plate = "plate_1",
  response_var = "mfi",
  independent_var = "concentration",
  is_log_response = TRUE,
  is_log_independent = TRUE
)

cr <- new_calibration_result(
  meta = meta,
  selection = list(best_model_name = "logistic4", criterion = "AIC")
)
cr
#> ── calibration_result (frequentist) ──
#>   Antigen : alpha
#>   Plate   : plate_1
#>   Feature : NA
#>   Package : curveRfreq v0.1.0
#>   Models  : 0 attempted, 0 converged
#>   Best    : logistic4 (by AIC)
#>   Grid    : 0 points

Data Preprocessing Pipeline

The preprocess_standards() function applies all transforms in the canonical order used by both packages:

Compute concentration from dilution x undiluted standard
Prozone (hook effect) correction
Blank handling (ignored / included / subtracted)
Log10 response transform (with adaptive floor for non-positive values)

# Example using the bead_assay_example dataset:
data(bead_assay_example)
curve_data <- filter_by_curve_id(bead_assay_example, curve_id = 1)

ac_bead <- new_antigen_constraints(
  antigen = "alpha", std_curve_conc = 10000, pcov_threshold = 15
)

result <- preprocess_standards(
  data                 = curve_data$standards,
  antigen_settings     = ac_bead,
  response_variable    = "mfi",
  independent_variable = "concentration",
  is_log_response      = TRUE,
  is_log_independent   = TRUE,
  apply_prozone        = TRUE,
  blank_option         = "ignored"
)

str(result$data[, c("concentration", "mfi")], vec.len = 3)

Next Steps

curveRfreq: Uses curveRcore model math + transforms to fit the NLS ensemble and select by AIC.
curveRbayes: Uses the same model math in Stan models with hierarchical priors across plates.
Cross-method comparison: compare_calibrations() merges freq and bayes calibration_result objects for paired analysis.