21 Introduction to Algorithm Principles
This chapter explains how to choose statistical units, connect the main ivdtools workflows, and understand what key arguments change mathematically. It is intended as a bridge between function documentation and audit of the implementation.
21.1 Replicates, Multiple Samples, and Statistical Units
21.1.1 First Decide What One Row Represents
Before choosing a function, identify whether one row is a patient/sample result, a replicate measurement, a run, a bottle, a time point, or a paired method result. Replicates are not automatically independent subjects. If a function expects sample-level means, aggregate according to the protocol and preserve replicate count and SD.
Functions that can receive replicate results directly include precision(), linearity(), interference_*(), dilution_recovery(), spike_recovery(), sensitivity_*(), and c5_c95() when their grouping arguments are supplied. Functions that require paired or sample-level results include mcr(), roc(), raw_to_table(), reference_interval(), and many sample-size calculations. Check the help page and define the analysis unit before reshaping data.
21.1.2 Joint or Separate Analysis for Multiple Samples
Use by when the function is designed to estimate every sample or concentration separately, such as precision(..., by="sample"). Use a common model only when the scientific question and variance structure require pooling. split_by_factors() and split_by_cutpoints() create named data lists for independent, prespecified analyses; they must not be used repeatedly to search for a favorable result.
21.2 Outliers and Normality
21.2.1 Grubbs Test
For a finite sample, the Grubbs statistic is
\[G=\frac{\max_i|x_i-\bar x|}{s}.\]
The test generally identifies at most one point per call and assumes independent, approximately normal observations. A flagged observation is a candidate for investigation, not an instruction to delete it.
21.2.2 Generalized ESD
Generalized ESD removes the most extreme standardized observation sequentially and compares each \(R_i\) with a critical value. r is the maximum number of candidate points; it does not mean that r points will be removed.
21.2.3 Dixon and IQR
Dixon uses an endpoint spacing divided by the range and is limited to small samples, currently \(3\le n\le30\) with alpha 0.01 or 0.05. IQR flags values below \(Q_1-kIQR\) or above \(Q_3+kIQR\) and is descriptive rather than a probability test.
21.2.4 Automatic Normality Branch
normal_test(method="auto") uses a deterministic sample-size rule:
| Sample size | Test |
|---|---|
| \(n<3\) | No test |
| \(3\le n\le50\) | Shapiro–Wilk |
| \(51\le n\le5000\) | Anderson–Darling |
| \(n>5000\) | Lilliefors |
P values depend strongly on sample size. Interpret them with Q-Q plots, histograms, experimental records, and the assumptions of the downstream model.
21.3 Method Comparison and Reference-Material Bias
21.3.1 Direction, Units, and Preprocessing
In mcr(), the candidate is Y and the reference is X; candidate-minus-reference differences and percent denominators must be recorded. A log transformation, unit conversion, or prespecified cutpoint changes the estimand and should be treated as a separate analysis.
21.3.2 Correlation Options
Pearson measures linear co-variation, Spearman measures monotone rank association, and Kendall measures concordance of ranks. None is an agreement measure. Use regression and Bland–Altman methods for bias and agreement.
21.3.3 Regression Methods
| Method | X error | Y variance | Point-estimate principle | Typical use |
|---|---|---|---|---|
ols |
Negligible | Constant or optional weights | Minimize vertical squared residuals | Reference X has negligible error |
wls |
Negligible | Changes with concentration | Weighted vertical residuals | Difference SD/CV changes with concentration |
deming |
Included | Constant error ratio \(λ\) | Errors-in-variables regression | Both procedures have measurement error |
wdeming |
Included | Common weights, constant CV, or precision profiles | Weighted/iterative errors-in-variables fit | Error changes with concentration |
pb |
Included | No normal-error assumption | Rank estimate of pairwise slopes | Linear relation with outliers or mixed variation |
OLS minimizes \(\sum_i(y_i-a-bx_i)^2\). Weighted regression minimizes a weighted version. Deming uses the relative error variance \(λ\) and should estimate it from independent precision evidence. Passing–Bablok is robust to some outliers but is not a guarantee against every data problem.
Common WLS weights include 1/y, 1/y^2, 1/x, 1/x^2, constant_cv, residual, and a user-supplied nonnegative finite column. 1/x^2 represents constant CV on the reference scale; residual updates a concentration-dependent residual SD and is iterative.
21.3.4 Confidence Intervals and Bland–Altman
ci_method="auto" uses analytic covariance for OLS/WLS, jackknife refitting for Deming variants, and rank intervals for Passing–Bablok. Bootstrap methods require enough replicates and a fixed seed for reproducibility. Bland–Altman center, SD, and LoA describe the distribution of paired differences; agree.level controls the proportion covered by LoA, whereas conf.level controls uncertainty in the estimated LoA.
21.3.5 Reference-Material Bias
Reference-material bias is not the same as clinical-sample method-comparison bias. reference_bias() combines the measurement-mean uncertainty and assigned-value uncertainty when supplied. State the assigned-value source, its uncertainty, the bias direction, and whether the material is commutable.
21.4 Reference-Material Commutability
EP14/EP30 methods compare materials with the clinical-sample procedure relationship, usually using Deming prediction intervals. IFCC Part 2 compares reference-material bias with clinical-sample bias and combines uncertainty. IFCC Part 3 evaluates before/after calibration. Different paths can produce different conclusions because they target different estimands.
21.5 Qualitative Performance, C5/C95, and ROC
For a four-fold table \(a=TP,b=FP,c=FN,d=TN\),
\[Se=\frac{a}{a+c},\quad Sp=\frac{d}{b+d},\quad PPV=\frac{a}{a+b},\quad NPV=\frac{d}{c+d}.\]
Kappa adjusts observed agreement for marginal agreement; PABAK uses \(2P_o-1\) and is not the same estimand. McNemar uses discordant pairs. C5/C95 solve for mean responses whose modeled positive probability is 0.05 or 0.95, using a precision profile. ROC AUC is the probability of correctly ranking a random positive above a random negative under the selected direction; it is not accuracy at a cutoff.
21.6 Precision, Bottle Variation, and Reference Intervals
For nested precision designs, total variance is the sum of variance components. variance() estimates components, ci() estimates intervals, and profile() models concentration-related SD/CV. Bottle ANOVA tests group means; it does not quantify whether a difference is practically acceptable. uncertainty_homogeneity() produces a standard uncertainty component.
Reference intervals have two independent settings: central coverage interval and endpoint confidence ci. Nonparametric, biweight, type-6, parametric, and Huber methods make different distributional and robust-estimation assumptions. A normality or outlier test is not an automatic method-selection or deletion rule.
21.7 Response-Equation Fitting and Linearity
fit_equation() fits registered response equations such as linear, quadratic, exponential, power, 4PLC, 5PLC, Gaussian, Michaelis–Menten, sinusoidal, and cubic forms. Equation fitting describes calibration response and should not be confused with EP06 method linearity. linearity() specifically evaluates level means, residuals, weights, and ADL. A high calibration \(R^2\) does not establish method linearity.
21.8 Interference, Dilution, Spike Recovery, and Hook Effect
Paired interference estimates a test-minus-control effect; dose response estimates how that effect changes with interferent concentration; patient-sample comparison assesses method differences in the presence of a factor. Dilution recovery corrects results by the dilution factor, spike recovery compares an observed increment with a theoretical mass-balance increment, and hook-effect analysis uses raw response at known high concentrations.
21.9 Analytical Sensitivity
LoB describes the blank upper tail, LoD describes distinguishability from blank, and LoQ adds a precision or total-error criterion. sensitivity_lob(), sensitivity_lod(), sensitivity_loq(), and sensitivity_verify() keep establishment and verification separate. A reported LoQ is constrained not to be below LoD when LoD is supplied.
21.10 Quality Control and Stability
For QC, \(z_i=(x_i-\mu_0)/\sigma_0\). The current rules use strict \(>2\), \(>3\), and \(>1\) thresholds for 1-2s, 1-3s, and 4-1s, while R-4s uses \(|z|\ge2\) and a range \(|\Delta z|\ge4\). qc_chart() scans one sequence in row order; it does not implement cross-level multirule decisions. Stability functions estimate node bias, time trends, time to a limit, Arrhenius acceleration, and MKT; they do not turn a model extrapolation into observed evidence.
21.11 Measurement Uncertainty and Traceability
Type A components come from statistical repeatability or replication; Type B components come from certificates, limits, distributions, or prior information. uncertainty_combine() combines standard uncertainties, uncertainty_propagate() applies a measurement model, and uncertainty_traceability() accumulates components across a traceability chain. Record the input quantity, standard uncertainty, sensitivity coefficient, degrees of freedom, covariance assumptions, and coverage factor.
21.12 Parameter Choices Change the Estimand
| Parameter | Mathematical object changed | Not changed | Selection principle |
|---|---|---|---|
method |
Loss function or rank estimator | X/Y direction | Match error and difference structure |
lambda |
Relative Deming penalty | Raw observations | Estimate independently, never tune to slope |
weights |
Contribution or concentration-dependent variance | Difference direction | Support with precision evidence |
ci_method |
Parameter uncertainty estimator | Point estimates | Prespecify the default or sensitivity path |
conf.level |
Coverage of parameter/center/LoA interval | agree.level |
Usually 0.95 |
agree.level |
Fraction of individual differences within LoA | CI confidence | Define agreement clinically |
percent_denominator |
Percent-difference denominator | Absolute difference | Choose reference or symmetric mean deliberately |
bootstrap_replicates |
Monte Carlo error | Target distribution | Use enough replicates and a fixed seed |
| cutpoints | Data ranges analyzed | Downstream estimator | Define clinical intervals before results |
The safest audit trail records the data snapshot, analysis unit, transformations, function arguments, package version, random seed, complete output, and any deviations from the standard algorithm.
21.13 Detailed Algorithm Specification
21.13.1 OLS and WLS Point Estimates
Ordinary least squares minimizes
\[ Q(a,b)=\sum_{i=1}^n(y_i-a-bx_i)^2. \]
Weighted least squares minimizes
\[ Q_w(a,b)=\sum_{i=1}^n w_i(y_i-a-bx_i)^2. \]
Let
\[ W=\sum_iw_i,\qquad \bar x_w=\frac{\sum_iw_ix_i}{W},\qquad \bar y_w=\frac{\sum_iw_iy_i}{W}, \]
\[ S_{xx,w}=\sum_iw_i(x_i-\bar x_w)^2,\qquad S_{xy,w}=\sum_iw_i(x_i-\bar x_w)(y_i-\bar y_w). \]
Then
\[ \hat b=\frac{S_{xy,w}}{S_{xx,w}},\qquad \hat a=\bar y_w-\hat b\bar x_w. \]
OLS is the special case w_i = 1. The weight choices change the assumed concentration dependence of the error variance:
| weights | Actual weight | Implied variance shape | Main distinction |
|---|---|---|---|
| NULL | unweighted | Var(Y | X) is constant under OLS |
| equal | 1 | Same as OLS | Explicit entry to the WLS interface |
| 1/y | 1/abs(y_i) | Var(Y_i) proportional to abs(y_i) | Depends on candidate observations |
| 1/y^2 | 1/y_i^2 | SD(Y_i) proportional to abs(y_i) | Empirical constant-CV candidate weighting |
| 1/x | 1/abs(x_i) | Var(Y_i) proportional to abs(x_i) | Depends on reference observations |
| 1/x^2 | 1/x_i^2 | SD(Y_i) proportional to abs(x_i) | Empirical constant-CV reference weighting |
| constant_cv | 1/x_i^2 | Same as above | Requires x_i > 0 in WLS |
| residual | 1 / estimated SD_i^2 | SD described by a linear function of concentration | WLS only; iteratively updated |
| user column | user-supplied w_i | Defined by the user | Must be finite and nonnegative; zero excludes the point from the objective |
Generic built-in weights such as 1/x protect division by zero with a machine-precision floor; this prevents a computational failure but does not make a relative-error model scientifically meaningful at zero concentration. constant_cv therefore rejects nonpositive concentrations directly.
Residual weighting first fits OLS. With residual e_i^(k) at iteration k, it fits
\[ |e_i^{(k)}|=\gamma_0^{(k)}+\gamma_1^{(k)}x_i+\varepsilon_i, \qquad w_i^{(k+1)}=\frac{1}{\left(\hat\gamma_0^{(k)}+ \hat\gamma_1^{(k)}x_i\right)^2}. \]
It then refits WLS. weight_iterations is the number of updates and defaults to four. Changing it can change the final weights and regression coefficients, but should not be used to match a published interval. A nonpositive or nonfinite predicted SD stops the procedure.
The analytic OLS/WLS covariance is
\[ \hat\sigma_w^2=\frac{\sum_iw_ie_i^2}{n-2},\qquad \widehat{\operatorname{Cov}}(\hat a,\hat b) =\hat\sigma_w^2(\mathbf X^\mathsf T\mathbf W\mathbf X)^{-1}, \]
and coefficient intervals are
\[ \hat\theta\pm t_{1-\alpha/2,n-2}SE(\hat\theta). \]
This is equivalent to the EP09c Appendix H calculation when the same final weights are used. YY/T 1789.2 Equation (8) does not multiply the medical-decision-level bias interval by the weighted residual scale; the two expressions need no residual-scale adjustment only when the weights are known absolute inverse variances.
21.13.2 Deming Regression and lambda
Deming regression defines
\[ S_{xx}=\sum_i(x_i-\bar x)^2,\quad S_{yy}=\sum_i(y_i-\bar y)^2,\quad S_{xy}=\sum_i(x_i-\bar x)(y_i-\bar y). \]
In the source, lambda is the ratio of Y-error variance to X-error variance. The slope and intercept are
\[ U=S_{yy}-\lambda S_{xx},\qquad b=\frac{U+\sqrt{U^2+4\lambda S_{xy}^2}}{2S_{xy}}, \qquad a=\bar y-b\bar x. \]
Thus
\[ \lambda=\frac{\sigma_Y^2}{\sigma_X^2}. \]
lambda = 1 means equal error variances; lambda > 1 means the candidate Y has greater error; lambda < 1 means the reference X has greater error. Reversing the ratio changes the Deming slope, so always enter candidate variance divided by reference variance. The reciprocal required by the lower-level algorithm is handled inside the package and does not change the user interface.
lambda = “replicates” estimates the variance ratio from summarized replicate SDs. If sample i has n_Xi and n_Yi replicates, first calculate degrees-of-freedom-pooled single-measurement variances:
\[ s_{p,X}^2= \frac{\sum_i(n_{Xi}-1)s_{Xi}^2}{\sum_i(n_{Xi}-1)},\qquad s_{p,Y}^2= \frac{\sum_i(n_{Yi}-1)s_{Yi}^2}{\sum_i(n_{Yi}-1)}. \]
Because the regression uses sample means, the input variance ratio is
\[ \hat\lambda_{\bar X,\bar Y} =\frac{s_{p,Y}^2/n_Y}{s_{p,X}^2/n_X} =\frac{s_{p,Y}^2}{s_{p,X}^2}\frac{n_X}{n_Y}. \]
The two methods may have different replicate counts, but each method’s count must be constant across samples. Otherwise one scalar lambda cannot describe different sample-mean variances and a precision profile should be used. Without replicate counts, the source assumes equal counts and n_X/n_Y cancels. lambda is estimated once; jackknife and bootstrap refits treat it as fixed external input rather than re-estimating it in each subsample.
21.13.3 Three Weighted-Deming Variants
For fixed observation weights, method = “wdeming” first sets w_tilde_i = w_i / sum_j w_j and replaces the ordinary central moments with
\[ \bar x_w=\sum_i\tilde w_ix_i,\quad \bar y_w=\sum_i\tilde w_iy_i, \]
\[ S_{xx,w}=\sum_i\tilde w_i(x_i-\bar x_w)^2,\quad S_{yy,w}=\sum_i\tilde w_i(y_i-\bar y_w)^2,\quad S_{xy,w}=\sum_i\tilde w_i(x_i-\bar x_w)(y_i-\bar y_w), \]
then applies the closed-form Deming solution. This gives one common weight to each pair and cannot represent separate concentration-dependent variance functions for X and Y.
constant_cv in WLS is simply 1/x_i^2; in weighted Deming it calls Linnet’s iterative constant-CV algorithm. At iteration k, the weight is formed from the mean of the two latent readings:
\[ w_i^{(k)}=\left(\frac{x_i^{*(k)}+y_i^{*(k)}}{2}\right)^{-2}. \]
The Deming line is updated with this weight. If d_i = y_i - a - b x_i and lambda = sigma_Y^2 / sigma_X^2, the latent updates are
\[ x_i^*=x_i+\frac{b\,d_i}{\lambda+b^2},\qquad y_i^*=y_i-\frac{\lambda d_i}{\lambda+b^2}. \]
Weights, line, and latent values are updated until relative change is below epsilon. Both concentrations must be strictly positive. After deleting a sample or drawing a bootstrap sample, the complete iteration starts again rather than reusing full-data weights.
precision_profile permits separate variance functions
\[ g(\mu)=\operatorname{Var}(X\mid\mu),\qquad h(\eta)=\operatorname{Var}(Y\mid\eta),\quad \eta=a+b\mu. \]
For each latent true value mu_i, the algorithm minimizes
\[ L(a,b,\boldsymbol\mu)= \sum_i\left[ \frac{(x_i-\mu_i)^2}{g(\mu_i)}+ \frac{(y_i-a-b\mu_i)^2}{h(a+b\mu_i)}+ \log\{g(\mu_i)h(a+b\mu_i)\} \right]. \]
Given a and b and current variances, the latent value update is
\[ \mu_i= \frac{h_i x_i+g_i b(y_i-a)}{h_i+g_i b^2}. \]
variance_models can contain variance-returning functions or precision objects that have already been processed by variance() and profile(). profile_components selects the variance component on each side; the default is “total” for both sides, and the predicted SD from the VFP is squared to obtain variance. Choosing “error”, “day”, or another component changes g or h directly; the available names depend on the precision() formula and must match the intended error layer. Values outside the precision-study range generate an extrapolation warning, while nonpositive or nonfinite predicted variances stop the fit. External precision profiles remain fixed during resampling.
These three variants are not interchangeable: fixed weights change each point’s empirical influence, constant-CV Deming models relative error on both axes, and precision-profile Deming allows different nonconstant variance functions on the two axes.
21.13.4 Regression-Parameter Confidence Intervals
ci_method determines the uncertainty algorithm for the intercept and slope. The auto mapping is:
| Regression method | auto method | Main basis |
|---|---|---|
| OLS, WLS | analytic | Analytic covariance from the regression design |
| Deming, constant-CV Deming, precision-profile Deming | jackknife | EP09c Appendix K1 |
| Passing–Bablok | rank | EP09c Appendix I |
Jackknife deletes observation i and completely refits, producing theta_(−i) = (a_(−i), b_(−i))^T. Its covariance is
\[ \widehat{\operatorname{Cov}}_{\mathrm{jack}}(\hat{\boldsymbol\theta}) =\frac{n-1}{n}\sum_{i=1}^n (\hat{\boldsymbol\theta}_{(-i)}-\bar{\boldsymbol\theta}_{(-\cdot)}) (\hat{\boldsymbol\theta}_{(-i)}-\bar{\boldsymbol\theta}_{(-\cdot)})^\mathsf T. \]
The interval remains centered on the full-data estimate:
\[ \hat\theta\pm t_{1-\alpha/2,n-2}SE_{\mathrm{jack}}(\hat\theta). \]
Bootstrap always samples complete (x_i, y_i) pairs with replacement. Because mcr() accepts one row per ID, there is no splitting of replicate rows from the same sample. Every bootstrap sample completely refits the regression. For bootstrap parameter estimates theta_star^(r),
\[ CI_{\mathrm{percentile}} =\left[q_{\alpha/2}(\hat\theta^*), q_{1-\alpha/2}(\hat\theta^*)\right], \]
\[ SE_{\mathrm{boot}}(\hat\theta)=SD(\hat\theta^*),\qquad CI_{\mathrm{standard}} =\hat\theta\pm t_{1-\alpha/2,n-2}SE_{\mathrm{boot}}(\hat\theta). \]
Both methods also use the joint bootstrap intercept–slope draws to calculate bias at any medical decision level. The source requires at least 5,000 replicates; seed controls reproducibility and is not a tuning parameter. Failed fits are removed and counted; the process stops if fewer than half or fewer than 100 fits are valid.
analytic is allowed only for OLS/WLS, and rank only for Passing–Bablok. Although other methods can explicitly request jackknife or bootstrap, that changes the interval question and should not be selected merely because it is closer to a published value.
21.13.5 Passing–Bablok
The source forms all pairwise slopes with x_j not equal to x_i:
\[ b_{ij}=\frac{y_j-y_i}{x_j-x_i}. \]
Exactly identical pairs contribute no slope; vertical pairs are retained as signed infinities. Algorithm I first counts slopes below −1 as K, removes slopes exactly equal to −1, and cyclically shifts the ranks of the point estimate and interval by K. The shifted median is b_hat, and
\[ \hat a=\operatorname{median}_i(y_i-\hat b x_i). \]
If m valid slopes remain, define
\[ C_\gamma=z_{1-\alpha/2} \sqrt{\frac{n(n-1)(2n+5)}{18}},\qquad M_1=\max\left\{1,\operatorname{round}\left(\frac{m-C_\gamma}{2}\right)\right\}, \]
and M_2 = m - M_1 + 1. The slope limits are the shifted M_1-th and M_2-th ordered slopes; intercept limits are
\[ a_L=\operatorname{median}_i(y_i-b_Ux_i),\qquad a_U=\operatorname{median}_i(y_i-b_Lx_i). \]
This is ci_method = “rank” for regression parameters. A rank interval does not directly provide a joint interval for bias at a medical decision level. Therefore bias() and predict() use paired percentile bootstrap for Passing–Bablok under ci_method = “auto”, rather than the jackknife excluded by EP09c Appendix K1. The object may retain a jackknife parameter covariance for explicitly requested extended prediction calculations, but it is not the default Passing–Bablok bias interval.
21.13.6 Differences, Percentage Differences, and Outlier Analysis
Absolute differences always use candidate minus reference:
\[ d_i=y_i-x_i. \]
The three percent_denominator choices target different estimands:
\[ p_{i,\mathrm{reference}}=100\frac{y_i-x_i}{x_i},\qquad p_{i,\mathrm{mean}}=100\frac{y_i-x_i}{(x_i+y_i)/2},\qquad p_{i,\mathrm{candidate}}=100\frac{y_i-x_i}{y_i}. \]
reference answers how far the candidate is above or below the reference and is the default. mean changes only sign when the two methods are exchanged and is useful without a natural reference. candidate uses the candidate as the denominator. The denominator changes the value, distribution, and outlier status, not just the plot label. outlier(type = “percent”) and bland_altman(type = “percent”) use the same options so the detected quantity and difference plot remain aligned.
outlier(type = “difference”) tests d_i; type = “percent” tests the selected p_i. method calls the Grubbs, generalized ESD, Dixon, or IQR rule described earlier. The first three use alpha; IQR uses coef for the fence width. Outlier tests flag observations and do not automatically remove them from later regression or Bland–Altman analyses.
21.13.7 Bland–Altman Axis, Center, and Limits of Agreement
x_axis changes only the horizontal coordinate:
\[ m_i=\frac{x_i+y_i}{2},\qquad x_i,\qquad y_i \]
for mean, reference, and candidate, respectively; it does not change the vertical difference. type selects
\[ \text{difference}:y_i-x_i,\qquad \text{ratio}:\frac{y_i}{x_i},\qquad \text{percent}:p_i. \]
The parametric method centers on the mean:
\[ \hat\delta=\bar d,\qquad CI_\delta=\bar d\pm t_{1-\alpha/2,n-1}\frac{s_d}{\sqrt n}. \]
For agreement level p, let z_p = z_(1+p)/2. The limits of agreement are
\[ L_-=\bar d-z_ps_d,\qquad L_+=\bar d+z_ps_d. \]
The standard error and interval for each limit are
\[ SE_L=s_d\sqrt{\frac1n+\frac{z_p^2}{2(n-1)}},\qquad CI(L_\pm)=L_\pm\pm t_{1-\alpha/2,n-1}SE_L. \]
agree.level controls the proportion of individual differences intended to lie within the LoA; conf.level controls statistical confidence in the center and LoA estimates. They have different meanings: changing agree.level from 0.95 to 0.99 widens the LoA, while changing conf.level from 0.95 to 0.99 primarily widens the confidence bands around the LoA.
The nonparametric method uses the ordinary sample median d_(ceil(n/2)), averaging the two middle values for even n. median_ci_method = “order_statistic” finds the largest k whose binomial coverage
\[ C(k)=1-2P\{\operatorname{Bin}(n,0.5)\le k-1\} \]
is at least conf.level, and reports
\[ [d_{(k)},d_{(n-k+1)}]. \]
This reaches at least the requested coverage and is usually conservative. interpolated uses the EP09c Appendix F continuous-rank normal approximation
\[ r_L=\frac{n+1-z_{1-\alpha/2}\sqrt n}{2},\qquad r_U=\frac{n+1+z_{1-\alpha/2}\sqrt n}{2}. \]
If r = j + q with j = floor(r) and 0 <= q < 1, linear interpolation gives
\[ d_{(r)}=d_{(j)}+q\{d_{(j+1)}-d_{(j)}\}. \]
The interpolated interval has nominal coverage and its endpoint need not equal an observed value; the order-statistic interval is discrete and conservative.
method = “hodges_lehmann” uses the median of all Walsh averages rather than the ordinary median:
\[ \hat\delta_{HL}=\operatorname{median}_{i\le j}\left(\frac{d_i+d_j}{2}\right). \]
Its center interval uses a Wilcoxon/Tukey normal approximation without continuity correction. The HL estimator is a pseudomedian of the difference distribution; it has an interpretation similar to the ordinary median only under conditions such as symmetry, and hl_estimate should not be labeled median_diff.
Nonparametric and HL LoA use empirical quantiles
\[ L_-=Q_{(1-p)/2}(d),\qquad L_+=Q_{1-(1-p)/2}(d), \]
with R type-7 linear quantiles. LoA confidence intervals are obtained by resampling d_i with replacement, recalculating both empirical quantiles, and taking bootstrap percentiles. bootstrap_replicates defaults to 5,000 with a minimum of 100; seed is only for reproducibility. These empirical LoA and bootstrap intervals are ivdtools extensions and should not be confused with the center-estimation intervals in EP09c Appendix F.
21.13.8 Bias at Medical Decision Levels and Prediction Intervals
For x_0, the fitted response is y_hat_0 = a_hat + b_hat x_0, with fitted-mean variance
\[ V_{\mathrm{fit}}(x_0)= SE_a^2+x_0^2SE_b^2+2x_0\operatorname{Cov}(a,b). \]
The mean interval uses sqrt(V_fit), while a single-result prediction interval uses sqrt(V_fit + s_e^2). Candidate-minus-reference bias is
\[ \widehat{\operatorname{bias}}(x_0) =(\hat a+\hat b x_0)-x_0 =\hat a+(\hat b-1)x_0. \]
Subtracting fixed x_0 adds no variance, so under analytic or jackknife parameter covariance
\[ SE\{\widehat{\operatorname{bias}}(x_0)\} =\sqrt{SE_a^2+x_0^2SE_b^2+2x_0\operatorname{Cov}(a,b)}. \]
Percentile bootstrap calculates B_star^(r)(x_0) = a_star^(r) + [b_star^(r) - 1] x_0 for each parameter draw and takes its quantiles. Standard bootstrap uses the SD of these bias draws and t_(n-2). With ci_method = “auto”, bias() and predict() use percentile bootstrap for Passing–Bablok and reuse the regression object’s selected analytic, jackknife, or bootstrap uncertainty for other methods. Existing compatible bootstrap samples are reused.
interval = “confidence” describes uncertainty in the mean regression relation or mean bias at x_0; interval = “prediction” adds s_e^2 and describes the scatter of one future candidate result; “both” returns both. Prediction intervals are an ivdtools extension and are not EP09c regression-parameter or bias-confidence intervals.
21.13.9 Reference-Material Bias
At each level, bias is the measurement mean minus the assigned reference value:
\[ b=\bar x-x_{\mathrm{ref}},\qquad u_{\bar x}=\frac{s}{\sqrt n},\qquad U_{\bar x}=k u_{\bar x}. \]
The source combines the supplied reference-value uncertainty directly with U_xbar:
\[ h=\sqrt{U_{\bar x}^2+U_{\mathrm{ref}}^2},\qquad [b-h,b+h]. \]
It does not divide reference uncertainty by an additional coverage factor. Relative bias is 100 b / x_ref; a zero reference value returns missing.
21.14 Reference-Material Commutability
21.14.1 EP14-Style Prediction Interval
For each clinical sample, first calculate replicate means and SDs in both procedures. The pooled repeatability variance for a procedure is
\[ s_p^2= \frac{\sum_j(n_j-1)s_j^2}{\sum_j(n_j-1)}. \]
The source sets lambda = s_p,Y^2 / s_p,X^2 and fits a Deming line to clinical-sample means, using population central moments divided by n. The slope variance is
\[ \operatorname{Var}(\hat b)= \frac{\hat b^2\max(0,S_{xx}S_{yy}-S_{xy}^2)}{nS_{xy}^2}. \]
At reference material x_R, the prediction variance is
\[ V_R=(x_R-\bar x)^2\operatorname{Var}(\hat b) +\left(1+\frac1n\right) \left(\frac{\hat b^2s_{p,X}^2}{n_{R,X}} +\frac{s_{p,Y}^2}{n_{R,Y}}\right). \]
Degrees of freedom are the smaller integer among the two procedures’ pooled repeatability degrees of freedom. The material is labeled commutable when its Y mean lies within a_hat + b_hat x_R +/- t sqrt(V_R).
21.14.2 IFCC-Style Difference Evaluation
Clinical-sample bias is b_i = y_i - x_i. When a reference measurement procedure is specified, the sorting coordinate is x_i; otherwise it is (x_i + y_i) / 2. After sorting, calculate
\[ s_{\mathrm{MSSD}} =\sqrt{\frac{\sum_{i=1}^{n-1}(b_{i+1}-b_i)^2}{2(n-1)}}. \]
The difference between reference-material bias and the selected clinical-sample mean bias is
\[ d=(y_R-x_R)-\bar b_C. \]
The clinical subset may be all samples, the q closest samples, or half on either side of the target. Its standard uncertainty is
\[ u_d=\sqrt{\frac{s_{X,R}^2}{p_X}+\frac{s_{Y,R}^2}{p_Y} +\frac{s_C^2}{q}}. \]
The expanded interval is d +/- k u_d, with k = 1.9 by default. If the complete interval is within [-Delta, Delta], the material is commutable; if it is entirely on one side outside the allowable band, it is noncommutable; otherwise the result is indeterminate.
21.15 Qualitative Performance, C5/C95, and ROC
21.15.1 Four-Fold Proportions and Intervals
The source treats the first row and first column as positive in the four-fold table and calculates
\[ Se=\frac{TP}{TP+FN},\quad Sp=\frac{TN}{TN+FP},\quad PPV=\frac{TP}{TP+FP},\quad NPV=\frac{TN}{TN+FN}, \]
\[ Acc=\frac{TP+TN}{n},\qquad Prev=\frac{TP+FN}{n}. \]
For x successes, the Wald interval is
\[ \hat p\pm z_{1-\alpha/2}\sqrt{\frac{\hat p(1-\hat p)}n}; \]
the continuity-corrected version adds 1/(2n) to the half-width. The Wilson center and half-width are
\[ c=\frac{\hat p+z^2/(2n)}{1+z^2/n},\qquad h=\frac{z}{1+z^2/n} \sqrt{\frac{\hat p(1-\hat p)}n+\frac{z^2}{4n^2}}. \]
Wilson continuity correction subtracts and adds 1/(2n) outside c +/- h. Agresti–Coull uses
\[ \tilde n=n+z^2,\qquad \tilde p=\frac{x+z^2/2}{\tilde n},\qquad \tilde p\pm z\sqrt{\frac{\tilde p(1-\tilde p)}{\tilde n}}. \]
Jeffreys uses Beta(x + 1/2, n - x + 1/2) quantiles. Clopper–Pearson uses Beta(x, n - x + 1) and Beta(x + 1, n - x) quantiles, with separate handling at x = 0 and x = n. All proportion intervals are finally clipped to [0, 1].
Given prevalence pi, predictive values are recalculated as
\[ PPV_\pi=\frac{Se\,\pi}{Se\,\pi+(1-Sp)(1-\pi)},\qquad NPV_\pi=\frac{Sp(1-\pi)}{(1-Se)\pi+Sp(1-\pi)}. \]
Intervals simply substitute sensitivity and specificity endpoints and do not propagate their covariance.
21.15.2 Likelihood Ratios, Agreement, and McNemar
\[LR^+=\frac{Se}{1-Sp},\qquad LR^-=\frac{1-Se}{Sp}.\]
Intervals are calculated on the log scale, for example
\[ SE\{\log(LR^+)\}= \sqrt{\frac1{TP}-\frac1{TP+FN}+\frac1{FP}-\frac1{FP+TN}}. \]
The odds ratio is OR = TP * TN / (FP * FN), with log-scale SE equal to the square root of the sum of the four cell reciprocals. A zero cell returns a missing interval; no continuity correction is added automatically.
Observed and chance agreement are
\[ p_o=\frac{TP+TN}{n},\qquad p_e=\frac{(TP+FP)(TP+FN)+(FN+TN)(FP+TN)}{n^2}, \]
\[\kappa=\frac{p_o-p_e}{1-p_e}.\]
Kappa intervals use an explicit asymptotic variance of the four-cell joint proportions and are clipped to [-1, 1]. The PABAK branch uses
\[PABAK=2p_o-1,\qquad \operatorname{Var}(PABAK)=\frac{4p_o(1-p_o)}n.\]
This branch is triggered by supplying prevalence, but the prevalence value does not enter the formula.
When b + c < 25, McNemar uses a two-sided exact test with B ~ Binomial(b + c, 0.5); otherwise it calculates
\[\chi^2=\frac{(|b-c|-1)^2}{b+c}\]
with one degree of freedom. When b + c = 0, p = 1.
21.15.3 C5/C95
Raw replicates are first summarized by sample-level mean mu_j, SD s_j, and n_j. At least three levels with n_j >= 2 and s_j > 0 are required. A concentration-dependent SD s(mu) is then fitted.
If a positive event is defined as result >= cutoff c,
\[P(Y\ge c\mid\mu)=1-\Phi\!\left(\frac{c-\mu}{s(\mu)}\right);\]
the reverse event uses the corresponding Phi expression. The source searches for sign changes between the target probability and the model on a 1,001-point grid, then refines the root by one-dimensional solving; if no root is found, it expands the search range up to 16 times. If multiple roots are found, it selects the root closest to the cutoff. This is inversion of a normal response distribution with a precision profile, not binomial probit regression.
21.15.4 Empirical ROC, AUC, and Multivariable Models
Candidate thresholds are -infinity, all unique observations, and +infinity. After sorting ROC points by increasing false-positive rate, AUC is trapezoidal area:
\[ A=\sum_i\frac{TPR_i+TPR_{i+1}}2(FPR_{i+1}-FPR_i). \]
For 0 < A < 1, the source uses the Hanley–McNeil approximation
\[Q_1=\frac{A}{2-A},\qquad Q_2=\frac{2A^2}{1+A},\]
\[ SE(A)=\sqrt{ \frac{A(1-A)+(n_+-1)(Q_1-A^2)+(n_--1)(Q_2-A^2)}{n_+n_-}}. \]
The AUC interval is A +/- 1.96 SE(A), clipped to [0, 1]; no SE is computed when A = 0 or 1. The Youden threshold maximizes Se + Sp - 1; the upper-left-corner rule minimizes sqrt((1 - Se)^2 + (1 - Sp)^2). With ties, the first point in the current order is returned.
Two AUCs from the same subjects cannot be compared as independent samples. auc_compare(method = “ep24”) uses the EP24-A2 Hanley–McNeil paired approximation: calculate Pearson correlations between the two scores within positive and negative subjects (Kendall can be selected for ordinal scores), average them, linearly interpolate the area correlation r from EP24 Table 5 using the mean AUC, and calculate
\[ SE(A_1-A_2)=\sqrt{SE_1^2+SE_2^2-2rSE_1SE_2}. \]
auc_compare(method = “delong”) builds a covariance matrix Sigma from rank contributions of both scores for each positive and negative subject and sets c = (1, -1)^T:
\[ SE(A_1-A_2)=\sqrt{\mathbf c^\mathsf T\widehat\Sigma\mathbf c}. \]
Both methods use a normal-approximation test of H0: A1 - A2 = 0. When scores are fitted probabilities from a training-set logistic model, they are treated as fixed scores; model-fit uncertainty and apparent-performance optimism are not corrected.
A multivariable model fits
\[ \operatorname{logit}(p_i)=\beta_0+\sum_{j=1}^q\beta_jx_{ij} \]
by binomial logistic regression and uses training predictions to construct the same empirical ROC. The linear-predictor interval is eta_hat +/- 1.96 SE(eta_hat), followed by the logistic inverse. Probability SE is SE(p_hat) = SE(eta_hat) p_hat (1 - p_hat), and the classification threshold is fixed at 0.5.
21.16 Precision, Bottle Variation, and Reference Intervals
21.16.1 Variance Components and Precision Profiles
The source first checks balance: balanced designs use ANOVA variance components and unbalanced designs use REML. Each random level returns SD_j = sqrt(VC_j) and CV_j = 100 SD_j / y_bar; total SD is
\[SD_T=\sqrt{\sum_j\widehat{VC}_j}.\]
Component intervals use the two-sided SD/CV intervals from the VCA inference object; the source does not implement a separate Satterthwaite–chi-square formula.
Concentration-related variance profiles use VFP. The ten registered models are:
| Model | Variance or CV form |
|---|---|
| 1 | sigma^2 = beta_1 |
| 2 | sigma^2 = beta_1 mu^2 |
| 3 | sigma^2 = beta_1 + beta_2 mu^2 |
| 4 | sigma^2 = (beta_1 + beta_2 mu)^2 |
| 5 | sigma^2 = beta_1 + beta_2 mu^K |
| 6 | sigma^2 = beta_1 + beta_2 mu + beta_3 mu^J |
| 7 | sigma^2 = beta_1 + beta_2 mu^J |
| 8 | sigma^2 = (beta_1 + beta_2 mu)^J |
| 9 | sigma^2 = beta_1 mu^J |
| 10 | CV = beta_1 mu^J |
When more than one candidate fits, the model with minimum AIC is selected. For CV input, SD is first reconstructed from the mean and then squared as VC for fitting.
21.16.2 Bottle Variation
Bottle, lot, or other grouping variables enter a fixed-effect analysis of variance. Within-group descriptive intervals are
\[ \bar x_j\pm t_{1-\alpha/2,n_j-1}\frac{s_j}{\sqrt{n_j}}. \]
Residuals receive a Shapiro–Wilk test, between-group variances receive Bartlett’s test, and pairwise comparisons use Tukey HSD. Bottle mean squares are not automatically converted into a homogeneity uncertainty; that calculation exists in the random-effects ANOVA branch of the uncertainty module.
21.16.3 Reference Intervals
For the default nonparametric method, alpha = (1 - p) / 2 and EP28 continuous ranks are
\[ h_q=q(n+1), \]
with linear interpolation between the order statistics at floor(h_q) and ceil(h_q). Reference-limit confidence intervals use binomial rank inversion. With delta = (1 - gamma)/2 and B ~ Bin(n, q), endpoint ranks are
\[ a=Q_B(\delta),\qquad b=Q_B(1-\delta)+1, \]
then clipped to [1, n]. Actual coverage is P(a <= B <= b - 1). If the unclipped rank lies outside the sample, the function marks that the requested coverage cannot be fully attained. type = “type6” retains the R type-6 point estimate but its CI still uses the existing binomial-rank normal approximation; it should not be called the exact-rank CI of the default nonparametric method.
The parametric method first tests the raw values. If p < 0.05 and all values are positive, it tests log values and uses the log scale only when the log-test p is at least 0.05. On the selected scale, limits are
\[\bar x\pm z_{1-\alpha}s.\]
Endpoint SE is
\[SE_L=s\sqrt{\frac1n+\frac{z_{1-\alpha}^2}{2n}}, \]
with t_(1 - gamma/2, max(n - 2, 1)) used for the endpoint interval; log-scale limits are exponentiated at the end.
biweight implements EP28 Appendix B. Initialize T_bi with the median and s = MAD / 0.6745. At each iteration use
\[ u_i=\frac{x_i-T_{bi}}{3.7s},\qquad w_i=(1-u_i^2)^2 I(|u_i|<1),\qquad T_{bi}\leftarrow\frac{\sum_iw_ix_i}{\sum_iw_i}. \]
Stop when relative change is below 0.001%. For tuning constant c,
\[ A=\sum_{|u_i|<1}(1-u_i^2)(1-5u_i^2),\qquad B=\sum_{|u_i|<1}u_i^2(1-u_i^2)^4, \]
\[ s_{bi}[c]=cs\sqrt{\frac{nB}{A\max(1,-1+A)}}. \]
Calculate s_bi[205.6] and s_bi[3.7] separately. S_T[3.7] uses u_i = (x_i - T_bi)/(3.7 s_bi[3.7]) and removes n from the numerator under the square root. The final limits are
\[ T_{bi}\pm t_{n-1,1-\alpha} \sqrt{s_{bi}^2[205.6]+S_T^2[3.7]}. \]
Endpoint intervals use percentile bootstrap from the same biweight estimator.
huber retains the package extension. Initialize location T with the median and scale S with MAD / 0.6745, then iterate
\[ u_i=\frac{x_i-T}{S},\qquad w_i=\begin{cases} 1,&|u_i|\le1.5,\\ 1.5/|u_i|,&|u_i|>1.5, \end{cases} \]
\[ T\leftarrow\frac{\sum_iw_ix_i}{\sum_iw_i},\qquad S\leftarrow\sqrt{\frac{\sum_iw_i(x_i-T)^2}{n-1}}. \]
The tolerance is 1e-8 with at most 100 iterations; limits are T +/- z_(1 - alpha) S. biweight and Huber endpoint intervals use 1,999 bootstrap draws by default; seed fixes the calculation without changing the caller’s existing random-number state.
21.17 Response-Equation Fitting and Linearity
21.17.1 Twelve Registered Response Equations
| ID | Equation in the source |
|---|---|
| E01 | a + b x |
| E02 | a + b x + c x^2 |
| E03 | a exp(b x) |
| E04 | a exp(-b x) |
| E05 | a x^b |
| E06 | a + b log x |
| E07 | A + (B - A) / [1 + (x/C)^D] |
| E08 | A + (B - A) / [1 + (x/C)^D]^E |
| E09 | a exp[-(x - b)^2 / (2 c^2)] |
| E10 | V_max x / (K_m + x) |
| E11 | a sin(bx + c) + d |
| E12 | a + b x + c x^2 + d x^3 |
E01, E02, E06, and E12 use a linear-model matrix; the remaining equations use Levenberg–Marquardt nonlinear least squares. The objective is
\[RSS_w=\sum_iw_i[y_i-f(x_i,\theta)]^2.\]
When only parameter bounds are supplied, the source first attempts bounded nonlinear least squares. General inequalities use COBYLA and equalities use an augmented Lagrangian. A failed constrained optimization receives the penalty objective
\[RSS_w+10^8\sum_r\{\text{constraint violation}_r\}^2.\]
For constrained fits, covariance uses the central-difference residual Jacobian J:
\[ \widehat V(\hat\theta) \approx(J^\mathsf TJ)^{-1}\frac{RSS}{n-p}. \]
This Jacobian covariance does not multiply J by the fitting weights. Constrained-model AIC is n log(RSS/n) + 2p, and R^2 = 1 - RSS/TSS.
Nonlinear prediction uses a central-difference gradient g(x):
\[V\{\hat f(x)\}=g(x)^\mathsf T\widehat V(\hat\theta)g(x).\]
A single-result prediction interval adds RSS/(n - p) to this variance.
The source uses analytic inverses for E01–E07. E09 takes the left Gaussian branch, E11 uses the principal asin branch, and E12 takes real polynomial roots and selects the root closest to the observed X mean. E08 currently uses
\[ x=C\left(\frac{B-A}{y-A}-1\right)^{1/(DE)}, \]
which is not a strict algebraic inverse of the current E08 equation when E is not 1. Other cases search an extended grid for sign changes, solve one-dimensional roots, and select the root closest to the observed X center.
21.17.2 Linearity Design Quantities and Regression Comparison
Expected concentration in a dilution panel is C(p) = p C_H + (1 - p) C_L. Replicate planning first sets p_L = (1 - alpha_F)^(1/L) and z = z_(1 + p_L)/2:
\[ n=\left\lceil \left(\frac{z\sigma}{ADL-d_{\mathrm{true}}}\right)^2 \right\rceil. \]
Linearity first summarizes each level with its mean, SD, and n_j. Equal weighting uses 1; level-variance weighting uses 1/s_j^2 for raw observations and n_j/s_j^2 for level means. pooled requires variance_group to define contiguous groups in x. For group g,
\[ \widehat\sigma_g^2= \frac{\sum_{j\in g}(n_j-1)s_j^2}{\sum_{j\in g}(n_j-1)}, \qquad w_j=\frac{n_j}{\widehat\sigma_g^2}. \]
The function does not form groups automatically because EP06 grouping also depends on the repeatability-SD versus concentration relationship. A precision-profile weight first fits SD against mean and then uses 1 / estimated SD_j^2.
The line fits level means; quadratic and cubic polynomials retain all observations. If any nonlinear term is significant, the model with the smallest s_y|x among models containing significant nonlinear terms is selected; otherwise the line is selected. Nonlinear deviation is
\[DL_j=\hat y_{\mathrm{poly},j}-\hat y_{\mathrm{linear},j},\]
and linear residual is y_bar_j - y_hat_linear,j. Percentage metrics can use the linear prediction, X, or the observed mean as denominator.
Each level’s point interval uses family coverage gamma^(1/L):
\[ h_j=z_{[1+\gamma^{1/L}]/2}\frac{\hat\sigma_j}{\sqrt{n_j}}. \]
The source centers the interval on the deviation. Allowable deviation is the maximum of direct absolute limits and absolute limits converted from relative limits; the final comparison is |metric_j| > ADL_j.
21.18 Interference, Dilution, Spike Recovery, and Hook Effect
21.18.1 Two-Group Interference and Difference-in-Differences
For a two-group interference study,
\[D=\bar y_T-\bar y_C,\qquad SE(D)=\sqrt{\frac{s_T^2}{n_T}+\frac{s_C^2}{n_C}}.\]
The interval uses a t critical value with n_T + n_C - 2 degrees of freedom, even though the standard error is unpooled. Percentage effect is 100 D / y_bar_C, with the SE and interval endpoints scaled directly.
For an analyte-spike design,
\[ D=(\bar y_{I,S}-\bar y_{I,U}) -(\bar y_{N,S}-\bar y_{N,U}), \]
with variance equal to the sum of s_g^2 / n_g across four cells and Welch–Satterthwaite degrees of freedom. Percentage effect is 100(I/N - 1), where I = y_bar_I,S - y_bar_I,U and N = y_bar_N,S - y_bar_N,U. The source calculates the variance by differentiating this expression with respect to the four means and applying the delta method.
When at least two control results exist, blank interference still uses the two-group difference; otherwise the theoretical control mean is zero and a one-sample standard error is used. Dose response calculates effects at each concentration relative to the observed baseline and can additionally fit an OLS line of raw result versus concentration. Forward and reverse conversions use linear interpolation between points and do not extrapolate; reverse interpolation requires a monotone response.
21.18.2 Dilution Recovery
For dilution factor D, the source sets f = 1 / D. With original assignment C_0 and diluent background C_d,
\[ C_{\mathrm{target}}=fC_0+(1-f)C_d, \]
\[ C_{\mathrm{corrected}}=\frac{\bar y-(1-f)C_d}{f}, \qquad Recovery=100\frac{\bar y}{C_{\mathrm{target}}}. \]
Recovery is measured diluted result divided by target; it is not corrected concentration divided by original assignment. For multiple samples, the overall interval is the t interval of sample-level recovery means.
21.18.3 Spike Recovery and Hook Effect
\[ C_{\mathrm{added}}=C_s\frac{V_s}{V_0+V_s}, \qquad C_{\mathrm{final}}=\frac{C_0V_0+C_sV_s}{V_0+V_s}, \]
\[ Recovery=100\frac{\bar y_{\mathrm{spike}}-\bar y_{\mathrm{solvent}}} {C_{\mathrm{added}}}. \]
Hook-effect analysis calculates a mean and t interval at each concentration and identifies the concentration with the maximum mean as the peak. After the peak, when the first level has concentration above ULoQ and a mean response no higher than the threshold, the source reports the preceding measured concentration as hook concentration; it does not interpolate between levels. If no decline is observed, it reports the highest tested concentration as a lower bound.
21.19 Analytical Sensitivity
21.19.1 Limits of Blank and Detection
The nonparametric LoB uses rank r = 0.5 + B(1 - alpha) and linearly interpolates between adjacent integer ranks. The parametric method sets nu = B - K and calculates
\[ c_p=\frac{z_{1-\alpha}}{\sqrt{1-1/(4\nu)}}, \qquad LoB=\bar y_B+c_ps_B. \]
The zero-blank branch sets LoB = 0 directly. With two or three reagent lots, the default calculates each lot and takes the maximum LoB; with other lot counts it pools the lots.
Classical LoD pools low-level SDs by degrees of freedom:
\[ s_p=\sqrt{\frac{\sum_j(n_j-1)s_j^2}{\sum_j(n_j-1)}}, \qquad LoD=LoB+c_ps_p. \]
The nonparametric validation method uses the median of all low-level results as a candidate LoD and retains it only when the fraction below LoB is no more than beta. The precision-profile method fits s(x) and solves x - LoB - c_p s(x) = 0, using the smallest qualifying root found.
The probit method fits
\[ \Phi^{-1}\{P(\mathrm{detected}\mid x)\}=\beta_0+\beta_1g(x), \]
where g(x) is log10(x) or x. For target detection probability 1 - beta,
\[g(LoD)=\frac{z_{1-\beta}-\hat\beta_0}{\hat\beta_1}.\]
When fitting by lot or stratum, the final LoD is the maximum group-specific LoD.
21.19.2 Limit of Quantitation and Limit Verification
For a low-level sample, bias is b = y_bar - x_a. The implemented error quantities are
\[ E_W=|b|+ks,\qquad E_{RMS}=\sqrt{s^2+b^2},\qquad E_P=s. \]
For a relative target, divide by abs(x_a). The observed candidate LoQ is the smallest measured mean whose error reaches the target. With a profile enabled, the candidate is the smallest root of a linear or quadratic regression of error against assigned value and the target line. The maximum candidate across lots is reported, and when both are finite
\[LoQ_{\mathrm{reported}}=\max(LoQ_{\mathrm{candidate}},LoD).\]
LoB verification calls y <= LoB consistent; LoD verification calls y >= LoB consistent; LoQ verification calls a result consistent when it lies within the allowed window around its assigned value. The source calculates one-sided Wilson intervals and uses an internal sample-size boundary table; a sample size not in the table takes the larger boundary from neighboring rows. It returns proportions, intervals, and differences from boundaries but does not generate a pass/fail conclusion.
21.20 Quality Control and Stability
21.20.1 Levey–Jennings and Westgard
QC results are standardized as z_i = (x_i - mu_0) / sigma_0. If target values are not supplied, the source uses the sample mean and sample SD of all finite results. The rules are:
| Rule | Source decision |
|---|---|
| 1-2s | abs(z_i) > 2 |
| 1-3s | abs(z_i) > 3 |
| 2-2s | two adjacent points on the same side, both strictly beyond 2s |
| R-4s | adjacent points each with abs(z) >= 2 and abs(delta z) >= 4 |
| 4-1s | four consecutive points on the same side, all strictly beyond 1s |
| 8x/10x | eight/ten consecutive points strictly on the same side of the mean |
The grouping variable is used only for plot facets; target estimation and rule scanning run across all finite results in row order. 1-3s is labeled out of control and other violations are warnings. Youden uses abs(z_1i) <= 2 and abs(z_2i) <= 2 to test membership in the +/-2s rectangle and calculates Pearson correlation; it does not fit an error ellipse.
21.20.2 Node Bias, Time Regression, and Time to Limit
Each condition uses its earliest time-point mean as baseline:
\[b_t=y_t-y_0,\qquad b_{t,\%}=100\frac{y_t-y_0}{y_0}.\]
Self-comparison first averages by time point. If all time-point SEs are finite and positive, it fits with 1 / SE_t^2 weights; otherwise it falls back to OLS. Absolute drift rate is b_1, relative drift rate is b_1 / b_0, with variance
\[ V\left(\frac{b_1}{b_0}\right)= \frac{V(b_1)}{b_0^2}+\frac{b_1^2V(b_0)}{b_0^4} -\frac{2b_1\operatorname{Cov}(b_0,b_1)}{b_0^3}. \]
The source uses a one-sided t critical value t_(gamma, nu). For a decline it takes the lower bound and for an increase the upper bound; if direction is uncertain, it takes the larger absolute side.
Condition comparison uses test-condition minus reference-condition bias at matched time points. Absolute-bias SE is sqrt(SE_T^2 + SE_R^2); relative-bias SE uses the ratio delta method. WLS or OLS is then selected according to SE availability.
Time to limit searches a 501-point grid for the first sign change of abs(b_hat(t)) - L or of a one-sided confidence bound versus L, then refines it with one-dimensional root finding. The default upper search is max(100 t_max, t_max + 100); a root beyond the largest observed time is marked extrapolated.
21.20.3 Mean Kinetic Temperature and Arrhenius
Let q_i = exp[-E_a / (R T_i)]. For equally spaced times,
\[ \bar q=\frac1n\sum_iq_i,\qquad T_{MKT}=-\frac{E_a}{R\log\bar q}. \]
With timestamps, trapezoidal integration is used:
\[ \bar q=\frac{\sum_i\Delta t_i(q_i+q_{i+1})/2}{\sum_i\Delta t_i}. \]
Multiple probes combine q by effective observation count or time duration, not by directly averaging probe-specific MKT values.
The Arrhenius branch fits zero-order y = a + b t or first-order log(y) = a + b t at each temperature and converts the slope to a positive rate according to the degradation direction. Automatic order selects the smaller root-mean-square normalized RMSE across temperatures; exact ties select first order. Positive rates then fit
\[\log k=\alpha+\beta/T,\qquad E_a=-\hat\beta R.\]
Zero-order shelf life is allowable change divided by the target-temperature rate. First-order relative decreases and increases are
\[t=-\frac{\log(1-r)}{k_*},\qquad t=\frac{\log(1+r)}{k_*}.\]
Shelf-life intervals convert both endpoints of the target-rate interval and sort them.
The stability design is not a continuous power formula; it uses two lookup matrices embedded in the source. Repeatability is reduced by sigma / sqrt(r), or combined as sqrt(sum_j sigma_j^2 / L_j). The resulting two ratios are rounded to two significant figures and looked up in the 80% or 90% power table.
21.21 Measurement Uncertainty and Traceability
21.21.1 Type A and Type B Standardization
The standard uncertainty of one repeat result is s; the standard uncertainty of a mean is s / sqrt(n), with n - 1 degrees of freedom. Relative components are stored internally as u / abs(x_bar).
For interval-based Type B components, the estimate must be the midpoint of the limits. With half-width a = (U - L) / 2:
| Distribution | Standard uncertainty |
|---|---|
| rectangular | a / sqrt(3) |
| triangular | a / sqrt(6) |
| U-shaped | a / sqrt(2) |
| normal with central probability p | a / z_(1+p)/2 |
| t with central probability p | a / t_(1+p)/2,nu |
A directly supplied uncertainty is treated as a standard uncertainty. Only when coverage factor k is separately supplied does the source calculate u = U / k. This path accepts standard or normal labels and does not divide by sqrt(3), sqrt(6), and so on according to a distribution label. Distribution divisors occur only on the limit-input path.
21.21.2 IQC, Homogeneity, Stability, and Characterization
For a one-way random-effects ANOVA,
\[ MS_W=SS_W/(N-m),\qquad MS_B=SS_B/(m-1), \]
\[ n_{\mathrm{eff}}= \frac{N-\sum_jn_j^2/N}{m-1},\qquad s_B^2=\max\left(0,\frac{MS_B-MS_W}{n_{\mathrm{eff}}}\right). \]
When IQC has within-run replicates, single-result uncertainty is
\[u_{\mathrm{single}}=\sqrt{s_W^2+s_B^2}.\]
If s_B > 0, mean uncertainty is sqrt(MS_B / N); otherwise the source falls back to the SD of all results divided by sqrt(N). Without a run column or with one result per run, it uses long-term SD and its mean standard error directly.
The homogeneity branch also calculates
\[ u_{\mathrm{hidden}}= \frac{s_W}{\sqrt{\bar n}} \left(\frac2{\nu_W}\right)^{1/4}, \qquad u_{bb}=\max(s_B,u_{\mathrm{hidden}}). \]
Long-term stability uncertainty fits y = a + b t for each condition and takes u_lts = SE(b_hat) T_s at shelf life T_s. Characterization first averages each laboratory, then assigns the unweighted mean of laboratory means:
\[u_{char}=\frac{SD(\bar x_1,\ldots,\bar x_L)}{\sqrt L}.\]
21.21.3 First-Order Combination and Effective Degrees of Freedom
For standard-uncertainty vector u and correlation matrix R,
\[ \Sigma=\operatorname{diag}(u)R\operatorname{diag}(u),\qquad u_c^2=c^\mathsf T\Sigma c. \]
Signed marginal contribution from component i is m_i = c_i (Sigma c)_i, and the contribution percentage is 100 m_i / u_c^2; with correlation, a contribution can be negative.
Only when R is the identity does the source calculate
\[ \nu_{\mathrm{eff}}= \frac{u_c^4}{\sum_i(c_iu_i)^4/\nu_i}. \]
Without a supplied coverage factor, finite effective degrees of freedom use a t quantile; otherwise a normal quantile is used. Expanded uncertainty is U = k u_c. When relative and absolute components are mixed, an output value must be supplied so relative components can be converted to absolute units.
21.21.4 General Measurement Models and Monte Carlo
For y = f(x_1, …, x_q), sensitivity coefficients use central differences:
\[ c_i\approx \frac{f(\ldots,x_i+h_i,\ldots)-f(\ldots,x_i-h_i,\ldots)}{2h_i}, \]
\[h_i=\epsilon^{1/3}\max\{1,|x_i|,u(x_i)\}.\]
The delta branch propagates coefficients through the covariance matrix and returns a symmetric interval y_hat +/- U.
The Monte Carlo branch samples from stored distributions: normal and standard use N(x_i, u_i^2); rectangular uses a uniform distribution; triangular uses the limits and estimate as mode; U-shaped uses x_i + a_i cos(U(0, pi)); t uses x_i + u_i T_(nu_i), so u_i is the t random-variable scale. Non-independent Monte Carlo accepts only normal or standard labels and generates correlated normal draws.
Simulated-output SD is the standard uncertainty, and the coverage interval is the central interval from R type-8 quantiles.
21.21.5 Stepwise Traceability
At stage r, all components from stages 1 through r are recombined. Sensitivity coefficients default to 1 and this interface has no between-stage correlation input, so
\[u_{c,r}=\sqrt{\sum_{i\in1:r}u_i^2}.\]
Effective degrees of freedom, coverage factor, and expanded uncertainty are recalculated at every stage; the final stage is the budget for all accumulated components.
21.22 Sample Size and Power
21.22.1 Bland–Altman Agreement
For difference mean mu, SD sigma, allowable limit Delta, and agreement quantile z_gamma,
\[ SE_L=\sigma\sqrt{\frac1n+\frac{z_\gamma^2}{2(n-1)}}, \]
\[ \tau_1=\frac{\Delta-\mu-z_\gamma\sigma}{SE_L},\qquad \tau_2=\frac{\Delta+\mu-z_\gamma\sigma}{SE_L}. \]
With t_c = t_(1 - (1 - C)/2, n - 1), power is
\[ \Pr\{T_{n-1}(\tau_1)>t_c\} +\Pr\{T_{n-1}(\tau_2)<-t_c\}. \]
Sample-size search uses a binary search within the stated bounds for the smallest integer reaching the target power.
21.22.2 Proportion Intervals and Target Tests
The closed-form Wald half-width sample size is
\[n=\left\lceil\frac{z_{1-\alpha/2}^2p(1-p)}{d^2}\right\rceil.\]
Other implemented intervals calculate the maximum actual endpoint half-width, then use exponential and binary search for the smallest n with half-width <= d. Methods involving x = round(np) can have nonmonotone jumps; the source scans downward to correct the result.
For a one-sample proportion target test, first find the critical success count. For an upper alternative it is the smallest integer whose selected one-sided interval lower limit L(x, n) is greater than p0; for a lower alternative it is the largest integer whose upper limit U(x, n) is less than p0. Therefore
\[ Power_{>}=1-F_{\mathrm{Bin}(n,p_1)}(x_c-1),\qquad Power_{<}=F_{\mathrm{Bin}(n,p_1)}(x_c). \]
A two-sided test allocates alpha / 2 to each tail. For sample-size search, a normal-approximation starting point is followed by upward search and downward backtracking; for alpha search, the discrete critical count is fixed first and alpha is solved so that the interval endpoint equals p0. Power is therefore a step function determined by the binomial critical count.