1.1 The curveR Package Ecosystem

The ecosystem consists of three packages with a strict dependency hierarchy:

curveRcore is the shared foundation. It defines the five canonical sigmoidal forward models and their analytical inverses, first and second derivatives, analytical gradients for delta-method error propagation, NLS formula constructors, prediction grid utilities, preprocessing transforms, eligibility gating logic, and the calibration_result S3 output class used by both fitting packages. It contains two example datasets (bead_assay_example and elisa_assay_example).

curveRfreq implements frequentist calibration via multi-start Levenberg–Marquardt NLS. It receives preprocessed data from curveRcore, fits an ensemble of models, computes per-model precision profiles, applies eligibility gates, selects the best eligible model by AIC, and returns a calibration_result_multiplate object.

curveRbayes implements Bayesian hierarchical calibration via Hamiltonian Monte Carlo (HMC) as implemented in Stan (Carpenter et al. 2017), accessed through cmdstanr. It fits all curve IDs simultaneously in a single hierarchical model per model family, applies the same eligibility gates as curveRfreq, selects the best eligible model by LOO-CV expected log predictive density (Vehtari, Gelman, and Gabry 2017), and returns an identical calibration_result_multiplate structure.

The symmetric output structure means that downstream comparison, agreement assessment, and reporting code is identical regardless of which fitting engine produced the result.

2.1 Dataset Description

The bead_assay_example dataset bundled with curveRcore is a fully synthetic, reproducible bead-based immunoassay experiment designed to reflect the characteristics of real Luminex multiplex data. It contains measurements for two analytes — alpha and beta — each run on three replicate plates, giving six curve_id values in total.

Standard curves. Each plate contains ten standard concentrations spanning four orders of magnitude: \(\{0.001, 0.003, 0.01, 0.03, 0.1, 0.3, 1, 3, 10, 30\}\) AU/mL, using a 1:3 serial dilution scheme from an undiluted stock of 10,000 AU/mL. One standard point per concentration per plate is measured (no within-plate replication of standards).

True curve shapes. The alpha antigen was simulated from a four-parameter Gompertz model with population parameters \(a = 18\), \(d = 20{,}000\), \(b = 1.20\), \(c = 0.10\) AU/mL. The beta antigen was simulated from a five-parameter logistic (5PL) model with \(a = 25\), \(d = 28{,}000\), \(b = 1.80\), \(c = 1.00\) AU/mL, \(g = 0.80\). Plate-to-plate variability was introduced by drawing per-plate values of \(d\) and \(c\) from log-normal distributions centred on the population values with multiplicative SDs of \(\exp(\mathcal{N}(0, 0.04))\) and \(\exp(\mathcal{N}(0, 0.05))\) respectively.

Noise model. Measurement noise is heteroscedastic and proportional to signal: \(\text{MFI}_{\text{obs}} \sim \mathcal{N}(\mu, \sigma(\mu))\) where \(\sigma(\mu) = 5 + 0.04 \lvert \mu \rvert\), corresponding approximately to a 4% coefficient of variation at high signal levels, floored at MFI = 1.

Test samples. Each plate contains 20 patient samples at a fixed serum dilution of 1:2000, with effective concentrations drawn uniformly on the log scale over the interval \([0.5 \times \min(\text{standards}),\; 1.5 \times \max(\text{standards})]\). True serum concentrations (before dilution) are retained internally for validation but are not exposed in the package data object.

Blanks. Four blank wells per plate are measured and used during preprocessing for optional blank subtraction or inclusion as an extra zero-concentration standard point.

In this vignette we use the three curve IDs (1, 2, 3) corresponding to the alpha antigen across its three replicate plates.

4.1 Four-Parameter Logistic (4PL)

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

The 4PL is the default model in most regulatory immunoassay guidance (Findlay and Dillard 2007). It is a symmetric sigmoid in \(x\): the inflection point falls at \((c,\,(a+d)/2)\), the exact midpoint of the response range, and the response is equidistant from \(a\) and \(d\) at equal distances above and below \(c\). The log-odds transform \(\ln\!\tfrac{y-a}{d-y} = (x-c)/b\) is linear in \(x\), which gives the 4PL a particularly transparent interpretation: \(b\) is the reciprocal of the slope of the linearised curve and \(c\) is the EC\(_{50}\) on the log scale.

The analytical inverse is \[ x = c + b \ln\!\frac{y - a}{d - y}, \] which exists for \(y \in (a, d)\).

Strengths: Numerically stable; well-conditioned covariance matrix in typical immunoassay ranges; theoretically justified as a logistic growth curve.

Limitations: Enforces strict symmetry; cannot accommodate assays whose saturation plateau is approached more gradually than the lower baseline, or vice versa.

4.2 Four-Parameter Log-Logistic (Hill Equation)

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

This is the classical Hill dose-response model, where \(c\) is the EC\(_{50}\) on the raw (not log-transformed) concentration scale and \(b\) is the Hill coefficient. It requires \(x > 0\).

Critical note: When the independent variable has been \(\log_{10}\)-transformed (as in all standard curveR workflows with is_log_independent = TRUE), the log-logistic 4PL is algebraically identical to the 4PL above. curveRcore automatically detects this redundancy and removes loglogistic4 from the fitting ensemble when is_log_independent = TRUE. It remains available for workflows on the raw concentration scale (e.g., certain flow cytometry or raw-scale ELISA applications).

4.3 Four-Parameter Gompertz

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

The Gompertz function is intrinsically asymmetric without requiring a fifth parameter. The inflection response is \(a + (d-a)e^{-1} \approx a + 0.632(d-a)\), which sits above the midpoint of the response range. Mechanistically, the Gompertz model describes a signal that rises steeply at low concentrations and then levels off gradually — a pattern common in assays with strong analyte binding at low occupancy and steric hindrance at high occupancy.

The analytical inverse is \[ x = c - \frac{1}{b}\ln\!\bigl(-\ln\tfrac{y-a}{d-a}\bigr), \] and both \(\tfrac{y-a}{d-a} \in (0, 1)\) and the resulting argument of the outer logarithm must be positive, restricting \(y\) to \((a, d)\) just as for the 4PL.

Strengths: Captures inherent asymmetry without the identifiability problems of 5-parameter models; often the correct model form for Luminex assays at high signal levels.

Limitations: Numerically sensitive due to the double exponential (nested exp); the asymmetry direction is fixed (always skewed toward \(d\), never toward \(a\)).

4.4 Five-Parameter Logistic (5PL)

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

The 5PL extends the 4PL with an asymmetry parameter \(g > 0\). When \(g = 1\) the model reduces exactly to the 4PL. The inflection point is displaced from \(c\) by \(b \ln g\): \[ x_{\text{infl}} = c + b \ln g, \] and the inflection response is \(a + (d-a)\bigl(\tfrac{g}{g+1}\bigr)^g\). Values \(g > 1\) shift the inflection toward higher concentrations and compress the lower part of the curve; \(g < 1\) has the opposite effect.

The analytical inverse is \[ x = c - b \ln\!\Bigl[\Bigl(\tfrac{d-a}{y-a}\Bigr)^{1/g} - 1\Bigr], \] valid for \(y \in (a, d)\).

Strengths: Most general of the logistic family; can model a wide range of asymmetric sigmoidal shapes; widely implemented in regulatory guidance as an extension of 4PL.

Limitations: Introduces strong correlation between \(b\), \(c\), and \(g\); the asymmetry parameter is often unidentifiable from typical 10-point standard curves, leading to boundary estimates (\(g\) at its lower constraint) and inflated covariance matrices. This is the primary failure mode that the eligibility gating system in curveR is designed to detect and exclude.

4.5 Five-Parameter Generalised Logistic (Richards / log-logistic 5PL)

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

This is the Richards generalisation of the logistic, also known as the log-logistic 5PL. It uses a different asymmetry mechanism from the 5PL: rather than raising the sigmoid to a power \(g\), it modulates the steepness of approach to each asymptote differently. The inflection point is at \(x_{\text{infl}} = c + \tfrac{\ln g}{b}\), and the model reduces to the 4PL in the limit \(g \to 1\).

The analytical inverse is \[ x = c - \frac{1}{b}\left[\ln\!\Bigl(\Bigl(\tfrac{y-a}{d-a}\Bigr)^{-g} - 1\Bigr) - \ln g\right], \] valid for \(y \in (a, d)\).

Strengths: Can approximate logistic, Gompertz, and other sigmoidal shapes; the most flexible model in the ensemble.

Limitations: Highest risk of overfitting and identifiability problems; \(g\) near zero leads to numerical instability in the outer logarithm.

4.6 Model Summary Table

On the \(\log_{10}\)-transformed concentration scale, four model families are fitted in the standard workflow (log-logistic 4PL is dropped when is_log_independent = TRUE):

5.1 Statistical Model

Let \(y_{ij}\) denote the (log\(_{10}\)-transformed) response for standard point \(j = 1, \ldots, n_i\) on plate \(i\), and let \(x_{ij}\) denote the corresponding (log\(_{10}\)-transformed) concentration. The frequentist model is

\[ y_{ij} = f(x_{ij};\, \boldsymbol{\theta}_i) + \varepsilon_{ij}, \quad \varepsilon_{ij} \overset{\text{iid}}{\sim} \mathcal{N}(0,\, \sigma_i^2), \]

where \(f\) is one of the five sigmoidal models above and \(\boldsymbol{\theta}_i\) is the plate-specific parameter vector. Plates are fitted independently; no cross-plate pooling is performed in the frequentist setting. The residual variance \(\sigma_i^2\) is estimated from the data as the mean squared residual \(\hat{\sigma}_i = \sqrt{\text{RSS}/(n_i - p)}\) where \(p\) is the number of free parameters.

5.2 Parameter Constraints and Bounds

Before fitting, curveRcore::adaptive_constraint_profile() inspects the observed response range and classifies the assay signal into one of three scale classes: high (log-max MFI \(> 2.5\) or raw max \(> 1000\), typical of Luminex), medium, or low (OD/absorbance-like, narrow dynamic range). From this classification, curveRfreq::compute_model_constraints() builds data-adaptive parameter bounds for all five free parameters:

where \(\Delta_y = y_{\max} - y_{\min}\) is the observed response range, \(\Delta_x\) is the observed concentration range, \(\bar{x}\) is the midpoint, \(m\) is a margin fraction (0.5 for high-signal assays), and \(p\) is a concentration padding fraction. When fixed_a is non-NULL, the lower asymptote is treated as a known constant and excluded from the optimisation.

5.3 Multi-Start Levenberg–Marquardt Optimisation

For each model family, fit_ensemble_nls() fits via multi-start Levenberg–Marquardt NLS using minpack.lm::nlsLM(). Starting values are generated as a deterministic grid: \(n_{\text{starts}}\) evenly-spaced points across the parameter bounds (3 starts for high/medium signal; 5 for low signal, with \(b\) biased small to avoid steep-slope initialisation failures).

For each start, nlsLM() is called with convergence tolerances of \(10^{-6}\) (parameter and function relative tolerances) and a maximum of 100 iterations (200 for low-signal data). The fit with the lowest AIC among all converged starts is retained.

If all nlsLM starts fail for low-signal data, two fallbacks are attempted: (1) nlsLM with bounds relaxed by 50%, and (2) base-R nls() with the "port" (box-constrained) algorithm and relaxed tolerances of \(10^{-10}\).

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: best-eligible model parameters by curve"
)

6.1 Statistical Model

curveRbayes fits all curve_id values simultaneously in a single hierarchical Stan model per model family. Let \(y_{ij}\) be the (log\(_{10}\)-transformed) response for standard point \(j\) on plate \(i\), with \(x_{ij}\) the corresponding log\(_{10}\) concentration. The likelihood is

\[ y_{ij} \sim \text{Student-}t\!\bigl(\nu,\; f(x_{ij};\, \boldsymbol{\theta}_i),\; \sigma_{\text{obs}}\bigr), \]

where \(f\) is the chosen forward model. The Student-\(t\) likelihood with estimated degrees of freedom \(\nu\) provides automatic robustness to outlying standard points compared with a Gaussian assumption (Y. Fong et al. 2012). Both \(\nu\) and \(\sigma_{\text{obs}}\) are shared across all plates.

6.2 Hierarchical Structure (Partial Pooling)

Plate-level parameters are drawn from plate-population distributions using a non-centred parameterisation to improve HMC geometry and reduce divergent transitions:

\[ \begin{aligned} \mu_a,\; \sigma_a &\quad\text{(population lower asymptote)}\\ a_i &= \mu_a + \sigma_a \cdot \tilde{a}_i, \quad \tilde{a}_i \sim \mathcal{N}(0, 1) \end{aligned} \]

and analogously for \(d_i\), \(\log b_i\), \(c_i\), and (for 5-parameter models) \(\log g_i\). Taking \(b_i = \exp(\log b_i)\) and \(g_i = \exp(\log g_i)\) enforces \(b_i > 0\) and \(g_i > 0\) automatically.

This partial-pooling structure regularises plate-level estimates toward the population mean, borrowing strength across replicate plates while still allowing plate-specific curves.

6.3 Data-Adaptive Priors

All prior hyperparameters are computed from the preprocessed data by compute_dynamic_priors() and passed to Stan as data, making the models scale-invariant:

\[ \begin{aligned} \mu_a &\sim \mathcal{N}(y_{\min},\; 0.3\Delta_y) \qquad \sigma_a \sim \mathcal{N}^+(0,\; 0.15\Delta_y)\\ \mu_d &\sim \mathcal{N}(y_{\max} + 0.1\Delta_y,\; 0.3\Delta_y) \qquad \sigma_d \sim \mathcal{N}^+(0,\; 0.15\Delta_y)\\ \mu_{\log b} &\sim \mathcal{N}(0,\; 0.7) \qquad \sigma_{\log b} \sim \mathcal{N}^+(0,\; 0.5)\\ \mu_c &\sim \mathcal{N}(\bar{x},\; 0.5\Delta_x) \qquad \sigma_c \sim \mathcal{N}^+(0,\; 0.25\Delta_x)\\ \sigma_{\text{obs}} &\sim \mathcal{N}^+(0,\; 0.3\Delta_y) \qquad \nu \sim \text{Gamma}(2,\; 0.1) \end{aligned} \]

For 5-parameter models, the asymmetry parameter \(g\) receives a log-scale prior centred at zero (i.e., \(g = 1\), the 4PL), regularising toward the simpler model: \[ \mu_{\log g} \sim \mathcal{N}(0,\; 0.5), \qquad \sigma_{\log g} \sim \mathcal{N}^+(0,\; 0.3). \]

This is a substantive prior choice: it encodes the belief that the 4PL is a plausible description a priori, and that asymmetry must be strongly evidenced in the data to move the posterior for \(g\) away from 1. This regularisation is a key reason why the Bayesian framework may validly select a 5-parameter model in cases where the frequentist framework correctly rejects it due to boundary estimation.

When fixed_a is supplied, the prior on \(\mu_a\) is tightened to a soft constraint centred on the fixed value: \(\mu_a \sim \mathcal{N}(\texttt{fixed\_a},\; \max(0.01\Delta_y,\; 10^{-4}))\), effectively pinning the lower asymptote while still allowing small numerical deviations.

6.4 MCMC Sampling

HMC sampling is performed via cmdstanr with default settings of 4 chains, 1000 warm-up iterations, and 1000 sampling iterations, using the NUTS sampler with adapt_delta = 0.9 and max_treedepth = 12. Convergence is monitored via E-BFMI and the number of divergent transitions stored in the result metadata.

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,
  run_loo            = TRUE,
  verbose            = FALSE
)

bayes_tbl <- summary_table_bayes(bayes_result)
knitr::kable(bayes_tbl, digits = 3,
             caption = "Bayesian: posterior means and SDs by curve")

7.1 The Limitation of AIC and LOO Alone

The purpose of a calibration curve is back-calculation: given an observed instrument response \(y^*\), invert the fitted model to recover the analyte concentration \(\hat{x} = f^{-1}(y^*;\, \hat{\boldsymbol{\theta}})\). AIC and LOO expected log predictive density (ELPD) measure forward-fit quality — how well the model predicts the observed standard responses — but are entirely agnostic to whether the fitted model can reliably invert.

A 5-parameter model may achieve a lower AIC than a 4-parameter model by capturing a small residual structure with its extra flexibility parameter. But if that parameter is boundary-constrained (e.g., \(\hat{g} = g_{\min}\) from nlsLM), the resulting NLS covariance matrix is singular or near-singular, the delta-method standard errors for back-calculated concentrations are inflated for every point on the grid, and the precision profile (pcov) is pinned at its ceiling of 150% throughout the entire working range. The model wins AIC selection but is useless for quantification.

AIC and LOO never see the precision profile. They cannot detect this failure mode.

7.2 The Four-Step Eligibility-Gated Selection Pipeline

Both curveRfreq and curveRbayes implement an identical four-step pipeline defined in curveRcore::eligibility.R:

# Frequentist AIC weights (always computed, regardless of gating)
aic_rows <- lapply(names(freq_result$plates), function(cid) {
  cr <- freq_result$plates[[cid]]
  if (is.null(cr) || is.null(cr$selection$weights)) return(NULL)
  w          <- cr$selection$weights
  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, delta-AIC, and Akaike weights for all converged models"
)

# Bayesian LOO-CV comparison and selected model
bayes_c1 <- bayes_result$plates[["1"]]
cat("Best Bayesian model (curve 1):", bayes_c1$selection$best_model_name, "\n")
#> Best Bayesian model (curve 1): logistic5
cat("Selection criterion          :", bayes_c1$selection$criterion, "\n")
#> Selection criterion          : LOO+eligibility
cat("Fallback triggered           :", bayes_c1$selection$fallback, "\n")
#> Fallback triggered           : FALSE
if (!is.null(bayes_c1$selection$loo_comparison))
  print(bayes_c1$selection$loo_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

# Side-by-side selection comparison across all curves
sel_rows <- list()
for (cid in names(freq_result$plates)) {
  freq_cr  <- freq_result$plates[[cid]]
  bayes_cr <- bayes_result$plates[[cid]]
  if (is.null(freq_cr) || is.null(bayes_cr)) next
  sel_rows[[length(sel_rows) + 1]] <- data.frame(
    curve_id        = cid,
    freq_best       = freq_cr$selection$best_model_name,
    freq_criterion  = freq_cr$selection$criterion,
    freq_aic_winner = freq_cr$selection$aic_best,
    bayes_best      = bayes_cr$selection$best_model_name,
    bayes_criterion = bayes_cr$selection$criterion,
    agree           = freq_cr$selection$best_model_name ==
                      bayes_cr$selection$best_model_name,
    stringsAsFactors = FALSE)
}
sel_df <- do.call(rbind, sel_rows)
knitr::kable(sel_df, caption = "Model selection agreement across methods and curves")

8.1 Definition

The prediction coefficient of variation (pcov) quantifies the relative uncertainty of a back-calculated concentration estimate at a given true concentration \(x_0\). It is the key output used to define the usable dynamic range of the assay. Both methods report pcov in percent CV units, where 20% is the conventional threshold for defining LLOQ and ULOQ (Ekins 1983).

8.2 Frequentist: The Delta Method

The frequentist pcov is derived from the delta method applied to the inverse calibration function \(\hat{x} = f^{-1}(y;\, \hat{\boldsymbol{\theta}})\). Let \(\Sigma\) be the estimated parameter covariance matrix (vcov(fit)), \(\hat{\sigma}\) be the residual standard deviation (summary(fit)$sigma), and define the two gradient vectors: \[ \nabla_{\theta} f^{-1}(y;\, \hat{\boldsymbol{\theta}}) = \frac{\partial f^{-1}}{\partial \boldsymbol{\theta}} \bigg|_{y = f(x_0;\, \hat{\boldsymbol{\theta}})} \qquad \frac{\partial f^{-1}}{\partial y} \bigg|_{y = f(x_0;\, \hat{\boldsymbol{\theta}})}. \]

The variance of the back-calculated concentration is (O’Connell, Belanger, and Harland 1993, eq. 3.3):

\[ \widehat{\text{Var}}(\hat{x}) = \underbrace{ \nabla_{\theta}^{\top} \Sigma\, \nabla_{\theta} }_{\text{parameter uncertainty}} + \underbrace{ \left(\frac{\partial f^{-1}}{\partial y}\right)^{\!2} \hat{\sigma}^2 }_{\text{observation noise}}. \]

The two terms represent, respectively, uncertainty propagated from estimation of \(\boldsymbol{\theta}\) (parameter covariance) and uncertainty from measurement noise in the response (observation noise, scaled by the Jacobian of the inverse). Including the observation noise term — via se_response = summary(fit)$sigma in predict_grid_freq() — is essential for aligning the frequentist precision profile with the Bayesian CDAN profile, which by construction includes both sources.

When is_log_x = TRUE (concentration on log\(_{10}\) scale), the pcov in percent CV units is: \[ \text{pcov}(x_0) = 100 \times \ln(10) \times \widehat{\text{SE}}(\hat{x}), \] reflecting the fact that a standard deviation of \(\delta\) on the log\(_{10}\) scale corresponds to approximately \(100\delta\ln(10)\%\) relative uncertainty on the original scale.

All analytical gradients are implemented in curveRcore::gradients.R for all five model families, in both the free-\(a\) and fixed-\(a\) variants. The factory function make_inv_and_grad_fixed(model, fixed_a) is constructed once per model and called inside the grid loop, avoiding repeated closure construction and achieving approximately 200× speedup over the prior per-grid-point implementation.

In addition to pcov, curveRfreq reports pcov_rmse — the relative root mean squared error, which incorporates the squared bias between the back-calculated estimate and the true grid concentration: \[ \text{pcov\_rmse}(x_0) = 100 \times \ln(10) \times \sqrt{\widehat{\text{Var}}(\hat{x}) + (\hat{x} - x_0)^2}. \]

Both pcov and pcov_rmse are capped at cv_x_max = 150%.

8.3 Bayesian: The CDAN Approach

The Bayesian precision profile uses the Concentration Distribution of Assay Noise (CDAN) method introduced by O’Malley (2008). For each of \(S\) posterior draws \(s = 1, \ldots, S\), the procedure is:

1. Forward prediction at the true grid concentration \(x_0\): \[ \mu_s = f(x_0;\, \boldsymbol{\theta}_s) \]

2. Inject observation noise from the posterior predictive distribution: \[ y_s^* = \mu_s + \sigma_{\text{obs},s} \cdot t_s, \quad t_s \sim \text{Student-}t(\nu_s) \]

3. Back-calculate concentration from the noisy response: \[ x_s^* = f^{-1}(y_s^*;\, \boldsymbol{\theta}_s) \]

The pcov is then the posterior standard deviation of the back-calculated concentrations, converted to percent CV: \[ \text{pcov}(x_0) = 100 \times \ln(10) \times \text{SD}(x_1^*, \ldots, x_S^*). \]

The pcov_rmse in the Bayesian setting is the relative RMSE of the back-calculated draws about the true concentration: \[ \text{pcov\_rmse}(x_0) = 100 \times \ln(10) \times \sqrt{\frac{1}{S}\sum_{s=1}^S (x_s^* - x_0)^2}. \]

Key distinction for test samples: When back-calculating concentrations for test samples (as opposed to the standard grid), the observed response \(y^*\) is already a noisy measurement. Therefore, predict_samples_bayes() does not inject additional noise: it applies the inverse function \(f^{-1}(y^*;\, \boldsymbol{\theta}_s)\) directly to the observed response at each draw. The result is a posterior distribution of plausible concentrations consistent with that observed response under all credible curve shapes.

# --- Model colour / line-type conventions (used throughout) ---
model_cols <- c(logistic4   = "steelblue",   logistic5  = "darkorange",
                gompertz4   = "forestgreen", loglogistic5 = "purple")
model_ltys <- c(logistic4   = 1, logistic5 = 2,
                gompertz4   = 3, loglogistic5 = 4)
model_pchs <- c(logistic4   = 16, logistic5 = 17,
                gompertz4   = 15, loglogistic5 = 18)

# --- LOQ from a pcov profile (linear interpolation at threshold crossings) ---
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))
  transitions <- diff(below)
  enter <- which(transitions ==  1)
  leave <- which(transitions == -1)
  lloq <- if (below[1]) x[1]
          else { idx <- enter[1];
            x[idx] + (threshold - pcov[idx]) /
              (pcov[idx+1] - pcov[idx]) * (x[idx+1] - x[idx]) }
  uloq <- if (below[length(below)]) x[length(x)]
          else { idx <- leave[length(leave)];
            x[idx] + (threshold - pcov[idx]) /
              (pcov[idx+1] - pcov[idx]) * (x[idx+1] - x[idx]) }
  list(lloq = lloq, uloq = uloq,
       dynamic_range_log10 = uloq - lloq,
       dynamic_range_fold  = 10^(uloq - lloq))
}

# --- Extract named parameter vector from an ensemble entry ---
get_params_freq  <- function(ens)
  stats::setNames(ens$parameters$estimate, ens$parameters$term)
get_params_bayes <- function(ens)
  stats::setNames(ens$parameters$mean, ens$parameters$term)

# --- Inflection point location and slope at inflection for any model ---
compute_sensitivity <- function(model_name, params) {
  a <- params["a"]; b <- params["b"]
  c <- params["c"]; d <- params["d"]; g <- params["g"]
  x_infl <- switch(model_name,
    logistic4    = c,
    logistic5    = c + b * log(g),
    gompertz4    = c,
    loglogistic4 = c,
    loglogistic5 = c + log(g) / b,
    c)
  slope <- switch(model_name,
    logistic4    = dydx_logistic4(x_infl, a, b, c, d),
    logistic5    = dydx_logistic5(x_infl, a, b, c, d, g),
    gompertz4    = dydx_gompertz4(x_infl, a, b, c, d),
    loglogistic4 = dydx_loglogistic4(x_infl, a, b, c, d),
    loglogistic5 = dydx_loglogistic5(x_infl, a, b, c, d, g),
    NA_real_)
  data.frame(inflection_x       = unname(x_infl),
             slope_at_inflection = unname(slope),
             stringsAsFactors   = FALSE)
}

8.4 Calibration Curve and Confidence Band (Best Model)

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

par(mar = c(4, 4, 2.5, 1))
plot(grid_comp$log10_concentration,
     grid_comp$frequentist_predicted_response,
     type = "l", col = "steelblue", lwd = 2.5,
     xlab = expression(log[10](concentration ~ (AU/mL))),
     ylab = expression(log[10](MFI)),
     main = "Calibration Curves: curve_id = 1 (best eligible model)")
lines(grid_comp$log10_concentration,
      grid_comp$bayesian_predicted_response,
      col = "firebrick", lwd = 2.5, 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.9)
legend("bottomright",
       legend = c(
         paste0("Frequentist (", freq_c1$selection$best_model_name, ")"),
         paste0("Bayesian (", bayes_c1$selection$best_model_name, ")"),
         "Standards"),
       col  = c("steelblue", "firebrick", "black"),
       lty  = c(1, 2, NA), pch = c(NA, NA, 16),
       lwd  = c(2.5, 2.5, NA), cex = 0.8, bg = "white")

Calibration curves for curve_id = 1 from both methods (best eligible model). Shaded bands are 95% CIs. Standard points shown as filled circles.

8.5 Precision Profiles: Best Model, Both Methods

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

# Panel 1: pcov (CV)
plot(grid_comp$log10_concentration, grid_comp$frequentist_pcov,
     type = "l", col = "steelblue", lwd = 2.5,
     xlab = expression(log[10](concentration ~ (AU/mL))),
     ylab = "pcov  (% CV)",  ylim = c(0, 100),
     main = "Precision Profile: CV")
lines(grid_comp$log10_concentration, grid_comp$bayesian_pcov,
      col = "firebrick", lwd = 2.5, lty = 2)
abline(h = pcov_threshold, col = "grey40", lty = 3, lwd = 1.5)
text(min(grid_comp$log10_concentration) + 0.4, pcov_threshold + 4,
     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.85, bg = "white")

# Panel 2: pcov_rmse (RRMSE)
plot(grid_comp$log10_concentration, grid_comp$frequentist_pcov_rmse,
     type = "l", col = "steelblue", lwd = 2.5,
     xlab = expression(log[10](concentration ~ (AU/mL))),
     ylab = "pcov_rmse  (% RRMSE)", ylim = c(0, 100),
     main = "Precision Profile: RRMSE")
lines(grid_comp$log10_concentration, grid_comp$bayesian_pcov_rmse,
      col = "firebrick", lwd = 2.5, lty = 2)
abline(h = pcov_threshold, col = "grey40", lty = 3, lwd = 1.5)
legend("top", legend = c("Frequentist", "Bayesian"),
       col = c("steelblue", "firebrick"),
       lty = c(1, 2), lwd = 2, cex = 0.85, bg = "white")

Precision profiles (pcov) for the best eligible model under each method. The horizontal dashed line marks the 20% CV threshold used to define LLOQ and ULOQ.

8.6 Precision Profiles by Model Form

Because grids are now computed for every converged model during eligibility assessment, per-model precision profiles are directly available from ensemble[[model_name]]$grid — no recomputation is required.

freq_c1 <- freq_result$plates[["1"]]

# Collect per-model grids directly from the ensemble
freq_model_grids <- list()
for (nm in names(freq_c1$ensemble)) {
  ens <- freq_c1$ensemble[[nm]]
  if (!isTRUE(ens$converged) || is.null(ens$grid)) next
  freq_model_grids[[nm]] <- ens$grid
}

par(mfrow = c(1, 1), mar = c(4, 4, 3, 1))
first <- TRUE
for (nm in names(freq_model_grids)) {
  g <- freq_model_grids[[nm]]
  if (first) {
    plot(g$log10_concentration, g$pcov,
         type = "l", col = model_cols[nm], lwd = 2, lty = model_ltys[nm],
         xlab = expression(log[10](concentration ~ (AU/mL))),
         ylab = "pcov (% CV)", ylim = c(0, 100),
         main = "Frequentist Precision by Model (curve_id = 1)")
    first <- FALSE
  } else {
    lines(g$log10_concentration, g$pcov,
          col = model_cols[nm], lwd = 2, lty = model_ltys[nm])
  }
}
abline(h = pcov_threshold, col = "grey40", lty = 3, lwd = 1.5)
legend("top", legend = names(freq_model_grids),
       col = model_cols[names(freq_model_grids)],
       lty = model_ltys[names(freq_model_grids)],
       lwd = 2, cex = 0.85, bg = "white")

Frequentist pcov profiles for all converged models (curve_id = 1). Profiles are taken directly from ensemble[[model]]$grid, computed during eligibility assessment.

par(mfrow = c(1, 1), mar = c(4, 4, 3, 1))
first <- TRUE
for (nm in names(bayes_c1$ensemble)) {
  ens <- bayes_c1$ensemble[[nm]]
  if (!isTRUE(ens$converged) || is.null(ens$grid)) next
  g <- ens$grid
  if (first) {
    plot(g$log10_concentration, g$pcov,
         type = "l", col = model_cols[nm], lwd = 2, lty = model_ltys[nm],
         xlab = expression(log[10](concentration ~ (AU/mL))),
         ylab = "pcov (% CV)", ylim = c(0, 100),
         main = "Bayesian CDAN Precision by Model (curve_id = 1)")
    first <- FALSE
  } else {
    lines(g$log10_concentration, g$pcov,
          col = model_cols[nm], lwd = 2, lty = model_ltys[nm])
  }
}
abline(h = pcov_threshold, col = "grey40", lty = 3, lwd = 1.5)
legend("top", legend = names(bayes_c1$ensemble),
       col = model_cols[names(bayes_c1$ensemble)],
       lty = model_ltys[names(bayes_c1$ensemble)],
       lwd = 2, cex = 0.85, bg = "white")

Bayesian CDAN pcov profiles for all converged models (curve_id = 1), taken directly from ensemble[[model]]$grid.

common_models <- intersect(names(freq_model_grids), names(bayes_c1$ensemble))
common_models <- common_models[vapply(common_models, function(nm)
  !is.null(bayes_c1$ensemble[[nm]]$grid), logical(1))]

n_mod <- length(common_models)
if (n_mod > 0) {
  par(mfrow = c(ceiling(n_mod / 2), 2), mar = c(4, 4, 3, 1))
  for (nm in common_models) {
    fg <- freq_model_grids[[nm]]
    bg <- bayes_c1$ensemble[[nm]]$grid
    plot(fg$log10_concentration, fg$pcov,
         type = "l", col = "steelblue", lwd = 2,
         xlab = expression(log[10](concentration)),
         ylab = "pcov (% CV)", ylim = c(0, 100),
         main = paste0(nm, ": Frequentist vs Bayesian"))
    lines(bg$log10_concentration, bg$pcov,
          col = "firebrick", lwd = 2, lty = 2)
    abline(h = pcov_threshold, col = "grey40", lty = 3, lwd = 1.5)
    legend("top", legend = c("Frequentist", "Bayesian"),
           col = c("steelblue", "firebrick"),
           lty = c(1, 2), lwd = 2, cex = 0.75, bg = "white")
  }
}

Side-by-side comparison of frequentist (delta method, blue) and Bayesian (CDAN, red) pcov profiles for each model family. Agreement between methods varies by model form.

10.1 Frequentist Sample Prediction

For test samples, predict_samples_freq() applies the inverse function and delta method to the observed sample response \(y^*\). The response is first log\(_{10}\)-transformed if is_log_response = TRUE. The back-calculated log\(_{10}\)-concentration is: \[ \hat{x}^* = f^{-1}(y^*;\, \hat{\boldsymbol{\theta}}) \]

with variance \[ \widehat{\text{Var}}(\hat{x}^*) = \nabla_{\theta}^{\top} \Sigma\, \nabla_{\theta} + \left(\frac{\partial f^{-1}}{\partial y}\right)^{\!2} \hat{\sigma}^2. \]

The final reported concentration applies the sample dilution: \(\hat{C} = 10^{\hat{x}^*} \times \text{dilution}\). For test samples, pcov_rmse = pcov because the true concentration is unknown; the bias term is undefined.

10.2 Bayesian Sample Prediction

For test samples, predict_samples_bayes() passes the observed log-response \(y^*\) directly through the inverse at each posterior draw — without injecting additional noise, since the measurement noise is already embodied in \(y^*\): \[ x_s^* = f^{-1}(y^*;\, \boldsymbol{\theta}_s), \quad s = 1, \ldots, S. \]

The reported concentration is the posterior median \(10^{\tilde{x}^*} \times \text{dilution}\), and the pcov is \(100 \times \ln(10) \times \text{SD}(x_1^*, \ldots, x_S^*)\).

sample_comp <- compare_samples(freq_c1, bayes_c1,
                               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 (log10)\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 = "Bayesian − Frequentist",
       main = sprintf("Bland–Altman\nBias = %.4f log10 units",
                      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
  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 ~ AU/mL)),
         ylab = expression(Bayes ~ log[10](final ~ conc ~ AU/mL)),
         main = "Final Concentration (dilution-adjusted)")
    abline(0, 1, col = "red", lty = 2, lwd = 2)
  }
}

Comparison of frequentist and Bayesian sample concentration predictions (curve_id = 1). Top row: scatter and Bland-Altman plots on the log10 scale. Bottom row: sample pcov agreement and final concentration comparison.

14.1 Model Forms and Their Suitability

The four-parameter logistic remains the most widely deployed calibration model for immunoassays, and in the simulated data presented here it performs well, producing well-conditioned parameter estimates and meaningful dynamic ranges under both methods. The Gompertz model is appropriate when the data exhibit inherent asymmetry — a pattern that arises naturally in bead-based assays where the signal rises steeply from background but saturates gradually — and offers this flexibility without the identifiability burden of a fifth parameter.

The five-parameter models (5PL and Richards logistic) offer the greatest flexibility but carry a corresponding cost. On a 10-point standard curve, the asymmetry parameter \(g\) is frequently non-identifiable: the optimiser converges with \(g\) at its lower bound, the covariance matrix becomes ill-conditioned, and the delta-method precision profile is pinned at its ceiling for the entire working range. The eligibility gating system in curveR is specifically designed to detect and exclude this failure mode before model selection proceeds.

14.2 Frequentist vs. Bayesian Precision Profiles

The frequentist delta-method and Bayesian CDAN precision profiles are now on closely comparable scales (Pearson correlation \(r \approx 0.96\)) following implementation of the two-term delta-method formula that includes both parameter covariance and residual observation noise. Prior to this, single-term implementations that omitted the observation noise term systematically underestimated pcov relative to CDAN.

The Bayesian precision estimates are typically slightly more conservative (higher pcov) than the frequentist estimates, consistent with the argument of O’Malley (2008) that the CDAN profile captures the full posterior predictive uncertainty. The frequentist delta method is a first-order Taylor approximation and will underestimate uncertainty in regions of high curvature of the inverse function.

14.3 The Role of Hierarchical Regularisation

The Bayesian framework’s partial pooling across replicate plates provides meaningful regularisation of the asymmetry parameter \(g\) through the log-scale prior centred at \(g = 1\). This regularisation is responsible for the key empirical observation noted in the architecture: the Bayesian framework may validly select a 5-parameter model on the same dataset where the frequentist framework correctly rejects it. The two frameworks are applying the same eligibility gates to precision profiles estimated by different (and differently regularised) methods.

Here is the continuation of the Discussion section, written at Biometrics manuscript quality and grounded precisely in the source code.

14.4 Identifiability, Boundary Estimates, and the Eligibility Gate System in Detail

14.5 Precision Profile Agreement Between Methods

14.6 Limits of Quantification and Dynamic Range in Practice

14.7 Computational Considerations

14.8 Relationship to Existing Software

The curveR ecosystem occupies a distinct niche relative to existing immunoassay calibration software.

nCal (Youyi Fong et al. 2013) implements a Bayesian hierarchical 5PL model for Luminex data and was one of the first open-source tools to provide posterior uncertainty quantification for immunoassay calibration. curveR extends this in three directions: it supports five model families rather than one; it implements the full CDAN precision profile rather than reporting concentration point estimates alone; and it provides a symmetric frequentist alternative on an identical interface.

drc (Ritz et al., 2015) provides a comprehensive frequentist dose-response framework covering many model families. It does not compute precision profiles or apply quantification-aware model selection.

The lm / nls workflow used in many laboratory pipelines fits a single user-chosen model without model selection, eligibility gating, or precision profiling. curveR automates and formalises this workflow while adding all three.

The key architectural innovations of curveR are: (1) the symmetric calibration_result output class enabling direct frequentist–Bayesian comparison; (2) the eligibility-gated model selection pipeline that conditions on quantification performance rather than forward-fit quality; and (3) the standardised precision profile (pcov and pcov_rmse) computed for every converged model in the ensemble, making model-form comparison of dynamic range a built-in feature rather than an afterthought.

14.9 Limitations and Future Directions

Replication structure. The current implementation treats each plate as an independent curve (frequentist) or as one level of a two-level hierarchy (Bayesian: plate within population). When multiple replicates of the same standard are measured within a plate — common in clinical laboratory settings — the within-plate replication structure is not explicitly modelled; replicates contribute to the residual variance \(\hat{\sigma}\) or \(\sigma_{\text{obs}}\) as ordinary residuals. A three-level hierarchy (replicate within standard concentration within plate) would be straightforward to add to the Stan models and would provide better-calibrated uncertainty estimates for high-replication designs.

Heteroscedastic noise. The current models assume homoscedastic residuals on the log-transformed response scale, which is approximately valid when the noise coefficient of variation is roughly constant in MFI units. Some assay platforms exhibit noise that is better described as a power-law function of the mean (\(\sigma \propto \mu^\alpha\)) or as an additive-multiplicative combination. Belanger, Davidian, and Giltinan (1996) showed that variance function misspecification can substantially inflate calibration intervals; extending curveR to support variance modelling (e.g., via gnls() in the frequentist setting or a heteroscedastic Student-t likelihood in Stan) is a natural next step.

Blank integration in the Bayesian model. The Stan models support blank data via an optional second likelihood: \(y_{\text{blank}} \sim \text{Student-}t(\nu, a_i, \sigma_{\text{blank}})\), where \(\sigma_{\text{blank}}\) is a separate noise parameter for blank wells. This anchors the lower asymptote \(a_i\) more precisely when blank data are available and is passed through build_stan_data(). The current vignette does not demonstrate this feature explicitly; a dedicated section showing the effect of blank integration on the precision profile near the LLOQ would be a valuable addition to a data analysis companion vignette.

Prozone correction robustness. The current hook-effect correction in correct_prozone() uses a simple heuristic: detect the peak response and apply a dampening factor to post-peak observations. This is appropriate when the hook effect is mild (a common pattern in Luminex assays at the highest standard concentrations) but may be insufficient for severe hook effects. A principled statistical treatment — for example, fitting a bi-phasic model or using a monotone spline with the hook region downweighted — would improve robustness for assays with strong antigen excess.

Ensemble model averaging. The current selection pipeline identifies a single best eligible model. Bayesian model averaging (BMA) or predictive stacking (Yao et al. 2018) over the eligible model ensemble would produce concentration estimates and precision profiles that are appropriately averaged over model-form uncertainty. The LOO stacking weights are already computed and stored in selection$loo_weights for the Bayesian pathway; using these weights to produce a stacked precision profile is a natural and methodologically principled extension.

Computational scaling. For very large multiplex panels (hundreds of analytes across many plates), the sequential per-curve-id fitting in curveRfreq and the per-model-family hierarchical fitting in curveRbayes can become a bottleneck. Parallelisation — across curve IDs in the frequentist setting, and across model families in both settings — is straightforward using parallel::mclapply() or future.apply::future_lapply() and is a planned enhancement.

14.10 Conclusion

The curveR ecosystem provides a complete, symmetric, and statistically rigorous framework for immunoassay calibration curve fitting and concentration back-calculation. The key methodological contributions are:

1. A shared model library covering four non-redundant sigmoidal families on the log-concentration scale, with analytical inverses, first derivatives, and full analytical gradient vectors for all five models.

2. A two-term delta-method precision profile (parameter covariance plus observation noise) for the frequentist pathway, aligned in scale with the Bayesian CDAN profile, enabling direct comparison of quantification uncertainty between methods.

3. A four-gate eligibility screening system that conditions model selection on quantification performance rather than forward-fit quality, eliminating the systematic failure mode in which 5-parameter models win AIC selection despite having zero usable dynamic range due to boundary- constrained asymmetry parameters.

4. Per-model precision grids, computed during eligibility assessment for every converged model in the ensemble and stored at ensemble[[model]]$grid, enabling transparent comparison of dynamic range across model forms without additional computation.

5. A symmetric output class (calibration_result_multiplate) returned by both fitting packages, making cross-method comparison, agreement quantification, and downstream reporting identical regardless of which inferential framework was used.

Together these features support a principled, reproducible, and fully auditable immunoassay calibration workflow suitable for both exploratory research and regulated analytical contexts.

Here is the complete continuation of the vignette from the end of the Discussion section through to the final closing sections.

---
# Appendix A: The curveRcore Mathematical Toolkit

This appendix documents the analytical functions in `curveRcore`
that underpin both fitting packages, providing the precise
mathematical expressions implemented in the source code.
These are intended as a reference for users who wish to extend
the ecosystem or verify the implemented formulas against an
independent source.

## A.1 Analytical Inverse Functions

All five models have closed-form analytical inverses, implemented
in `curveRcore::inverses.R`.  Each inverse is valid only for
$y \in (a, d)$; responses outside this interval return `NA`.

**4PL inverse:**
$$
x = c + b \ln\!\frac{y - a}{d - y}
$$

**Log-logistic 4PL inverse** (requires $x > 0$; redundant with
4PL on log scale):
$$
x = \frac{c}{\left(\dfrac{d - y}{y - a}\right)^{1/b}}
$$

**Gompertz inverse:**
$$
x = c - \frac{1}{b}\ln\!\left(-\ln\frac{y - a}{d - a}\right)
$$

**5PL inverse:**
$$
x = c - b\ln\!\left[\left(\frac{d - a}{y - a}\right)^{1/g} - 1\right]
$$

**Richards (loglogistic5) inverse:**
$$
x = c - \frac{1}{b}\left[\ln\!\left(\left(\frac{y - a}{d - a}
\right)^{-g} - 1\right) - \ln g\right]
$$

Each function has two variants: a *free-a* form
(`inv_<model>(y, a, b, c, d [, g])`) where $a$ is estimated
from data, and a *fixed-a* form
(`inv_<model>_fixed(y, fixed_a, b, c, d [, g])`) where $a$ is
a pre-specified constant.  The fixed-a forms are used when
`fixed_a` is non-`NULL`, removing $a$ from the gradient
calculation and reducing the parameter covariance matrix
by one row and one column.

## A.2 Analytical First Derivatives (Forward Models)

The first derivative $dy/dx$ is used for assay sensitivity
calculations and for the second-order correction in the
delta method.  All derivatives are implemented in
`curveRcore::derivatives.R`.

**4PL:**
$$
\frac{dy}{dx} = \frac{(d-a)\,u}{b(1+u)^2}, \quad
u = \exp\!\left(-\frac{x-c}{b}\right)
$$

**Gompertz:**
$$
\frac{dy}{dx} = b(d-a)\,v\exp(-v), \quad
v = \exp(-b(x-c))
$$

**5PL:**
$$
\frac{dy}{dx} = \frac{g(d-a)\,u}{b(1+u)^{g+1}}, \quad
u = \exp\!\left(-\frac{x-c}{b}\right)
$$

**Richards (loglogistic5):**
$$
\frac{dy}{dx} = b(d-a)\exp(-b(x-c))\,
\bigl(1 + g\exp(-b(x-c))\bigr)^{-1/g - 1}
$$

Second derivatives $d^2y/dx^2$ are also implemented: analytically
for the 4PL and Gompertz, and via central-difference numerical
approximation ($h = 10^{-5}$) for the three remaining models
where closed-form expressions become unwieldy.

## A.3 Analytical Gradients of the Inverse Functions

The delta method requires two gradient vectors evaluated at
$y_0 = f(x_0;\, \hat{\boldsymbol{\theta}})$:

- $\nabla_{\boldsymbol{\theta}} f^{-1}(y_0;\, \hat{\boldsymbol{\theta}})$:
  the gradient with respect to the parameter vector (used for
  the parameter-uncertainty variance component `var_par`)
- $\partial f^{-1}/\partial y \big|_{y=y_0}$: the scalar
  response-Jacobian (used for the observation-noise variance
  component `var_y`)

These are implemented in `curveRcore::gradients.R`.  We give
the full expressions for the two most commonly selected models:

### 4PL gradients (free $a$)

Let $y_0 \in (a, d)$.  Then:

$$
\frac{\partial x}{\partial a} = -\frac{b}{y_0 - a}, \quad
\frac{\partial x}{\partial b} = \ln\frac{y_0 - a}{d - y_0}, \quad
\frac{\partial x}{\partial c} = 1, \quad
\frac{\partial x}{\partial d} = -\frac{b}{d - y_0}
$$

$$
\frac{\partial x}{\partial y} = \frac{b(d-a)}{(y_0-a)(d-y_0)}
$$

### Gompertz gradients (free $a$)

Let $R = (y_0 - a)/(d - a)$ and $L = -\ln R > 0$.  Then:

$$
\frac{\partial x}{\partial a} = \frac{1}{b L (y_0 - a)}, \quad
\frac{\partial x}{\partial b} = \frac{\ln L}{b^2}, \quad
\frac{\partial x}{\partial c} = 1, \quad
\frac{\partial x}{\partial d} = -\frac{1}{b L (d - a)}
$$

$$
\frac{\partial x}{\partial y} = \frac{1}{b L (y_0 - a)}
$$

Note that $\partial x/\partial a = \partial x/\partial y$ for
the Gompertz, a consequence of the symmetric role of $a$ and
$y_0$ in the expression for $R$.

### Fixed-$a$ variants

When $a$ is treated as a known constant, the $\partial x/\partial a$
term is excluded from `grad_theta`, reducing the gradient vector
by one element.  The response-Jacobian $\partial x/\partial y$
is unchanged.  All five fixed-$a$ gradient functions are
implemented as separate exported functions
(`grad_inv_<model>_fixed()`) and as separate `grad_y_<model>_fixed()`
response-Jacobian functions, enabling reuse in both the grid
and sample prediction loops.

### The factory pattern

The `make_inv_and_grad_fixed(model, fixed_a)` factory constructs
three closures — `$inv(y, p)`, `$grad(y, p)`, `$grad_y(y, p)` —
from the appropriate analytical functions based on the model name
and fixed-$a$ status.  Parameters are extracted by name
(`p[["b"]]`, `p[["g"]]`, etc.) inside every closure, so the
order of elements in the parameter vector `p` does not affect
the result.


``` r
# Demonstrate the factory: build closures for logistic4, a free
fns <- curveRcore::make_inv_and_grad_fixed("logistic4", fixed_a = NULL)

# Hypothetical parameter vector
theta <- c(a = 0.5, b = 0.8, c = 2.0, d = 4.5)

# Evaluate at the inflection response (a+d)/2
y_infl <- (theta["a"] + theta["d"]) / 2

cat("Back-calculated x at inflection response:",
    round(fns$inv(y_infl, theta), 6), "\n")
#> Back-calculated x at inflection response: 2
cat("Should equal c =", theta["c"], "\n\n")
#> Should equal c = 2

cat("Gradient w.r.t. theta:\n")
#> Gradient w.r.t. theta:
print(round(fns$grad(y_infl, theta), 6))
#>  a.a  b.a    c  d.a 
#> -0.4  0.0  1.0 -0.4

cat("\nGradient w.r.t. y (response Jacobian):\n")
#> 
#> Gradient w.r.t. y (response Jacobian):
print(round(fns$grad_y(y_infl, theta), 6))
#>   a 
#> 0.8

15.1 B.1 Top-Level Structure

calibration_result_multiplate
├── $meta          Named list — multiplate metadata
│   ├── method          "frequentist" or "bayesian"
│   ├── package         "curveRfreq" or "curveRbayes"
│   ├── curve_ids       Integer vector of curve IDs fitted
│   ├── n_curves        Integer
│   ├── response_var    Character — response column name
│   ├── is_log_response Logical
│   ├── is_log_independent Logical
│   └── timestamp       POSIXct
│
└── $plates        Named list — one entry per curve_id
    ├── "1" → calibration_result   (see B.2)
    ├── "2" → calibration_result
    └── ...

15.2 B.2 Per-Plate calibration_result Structure

calibration_result          (class: c("calibration_result", "list"))
│
├── $meta          Named list — plate metadata
│   ├── method, package, version, timestamp
│   ├── curve_id            Character
│   ├── response_var, independent_var
│   ├── is_log_response, is_log_independent
│   ├── n_standards, n_samples
│   └── pcov_threshold      Numeric — threshold used for LOQ gating
│
├── $ensemble      Named list — one entry per model attempted
│   └── "logistic4" (example)
│       ├── $model_name     Character
│       ├── $converged      Logical
│       ├── $parameters     Data frame
│       │     Frequentist: term | estimate | std_error
│       │     Bayesian:    term | mean | sd | q2.5 | q50 | q97.5
│       ├── $fit_stats      List
│       │     Frequentist: n_obs, n_params, df_resid, rss, aic, bic, log_lik
│       │     Bayesian:    n_divergent, n_max_treedepth, ebfmi
│       ├── $raw_fit        nlsLM object (freq) or fit_bayes_single() output (Bayes)
│       ├── $grid           Data frame — per-model pcov profile (ALWAYS present
│       │                   for converged models; computed during Step 2 of the
│       │                   eligibility pipeline and stored here directly)
│       │     Columns: log10_concentration | concentration | x_fit |
│       │              predicted_response | ci_lower | ci_upper |
│       │              predicted_concentration | se_concentration |
│       │              pcov | pcov_rmse | pcov_pass
│       └── $eligibility    List — assess_model_eligibility() output
│             ├── $eligible             Logical
│             ├── $gates                Data frame: gate | passed | detail
│             ├── $dynamic_range_log10  Numeric
│             ├── $lloq                 Numeric (log10 scale)
│             └── $uloq                 Numeric (log10 scale)
│
├── $selection     Named list — model selection outcome
│   ├── $best_model_name    Character
│   ├── $criterion          "AIC+eligibility" or "LOO+eligibility"
│   ├── $fallback           Logical
│   ├── $fallback_reason    Character
│   ├── $eligible_models    Character vector
│   ├── $assessments        Named list of assess_model_eligibility() outputs
│   ├── $weights            Data frame — AIC weights (freq); always present
│   ├── $aic_best           Character — AIC winner before gating (freq only)
│   ├── $loo_comparison     loo_compare() matrix (Bayes only)
│   └── $loo_weights        Stacking weights (Bayes only)
│
├── $grid          Data frame — best eligible model grid (backwards compatible)
│   (identical to $ensemble[[best_model_name]]$grid)
│
└── $samples       Data frame or NULL — test sample predictions
      Columns: [original sample columns] |
               raw_assay_response | observed_response_fit |
               predicted_log10_concentration | predicted_concentration |
               final_concentration | se_concentration |
               pcov | pcov_rmse | pcov_pass

15.3 B.3 Accessing Key Results

# ── Best model name and selection criterion ──
cr <- freq_result$plates[["1"]]
cat("Best model   :", cr$selection$best_model_name, "\n")
#> Best model   : logistic4
cat("Criterion    :", cr$selection$criterion, "\n")
#> Criterion    : AIC+eligibility
cat("Fallback     :", cr$selection$fallback, "\n")
#> Fallback     : FALSE

# ── AIC weights table ──
cat("\nAIC weights:\n")
#> 
#> AIC weights:
print(cr$selection$weights)
#>     model_name converged        aic delta_aic       weight
#> 1    gompertz4      TRUE -11.218252 18.975924 6.070419e-05
#> 2    logistic4      TRUE -27.404848  2.789328 1.986522e-01
#> 3    logistic5      TRUE -30.194176  0.000000 8.012871e-01
#> 4 loglogistic5      TRUE   4.416831 34.611007 2.444013e-08

# ── Eligibility gate outcomes for each model ──
cat("\nEligibility gate summary:\n")
#> 
#> Eligibility gate summary:
for (nm in names(cr$selection$assessments)) {
  a      <- cr$selection$assessments[[nm]]
  status <- if (isTRUE(a$eligible)) "PASS" else "FAIL"
  cat(sprintf("  %-14s  %s  (dynamic range = %.2f log10)\n",
              nm, status, a$dynamic_range_log10 %||% NA))
}
#>   logistic4       PASS  (dynamic range = 1.67 log10)
#>   logistic5       FAIL  (dynamic range = 0.00 log10)
#>   gompertz4       FAIL  (dynamic range = 0.00 log10)
#>   loglogistic5    FAIL  (dynamic range = 0.00 log10)

# ── Per-model grids — now always present for all converged models ──
cat("\nPer-model grid availability:\n")
#> 
#> Per-model grid availability:
for (nm in names(cr$ensemble)) {
  ens      <- cr$ensemble[[nm]]
  has_grid <- !is.null(ens$grid)
  n_pts    <- if (has_grid) nrow(ens$grid) else 0L
  cat(sprintf("  %-14s  converged = %-5s  grid rows = %d\n",
              nm, ens$converged, n_pts))
}
#>   logistic4       converged = TRUE   grid rows = 200
#>   logistic5       converged = TRUE   grid rows = 200
#>   gompertz4       converged = TRUE   grid rows = 200
#>   loglogistic5    converged = TRUE   grid rows = 200

# ── Best-model parameters ──
cat("\nBest model parameters:\n")
#> 
#> Best model parameters:
best_ens <- cr$ensemble[[cr$selection$best_model_name]]
print(best_ens$parameters)
#>   term  estimate  std_error
#> 1    a 1.7740750 0.12594608
#> 2    d 4.3391150 0.02610856
#> 3    c 1.9835127 0.06453522
#> 4    b 0.4836332 0.04526522

# ── Sample predictions (first 5 rows) ──
if (!is.null(cr$samples)) {
  cat("\nSample predictions (first 5):\n")
  print(head(cr$samples[, c("sampleid", "predicted_log10_concentration",
                             "final_concentration", "pcov")], 5))
}
#> 
#> Sample predictions (first 5):
#>   sampleid predicted_log10_concentration final_concentration     pcov
#> 1     a001                      3.670156             9358066  78.7173
#> 2     a002                      3.868780            14784624 116.1580
#> 3     a003                      4.040622            21960973 150.0000
#> 4     a004                      3.900223            15894714 123.7074
#> 5     a005                      4.064427            23198354 150.0000

16.1 C.1 The Full Preprocessing Pipeline

curveRcore::preprocess_standards() applies the following transformations in strict sequence. All steps operate on the data frame in place; the function returns both the transformed data frame and an antigen_fit_options record.

# Step 1 — concentration computation
data <- compute_concentration(
  data,
  undiluted_sc_concentration = antigen_settings$standard_curve_concentration,
  independent_variable       = independent_variable,
  is_log_concentration       = is_log_independent
)

# Step 2 — prozone correction (if apply_prozone = TRUE)
data <- correct_prozone(
  stdframe             = data,
  prop_diff            = 0.1,   # dampening factor
  dil_scale            = 2,     # dilution scale
  response_variable    = response_variable,
  independent_variable = independent_variable
)

# Step 3 — blank handling (if blank_option != "ignored")
data <- perform_blank_operation(
  blank_data           = blank_data,
  data                 = data,
  response_variable    = response_variable,
  independent_variable = independent_variable,
  is_log_response      = is_log_response,
  blank_option         = blank_option
)

# Step 4 — log10 response transform (if is_log_response = TRUE)
data <- compute_log_response(
  data              = data,
  response_variable = response_variable,
  is_log_response   = is_log_response,
  floor_method      = "adaptive"   # floors to 1% of min positive value
)

16.2 C.2 Blank Handling Options

The five blank_option values control how blank well measurements are incorporated into the standard curve data before fitting:

The multiplication factors 3× and 10× are conventional in immunoassay practice as conservative estimates of the blank distribution tail, analogous to the 2SD and 3SD rules for limit-of-detection determination (Rodbard 1974).

16.3 C.3 The fixed_a Lower Asymptote Constraint

In some assay configurations — particularly when blank measurements are plentiful and reliable — the lower asymptote \(a\) can be fixed at a pre-determined value rather than estimated from the standard curve data. This reduces the degrees of freedom in the fit and can improve identifiability of the remaining parameters when the lowest standard concentrations lie well above background.

The constraint is specified through the l_asy_constraint_method argument of new_antigen_constraints(). Four methods are supported:

In the frequentist setting, fixed_a is baked directly into the NLS formula as a numeric constant — the parameter does not appear in the optimisation at all. In the Bayesian setting, a tight normal prior centred on fixed_a with scale \(\max(0.01\Delta_y, 10^{-4})\) acts as a soft constraint.

17.1 D.1 Reading the Eligibility Audit Trail

The full gate-level audit trail for every model is stored in selection$assessments. This is the primary tool for diagnosing why a particular model was excluded.

cr <- freq_result$plates[["1"]]

cat("=== Eligibility audit: curve_id = 1 ===\n\n")
#> === Eligibility audit: curve_id = 1 ===

for (nm in names(cr$selection$assessments)) {
  a <- cr$selection$assessments[[nm]]
  cat(sprintf("Model: %s  |  Eligible: %s\n", nm,
              if (isTRUE(a$eligible)) "YES" else "NO"))

  if (!is.null(a$gates) && nrow(a$gates) > 0) {
    for (i in seq_len(nrow(a$gates))) {
      g <- a$gates[i, ]
      status <- if (g$passed) "  PASS" else "  FAIL"
      detail <- if (nchar(g$detail) > 0)
        paste0("  [", g$detail, "]") else ""
      cat(sprintf("  gate %-18s %s%s\n", g$gate, status, detail))
    }
  }

  if (is.finite(a$dynamic_range_log10)) {
    cat(sprintf("  dynamic range: %.3f log10 units  ",
                a$dynamic_range_log10))
    cat(sprintf("LLOQ = %.3f  ULOQ = %.3f (log10 AU/mL)\n",
                a$lloq, a$uloq))
  }
  cat("\n")
}
#> Model: logistic4  |  Eligible: YES
#>   gate at_bound             PASS
#>   gate vcov_condition       PASS
#>   gate rel_se               PASS
#>   gate dynamic_range        PASS  [dynamic range = 1.67 log10]
#>   dynamic range: 1.665 log10 units  LLOQ = 1.163  ULOQ = 2.829 (log10 AU/mL)
#> 
#> Model: logistic5  |  Eligible: NO
#>   gate at_bound             FAIL  [g at lower bound]
#>   gate vcov_condition       PASS
#>   gate rel_se               PASS
#>   gate dynamic_range        FAIL  [dynamic range = 0 log10 (need >= 0.5)]
#>   dynamic range: 0.000 log10 units  LLOQ = NA  ULOQ = NA (log10 AU/mL)
#> 
#> Model: gompertz4  |  Eligible: NO
#>   gate at_bound             FAIL  [a at upper bound]
#>   gate vcov_condition       PASS
#>   gate rel_se               PASS
#>   gate dynamic_range        FAIL  [dynamic range = 0 log10 (need >= 0.5)]
#>   dynamic range: 0.000 log10 units  LLOQ = NA  ULOQ = NA (log10 AU/mL)
#> 
#> Model: loglogistic5  |  Eligible: NO
#>   gate at_bound             FAIL  [a at lower bound; g at upper bound; b at upper bound]
#>   gate vcov_condition       PASS
#>   gate rel_se               FAIL  [a: rel_se=13.37; g: rel_se=8.58]
#>   gate dynamic_range        FAIL  [dynamic range = 0 log10 (need >= 0.5)]
#>   dynamic range: 0.000 log10 units  LLOQ = NA  ULOQ = NA (log10 AU/mL)

cat("Selected model   :", cr$selection$best_model_name, "\n")
#> Selected model   : logistic4
cat("Selection notes  :", cr$selection$fallback_reason %||%
      "All eligible; ranked by AIC", "\n")
#> Selection notes  :

17.2 D.2 Comparing AIC Ranking vs. Eligibility-Gated Ranking

A key diagnostic question is: how often does eligibility gating change the selected model relative to pure AIC? The selection$aic_best field (frequentist) records the AIC winner before gating, enabling direct comparison.

gating_changed <- vapply(names(freq_result$plates), function(cid) {
  cr <- freq_result$plates[[cid]]
  if (is.null(cr)) return(FALSE)
  aic_winner   <- cr$selection$aic_best
  gated_winner <- cr$selection$best_model_name
  !is.na(aic_winner) && !is.na(gated_winner) &&
    aic_winner != gated_winner
}, logical(1))

names(gating_changed) <- names(freq_result$plates)

comparison_df <- do.call(rbind, lapply(names(freq_result$plates), function(cid) {
  cr <- freq_result$plates[[cid]]
  if (is.null(cr)) return(NULL)
  data.frame(
    curve_id      = cid,
    aic_winner    = cr$selection$aic_best %||% NA_character_,
    gated_winner  = cr$selection$best_model_name,
    changed       = gating_changed[cid],
    fallback      = isTRUE(cr$selection$fallback),
    stringsAsFactors = FALSE
  )
}))

knitr::kable(comparison_df,
             caption = paste0(
               "AIC winner vs. eligibility-gated winner. ",
               "When 'changed = TRUE', the gating system overrode ",
               "the pure-AIC selection."))

17.3 D.3 Visualising Gate Outcomes Across All Models and Curves

# Collect gate outcomes into a flat data frame
gate_rows <- list()
for (cid in names(freq_result$plates)) {
  cr <- freq_result$plates[[cid]]
  if (is.null(cr)) next
  for (nm in names(cr$selection$assessments)) {
    a <- cr$selection$assessments[[nm]]
    if (is.null(a$gates) || nrow(a$gates) == 0) next
    for (i in seq_len(nrow(a$gates))) {
      gate_rows[[length(gate_rows) + 1]] <- data.frame(
        curve_id = cid, model = nm,
        gate     = a$gates$gate[i],
        passed   = a$gates$passed[i],
        detail   = a$gates$detail[i],
        stringsAsFactors = FALSE)
    }
    # Add "not applicable" rows for Bayesian-only or missing gates
    all_gates  <- c("at_bound", "vcov_condition", "rel_se", "dynamic_range")
    seen_gates <- a$gates$gate
    for (g in setdiff(all_gates, seen_gates)) {
      gate_rows[[length(gate_rows) + 1]] <- data.frame(
        curve_id = cid, model = nm,
        gate     = g, passed = NA, detail = "N/A",
        stringsAsFactors = FALSE)
    }
  }
}
gate_df <- do.call(rbind, gate_rows)

# Plot
gate_names  <- c("at_bound", "vcov_condition", "rel_se", "dynamic_range")
model_names <- unique(gate_df$model)
curve_ids   <- unique(gate_df$curve_id)

pass_col <- c("TRUE" = "#4dac26", "FALSE" = "#d01c8b", "NA" = "grey85")

par(mar = c(5, 11, 3, 8))
n_rows <- length(model_names) * length(curve_ids)
n_cols <- length(gate_names)

plot(NA, xlim = c(0.5, n_cols + 0.5), ylim = c(0.5, n_rows + 0.5),
     xaxt = "n", yaxt = "n", xlab = "", ylab = "",
     main = "Eligibility Gate Outcomes (Frequentist)")

row_idx <- 0
y_labels <- character(0)
for (cid in curve_ids) {
  for (nm in model_names) {
    row_idx  <- row_idx + 1
    y_labels <- c(y_labels, paste0(nm, " [c", cid, "]"))
    for (j in seq_along(gate_names)) {
      sub <- gate_df[gate_df$curve_id == cid &
                       gate_df$model == nm &
                       gate_df$gate == gate_names[j], ]
      col <- if (nrow(sub) == 0) "grey85"
             else pass_col[as.character(sub$passed[1])]
      rect(j - 0.45, row_idx - 0.45, j + 0.45, row_idx + 0.45,
           col = col, border = "white", lwd = 1.5)
    }
  }
}

axis(1, at = seq_along(gate_names), labels = gate_names,
     las = 2, cex.axis = 0.85)
axis(2, at = seq_len(n_rows), labels = y_labels,
     las = 1, cex.axis = 0.70)
legend("topright", inset = c(-0.22, 0),
       legend = c("Pass", "Fail", "N/A"),
       fill   = c("#4dac26", "#d01c8b", "grey85"),
       border = "white", cex = 0.85, xpd = TRUE, bty = "n")

Gate outcome matrix for all models across all curves (frequentist). Green = passed, red = failed, grey = not applicable.

18.1 E.1 Adding a New Model Family

To add a sixth model family to the curveR ecosystem, the following components must be implemented in curveRcore and the Stan models must be added to curveRbayes:

In curveRcore:

1. models.R — add the forward function newmodel(x, a, b, c, d, ...).
2. inverses.R — add inv_newmodel() and inv_newmodel_fixed().
3. derivatives.R — add dydx_newmodel() and d2x_newmodel().
4. gradients.R — add grad_newmodel(), grad_inv_newmodel_fixed(), and grad_y_newmodel_fixed(); register the model in the make_inv_and_grad_fixed() dispatcher.
5. utils.R — add the model name to available_models() and the parameter names to model_params().
6. formulas.R — add the NLS formula template to build_nls_formulas().

In curveRbayes:

7. Add inst/stan/hierarchical_newmodel.stan.
8. Register the .stan file path in stan_model_path().
9. Add the parameter extraction pattern in extract_curve_params().
10. Add data-adaptive prior computation in compute_dynamic_priors().
11. Add Stan data construction in build_stan_data().

No changes to curveRfreq are required beyond ensuring the model name is passed in model_names.

18.2 E.2 Using a Fixed Lower Asymptote from Blank Data

# Compute fixed_a from blank data (e.g. geometric mean of blanks)
blanks <- bead_assay_example$blanks[
  bead_assay_example$blanks$curve_id %in% c(1, 2, 3), ]

blank_geomean <- exp(mean(log(blanks$mfi[blanks$mfi > 0])))
cat("Geometric mean of blanks (raw MFI):", round(blank_geomean, 2), "\n")

# Validate and log10-transform for use with is_log_response = TRUE
fixed_a_raw <- curveRcore::validate_fixed_lower_asymptote(
  blank_geomean, verbose = TRUE)
fixed_a_log <- if (!is.null(fixed_a_raw)) log10(fixed_a_raw) else NULL

cat("Fixed a (log10 scale):", round(fixed_a_log, 4), "\n")

# Pass to fitting functions
freq_result_fixed <- fit_calibration_freq_multiplate(
  standards          = standards,
  samples            = samples,
  response_var       = "mfi",
  model_names        = c("logistic4", "gompertz4"),
  is_log_response    = TRUE,
  is_log_independent = TRUE,
  std_curve_conc     = 10000,
  fixed_a            = fixed_a_log,   # <-- fixed lower asymptote
  verbose            = FALSE
)

When fixed_a is supplied:

• Frequentist: the parameter \(a\) is substituted as a numeric constant in the NLS formula and does not appear in the optimisation or the covariance matrix. The gradient vector returned by make_inv_and_grad_fixed() omits the \(\partial x / \partial a\) term.
• Bayesian: a tight normal prior \(\mu_a \sim \mathcal{N}(\texttt{fixed\_a},\; \max(0.01\Delta_y, 10^{-4}))\) anchors the population-level lower asymptote at the specified value while allowing small deviations.

18.3 E.3 Adjusting Gate Thresholds for Specific Assay Contexts

All four gate thresholds are exposed as named arguments on both top-level fitting functions. The table below gives guidance on when to deviate from the defaults:

# Tighter thresholds for a regulatory submission context
freq_result_strict <- 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,
  pcov_threshold          = 15,     # tighter: 15% CV defines LOQ
  min_dynamic_range_log10 = 1.0,    # wider: require 10-fold range
  max_rel_se              = 2.0,    # stricter: SE < 2x estimate
  bound_tol               = 1e-3,   # slightly more permissive bound detection
  verbose                 = FALSE
)