我很困惑在调整多个测试的 p 值时是否应该排除缺失的观察结果。R 函数之间似乎没有就是否这样做达成共识。stats::p.adjust(x)如果n = length(x)显式指定默认值(计算 NA)与未显式指定默认值(不计算 NA),则行为不同。multtest::rawp2adjp(x)计算 NA。什么是正确的行为?
编辑:NA评论中要求澄清此数据的含义。在混合效应模型拟合以识别数据中的异常值后,按残差计算 p 值。实验过程复杂,由许多实验者并行进行,因此可能出现错误。显着性识别数据点,这些数据点出乎意料地远离拟合值,因为残差是居中且正态分布的 -> 可能是异常值[ref]。在拟合之前必须删除一些观测值,例如,因为某个组的观测值太少,这会导致模型拟合出现问题,或者由于实验失败而在分析之前已经丢失。
MWE
## Generate some p values and compare the three possibilities of
## adjusting for multiple testing
n <- 10000
x <- pmin(rexp(n,rate =1/0.01), 1) # Generate some p values
x[sample(c(F,T), n/10, TRUE)] <- NA # delete some observations
# Three different methods of p value calculation
x1 <- p.adjust(x, method = 'holm')
x2 <- p.adjust(x, method = 'holm', n = length(x))
x3 <- multtest::mt.rawp2adjp(x, proc = 'Holm')
x3 <- x3$adjp[order(x3$index),"Holm"]
# Compare p.adjust and mt.rawp2adjp
par(mfrow=c(1,2))
plot(x1, x3); title('p.adj(x) \n vs. mt.rawp2adjp')
plot(x2, x3); title('p.adj(x, len = length(x)) \n vs. mt.rawp2adjp')
附录:为什么stats::p.adjust会有这样的行为?
p.adjust源代码的开头,R 3.4.3:
在头部,n定义为length(p),然而,R 在需要参数之前不会评估它们
function (p, method = p.adjust.methods, n = length(p))
{
method <- match.arg(method)
if (method == "fdr")
method <- "BH"
nm <- names(p)
p <- as.numeric(p)
p0 <- setNames(p, nm)
if (all(nna <- !is.na(p)))
nna <- TRUE
在这里,p 被剥夺了所有NA的 s,n直到此时都不需要
p <- p[nna]
lp <- length(p)
现在,n第一次使用,这意味着length(p)现在才评估。因此,如果保留默认设置,length(p[!is.na(p)])则计算。
stopifnot(n >= lp)
[ remaining source code omitted ]
}
