代码
library(ivdtools)
library(readr)瓶间差(bottle-to-bottle variability)评价同一产品不同瓶(或不同批次)之间测量结果的差异。分析通常分两步:先用 bottle_anova() 检验各瓶(或批次) 均值是否有统计学差异,再用方差分量分析把瓶间变异与瓶内变异拆开,判断差异大小是否满足预设的精密度 限值。本文演示这两种场景:15 瓶 × 5 次测定的瓶间差,以及不同试剂批测定同一质控品的批间差。
示例数据为确定性的教学数据;统计显著不等于产品不可接受,只能依据方案预先规定的接受限作出判定。
library(ivdtools)
library(readr)| 函数 | 主要用途 | 关键输入或输出 |
|---|---|---|
bottle_anova() |
单因素或嵌套 ANOVA | value ~ batch 或 value ~ batch/vial |
tukey() |
Tukey HSD 两两比较与紧凑字母 | term 指定效应 |
precision() / variance() |
拆分瓶间/瓶内方差分量 | value ~ bottle 单因素 |
bottle <- read_csv("./data/bottle-difference.csv", show_col_types = FALSE)
bottle <- as.data.frame(bottle)
bottle$bottle <- factor(bottle$bottle)
str(bottle)'data.frame': 75 obs. of 3 variables:
$ bottle : Factor w/ 15 levels "B01","B02","B03",..: 1 1 1 1 1 2 2 2 2 2 ...
$ replicate: num 1 2 3 4 5 1 2 3 4 5 ...
$ value : num 97.5 100.5 98.8 97.5 97.8 ...
dim(bottle)[1] 75 3
单因素模型 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 正态性检验。
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)),字母不同 的瓶表示均值差异显著。本例瓶数较多,部分瓶间差异显著属正常现象,关键是差异大小是否可接受。
用 precision(value ~ bottle) 把总变异拆成瓶间分量(bottle)与瓶内分量(error):
bp <- precision(bottle, value ~ bottle)
bp <- variance(bp)Variance components
sample component VC %Total SD CV[%]
---------------------------------------------
all bottle 1.3236 55.1 1.1505 1.1485
all error 1.0768 44.9 1.0377 1.0359
瓶间 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
3 个试剂批(L1、L2、L3)分别对同一质控品测定 15 次:
lot_data <- read_csv("./data/reagent-lot.csv", show_col_types = FALSE)
lot_data <- as.data.frame(lot_data)
lot_data$lot <- factor(lot_data$lot)
str(lot_data)'data.frame': 45 obs. of 3 variables:
$ lot : Factor w/ 3 levels "L1","L2","L3": 1 1 1 1 1 1 1 1 1 1 ...
$ replicate: num 1 2 3 4 5 6 7 8 9 10 ...
$ value : num 100 102 101 102 101 ...
lot_summary <- aggregate(value ~ lot, lot_data, function(x) c(n = length(x), mean = mean(x), sd = sd(x)))
lot_summary <- data.frame(lot = lot_summary$lot, lot_summary$value)
lot_summaryba_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 值下结论。