18  瓶间差分析

本文用瓶间差和批次间比较判断物料的均匀性,并演示如何把均匀性结果转化为不确定度分量。显著性检验只说明是否观察到差异,不能单独说明差异是否具有实际影响;最终结论还要结合重复性、目标值和预设可接受标准。

重要使用边界

bottle_anova() 的 F 检验回答“组均值是否完全相同”,不回答差异是否满足产品均匀性限。tukey() 只在整体模型和假设适用时定位成对差异。均匀性不确定度或瓶间 SD 必须由与抽样、测量顺序和重复设计一致的方差模型估计,再与预设技术限比较;不能只凭 p 值放行或拒绝。

分析流程:从瓶间差异到均匀性不确定度

%%{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

瓶间差(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
YY/T 1709 附录 A:ivdtools 与标准原文结果对照
标准位置 统计量 标准结果 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