21  算法原理简介

本文只描述 ivdtools 当前源码实际执行的计算。涉及同一函数的多个算法参数时,本章同时列出其数学定义、默认选择和适用条件,以便判断改变参数究竟改变了估计对象、误差模型,还是只改变了区间算法;具体调用和报告流程仍见相应专题章节。

以下统一记样本量为 \(n\),样本均值为 \(\bar x\),样本标准差为 \(s\),正态分布和 t 分布的 \(p\) 分位数分别为 \(z_p\) 和 \(t_{p,\nu}\)。除特别说明外,源码使用 R 的双精度未舍入数值完成中间计算。

21.1 重复测量、多样本与统计单位

21.1.1 先确定“哪一行是一个观测”

ivdtools 的函数并不把“长表”自动等同于“支持重复测量”。在调用函数前先确定统计单位:

  • 单次测量:每行本身就是独立测量,例如精密度方差分量、线性各水平的重复结果、稀释回收的重复结果;
  • 样本级代表值:每个受试者或样本只保留一行,例如方法学比较、普通 ROC、参考区间和四格表;
  • 样本内重复结构:重复次数、日、运行、批次或中心具有模型含义,例如 precision(result ~ day/run);
  • 分组后的独立分析:不同产品、分析物、批次、参考人群分区或质控品分别形成多个数据集。

如果函数没有明确的重复结构参数,同一个 ID 的多行不会自动变成相关观测,也不会自动计算均值;这些行可能直接增加 n、进入回归或被计入四格表。因此,重复数据的处理顺序应为:

  1. 依据方案确认重复测量的层级、有效重复数和最终统计单位;
  2. 若目标函数需要样本级一行,按预先规定的均值、中位数、首次结果或最终判定规则汇总;
  3. 若目标函数本身估计重复性或方差分量,保留长表并提供日、运行、批次、中心等分层列;
  4. 若是多个独立项目,先按分类变量或连续界值拆分,再对每个数据集执行同一分析流程;
  5. 同时保存原始行、汇总规则、实际有效重复数和拆分记录,避免把测量次数误报为独立样本数。

21.1.2 可以直接接收重复结果的函数

precision() 将重复结果作为方差模型的观测,按公式中的嵌套层级估计残差、日内/日间、批内/批间或站点等方差分量;by 只负责在外层分组后分别建立模型,不改变组内方差公式。linearity() 先按浓度水平计算重复均值、SD 和有效重复数,再按 weights 形成水平权重;method = "linear" 的最终直线和偏离分析均使用这些水平汇总量,探索性多项式则按当前接口的观测层级拟合。c5_c95() 同样先按水平得到 \((mu_j,s_j,n_j)\),再拟合 SD—均值关系,不能把二分类阳性数伪装成连续响应。

interference_paired()、 interference_dose_response()、 interference_patient()、 dilution_recovery()、 spike_recovery()、 hook_effect()、 commutability()、 reference_bias() 和 sensitivity_*() 的共同原则是:在函数明确的样本、水平、方法、条件或批次键内先汇总重复结果,再计算效应、区间或研究级汇总。它们返回的样本级槽位与总体槽位不是同一个统计对象,报告时应分别说明。

21.1.3 必须先归并为样本级结果的函数

普通 mcr()、roc()、auc()、 cutoff()、 mlr()、 raw_to_table()、 kappa()、 mcnemar()、 reference_interval()、 outliers_test() 和 normal_test() 不会为普通数据框建立样本内相关结构。方法学比较应先把候选和参考方法的重复结果归并为每个样本一对结果;ROC 和四格表应先形成每名受试者一个结果或一个最终判定;参考区间和普通分布检验应先形成每名参考个体一个值。

对方法学比较,若第 \(i\) 个样本两侧分别有 \(r_{Xi}\)、\(r_{Yi}\) 次重复,至少应保存

\[ \bar x_i=\frac1{r_{Xi}}\sum_{k=1}^{r_{Xi}}x_{ik},\qquad \bar y_i=\frac1{r_{Yi}}\sum_{k=1}^{r_{Yi}}y_{ik}, \]

以及两侧 SD 和重复次数。回归和 bootstrap 以 \((\bar x_i,\bar y_i)\) 为抽样单位;重复 SD 只在 lambda = "replicates" 或精密度轮廓权重中用于误差模型,不把重复行再加入回归样本量。

21.1.4 多样本是“共同分析”还是“分别分析”

多份患者样本共同估计一条方法学回归线,多名受试者共同形成一条 ROC 曲线,多名参考个体共同估计参考区间,多瓶共同形成瓶间 ANOVA;这些场景不能按样本逐个拆开后再把结果平均。相反,不同分析物、产品、批次、参考人群分区、质控品、独立 panel 或方法组合若要分别报告,必须先形成多个数据集。

连续变量分段使用 split_by_cutpoints(),得到 [−Inf,c_1)、[c_1,c_2)、…、[c_K,Inf) 的具名 dataframe list;分类列使用 split_by_factors(),按 factor levels 或非缺失值首次出现顺序拆分。拆分函数只改变数据范围,不改变下游统计公式:

代码
segments <- split_by_cutpoints(
  method_data,
  variable = "reference",
  cutpoints = c(100, 500)
)

results <- lapply(
  segments,
  function(d) {
    mcr(d, id = "sample", candidate = "candidate", reference = "reference") |>
      describe() |>
      correlation() |>
      regression(method = "deming") |>
      bland_altman()
  }
)

这种预处理方式使每个结果对象的样本范围、缺失处理、异常值记录和置信区间都可独立审计;如果需要整体结论,应在原始总体数据上另建对象,而不是把各段结果简单加权平均。

21.2 异常值与正态性

21.2.1 Grubbs 检验

对有限观测计算

\[ G=\frac{\max_i|x_i-\bar x|}{s}. \]

令

\[ t=t_{1-\alpha/(2n),\,n-2},\qquad G_{\mathrm{crit}} =\frac{n-1}{\sqrt n} \sqrt{\frac{t^2}{n-2+t^2}}. \]

当 \(G>G_{\mathrm{crit}}\) 时,只标记本次绝对离均差最大的一个观测。源码没有在标记后自动重复 Grubbs 检验;\(n<3\) 或 \(s=0\) 时不产生异常值。

21.2.2 广义 ESD

第 \(i\) 轮在当前剩余样本中计算

\[ R_i=\frac{\max_j|x_j-\bar x_i|}{s_i}, \]

移除该轮最极端值,并以

\[ t_i=t_{1-\alpha/[2(n-i+1)],\,n-i-1}, \]

\[ \lambda_i= \frac{(n-i)t_i} {\sqrt{(n-i-1+t_i^2)(n-i+1)}} \]

作为临界值。最终异常值数是满足 \(R_i>\lambda_i\) 的最大轮次,而不是逐轮遇到不显著即停止。源码要求 \(n\ge5\);未指定最大异常数时取 \(\min(5,n-2)\)。

21.2.3 Dixon 与 IQR 规则

Dixon 计算两端间隙占极差的比例:

\[ Q_{\min}=\frac{x_{(2)}-x_{(1)}}{x_{(n)}-x_{(1)}},\qquad Q_{\max}=\frac{x_{(n)}-x_{(n-1)}}{x_{(n)}-x_{(1)}}. \]

临界值来自源码内置的 \(n=3,\ldots,30\)、\(\alpha\in\{0.01,0.05\}\) 表。双端检验时源码先检查低端;低端达到临界值后立即返回,不再同时标记高端。

IQR 规则使用 R 默认分位数算法得到 \(Q_1\)、\(Q_3\),并以

\[ [Q_1-k(Q_3-Q_1),\;Q_3+k(Q_3-Q_1)] \]

之外的观测作为异常值,默认 \(k=1.5\)。

21.2.4 正态性检验的自动分支

源码对有限观测要求 \(n\ge3\) 且 \(s>0\)。自动选择规则为:

样本量 实际检验
\(3\le n\le50\) Shapiro–Wilk
\(51\le n\le5000\) Anderson–Darling
\(n>5000\) Lilliefors

Cramér–von Mises 仅在明确选择时执行。后面三种检验直接使用 nortest 包实现;源码只保存检验统计量和 p 值,没有计算偏离正态的效应量,也没有自动进行数据变换。

21.3 方法比较与参考物质偏倚

21.3.1 分析方向、数据单位与预处理分段

以下统一以参考或比较方法为 \(X\)、候选或待评价方法为 \(Y\),并建立

\[Y=a+bX.\]

因此,正偏倚表示候选方法结果高于参考方法。源码只保留两列均为有限值的完整配对,并要求每个样本 ID 只有一行。原始重复测量必须先汇总为每个样本的两种方法均值;如果还要由重复测量估计 Deming 方差比,则同时保留两侧 SD 和重复次数。这样,回归和 bootstrap 的统计单位始终是“样本”,不会把同一样本的重复结果误当成相互独立的临床样本。

若预设界值为 \(c_1<\cdots<c_K\),split_by_cutpoints() 在建立 mcr 对象前把指定数值变量划为

\[ X<c_1,\quad c_1\le X<c_2,\quad\ldots,\quad X\ge c_K. \]

函数返回由 dataframe 构成的具名 list,每个区段分别建立一个 mcr 对象。由此,一个 mcr 对象只对应一个确定的分析数据集,描述、相关、离群点、回归、差值分析及绘图使用相同样本。如需合并区段,应在建立对象前显式合并相应 dataframe,而不是在分析函数内部选择区段。

分类变量预处理使用 split_by_factors()。factor 按声明的 levels 顺序,character 和 logical 按首次出现顺序形成组;每行只属于一个水平。默认排除分组缺失并删除未使用的 factor level,返回值保留原始分组列和各组行数。它只改变后续函数接收的数据集范围,不改变任何统计公式。

21.3.2 相关系数的三个选项

correlation(method=) 只描述共同变化强度,不估计一致性。Pearson 相关为

\[ r=\frac{\sum_i(x_i-\bar x)(y_i-\bar y)} {\sqrt{\sum_i(x_i-\bar x)^2\sum_i(y_i-\bar y)^2}}, \]

检验的是总体线性相关为零。Spearman 先分别把 \(x_i,y_i\) 转换为秩 \(R_i,S_i\),再计算 Pearson 秩相关

\[ \rho_s=\operatorname{Cor}(R,S). \]

没有并列秩时也可写成

\[ \rho_s=1-\frac{6\sum_i(R_i-S_i)^2}{n(n^2-1)}. \]

Kendall 依据所有可比较样本对的一致与不一致方向。无并列值时

\[ \tau=\frac{N_C-N_D}{\binom n2}; \]

实际由 stats::cor.test() 处理并列秩修正。Pearson 针对线性关系且更受极端值影响;Spearman 针对单调关系;Kendall 直接解释为样本对方向一致性的净比例。无论哪一种,高相关都不能证明两方法偏倚小或可以互换,因为 \(Y=a+bX\) 即使 \(a\ne0\) 或 \(b\ne1\) 也可能具有接近 1 的相关系数。

回归 method 的选择首先改变误差模型:

method X 轴误差 Y 轴误差/方差形态 点估计原则 典型用途
"ols" 视为可忽略 恒定或由可选权重描述 最小化纵向平方残差 参考 X 的误差相对很小
"wls" 视为可忽略 随浓度变化 最小化加权纵向平方残差 差值 SD/CV 随浓度变化
"deming" 不忽略 两轴误差比为常数 \(\lambda\) 误差变量回归 两方法都有测量误差且方差形态近似恒定
"wdeming" 不忽略 由共同权重、恒定 CV 或两侧方差轮廓描述 加权或迭代误差变量回归 两轴误差随浓度变化
"pb" 不忽略 不要求正态误差模型 两两斜率的秩估计 线性但有异常值、混合变异或不宜作正态假设

选择回归方法依据的是研究设计和差值图所揭示的误差结构,而不是哪个方法产生更接近 1 的斜率。

21.3.3 OLS 与 WLS 点估计

普通最小二乘最小化

\[ Q(a,b)=\sum_{i=1}^n(y_i-a-bx_i)^2. \]

WLS 最小化

\[ Q_w(a,b)=\sum_{i=1}^n w_i(y_i-a-bx_i)^2. \]

记

\[ 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), \]

则

\[ \hat b=\frac{S_{xy,w}}{S_{xx,w}},\qquad \hat a=\bar y_w-\hat b\bar x_w. \]

OLS 是 \(w_i=1\) 的特例。weights 各选项改变的是误差方差随浓度变化的假设:

weights 实际权重 隐含的相对方差形态 主要区别
NULL 无权重 OLS 假定 \(\operatorname{Var}(Y\mid X)\) 恒定 method="wls" 时不能为 NULL
"equal" \(1\) 与 OLS 相同 用于显式进入 WLS 接口
"1/y" \(1/\lvert y_i\rvert\) \(\operatorname{Var}(Y_i)\propto \lvert y_i\rvert\) 权重依赖候选方法观测值
"1/y^2" \(1/y_i^2\) \(SD(Y_i)\propto |y_i|\) 候选方法恒定 CV 的经验权重
"1/x" \(1/\lvert x_i\rvert\) \(\operatorname{Var}(Y_i)\propto \lvert x_i\rvert\) 权重依赖参考方法观测值
"1/x^2" \(1/x_i^2\) \(SD(Y_i)\propto |x_i|\) 参考尺度上的恒定 CV 经验权重
"constant_cv" \(1/x_i^2\) 与上一行相同 WLS 中严格要求 \(x_i>0\);名称明确其模型含义
"residual" \(1/\hat\sigma_i^2\) SD 由浓度的一次函数描述 仅用于 WLS,并迭代更新
数据列名 用户给定 \(w_i\) 由用户定义 要求非负、有限;零权重等于不让该点贡献目标函数

"1/x" 等通用内置权重遇到零分母时以机器精度下限保护数值;这只避免程序除零,并不使零浓度上的相对误差模型具有统计意义。"constant_cv" 因而直接拒绝非正浓度。

残差加权先拟合 OLS,令第 \(k\) 轮残差为

\[ e_i^{(k)}=y_i-\hat a^{(k)}-\hat b^{(k)}x_i, \]

再拟合绝对残差函数

\[ |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}. \]

随后重新拟合 WLS。weight_iterations 是上述更新次数,默认 4;改变它可能改变最终权重和回归系数,但不应为了贴合某个已公布区间而调节。任何一轮预测 SD 非正或非有限都会停止。

OLS/WLS 的解析协方差为

\[ \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}, \]

系数区间为 \(\hat\theta\pm t_{1-\alpha/2,n-2}SE(\hat\theta)\)。这与 EP09c 附录 H 在相同最终权重下的计算等价。YY/T 1789.2 式(8)在医学决定水平偏倚区间中没有再次乘 \(\hat\sigma_w\);只有当权重是已知绝对逆方差时,两种写法才不需要残差尺度校正。

21.3.4 Deming 回归与 lambda

Deming 回归定义

\[ 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), \]

其中 \(\lambda\) 在源码中表示 Y 误差方差与 X 误差方差之比。斜率和截距为

\[ 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. \]

因此

\[ \lambda=\frac{\sigma_Y^2}{\sigma_X^2}. \]

lambda = 1 表示两种方法的误差方差相同;\(\lambda>1\) 表示候选方法 Y 的误差较大,\(\lambda<1\) 表示参考方法 X 的误差较大。把比值方向写反会改变 Deming 斜率,因而必须始终按“候选方差除以参考方差”录入。调用底层算法时所需的倒数换算在包内部完成,不改变用户接口。

lambda = "replicates" 从已经汇总的重复测量 SD 估计方差比。若第 \(i\) 个样本在两种方法分别有 \(n_{Xi}\)、\(n_{Yi}\) 次重复,先计算自由度合并的单次测量方差

\[ 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)}. \]

回归使用的是样本均值,所以实际输入的方差比为

\[ \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}. \]

两种方法的重复次数可以彼此不同,但每种方法内部必须在样本间恒定;否则单个标量 \(\lambda\) 不能描述各样本均值的不同误差方差,应改用精密度轮廓。没有提供重复次数时,源码假定两侧重复次数相等,\(n_X/n_Y\) 抵消。\(\lambda\) 只估计一次;jackknife 和 bootstrap 重拟合时把它视为固定的外部输入,不在每个子样本内重新估计。

21.3.5 三类加权 Deming

固定观测权重的 method="wdeming" 先令 \(\tilde w_i=w_i/\sum_jw_j\),以

\[ \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) \]

替换普通中心矩,再代入 Deming 闭式。它只给每一对观测一个共同权重,不能表达 X、Y 两侧各自不同且随浓度变化的方差函数。

weights="constant_cv" 在 WLS 中只是 \(1/x_i^2\);在加权 Deming 中则调用 Linnet 的迭代恒定 CV 算法。第 \(k\) 轮以两侧隐变量读数的均值构造

\[ w_i^{(k)}=\left(\frac{x_i^{*(k)}+y_i^{*(k)}}{2}\right)^{-2}, \]

用该权重更新 Deming 直线。若 \(d_i=y_i-a-bx_i\),按 ivdtools 的 \(\lambda=\sigma_Y^2/\sigma_X^2\) 记法,隐变量更新可写为

\[ x_i^*=x_i+\frac{b\,d_i}{\lambda+b^2},\qquad y_i^*=y_i-\frac{\lambda d_i}{\lambda+b^2}. \]

权重、直线和隐变量反复更新至相对变化不超过 epsilon。两侧浓度必须严格为正;删除一个样本或抽取一份 bootstrap 样本后,整个迭代从头完成,而不是沿用全数据的最终权重。

weights="precision_profile" 允许参考方法和候选方法分别给出

\[ g(\mu)=\operatorname{Var}(X\mid\mu),\qquad h(\eta)=\operatorname{Var}(Y\mid\eta),\quad \eta=a+b\mu. \]

对每个样本的潜在真值 \(\mu_i\),算法最小化带方差尺度项的目标函数

\[ 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]. \]

给定 \(a,b\) 和当前方差时,潜在真值更新为

\[ \mu_i= \frac{h_i x_i+g_i b(y_i-a)}{h_i+g_i b^2}. \]

variance_models 的两项可以是返回方差的函数,也可以是已经完成 profile() 的 precision 对象。后者由 profile_components 分别选择参考侧和候选侧的方差组分,默认两侧均为 "total",并把 VFP 预测 SD 平方后作为方差。选择 "error"、"day" 或设计中其他实际可用的方差组分会直接改变 \(g\) 或 \(h\);具体名称取决于 precision() 的公式,且必须与 方法比较数据所需的误差层级一致。超出精密度研究浓度范围只警告外推,预测方差非正或非有限则停止。外部精密度轮廓在重抽样中保持固定。

三类加权 Deming 的含义不能互换:固定权重是经验性地改变每个点的影响,恒定 CV Deming 同时建模两轴的相对误差,精密度轮廓 Deming 则允许两轴拥有不同的非恒定方差函数。

21.3.6 回归参数置信区间

ci_method 决定截距和斜率的不确定度算法。"auto" 的映射为:

回归方法 auto 实际方法 主要依据
OLS、WLS analytic 回归设计矩阵的解析协方差
Deming、恒定 CV Deming、精密度轮廓 Deming jackknife EP09c 附录 K1
Passing–Bablok rank EP09c 附录 I

jackknife 每次删除第 \(i\) 个样本并完整重拟合,得到 \(\hat{\boldsymbol\theta}_{(-i)}=(\hat a_{(-i)},\hat b_{(-i)})^\mathsf T\),其协方差为

\[ \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. \]

区间仍以全数据点估计为中心:

\[ \hat\theta\pm t_{1-\alpha/2,n-2}SE_{\mathrm{jack}}(\hat\theta). \]

bootstrap 始终对完整 \((x_i,y_i)\) 样本对有放回抽样;由于 mcr() 只接受每个 ID 一行,不存在把同一样本的重复行拆开抽样的问题。每份 bootstrap 样本都完整重拟合回归。若第 \(r\) 次得到 \(\hat{\boldsymbol\theta}^{*(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). \]

两种方法还用 bootstrap 截距、斜率样本的联合协方差计算任意浓度的偏倚。源码要求至少 5000 次;seed 只控制可复现性,不应作为改变结论的调节参数。失败拟合被删除并计数;有效拟合不足一半或不足 100 次时停止。

"analytic" 只允许 OLS/WLS,"rank" 只允许 Passing–Bablok。虽然可以显式要求其他方法使用 jackknife 或 bootstrap,但这会改变区间问题,不能把“结果更接近某标准公布值”作为选择依据。

21.3.7 Passing–Bablok

源码形成全部 \(x_j\ne x_i\) 的两两斜率

\[ b_{ij}=\frac{y_j-y_i}{x_j-x_i}, \]

完全重复的点对不提供斜率;垂直点对保留为带符号无穷。Algorithm I 先计算小于 \(-1\) 的斜率数 \(K\),删除恰好等于 \(-1\) 的斜率,再将点估计和区间的秩循环平移 \(K\) 位。平移后的中位数为 \(\hat b\),并令

\[ \hat a=\operatorname{median}_i(y_i-\hat b x_i). \]

若剩余有效斜率数为 \(m\),令

\[ 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\}, \]

并令 \(M_2=m-M_1+1\)。斜率上下限取平移后的第 \(M_1\)、\(M_2\) 个顺序斜率;截距界限分别为

\[ a_L=\operatorname{median}_i(y_i-b_Ux_i),\qquad a_U=\operatorname{median}_i(y_i-b_Lx_i). \]

这就是回归参数的 ci_method="rank"。秩区间并不直接给出某个医学决定水平处偏倚的联合区间;因此 bias() 和 predict() 对 Passing–Bablok 的 ci_method="auto" 改用成对百分位 bootstrap,而不使用 EP09c K1 明确排除的 jackknife。对象中保留的 jackknife 参数协方差只服务于显式要求的扩展型预测计算,不是 Passing–Bablok 默认偏倚区间。

21.3.8 差值、百分比差与离群分析

绝对差始终按候选减参考:

\[ d_i=y_i-x_i. \]

百分比差的 percent_denominator 有三个不同估计对象:

\[ 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" 回答“相对于参考方法高或低多少”,是默认值;"mean" 对调两种方法时只改变符号,适合没有天然参考方法的对称比较;"candidate" 以候选结果为基数。分母选择会同时改变数值、分布和离群点,不能只改图轴标签。outlier(type="percent") 与 bland_altman(type="percent") 使用同一组选项,因而可保持检测量和差值图完全一致。

outlier(type="difference") 对 \(d_i\) 检验;type="percent" 对所选 \(p_i\) 检验。method 分别调用本章前述 Grubbs、广义 ESD、Dixon 或 IQR 规则:前三者由 alpha 控制显著性,IQR 由 coef 控制围栏宽度。离群检验只标记样本,不会自动从后续回归或 Bland–Altman 中删除。

21.3.9 Bland–Altman 横轴、中心与一致性限

x_axis 只改变图的横坐标:

\[ m_i=\frac{x_i+y_i}{2},\qquad x_i,\qquad y_i \]

分别对应 "mean"、"reference"、"candidate";它不改变纵轴差值。type 决定纵轴:

\[ \text{difference}:y_i-x_i,\qquad \text{ratio}:\frac{y_i}{x_i},\qquad \text{percent}:p_i. \]

参数法 method="parametric" 以均值为中心:

\[ \hat\delta=\bar d,\qquad CI_\delta=\bar d\pm t_{1-\alpha/2,n-1}\frac{s_d}{\sqrt n}. \]

若 agree.level 为 \(p\),令 \(z_p=z_{(1+p)/2}\),一致性限为

\[ L_-=\bar d-z_ps_d,\qquad L_+=\bar d+z_ps_d. \]

各一致性限的标准误和区间为

\[ 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 控制一致性限希望覆盖的个体差值比例;conf.level 控制中心和一致性限本身的统计置信度。两者含义不同,例如把 agree.level 从 0.95 改为 0.99 会扩大 LoA,而把 conf.level 从 0.95 改为 0.99 主要扩大 LoA 外面的置信区间。

非参数法 method="nonparametric" 以普通样本中位数 \(d_{(\lceil n/2\rceil)}\)(偶数时取中间两项均值)为中心。median_ci_method="order_statistic" 对候选下限秩 \(k\) 计算实际覆盖率

\[ C(k)=1-2P\{\operatorname{Bin}(n,0.5)\le k-1\}, \]

选择满足 \(C(k)\ge\text{conf.level}\) 的最大 \(k\),区间为

\[ [d_{(k)},d_{(n-k+1)}]. \]

它至少达到请求覆盖率,通常略保守,并报告实际覆盖率。"interpolated" 使用 EP09c 附录 F 的连续秩正态近似

\[ 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}, \]

若 \(r=j+q\)、\(j=\lfloor r\rfloor\)、\(0\le q<1\),则线性插值得到

\[ d_{(r)}=d_{(j)}+q\{d_{(j+1)}-d_{(j)}\}. \]

插值法给出名义覆盖率,边界不必等于某个实际观测;顺序统计量法给出离散而保守的区间。

method="hodges_lehmann" 的中心不是普通中位数,而是所有 Walsh 均值的中位数:

\[ \hat\delta_{HL}=\operatorname{median}_{i\le j}\left(\frac{d_i+d_j}{2}\right). \]

其中心区间采用无连续性校正的 Wilcoxon/Tukey 正态近似。HL 估计的是差值分布的伪中位数;只有在分布对称等条件下,它才与普通中位数具有相近解释,不能把 hl_estimate 标成 median_diff。

非参数和 HL 方法的一致性限均取经验分位数

\[ L_-=Q_{(1-p)/2}(d),\qquad L_+=Q_{1-(1-p)/2}(d), \]

使用 R 的 type 7 线性插值分位数。两端 LoA 的置信区间通过对 \(d_i\) 有放回抽样、每次重新计算两个经验分位数,再取 bootstrap 百分位数得到。bootstrap_replicates 默认 5000、最低 100;seed 只保证复现。这个经验 LoA 及其 bootstrap 区间是 ivdtools 保留的扩展,不应与 EP09c 附录 F 的中心估计区间混为一谈。

21.3.10 医学决定水平偏倚与预测区间

对给定 \(x_0\),回归预测为 \(\hat y_0=\hat a+\hat b x_0\),拟合均值方差为

\[ V_{\mathrm{fit}}(x_0)= SE_a^2+x_0^2SE_b^2+2x_0\operatorname{Cov}(a,b). \]

均值区间使用 \(\sqrt{V_{\mathrm{fit}}}\),单次结果预测区间使用 \(\sqrt{V_{\mathrm{fit}}+s_e^2}\)。医学决定水平上的候选方法减参考方法偏倚定义为

\[ \widehat{\operatorname{bias}}(x_0) =(\hat a+\hat b x_0)-x_0 =\hat a+(\hat b-1)x_0. \]

减去固定的 \(x_0\) 不增加方差,因此解析或 jackknife 参数协方差下

\[ SE\{\widehat{\operatorname{bias}}(x_0)\} =\sqrt{SE_a^2+x_0^2SE_b^2+2x_0\operatorname{Cov}(a,b)}. \]

百分位 bootstrap 则对每一份参数样本直接计算

\[ B^{*(r)}(x_0)=a^{*(r)}+[b^{*(r)}-1]x_0 \]

并取其分位数;standard bootstrap 使用这些偏倚样本的 SD 和 \(t_{n-2}\)。bias() 与 predict() 的 ci_method="auto" 对 Passing–Bablok 使用百分位 bootstrap,对其他方法复用回归对象已经选择的解析、jackknife 或 bootstrap 参数不确定度。若回归对象已经保存相同方法、足够次数且 seed 相容的 bootstrap 样本,就直接复用。

interval="confidence" 描述给定 \(x_0\) 处平均回归关系或平均偏倚的不确定度;"prediction" 再加入 \(s_e^2\),描述未来单次候选方法结果的散布;"both" 同时返回。预测区间是 ivdtools 扩展,并不是 EP09c 的回归参数或偏倚置信区间。

21.3.11 参数选择改变了什么

参数 改变的数学对象 不改变的内容 选择原则
method 直线的损失函数或秩估计量 X/Y 分析方向 按两轴误差及差值形态预先选择
lambda Deming 两轴误差的相对惩罚 原始观测值 用独立精密度或合规重复测量估计,不为贴合斜率调参
weights 各样本贡献或浓度相关方差模型 候选减参考的方向 恒定 SD、恒定 CV、残差 SD 或精密度轮廓应有实验依据
profile_components 精密度轮廓中的方差组分 \(g,h\) 方法比较均值 选择与研究目标一致的重复性、中间精密度或总精密度
ci_method 参数不确定度的估计方法 全数据点估计 auto 遵循方法对应的 EP09c 默认;其他方法用于预设的敏感性分析
conf.level 参数、中心或 LoA 的置信覆盖率 agree.level 通常预先设为 0.95
agree.level LoA 覆盖的个体差值比例 中心估计的置信度 按一致性定义设定,不能当作 CI 水平
median_ci_method 普通中位数 CI 的离散/连续秩算法 中位数点估计 保守实际覆盖选顺序统计量;复现附录 F 连续秩选插值法
percent_denominator 百分比差的分母和估计对象 绝对差 有参考方法选 reference;对称比较可选 mean
bootstrap_replicates Monte Carlo 误差 bootstrap 目标分布 回归至少 5000;更多次数提高分位数稳定性
seed 伪随机抽样序列 统计方法和目标参数 用于复现,不能用来挑选有利结果
split_by_cutpoints() 的 cutpoints 分析前各 dataframe 的样本范围 后续统计方法 依据预设临床区间,不按结果反复切分

21.3.12 参考物质偏倚

每个水平以测量均值 \(\bar x\) 与参考值 \(x_{\mathrm{ref}}\) 之差为偏倚:

\[ b=\bar x-x_{\mathrm{ref}},\qquad u_{\bar x}=\frac{s}{\sqrt n},\qquad U_{\bar x}=k u_{\bar x}. \]

源码将传入的参考值不确定度直接与 \(U_{\bar x}\) 合成:

\[ h=\sqrt{U_{\bar x}^2+U_{\mathrm{ref}}^2},\qquad [b-h,b+h]. \]

这里没有在内部把参考值不确定度再除以覆盖因子。相对偏倚为 \(100b/x_{\mathrm{ref}}\);参考值为零时返回缺失。

21.4 参考物质互换性

21.4.1 EP14 型预测区间

每个临床样本先在两个测量程序内求重复测量均值与标准差。各程序的临床样本合并重复性方差为

\[ s_{p}^2= \frac{\sum_j(n_j-1)s_j^2}{\sum_j(n_j-1)}. \]

源码令 \(\lambda=s_{p,Y}^2/s_{p,X}^2\),以临床样本均值拟合 Deming 直线。该分支使用除以 \(n\) 的总体中心矩。斜率方差为

\[ \operatorname{Var}(\hat b)= \frac{\hat b^2\max(0,S_{xx}S_{yy}-S_{xy}^2)}{nS_{xy}^2}. \]

参考物质在 \(x_R\) 处的预测方差为

\[ 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). \]

自由度取两个程序合并重复性自由度的较小值并向下取整。参考物质的 Y 均值落入 \(\hat a+\hat b x_R\pm t\sqrt{V_R}\) 时,源码标记为可互换。

21.4.2 IFCC 型差异评价

临床样本偏倚为 \(b_i=y_i-x_i\)。若明确参考测量程序,排序坐标用 \(x_i\);否则用 \((x_i+y_i)/2\)。按坐标排序后计算

\[ s_{\mathrm{MSSD}} =\sqrt{\frac{\sum_{i=1}^{n-1}(b_{i+1}-b_i)^2}{2(n-1)}}. \]

参考物质偏倚与所选临床样本偏倚均值之差记为

\[d=(y_R-x_R)-\bar b_C.\]

临床子集可为全部样本、距离最近的 \(q\) 个样本,或目标位置两侧各取一半。区间标准不确定度为

\[ u_d=\sqrt{\frac{s_{X,R}^2}{p_X}+\frac{s_{Y,R}^2}{p_Y} +\frac{s_C^2}{q}}. \]

扩展区间为 \(d\pm k u_d\),源码默认 \(k=1.9\)。若整个区间位于允许差异带 \([-\Delta,\Delta]\) 内,结果为可互换;若整个区间位于允许带同一侧之外,结果为不可互换;其余为不确定。

21.5 定性性能、C5/C95 与 ROC

21.5.1 四格表比例及区间

源码按第一行、第一列为阳性组织四格表,并计算

\[ 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}. \]

对成功数 \(x\),Wald 区间为

\[ \hat p\pm z_{1-\alpha/2}\sqrt{\frac{\hat p(1-\hat p)}n}; \]

连续性校正版在半宽上再加 \(1/(2n)\)。Wilson 区间的中心和半宽为

\[ 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 连续性校正是在 \(c\pm h\) 外再减、加 \(1/(2n)\)。Agresti–Coull 使用

\[ \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 区间使用 \(\operatorname{Beta}(x+1/2,n-x+1/2)\) 分位数;Clopper–Pearson 区间使用 \(\operatorname{Beta}(x,n-x+1)\) 与 \(\operatorname{Beta}(x+1,n-x)\) 分位数,边界 \(x=0,n\) 单独处理。所有比例区间最终截断到 \([0,1]\)。

给定患病率 \(\pi\) 时,预测值重新按

\[ PPV_\pi=\frac{Se\,\pi}{Se\,\pi+(1-Sp)(1-\pi)},\qquad NPV_\pi=\frac{Sp(1-\pi)}{(1-Se)\pi+Sp(1-\pi)} \]

计算。区间只是把灵敏度和特异度各自的上下限代入,没有传播协方差。

21.5.2 似然比、一致性与 McNemar

\[LR^+=\frac{Se}{1-Sp},\qquad LR^-=\frac{1-Se}{Sp}.\]

源码在对数尺度上计算区间,例如

\[ SE\{\log(LR^+)\}= \sqrt{\frac1{TP}-\frac1{TP+FN}+\frac1{FP}-\frac1{FP+TN}}. \]

优势比为 \(OR=TP\cdot TN/(FP\cdot FN)\),对数标准误为四个格子倒数之和的平方根。相关格子为零时返回缺失区间,不自动加连续性校正。

观察一致率和随机一致率为

\[ 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\) 区间,并截断到 \([-1,1]\)。PABAK 分支使用

\[PABAK=2p_o-1,qquad \operatorname{Var}(PABAK)=\frac{4p_o(1-p_o)}n.\]

该分支由是否传入患病率触发,但患病率数值不进入公式。

McNemar 在 \(b+c<25\) 时,以 \(B\sim\operatorname{Binomial}(b+c,0.5)\) 做双侧精确检验;否则计算

\[\chi^2=\frac{(|b-c|-1)^2}{b+c}\]

并取 1 个自由度。当 \(b+c=0\) 时,p 值置为 1。

21.5.3 C5/C95

原始重复测量先按样本水平形成均值 \(\mu_j\)、标准差 \(s_j\) 和样本量 \(n_j\)。源码要求至少三个 \(n_j\ge2\) 且 \(s_j>0\) 的水平,再拟合方差函数得到浓度相关标准差 \(s(\mu)\)。

若事件定义为结果不低于 cutoff \(c\),则

\[P(Y\ge c\mid\mu)=1-\Phi\!\left(\frac{c-\mu}{s(\mu)}\right);\]

反向事件改用对应的 \(\Phi\)。源码在 1001 点网格上搜索目标概率与该函数的符号变化,再用一维求根细化;无根时最多扩展搜索范围 16 次。若出现多个根,选择最接近 cutoff 的根。这是“正态响应分布 + 方差轮廓”的反演,不是二项 probit 回归。

21.5.4 经验 ROC、AUC 与多变量模型

候选阈值为 \(-\infty\)、全部唯一观测值和 \(+\infty\)。将 ROC 点按假阳性率递增排序后AUC 用梯形和计算:

\[ A=\sum_i\frac{TPR_i+TPR_{i+1}}2(FPR_{i+1}-FPR_i). \]

当 \(0<A<1\) 时,源码用 Hanley–McNeil 近似:

\[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_-}}. \]

AUC 区间固定为 \(A\pm1.96SE(A)\) 并截断到 \([0,1]\);\(A=0\) 或 \(1\) 时不计算标准误。Youden 截断值最大化 \(Se+Sp-1\);左上角法最小化 \(\sqrt{(1-Se)^2+(1-Sp)^2}\)。并列时返回当前顺序的第一个点。

同一批受试者的两条 AUC 不能按独立样本比较。auc_compare(method="ep24") 采用 EP24-A2 Hanley–McNeil 配对近似:分别求阳性组和阴性组内两项评分的 Pearson 相关(有序评分可选 Kendall),取平均后与平均 AUC 对 EP24 表 5 线性插值得到面积相关 \(r\),再计算

\[ SE(A_1-A_2)=\sqrt{SE_1^2+SE_2^2-2rSE_1SE_2}. \]

auc_compare(method="delong") 以两条评分对每名阳性和阴性受试者的秩贡献建立协方差矩阵 \(\widehat\Sigma\),并令 \(\mathbf c=(1,-1)^\mathsf T\):

\[ SE(A_1-A_2)=\sqrt{\mathbf c^\mathsf T\widehat\Sigma\mathbf c}. \]

两种方法都用正态近似检验 \(H_0:A_1-A_2=0\)。当曲线来自训练集 Logistic 拟合概率时,它们把评分当作既定值,不校正模型拟合不确定性和表观性能的乐观偏倚。

多指标模型实际拟合

\[ \operatorname{logit}(p_i)=\beta_0+\sum_{j=1}^q\beta_jx_{ij} \]

的二项 Logistic 回归,并用训练数据预测概率建立同一经验 ROC。新观测的线性预测量区间为 \(\hat\eta\pm1.96SE(\hat\eta)\),再经 Logistic 逆变换。概率标准误采用 \(SE(\hat p)=SE(\hat\eta)\hat p(1-\hat p)\),分类阈值固定为 0.5。

21.6 精密度、瓶间差与参考区间

21.6.1 方差分量与精密度轮廓

源码先判断设计是否平衡:平衡设计使用 ANOVA 方差分量,不平衡设计使用 REML。每个随机层级返回 \(SD_j=\sqrt{\widehat{VC}_j}\) 和 \(CV_j=100SD_j/\bar y\);总标准差为

\[SD_T=\sqrt{\sum_j\widehat{VC}_j}.\]

方差分量区间完全采用 VCA 推断对象中的双侧 SD/CV 区间,源码没有另写 Satterthwaite–卡方区间公式。

浓度相关方差轮廓由 VFP 拟合。源码注册的十个模型为:

编号 实际方差或 CV 模型
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\)

候选模型多于一个时按 AIC 最小选择。输入为 CV 时,源码先结合均值还原 SD,再以 \(VC=SD^2\) 进入拟合。

21.6.2 瓶间差

瓶、批或其他分组直接进入固定效应方差分析。组内描述区间为

\[ \bar x_j\pm t_{1-\alpha/2,n_j-1}\frac{s_j}{\sqrt{n_j}}. \]

残差执行 Shapiro–Wilk,组间方差执行 Bartlett;两两比较直接使用 Tukey HSD。源码没有将瓶间均方自动转换为均匀性不确定度;该计算只存在于不确定度模块的随机效应 ANOVA 分支。

21.6.3 参考区间

默认非参数法令 \(\alpha=(1-p)/2\),按 EP28 计算连续秩

\[ h_q=q(n+1), \]

并在 \(\lfloor h_q\rfloor\) 与 \(\lceil h_q\rceil\) 对应的顺序统计量间线性插值。参考限置信区间使用二项分布秩反演。令 \(\delta=(1-\gamma)/2\)、\(B\sim\mathrm{Bin}(n,q)\),秩端点为

\[ a=Q_B(\delta),\qquad b=Q_B(1-\delta)+1, \]

再限制到 \([1,n]\);实际覆盖率为 \(P(a\le B\le b-1)\)。若截断前秩超出样本范围,函数明确标记请求覆盖率不可完全达到。method="type6" 保留相同的 R type-6 点估计,但其 CI 仍使用原有二项秩正态近似,不能与默认非参数法的精确秩 CI 混称。

参数法先检验原始值;若 p 值小于 0.05 且所有值为正,再检验对数值,并只在对数检验 p 值不小于 0.05 时使用对数尺度。所用尺度上的参考限为

\[\bar x\pm z_{1-\alpha}s.\]

每个端点的标准误按

\[SE_L=s\sqrt{\frac1n+\frac{z_{1-\alpha}^2}{2n}}\]

计算,以 \(t_{1-\gamma/2,\max(n-2,1)}\) 构造端点区间;若采用对数尺度,最后整体取指数。

method="biweight" 实现 EP28 附录 B。以中位数初始化 \(T_{bi}\),令 \(s=MAD/0.6745\),每轮使用

\[ 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}, \]

直到相对变化小于 0.001%。对调节常数 \(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)}}. \]

按附录分别计算 \(s_{bi}[205.6]\)、\(s_{bi}[3.7]\);\(S_T[3.7]\) 改用 \(u_i=(x_i-T_{bi})/(3.7s_{bi}[3.7])\),并在根号分子中去掉 \(n\)。最终参考限为

\[ T_{bi}\pm t_{n-1,1-\alpha} \sqrt{s_{bi}^2[205.6]+S_T^2[3.7]}. \]

端点区间由同一 biweight 估计量的 percentile bootstrap 得到。

method="huber" 保留包原有扩展:以中位数初始化位置 \(T\),以 \(MAD/0.6745\) 初始化尺度 \(S\),每轮令

\[ 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}}. \]

迭代容差为 \(10^{-8}\),最多 100 轮,参考限为 \(T\pm z_{1-\alpha}S\)。biweight 与 Huber 的端点区间默认来自 1999 次有放回抽样;seed 可固定计算且不改变调用者原有随机数状态。

21.7 响应方程拟合与线性评价

21.7.1 已实现的十二个响应方程

编号 源码中的方程
E01 \(a+bx\)
E02 \(a+bx+cx^2\)
E03 \(a\exp(bx)\)
E04 \(a\exp(-bx)\)
E05 \(ax^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/(2c^2)]\)
E10 \(V_{\max}x/(K_m+x)\)
E11 \(a\sin(bx+c)+d\)
E12 \(a+bx+cx^2+dx^3\)

E01、E02、E06、E12 用线性模型矩阵估计,其余方程用 Levenberg–Marquardt 非线性最小二乘。目标函数为

\[RSS_w=\sum_iw_i[y_i-f(x_i,\theta)]^2.\]

仅有参数上下界时先尝试有界非线性最小二乘;一般不等式约束使用 COBYLA,等式约束使用增广 Lagrangian。优化失败后的惩罚目标为

\[RSS_w+10^8\sum_r\{\text{constraint violation}_r\}^2.\]

约束拟合的协方差由中心差分残差 Jacobian \(J\) 近似:

\[ \widehat V(\hat\theta) \approx(J^\mathsf TJ)^{-1}\frac{RSS}{n-p}. \]

该分支的 Jacobian 协方差没有把拟合权重乘入 \(J\)。约束模型的 AIC 为 \(n\log(RSS/n)+2p\),而 \(R^2=1-RSS/TSS\)。

非线性预测用中心差分求梯度 \(g(x)\),并令

\[V\{\hat f(x)\}=g(x)^\mathsf T\widehat V(\hat\theta)g(x).\]

单次预测区间在该方差上再加 \(RSS/(n-p)\)。

源码对 E01–E07 使用解析逆式;E09 固定取高斯峰左支,E11 取 asin 主值分支,E12 对多项式根取实部后选择最接近观测 X 均值的根。E08 当前逆式为

\[ x=C\left(\frac{B-A}{y-A}-1\right)^{1/(DE)}, \]

它不是 \(E\ne1\) 时当前 E08 方程的严格代数逆。其他情况在扩展网格上搜索符号变化,用一维求根求解并选择最接近观测 X 中心的根。

21.7.2 线性设计量与回归比较

稀释面板的期望浓度按 \(C(p)=pC_H+(1-p)C_L\) 计算。重复数规划先令每个水平的置信概率 \(p_L=(1-\alpha_F)^{1/L}\),再取 \(z=z_{(1+p_L)/2}\):

\[ n=\left\lceil \left(\frac{z\sigma}{ADL-d_{\mathrm{true}}}\right)^2 \right\rceil. \]

线性分析先按水平形成均值、SD 和 \(n_j\)。等权时水平权重为 1;逐水平方差法对原始观测用 \(1/s_j^2\),对水平平均值用 \(n_j/s_j^2\)。weights="pooled" 要求 variance_group 明确指定按 \(x\) 连续的水平组,对组 \(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}. \]

函数不根据水平数自动分组,因为 EP06 分组还要结合重复性 SD—浓度关系。精密度轮廓法先以普通直线拟合 SD 对均值,再用预测 \(1/\hat\sigma_j^2\) 作权重。

直线拟合水平均值;二次和三次多项式拟合全部保留观测。若任何非线性项显著,则在有显著非线性项的模型中选择 \(s_{y\cdot x}\) 最小者;否则选择直线。非线性偏离为

\[DL_j=\hat y_{\mathrm{poly},j}-\hat y_{\mathrm{linear},j},\]

线性残差为 \(\bar y_j-\hat y_{\mathrm{linear},j}\)。百分比可分别以线性预测值、X 值或实测均值作分母。

每个水平的点区间使用家族置信概率 \(\gamma^{1/L}\):

\[ h_j=z_{[1+\gamma^{1/L}]/2}\frac{\hat\sigma_j}{\sqrt{n_j}}. \]

源码以偏离量为区间中心。允许偏离限取直接绝对限、绝对限和相对限换算值中的最大值,最终仅比较 \(|\text{metric}_j|>ADL_j\)。

21.8 干扰、稀释、加标与钩状效应

21.8.1 双组干扰与差中之差

双组干扰效应为

\[D=\bar y_T-\bar y_C,qquad SE(D)=\sqrt{\frac{s_T^2}{n_T}+\frac{s_C^2}{n_C}}.\]

源码用 \(n_T+n_C-2\) 自由度的 t 临界值构造区间,虽然标准误采用未合并方差形式。百分比效应 为 \(100D/\bar y_C\),标准误和区间端点也直接缩放。

分析物加标设计的效应是

\[ D=(\bar y_{I,S}-\bar y_{I,U}) -(\bar y_{N,S}-\bar y_{N,U}), \]

其方差为四个单元 \(s_g^2/n_g\) 之和,自由度用 Welch–Satterthwaite 合成。百分比效应写为 \(100(I/N-1)\),其中 \(I=\bar y_{I,S}-\bar y_{I,U}\)、\(N=\bar y_{N,S}-\bar y_{N,U}\);源码按该表达式对四个均值求梯度并用 delta 法计算方差。

空白干扰在有至少两个对照结果时仍用双组差;否则把理论对照均值置为零,使用单样本标准误。剂量响应按各浓度均值相对观测基线计算差异,可另拟合原始结果对浓度的 OLS 直线。点间正向和反向换算使用线性插值,不外推;反向插值要求响应单调。

21.8.2 稀释回收

稀释倍数为 \(D\) 时,源码令 \(f=1/D\)。原样本赋值为 \(C_0\)、稀释液本底为 \(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}}}. \]

回收率使用实测稀释结果与目标值之比,不使用校正后浓度与原样本赋值之比。多个样本的总体回收率区间是样本级回收率均值的 t 区间。

21.8.3 加标回收与钩状效应

\[ 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}}}. \]

钩状效应先按浓度计算均值及 t 区间,取第一个最大均值所在浓度为峰值。峰后首次出现“浓度高于 ULoQ 且平均响应不高于阈值”的水平时,源码报告其前一个已测浓度为钩状效应浓度;不做水平间插值。未出现回落时,只报告最高测试浓度为下界。

21.9 分析灵敏度

21.9.1 空白限与检出限

非参数空白限以 \(r=0.5+B(1-\alpha)\) 确定排序秩,并在相邻整数秩间线性插值。参数法令 \(\nu=B-K\),计算

\[ c_p=\frac{z_{1-\alpha}}{\sqrt{1-1/(4\nu)}}, \qquad LoB=\bar y_B+c_ps_B. \]

零空白限分支直接令 \(LoB=0\)。存在 2 或 3 个试剂批时默认逐批计算并取最大 LoB;其他批次数默认合并。

经典 LoD 对低水平标准差作自由度加权合并:

\[ s_p=\sqrt{\frac{\sum_j(n_j-1)s_j^2}{\sum_j(n_j-1)}}, \qquad LoD=LoB+c_ps_p. \]

非参数试验法以全部低水平结果中位数为候选 LoD,仅当低于 LoB 的比例不超过 \(\beta\) 时保留。精密度轮廓法拟合 \(s(x)\) 后求解

\[x-LoB-c_ps(x)=0,\]

并取搜索到的最小合格根。

Probit 法拟合

\[ \Phi^{-1}\{P(\mathrm{detected}\mid x)\}=\beta_0+\beta_1g(x), \]

其中 \(g(x)=\log_{10}x\) 或 \(x\)。目标检出概率为 \(1-\beta\):

\[g(LoD)=\frac{z_{1-\beta}-\hat\beta_0}{\hat\beta_1}.\]

分批或分层拟合时,最终报告各组 LoD 的最大值。

21.9.2 定量限与限值验证

每个低水平样本的偏倚为 \(b=\bar y-x_a\)。源码实现的误差量为

\[ E_W=|b|+ks,qquad E_{RMS}=\sqrt{s^2+b^2},qquad E_P=s. \]

相对目标时再除以 \(|x_a|\)。观测候选 LoQ 是误差量达到目标的样本中最小实测均值;启用轮廓时,以误差量对赋值的一次或二次回归与目标线的最小根作为候选。多批取最大候选,最后在两者有限时取

\[LoQ_{\mathrm{reported}}=\max(LoQ_{\mathrm{candidate}},LoD).\]

LoB 验证把 \(y\le LoB\) 记为一致,LoD 验证把 \(y\ge LoB\) 记为一致,LoQ 验证把结果落入赋值的允许窗记为一致。源码计算单侧 Wilson 区间,并按内置样本量–边界表选择验证边界;样本量不在表中时取相邻表行中较大的边界。只返回比例、区间和与边界之差,不生成通过/不通过结论。

21.10 质量控制与稳定性

21.10.1 Levey–Jennings 与 Westgard

控制结果先标准化为 \(z_i=(x_i-\mu_0)/\sigma_0\)。未给目标量时,源码用全部有限结果的样本均值和样本 SD。实际规则为:

规则 源码判定
1-2s \(|z_i|>2\)
1-3s \(|z_i|>3\)
2-2s 相邻两点同侧且都严格超过 \(2s\)
R-4s 相邻两点各自 \(|z|\ge2\) 且 \(|z_i-z_{i+1}|\ge4\)
4-1s 连续四点同侧且都严格超过 \(1s\)
8x/10x 连续 8/10 点严格位于均值同一侧

分组变量只用于图形分面;目标量估计和规则扫描跨全部有限结果连续执行。源码把 1-3s 标为失控,其他违反标为警告。Youden 以 \(|z_{1i}|\le2\) 且 \(|z_{2i}|\le2\) 判断是否位于 \(\pm2s\) 矩形内,并计算 Pearson 相关;没有拟合误差椭圆。

21.10.2 节点偏倚、时间回归与越限时间

每个条件以自身最早时间点均值为基线:

\[b_t=y_t-y_0,qquad b_{t,\%}=100\frac{y_t-y_0}{y_0}.\]

自比较先按时间点求均值。若各时间点标准误均有限且为正,以 \(1/SE_t^2\) 加权拟合;否则退回 OLS。绝对漂移率为 \(b_1\),相对漂移率为 \(b_1/b_0\),其方差按

\[ 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} \]

计算。源码使用单侧临界值 \(t_{\gamma,\nu}\)。下降取下界、上升取上界;方向无法确定时取绝对值更大的一侧。

条件比较以匹配时间点的试验条件减参考条件为偏倚。绝对偏倚标准误为 \(\sqrt{SE_T^2+SE_R^2}\);相对偏倚标准误按比值 delta 法计算。随后仍按标准误可用性选择 WLS 或 OLS。

越限时间通过 501 点网格寻找 \(|\hat b(t)|-L=0\) 或单侧置信界与 \(L\) 的首次符号变化,再用一维求根细化。默认搜索上限为 \(\max(100t_{\max},t_{\max}+100)\);超过观测最大时间的根标为外推。

21.10.3 平均动力学温度与 Arrhenius

令 \(q_i=\exp[-E_a/(RT_i)]\)。等时间间隔时

\[ \bar q=\frac1n\sum_iq_i,qquad T_{MKT}=-\frac{E_a}{R\log\bar q}. \]

有时间戳时用梯形积分

\[ \bar q=\frac{\sum_i\Delta t_i(q_i+q_{i+1})/2}{\sum_i\Delta t_i}. \]

多个探头按各自有效观测数或时间长度合并 \(q\),不是直接平均各探头 MKT。

Arrhenius 分支在每个温度拟合零级 \(y=a+bt\) 或一级 \(\log y=a+bt\),按退化方向转换为正速率。自动级数以各温度原尺度归一化 RMSE 的均方根最小为准;完全相同时选择一级。正速率再拟合

\[\log k=\alpha+\beta/T,qquad E_a=-\hat\beta R.\]

零级货架期为允许变化量除以目标温度速率。一级相对下降和上升分别为

\[t=-\frac{\log(1-r)}{k_*},qquad t=\frac{\log(1+r)}{k_*}.\]

货架期区间由目标速率置信区间两端分别换算后排序得到。

稳定性设计不是连续功效公式,而是源码内置的两个查表矩阵。重复性按 \(\sigma/\sqrt r\) 缩减,或按 \(\sqrt{\sum_j\sigma_j^2/L_j}\) 合成;结果保留两位有效数字后形成两个比值,再查 80% 或 90% 功效表。

21.11 测量不确定度与溯源

21.11.1 A 类与 B 类标准化

重复测量的单次结果标准不确定度为 \(s\),均值标准不确定度为 \(s/\sqrt n\),自由度均记为 \(n-1\)。相对分量在内部保存为 \(u/|\bar x|\)。

区间型 B 类分量要求估计值恰为上下限中点。令半宽度 \(a=(U-L)/2\),源码换算为:

分布 标准不确定度
矩形 \(a/\sqrt3\)
三角 \(a/\sqrt6\)
U 形 \(a/\sqrt2\)
正态中央覆盖概率 \(p\) \(a/z_{(1+p)/2}\)
t 分布中央覆盖概率 \(p\) \(a/t_{(1+p)/2,\nu}\)

若直接提供一个不确定度数值,源码把它视为标准不确定度;只有另给覆盖因子 \(k\) 时才执行 \(u=U/k\)。此路径只允许 standard 或 normal 标签,不再根据分布名称自动除以 \(\sqrt3\)、\(\sqrt6\) 等。因此分布除数只在上下限路径中产生。

21.11.2 IQC、均匀性、稳定性和定值

一元随机效应 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). \]

IQC 有运行内重复时,单次结果不确定度为

\[u_{\mathrm{single}}=\sqrt{s_W^2+s_B^2}.\]

若 \(s_B>0\),总体均值不确定度按 \(\sqrt{MS_B/N}\) 计算;否则退回全部结果的 \(s/\sqrt N\)。没有运行列或每次运行只有一个结果时,直接使用长期 SD 及其均值标准误。

均匀性分支另计算

\[ 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}}). \]

长期稳定性不确定度对每个条件拟合 \(y=a+bt\),在货架期 \(T_s\) 处取 \(u_{lts}=SE(\hat b)T_s\)。定值分支先求每个实验室均值,再以实验室均值的非加权平均赋值:

\[u_{char}=\frac{SD(\bar x_1,\ldots,\bar x_L)}{\sqrt L}.\]

21.11.3 一阶合成与有效自由度

令标准不确定度向量为 \(u\),相关矩阵为 \(R\):

\[ \Sigma=\operatorname{diag}(u)R\operatorname{diag}(u),qquad u_c^2=c^\mathsf T\Sigma c. \]

第 \(i\) 个分量的有符号边际贡献为 \(m_i=c_i(\Sigma c)_i\),贡献率为 \(100m_i/u_c^2\);存在相关性时贡献率可以为负。

仅在 \(R\) 为单位阵时,源码计算

\[ \nu_{\mathrm{eff}}= \frac{u_c^4}{\sum_i(c_iu_i)^4/\nu_i}. \]

未指定覆盖因子时,有限有效自由度使用 t 分位数,否则使用正态分位数;扩展不确定度为 \(U=ku_c\)。混合相对和绝对分量时,源码要求给出输出量值,并以其绝对值把相对分量转换为绝对分量。

21.11.4 一般测量模型与 Monte Carlo

对 \(y=f(x_1,\ldots,x_q)\),灵敏度系数用中心差分:

\[ 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)\}.\]

delta 分支将系数代入协方差传播,并返回对称区间 \(\hat y\pm U\)。

Monte Carlo 分支按保存分布抽样:normal 和 standard 标签用 \(N(x_i,u_i^2)\);矩形用均匀分布;三角分布用上下限和估计值为众数;U 形分布用 \(x_i+a_i\cos\{U(0,\pi)\}\);t 分布用 \(x_i+u_iT_{\nu_i}\)。因此 t 分支把 \(u_i\) 作为 t 随机变量的尺度。非独立 Monte Carlo 只允许 normal 或 standard 标签,并通过相关正态矩阵生成。

模拟输出的标准不确定度为模拟结果 SD,覆盖区间采用 R 第 8 类分位数的中央区间。

21.11.5 逐级溯源

第 \(r\) 阶段对前 \(r\) 阶段全部分量重新合成。灵敏度系数默认为 1,且该接口没有阶段间相关矩阵入口, 所以

\[u_{c,r}=\sqrt{\sum_{i\in1:r}u_i^2}.\]

每一级重新计算有效自由度、覆盖因子和扩展不确定度;最终级等于全部累计分量的预算。

21.12 样本量与功效

21.12.1 Bland–Altman 一致性

令差值均值为 \(\mu\)、SD 为 \(\sigma\)、允许界为 \(\Delta\)、一致性覆盖分位数为 \(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}. \]

令 \(t_c=t_{1-(1-C)/2,n-1}\),功效为

\[ \Pr\{T_{n-1}(\tau_1)>t_c\} +\Pr\{T_{n-1}(\tau_2)<-t_c\}. \]

求样本量时,在给定上下限内二分搜索达到目标功效的最小整数。

21.12.2 单比例区间与目标值检验

Wald 半宽的闭式样本量为

\[n=\left\lceil\frac{z_{1-\alpha/2}^2p(1-p)}{d^2}\right\rceil.\]

其他已实现区间用实际端点计算最大半宽,再经指数搜索和二分搜索寻找满足半宽 \(\le d\) 的最小整数。涉及 \(x=\operatorname{round}(np)\) 的方法会出现非单调跳变,源码还向下回扫修正结果。

单组比例目标值检验先寻找临界成功数。上侧检验的临界数是满足所选单侧比例区间下限 \(L(x,n)>p_0\) 的最小整数;下侧检验是满足区间上限 \(U(x,n)<p_0\) 的最大整数。因此

\[ Power_{>}=1-F_{\mathrm{Bin}(n,p_1)}(x_c-1),qquad Power_{<}=F_{\mathrm{Bin}(n,p_1)}(x_c). \]

双侧检验把 \(\alpha/2\) 分配到两端。求 \(n\) 时先用正态近似产生起点,再逐个增加并向下回查;求 \(\alpha\) 时先确定离散临界成功数,再反解使比例区间边界等于 \(p_0\) 的 \(\alpha\)。因此功效是由二项临界数决定的阶梯函数。