方法比较与参考物质偏倚
分析方向、数据单位与预处理分段
以下统一以参考或比较方法为 \(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,返回值保留原始分组列和各组行数。它只改变后续函数接收的数据集范围,不改变任何统计公式。
相关系数的三个选项
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 的选择首先改变误差模型:
"ols" |
视为可忽略 |
恒定或由可选权重描述 |
最小化纵向平方残差 |
参考 X 的误差相对很小 |
"wls" |
视为可忽略 |
随浓度变化 |
最小化加权纵向平方残差 |
差值 SD/CV 随浓度变化 |
"deming" |
不忽略 |
两轴误差比为常数 \(\lambda\) |
误差变量回归 |
两方法都有测量误差且方差形态近似恒定 |
"wdeming" |
不忽略 |
由共同权重、恒定 CV 或两侧方差轮廓描述 |
加权或迭代误差变量回归 |
两轴误差随浓度变化 |
"pb" |
不忽略 |
不要求正态误差模型 |
两两斜率的秩估计 |
线性但有异常值、混合变异或不宜作正态假设 |
选择回归方法依据的是研究设计和差值图所揭示的误差结构,而不是哪个方法产生更接近 1 的斜率。
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 各选项改变的是误差方差随浓度变化的假设:
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\);只有当权重是已知绝对逆方差时,两种写法才不需要残差尺度校正。
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 重拟合时把它视为固定的外部输入,不在每个子样本内重新估计。
三类加权 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 则允许两轴拥有不同的非恒定方差函数。
回归参数置信区间
ci_method 决定截距和斜率的不确定度算法。"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,但这会改变区间问题,不能把“结果更接近某标准公布值”作为选择依据。
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 默认偏倚区间。
差值、百分比差与离群分析
绝对差始终按候选减参考:
\[
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 中删除。
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 的中心估计区间混为一谈。
医学决定水平偏倚与预测区间
对给定 \(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 的回归参数或偏倚置信区间。
参数选择改变了什么
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 的样本范围 |
后续统计方法 |
依据预设临床区间,不按结果反复切分 |
参考物质偏倚
每个水平以测量均值 \(\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}}\);参考值为零时返回缺失。
定性性能、C5/C95 与 ROC
四格表比例及区间
源码按第一行、第一列为阳性组织四格表,并计算
\[
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)}
\]
计算。区间只是把灵敏度和特异度各自的上下限代入,没有传播协方差。
似然比、一致性与 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。
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 回归。
经验 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。
精密度、瓶间差与参考区间
方差分量与精密度轮廓
源码先判断设计是否平衡:平衡设计使用 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 拟合。源码注册的十个模型为:
| 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\) 进入拟合。
瓶间差
瓶、批或其他分组直接进入固定效应方差分析。组内描述区间为
\[
\bar x_j\pm t_{1-\alpha/2,n_j-1}\frac{s_j}{\sqrt{n_j}}.
\]
残差执行 Shapiro–Wilk,组间方差执行 Bartlett;两两比较直接使用 Tukey HSD。源码没有将瓶间均方自动转换为均匀性不确定度;该计算只存在于不确定度模块的随机效应 ANOVA 分支。
参考区间
默认非参数法令 \(\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 可固定计算且不改变调用者原有随机数状态。
测量不确定度与溯源
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\) 等。因此分布除数只在上下限路径中产生。
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}.\]
一阶合成与有效自由度
令标准不确定度向量为 \(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\)。混合相对和绝对分量时,源码要求给出输出量值,并以其绝对值把相对分量转换为绝对分量。
一般测量模型与 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 类分位数的中央区间。
逐级溯源
第 \(r\) 阶段对前 \(r\) 阶段全部分量重新合成。灵敏度系数默认为 1,且该接口没有阶段间相关矩阵入口, 所以
\[u_{c,r}=\sqrt{\sum_{i\in1:r}u_i^2}.\]
每一级重新计算有效自由度、覆盖因子和扩展不确定度;最终级等于全部累计分量的预算。