16 Introduction to Algorithm Principles
16.1 Introduction
This document summarizes the algorithm principles underlying ivdtools, to support understanding of the statistical basis of each analysis module and the mapping with the package interface. It covers outliers/distribution, method comparison/qualitative evaluation, ROC, precision/ANOVA/reference intervals, and QC/stability/equation fitting. For ease of cross-reference, method comparison uniformly adopts the convention of candidate \(X\), reference \(Y\), and regression model \(Y=a+bX\); the specific parameter names and directions in the package are governed by the registries such as list_equation(), list_sadler(), and list_westgard() in ivdtools.
The formulas are used to illustrate the principles; they do not replace textbook derivations, nor do they replace the judgment about method selection in the protocol.
16.2 General Conventions and Notation
Ordered samples \(x_{(1)}\le\cdots\le x_{(n)}\), mean \(\bar x\), standard deviation \(s\), significance level \(\alpha=1-\text{level}\). Most tests require approximately independent samples; a detected “outlier” or “failure to follow some distribution” only indicates the strength of conflict with the assumed model, cannot quantify its actual effect on the final performance estimate, and does not automatically trigger deletion.
16.3 Outlier and Distribution Algorithms
16.3.1 Outliers
- Grubbs: \(G=\max_i|x_i-\bar x|/s\), judged against a critical value derived from the t distribution for the most extreme point; approximately normal and applies to one point at a time. Mapping:
outliers_test(method="grubbs", alpha=...). - Generalized ESD: successively removes the most extreme point, computes \(R_i\) and compares it with \(\lambda_i\) that changes with the remaining n; the maximum number of outliers is limited by an additional parameter.
- Dixon: uses the ratio of the endpoint gap to the sample range; the critical value depends on n and is suitable only for a small-sample range.
- IQR: marks values less than \(Q_1-k\,IQR\) or greater than \(Q_3+k\,IQR\); the default coefficient is determined by the interface implementation; this is a rule rather than a significance test.
outliers_test() first excludes missing values and retains the original row numbers, then returns the statistic, critical value, and outlier flag. A constant column makes statistics with s or the range as the denominator undefined; multiple testing increases false positives.
16.3.2 Normality
- Shapiro–Wilk: \(W=(\sum a_i x_{(i)})^2/\sum(x_i-\bar x)^2\), with weights from the normal order statistics.
- Anderson–Darling: the weighted squared distance between the empirical distribution and the fitted normal CDF, with greater weight on the tails.
- Lilliefors: a Kolmogorov–Smirnov-type statistic when parameters are estimated from the sample.
- Cramér–von Mises: the overall integrated squared distance between the empirical CDF and the theoretical CDF.
normal_test(method="auto") selects an available test based on the sample conditions (Shapiro–Wilk for n≤50, Anderson–Darling for larger samples), and level maps to the rejection threshold. Graphics provide Q-Q and histogram; the p-value cannot quantify the actual effect of the deviation on the performance estimate.
16.4 Method Comparison and Qualitative Evaluation Algorithms
16.4.1 Regression, Correlation, and Bias
The unified direction is candidate \(X\), reference \(Y\), model \(Y=a+bX\).
- OLS minimizes \(\sum(y_i-a-bx_i)^2\); assumes the X error is negligible.
- Deming minimizes the weighted orthogonal distance to the line; the error variance ratio maps to
lambda; the direction of \(\lambda=\sigma_Y^2/\sigma_X^2\) must be checked against the source code or protocol. - Passing–Bablok estimates b from the median of all pairwise slopes, then estimates a from the median of \(y_i-bx_i\); it is robust to outliers but still requires linearity and independent samples.
- Pearson measures linear correlation; Spearman/Kendall measure rank association;
correlation(method=...)answers “correlation” rather than “agreement” and cannot replace bias and agreement-limit evaluation.
The absolute Bland–Altman difference is recorded as \(d_i=X_i-Y_i\) (the report must re-check the actual direction in the package object), with mean difference \(\bar d\) and agreement limits \(\bar d\pm z_{(1+p)/2}s_d\); relative BA uses percentage differences. agree.level controls the coverage proportion and conf.level controls the estimation interval. bias() maps the regression equation to the specified level and outputs the predicted bias and confidence/prediction intervals.
16.4.2 Four-Fold Table
From TP, FP, TN, FN:
\[Se=\frac{TP}{TP+FN},\quad Sp=\frac{TN}{TN+FP},\quad PPV=\frac{TP}{TP+FP}.\]
The Wilson interval is inverted from a score test; the exact interval is inverted from binomial tail probabilities. Cohen \(\kappa=(P_o-P_e)/(1-P_e)\), where \(P_e\) is obtained from the marginal proportions. McNemar uses only the discordant cells b and c, with the asymptotic statistic approximately \((b-c)^2/(b+c)\) (with the continuity correction or exact version as recorded in the object). The package interfaces map to diagnostics(), kappa(), and mcnemar(); zero denominators, small cell counts, and extreme prevalence all require explicit explanation.
16.5 ROC Algorithms
For each threshold t, compute \(TPR(t)=Se(t)\) and \(FPR(t)=1-Sp(t)\) according to the marker direction. The empirical ROC treats all observations as candidate thresholds and connects them as a step/broken line; the AUC is summed by the trapezoidal rule over adjacent points:
\[AUC=\sum_i \frac{TPR_i+TPR_{i+1}}{2}(FPR_{i+1}-FPR_i).\]
auc(roc_obj, cols=...) saves the per-marker AUC table; the empirical AUC can be interpreted as the probability that a random positive score exceeds a random negative score (after direction adjustment). Sample-correlated multi-curve comparison, partial AUC, or complex sampling are not automatically resolved by this point estimate.
The Youden index \(J(t)=Se(t)+Sp(t)-1\); the closest-corner method minimizes \(\sqrt{(1-Se)^2+(1-Sp)^2}\). cutoff() outputs both types of optimal points. Tied thresholds, measurement resolution, misclassification costs, and clinical prevalence change the actual choice, so “mathematically optimal” is not an automatic clinical cutoff.
The multi-marker model is
\[\operatorname{logit}\{P(Y=1)\}=\beta_0+\sum_j\beta_jx_j,\]
fitted by binomial logistic regression via mlr(cols, name); predict() transforms the linear predictor into a probability through the logistic inverse, classifies at the default 0.5, and computes intervals using the model covariance. Perfect separation, collinearity, missingness, and small event counts lead to instability. Re-drawing the ROC on the training set is an apparent performance; development and validation data must be distinguished.
16.6 Precision, ANOVA, and Reference Interval Algorithms
16.6.1 Variance Components
A nested design can be written as \(y_{ijkl}=\mu+D_i+R_{j(i)}+\varepsilon_{k(ij)}\), with total variance \(\sigma_T^2=\sigma_D^2+\sigma_R^2+\sigma_e^2\). precision(form = y ~ day/run) passes the formula to VCA::anovaVCA(); variance()/vc() summarize the VC, \(SD=\sqrt{VC}\), \(CV=100SD/\bar y\), and proportions. Satterthwaite uses
\[\nu\approx\frac{(\sum c_jMS_j)^2}{\sum(c_jMS_j)^2/\nu_j}\]
to approximate the effective degrees of freedom, then constructs variance/SD intervals from chi-square quantiles. Negative components, low degrees of freedom, and imbalance degrade the approximation.
The Sadler precision profile expresses the variance or SD as a function of concentration; list_sadler() lists 10 candidate forms, and profile(model.no=...) calls the VFP fit and selects by AIC. Different models extrapolate very differently at low concentrations; the concentration coverage and residuals must be checked.
16.6.2 ANOVA and Reference Intervals
bottle_anova() uses stats::aov() to decompose the between-group/within-group sums of squares; Bartlett tests variance homogeneity, the residual Shapiro–Wilk checks normality; Tukey HSD uses the studentized range to control the family-wise error rate of all pairwise comparisons at once.
The parametric reference interval is approximated as \(\bar x\pm z_{1-\alpha/2}s\) under the normality assumption; the nonparametric percentile method uses R quantile type 6, with endpoints from the order statistics and endpoint intervals computed by binomial ordering or bootstrap paths. interval is the coverage and ci is the endpoint confidence; the two have different meanings. Sample selection bias cannot be repaired by any interval algorithm.
16.7 QC, Stability, and Equation Fitting Algorithms
16.7.1 QC and Stability
The Levey–Jennings standardization is \(z_i=(x_i-\mu_0)/\sigma_0\); the Westgard rules combine these z values across single points, consecutive same-side/same-direction runs, or multiple levels in the same batch. qc_chart() evaluates in run order, and youden_plot() uses the target means/SDs of two control materials to form quadrant and ellipse-style judgments; the order, target values, and lot boundaries are inputs to the algorithm.
The node relative bias is \(100(x_t-x_0)/x_0\) (relative) or \(x_t-x_0\) (absolute). stability_regression() fits a time trend, and stability_time() finds the intersection of the confidence boundary/trend with the allowable limit. MKT is the Arrhenius-weighted temperature history:
\[T_k=-\frac{E_a}{R}\left[\log\left\{\sum_i w_i\exp(-E_a/RT_i)\right\}\right]^{-1}.\]
The Arrhenius relationship is \(k=Ae^{-E_a/(RT)}\); zero-order is \(C_t=C_0-kt\), first-order is \(\log C_t=\log C_0-kt\). mkt(ea=83.144), arrhenius(order, target_temp, limit) map these quantities. Temperature must be converted to Kelvin for the formulas, and extrapolation is limited by the assumption of an unchanged reaction mechanism.
16.7.2 Equation Fitting
Linear/polynomial are estimated by least squares; exponential and 4PL/5PL use nlsLM/constrained optimization. A typical 4PL is
\[y=d+\frac{a-d}{1+(x/c)^b},\]
and 5PL adds an asymmetry exponent g, but the specific parameter names and directions in the package are governed by list_equation(). Weighted fitting minimizes \(\sum w_i(y_i-f(x_i,\theta))^2\). AIC approximates \(-2\log L+2k\); residuals reveal heteroscedasticity, curvature, and anomalous points.
The mean confidence interval reflects the uncertainty of \(E(Y|X)\), while the prediction interval also includes single-observation error. Nonlinear inverse prediction requires numerical root-finding; small signal errors near a flat plateau can translate into enormous concentration errors. Starting values, bounds, weights, equation versions, convergence codes, and failed models should all enter the audit record.
16.8 Module-to-Guideline Mapping
| Topic | Main functions | Corresponding CLSI guideline |
|---|---|---|
| Outliers and normality | outliers_test(), normal_test() |
General (pre-analysis steps of the EP guidelines) |
| Method comparison | mcr() series |
EP09 |
| Qualitative evaluation | fourfold_table series |
EP12 |
| ROC | roc() series |
EP24 (reference) |
| Precision | precision() series |
EP05 |
| Bottle/lot ANOVA | bottle_anova() |
EP15 |
| Reference intervals | reference_interval() |
EP28-A3c |
| Stability | stability_*(), arrhenius(), mkt() |
EP25 |
| QC | qc_chart(), youden_plot() |
Westgard rules |
| Equation fitting | fit_equation() series |
EP06 |
| Analytical specificity | lob_lod_loq() |
EP17 |
16.9 Scope of Use
- Formulas and statistics are valid only under reasonable assumptions; when assumptions are not met, the results cannot be applied mechanically.
- “Statistically significant” and “mathematically optimal” are not the same as “clinically/product acceptable”; acceptance criteria must be defined before analysis.
- Performance on training data (such as ROC AUC, regression goodness of fit) is apparent performance; it should be validated on independent data before generalization.
- All computations should be conducted under the guidance of the protocol, standards, and professional judgment, and reviewed by appropriate personnel.