%%{init: {"theme":"base","themeVariables":{"primaryColor":"#EAF2FB","primaryBorderColor":"#2780E3","primaryTextColor":"#212529","lineColor":"#6C757D","secondaryColor":"#EAF7E6","tertiaryColor":"#FFF4E8","fontSize":"13px"},"flowchart":{"nodeSpacing":18,"rankSpacing":24,"padding":6}}}%%
flowchart LR
A["瓶/批抽样与<br/>重复设计检查"] --> B["bottle_anova()<br/>检验组间差异"]
B --> C["tukey()<br/>必要时定位差异"]
B --> D["precision() / variance()<br/>分解瓶间与瓶内方差"]
D --> E["uncertainty_homogeneity()<br/>均匀性不确定度"]
C --> F["结合预设均匀性限"]
E --> F
F --> G["形成放行依据"]
classDef input fill:#EAF2FB,stroke:#2780E3,color:#212529
classDef process fill:#EAF7E6,stroke:#3FB618,color:#212529
classDef output fill:#F4EEFB,stroke:#9954BB,color:#212529
class A input
class B,C,D,E,F process
class G output
18 瓶间差分析
本文用瓶间差和批次间比较判断物料的均匀性,并演示如何把均匀性结果转化为不确定度分量。显著性检验只说明是否观察到差异,不能单独说明差异是否具有实际影响;最终结论还要结合重复性、目标值和预设可接受标准。
bottle_anova() 的 F 检验回答“组均值是否完全相同”,不回答差异是否满足产品均匀性限。tukey() 只在整体模型和假设适用时定位成对差异。均匀性不确定度或瓶间 SD 必须由与抽样、测量顺序和重复设计一致的方差模型估计,再与预设技术限比较;不能只凭 p 值放行或拒绝。
分析流程:从瓶间差异到均匀性不确定度
瓶间差(bottle-to-bottle variability)评价同一产品不同瓶(或不同批次)之间测量结果的差异。分析通常分两步:先用 bottle_anova() 检验各瓶(或批次)均值是否有统计学差异,再用方差分量分析把瓶间变异与瓶内变异拆开,判断差异大小是否满足预设的精密度限值。本文演示这两种场景:15 瓶 × 5 次测定的瓶间差,以及不同试剂批测定同一质控品的批间差。
示例数据为确定性的教学数据;统计显著不等于产品不可接受,只能依据方案预先规定的接受限作出判定。
代码
library(ivdtools)
library(readr)18.1 主要函数
| 函数 | 主要用途 | 关键输入或输出 |
|---|---|---|
bottle_anova() |
单因素或嵌套 ANOVA | value ~ batch 或 value ~ batch/vial |
tukey() |
Tukey HSD 两两比较与紧凑字母 | term 指定效应 |
precision() / variance() |
拆分瓶间/瓶内方差分量 | value ~ bottle 单因素 |
18.1.1 ANOVA、Tukey 与方差分量
代码
bottle_anova(data, formula, conf.level = 0.95)
tukey(x, term = NULL, conf.level = NULL)
precision(data, value ~ bottle)
variance(x, method = c("auto", "anova", "reml"))
uncertainty_homogeneity(data, unit = "bottle", result = "value", ...)单因素模型写为
\[ Y_{ij}=\mu+B_i+\varepsilon_{ij}, \]
整体零假设为全部瓶均值相等。经典 F 统计量为
\[ F=\frac{MS_{\mathrm{between}}}{MS_{\mathrm{within}}}. \]
bottle_anova() 删除响应或分组缺失/非有限的行,把右侧变量转为因子,并用 aov() 拟合;嵌套公式如 value ~ lot/vial 展开为 lot 和 lot:vial。对象还给出各组均值 t 区间、Bartlett 方差齐性检验和模型残差 Shapiro–Wilk 检验。这些假设检验同样不能自动选择模型。
tukey() 对指定效应做 Tukey HSD,返回差值、同时区间、调整后 p 值和紧凑字母;默认效应在复杂模型中应核对,必要时显式给出 term。方差分量分析则估计
\[ \sigma_{\mathrm{total}}^2=\sigma_B^2+\sigma_e^2, \]
瓶间 SD 为 \(\sqrt{\sigma_B^2}\),与 ANOVA p 值是不同尺度的结果。均匀性结论应优先围绕预设最大瓶间变异或其不确定度贡献表达。
bottle_anova() 返回 ANOVA 模型和组间诊断;tukey() 返回两两比较表;precision()/variance() 返回方差分量;uncertainty_homogeneity() 返回均匀性标准不确定度组件。
18.2 分析案例
18.2.1 示例一:15 瓶 × 5 次测定的瓶间差
18.2.1.1 读取和核验数据
代码
bottle <- read_csv("./data/018-瓶间差/bottle-difference.csv", show_col_types = FALSE)
bottle <- as.data.frame(bottle)
bottle$bottle <- factor(bottle$bottle)
head(bottle, 6) bottle replicate value
1 B01 1 97.52
2 B01 2 100.53
3 B01 3 98.81
4 B01 4 97.55
5 B01 5 97.85
6 B02 1 100.92
18.2.1.2 第一步:差异显著性检验
单因素模型 value ~ bottle 检验 15 个瓶均值是否相同:
代码
ba <- bottle_anova(bottle, value ~ bottle)
print(ba)
Bottle ANOVA Analysis
Response: value
Group variable: bottle (15 groups)
Complete cases: 75 / 75
Confidence level: 0.95
Descriptive Statistics:
------------------------------------------------------------
bottle n Mean SD
------------------------------------------------------------
B01 5 98.4520 1.2739
B02 5 101.7140 0.9360
B03 5 99.8060 0.8153
B04 5 102.7980 0.4992
B05 5 98.7520 0.8662
B06 5 100.0960 1.4971
B07 5 100.7180 0.8466
B08 5 100.5000 1.3179
B09 5 99.8160 1.0246
B10 5 101.6400 0.8034
B11 5 99.3300 0.9360
B12 5 98.4700 1.0884
B13 5 99.7440 1.1004
B14 5 100.9040 0.4803
B15 5 99.7920 1.4477
ANOVA Table:
------------------------------------------------------------
Term Df Sum Sq Mean Sq F-value p-value Sig
------------------------------------------------------------
bottle 14 107.7233 7.6945 7.1460 0.0000 ***
Residuals 60 64.6057 1.0768
------------------------------------------------------------
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Assumption Checks:
Bartlett test: K-squared = 9.9325, df = 14, p-value = 0.7671 (OK)
Shapiro-Wilk: W = 0.9828, p-value = 0.4011 (OK)
Post-hoc: use tukey() for pairwise comparisons
ANOVA 表给出瓶间 F 检验;assumptions 部分报告 Bartlett 方差齐性与 Shapiro-Wilk 正态性检验。
18.2.1.3 两两比较
代码
posthoc <- tukey(ba)
posthoc
Tukey HSD Post-hoc Test
Response: value
Complete cases: 75 / 75
Confidence level: 0.95
Term: bottle
------------------------------------------------------------
Comparison Diff lwr upr p-value Sig
------------------------------------------------------------
B02-B01 3.2620 0.9412 5.5828 0.0005 ***
B03-B01 1.3540 -0.9668 3.6748 0.7526
B04-B01 4.3460 2.0252 6.6668 0.0000 ***
B05-B01 0.3000 -2.0208 2.6208 1.0000
B06-B01 1.6440 -0.6768 3.9648 0.4526
B07-B01 2.2660 -0.0548 4.5868 0.0624
B08-B01 2.0480 -0.2728 4.3688 0.1408
B09-B01 1.3640 -0.9568 3.6848 0.7431
B10-B01 3.1880 0.8672 5.5088 0.0008 ***
B11-B01 0.8780 -1.4428 3.1988 0.9896
B12-B01 0.0180 -2.3028 2.3388 1.0000
B13-B01 1.2920 -1.0288 3.6128 0.8078
B14-B01 2.4520 0.1312 4.7728 0.0287 *
B15-B01 1.3400 -0.9808 3.6608 0.7656
B03-B02 -1.9080 -4.2288 0.4128 0.2229
B04-B02 1.0840 -1.2368 3.4048 0.9391
B05-B02 -2.9620 -5.2828 -0.6412 0.0026 **
B06-B02 -1.6180 -3.9388 0.7028 0.4794
B07-B02 -0.9960 -3.3168 1.3248 0.9686
B08-B02 -1.2140 -3.5348 1.1068 0.8674
B09-B02 -1.8980 -4.2188 0.4228 0.2298
B10-B02 -0.0740 -2.3948 2.2468 1.0000
B11-B02 -2.3840 -4.7048 -0.0632 0.0384 *
B12-B02 -3.2440 -5.5648 -0.9232 0.0006 ***
B13-B02 -1.9700 -4.2908 0.3508 0.1830
B14-B02 -0.8100 -3.1308 1.5108 0.9952
B15-B02 -1.9220 -4.2428 0.3988 0.2134
B04-B03 2.9920 0.6712 5.3128 0.0022 **
B05-B03 -1.0540 -3.3748 1.2668 0.9507
B06-B03 0.2900 -2.0308 2.6108 1.0000
B07-B03 0.9120 -1.4088 3.2328 0.9852
B08-B03 0.6940 -1.6268 3.0148 0.9990
B09-B03 0.0100 -2.3108 2.3308 1.0000
B10-B03 1.8340 -0.4868 4.1548 0.2778
B11-B03 -0.4760 -2.7968 1.8448 1.0000
B12-B03 -1.3360 -3.6568 0.9848 0.7693
B13-B03 -0.0620 -2.3828 2.2588 1.0000
B14-B03 1.0980 -1.2228 3.4188 0.9330
B15-B03 -0.0140 -2.3348 2.3068 1.0000
B05-B04 -4.0460 -6.3668 -1.7252 0.0000 ***
B06-B04 -2.7020 -5.0228 -0.3812 0.0092 **
B07-B04 -2.0800 -4.4008 0.2408 0.1258
B08-B04 -2.2980 -4.6188 0.0228 0.0549
B09-B04 -2.9820 -5.3028 -0.6612 0.0023 **
B10-B04 -1.1580 -3.4788 1.1628 0.9026
B11-B04 -3.4680 -5.7888 -1.1472 0.0002 ***
B12-B04 -4.3280 -6.6488 -2.0072 0.0000 ***
B13-B04 -3.0540 -5.3748 -0.7332 0.0016 **
B14-B04 -1.8940 -4.2148 0.4268 0.2326
B15-B04 -3.0060 -5.3268 -0.6852 0.0020 **
B06-B05 1.3440 -0.9768 3.6648 0.7619
B07-B05 1.9660 -0.3548 4.2868 0.1854
B08-B05 1.7480 -0.5728 4.0688 0.3515
B09-B05 1.0640 -1.2568 3.3848 0.9470
B10-B05 2.8880 0.5672 5.2088 0.0037 **
B11-B05 0.5780 -1.7428 2.8988 0.9999
B12-B05 -0.2820 -2.6028 2.0388 1.0000
B13-B05 0.9920 -1.3288 3.3128 0.9696
B14-B05 2.1520 -0.1688 4.4728 0.0968
B15-B05 1.0400 -1.2808 3.3608 0.9556
B07-B06 0.6220 -1.6988 2.9428 0.9997
B08-B06 0.4040 -1.9168 2.7248 1.0000
B09-B06 -0.2800 -2.6008 2.0408 1.0000
B10-B06 1.5440 -0.7768 3.8648 0.5576
B11-B06 -0.7660 -3.0868 1.5548 0.9972
B12-B06 -1.6260 -3.9468 0.6948 0.4711
B13-B06 -0.3520 -2.6728 1.9688 1.0000
B14-B06 0.8080 -1.5128 3.1288 0.9953
B15-B06 -0.3040 -2.6248 2.0168 1.0000
B08-B07 -0.2180 -2.5388 2.1028 1.0000
B09-B07 -0.9020 -3.2228 1.4188 0.9866
B10-B07 0.9220 -1.3988 3.2428 0.9837
B11-B07 -1.3880 -3.7088 0.9328 0.7199
B12-B07 -2.2480 -4.5688 0.0728 0.0670
B13-B07 -0.9740 -3.2948 1.3468 0.9739
B14-B07 0.1860 -2.1348 2.5068 1.0000
B15-B07 -0.9260 -3.2468 1.3948 0.9831
B09-B08 -0.6840 -3.0048 1.6368 0.9992
B10-B08 1.1400 -1.1808 3.4608 0.9125
B11-B08 -1.1700 -3.4908 1.1508 0.8956
B12-B08 -2.0300 -4.3508 0.2908 0.1498
B13-B08 -0.7560 -3.0768 1.5648 0.9976
B14-B08 0.4040 -1.9168 2.7248 1.0000
B15-B08 -0.7080 -3.0288 1.6128 0.9988
B10-B09 1.8240 -0.4968 4.1448 0.2859
B11-B09 -0.4860 -2.8068 1.8348 1.0000
B12-B09 -1.3460 -3.6668 0.9748 0.7601
B13-B09 -0.0720 -2.3928 2.2488 1.0000
B14-B09 1.0880 -1.2328 3.4088 0.9374
B15-B09 -0.0240 -2.3448 2.2968 1.0000
B11-B10 -2.3100 -4.6308 0.0108 0.0523
B12-B10 -3.1700 -5.4908 -0.8492 0.0009 ***
B13-B10 -1.8960 -4.2168 0.4248 0.2312
B14-B10 -0.7360 -3.0568 1.5848 0.9982
B15-B10 -1.8480 -4.1688 0.4728 0.2668
B12-B11 -0.8600 -3.1808 1.4608 0.9914
B13-B11 0.4140 -1.9068 2.7348 1.0000
B14-B11 1.5740 -0.7468 3.8948 0.5257
B15-B11 0.4620 -1.8588 2.7828 1.0000
B13-B12 1.2740 -1.0468 3.5948 0.8226
B14-B12 2.4340 0.1132 4.7548 0.0311 *
B15-B12 1.3220 -0.9988 3.6428 0.7819
B14-B13 1.1600 -1.1608 3.4808 0.9015
B15-B13 0.0480 -2.2728 2.3688 1.0000
B15-B14 -1.1120 -3.4328 1.2088 0.9266
------------------------------------------------------------
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Compact Letter Display:
------------------------------------------------------------
B04: a
B02: a
B10: a
B14: a
B07: ab
B08: ab
B06: b
B09: b
B03: b
B15: b
B13: b
B11: b
B05: b
B12: b
B01: b
Tukey HSD 给出全部两两比较与紧凑字母标记(紧凑字母显示(compact letter display)),字母不同的瓶表示均值差异显著。本例瓶数较多,部分瓶间差异显著属正常现象,关键是差异大小是否可接受。
18.2.1.4 第二步:差异大小是否满足预期
用 precision(value ~ bottle) 把总变异拆成瓶间分量(bottle)与瓶内分量(error):
代码
bp <- precision(bottle, value ~ bottle)
bp <- variance(bp)Variance components
sample component VC %Total SD CV[%] Method
----------------------------------------------------
all bottle 1.3236 55.1 1.1505 1.1485 anova
all error 1.0768 44.9 1.0377 1.0359 anova
瓶间 SD ≈ 1.15、CV ≈ 1.15%,瓶内 SD ≈ 1.04、CV ≈ 1.04%。假设方案规定总精密度限值为 2%(CV),则瓶间差异在限值内。注意: 这里的”限值”仅为示例,正式研究必须使用方案或标准规定的限值,且应结合总 CV(含瓶间与瓶内)与各组分分别比较。
代码
bp <- ci(bp)Confidence intervals (SD)
sample component estimate lower upper
----------------------------------------
all total 1.5493 1.2449 2.0519
all bottle 1.1505 0.4254 1.5704
all error 1.0377 0.8807 1.2633
Confidence intervals (%CV)
sample component estimate lower upper
----------------------------------------
all total 1.5467 1.2428 2.0485
all bottle 1.1485 0.4246 1.5678
all error 1.0359 0.8792 1.2612
18.2.2 示例二:不同试剂批测定同一质控品
本例的三个试剂批共同构成一次批间差模型,因此不能按 lot 拆开。如果原始长表还同时包含多个产品或分析物,应先用 split_by_factors(data, variable="product") 按独立评价项目拆分,再在每个子数据框内保留全部 lot,分别调用 bottle_anova() 或 precision()。
18.2.2.1 读取数据
3 个试剂批(L1、L2、L3)分别对同一质控品测定 15 次:
代码
lot_data <- read_csv("./data/018-瓶间差/reagent-lot.csv", show_col_types = FALSE)
lot_data <- as.data.frame(lot_data)
lot_data$lot <- factor(lot_data$lot)18.2.2.2 批次间差异检验
代码
ba_lot <- bottle_anova(lot_data, value ~ lot)
print(ba_lot)
Bottle ANOVA Analysis
Response: value
Group variable: lot (3 groups)
Complete cases: 45 / 45
Confidence level: 0.95
Descriptive Statistics:
------------------------------------------------------------
lot n Mean SD
------------------------------------------------------------
L1 15 100.8173 1.1626
L2 15 101.2760 1.4271
L3 15 100.5673 0.7829
ANOVA Table:
------------------------------------------------------------
Term Df Sum Sq Mean Sq F-value p-value Sig
------------------------------------------------------------
lot 2 3.8754 1.9377 1.4528 0.2454
Residuals 42 56.0183 1.3338
------------------------------------------------------------
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Assumption Checks:
Bartlett test: K-squared = 4.6245, df = 2, p-value = 0.0990 (OK)
Shapiro-Wilk: W = 0.9817, p-value = 0.6870 (OK)
Post-hoc: use tukey() for pairwise comparisons
代码
tukey(ba_lot)
Tukey HSD Post-hoc Test
Response: value
Complete cases: 45 / 45
Confidence level: 0.95
Term: lot
------------------------------------------------------------
Comparison Diff lwr upr p-value Sig
------------------------------------------------------------
L2-L1 0.4587 -0.5659 1.4832 0.5269
L3-L1 -0.2500 -1.2745 0.7745 0.8247
L3-L2 -0.7087 -1.7332 0.3159 0.2245
------------------------------------------------------------
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Compact Letter Display:
------------------------------------------------------------
L2: a
L1: a
L3: a
紧凑字母标记显示哪些批之间均值差异显著、差异方向如何。批间均值差异即使统计显著,仍需结合预设的可接受偏差(例如质控品允许的靶值范围)判断是否在可接受范围内,不能仅凭 p 值下结论。
18.2.3 标准案例:YY/T 1709—2020 校准品均匀性
附录 A 从一批 CRP 产品校准品中抽取 14 瓶,每瓶测量 3 次,以单因素 ANOVA 分离瓶间和瓶内变异。
代码
yyt1709_homogeneity <- read_csv(
"./data/018-瓶间差/std-yyt1709-homogeneity.csv",
show_col_types = FALSE
)
yyt1709_homogeneity$bottle <- factor(yyt1709_homogeneity$bottle)
yyt1709_fit <- ivdtools::bottle_anova(
yyt1709_homogeneity, result ~ bottle, conf.level = 0.95
)
yyt1709_fit
Bottle ANOVA Analysis
Response: result
Group variable: bottle (14 groups)
Complete cases: 42 / 42
Confidence level: 0.95
Descriptive Statistics:
------------------------------------------------------------
bottle n Mean SD
------------------------------------------------------------
1 3 7.1433 0.1415
2 3 7.2333 0.1301
3 3 7.0433 0.1210
4 3 7.1333 0.1106
5 3 7.0233 0.0723
6 3 7.1133 0.0987
7 3 7.2533 0.1069
8 3 7.2033 0.0462
9 3 7.2500 0.0872
10 3 7.1767 0.1115
11 3 7.0767 0.1234
12 3 7.1433 0.1115
13 3 7.1633 0.1563
14 3 7.0800 0.1153
ANOVA Table:
------------------------------------------------------------
Term Df Sum Sq Mean Sq F-value p-value Sig
------------------------------------------------------------
bottle 13 0.2114 0.0163 1.2807 0.2809
Residuals 28 0.3556 0.0127
------------------------------------------------------------
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Assumption Checks:
Bartlett test: K-squared = 3.1688, df = 13, p-value = 0.9973 (OK)
Shapiro-Wilk: W = 0.9532, p-value = 0.0839 (OK)
Post-hoc: use tukey() for pairwise comparisons
代码
yyt1709_u <- ivdtools::uncertainty_homogeneity(
yyt1709_homogeneity, unit = "bottle", result = "result"
)
yyt1709_u
Between-unit homogeneity uncertainty
Units : 14
Results : 42
Mean : 7.1454762
Effective degrees of freedom : 3
Average replicates : 3
Within-unit mean square : 0.0127
Between-unit mean square : 0.016264652
Within-unit degrees of freedom : 28
Between-unit degrees of freedom : 13
Repeatability SD : 0.11269428
Between-unit SD : 0.034470529
Hidden inhomogeneity : 0.033636378
Between-unit uncertainty : 0.034470529
Relative repeatability : 0.015771416
Relative between-unit : 0.0048241052
Relative hidden inhomogeneity : 0.0047073669
Relative between-unit uncertainty : 0.0048241052
| 标准位置 | 统计量 | 标准结果 | ivdtools结果 |
|---|---|---|---|
| YY/T 1709 Appendix A | SS_between | 0.211440 | 0.211440 |
| YY/T 1709 Appendix A | SS_within | 0.355600 | 0.355600 |
| YY/T 1709 Appendix A | MS_between | 0.016265 | 0.016265 |
| YY/T 1709 Appendix A | MS_within | 0.012700 | 0.012700 |
| YY/T 1709 Appendix A | F | 1.280681 | 1.280681 |
| YY/T 1709 Appendix A | u_bb | 0.034492 | 0.034471 |