当 Wilcoxon 测试 returns 一些 0 p 值时该怎么办？

Question

我在一个大列表中使用 R 执行 Wilcoxon 测试（包含 86 个数据帧，行数和值数可变）。我不明白为什么在 p 值 < 10^-307 之后它 returns 一些 p 值 = 0 然后它 returns 具有“正常”值。这种情况发生在具有非常多行数（大约 4000 或更多）的数据帧中，但不会出现在具有大约 1000 行的数据帧中。我的脚本是错误的还是有可能得到 0 个值？如果是这样，是否有特定的方式来解释它们？

这里我报告我的脚本：

for (i in 1 : length(List)) {   
    for(j in 2 : (nrow(List[[i]]) - 1)){
        divider <- (List[[i]][j,2])
        ValueInfe <- List[[i]][List[[i]][,2] < divider ,]
        ValueSupUgu <- List[[i]][List[[i]][,2] >= divider ,]
        if(j == 2){Num_ValoreInfe <- as.numeric(ValoreInfe[2])}
        if(j!= 2){Num_ValoreInfe <- as.numeric(ValoreInfe[,2])}
        Num_ValoreSupUgu <- as.numeric(ValoreSupUgu[,2])
        b <- wilcox.test(Num_ValoreInfe, Num_ValoreSupUgu)
        List2.0[[i]][j,3] <- b$p.value
    }

}

这是我的结果示例：

0.000000e+00
8.343024e-02
1.435822e-02
2.716505e-03
5.370877e-04
1.089895e-04
2.250558e-05
4.706192e-06
9.936437e-07
2.114061e-07
4.526195e-08
9.741929e-09
2.106339e-09
4.572291e-10
9.960156e-11
9.960156e-11
[...]
0
0
0
0
0
[...]
2.114061e-07
9.936437e-07
4.706192e-06
2.250558e-05
1.089895e-04
5.370877e-04
2.716505e-03
1.435822e-02
0.000000e+00

Answer 1

通常，R 可以区分零的最小数字约为 1e-308（即 10^(-308)） - 具体来说，.Machine$double.xmin=2.225074e-308。更准确地说，R 可以处理稍微较小的值：?Machine 说：

Note that on most platforms smaller positive values than ‘.Machine$double.xmin’ can occur. On a typical R platform the smallest positive double is about ‘5e-324’.

如果你想处理比这更小的数字，你必须做一些聪明的事情，比如记录它们的对数（log(.Machine$double.xmin) 是 -708，你可以很容易地记录这样小得多的数字）。 R 中的一些 p 值计算允许您检索 log-p 值而不是 p 值，但 Wilcoxon 检验没有这样的功能。

虽然如果您非常需要这种能力，也许可以从头开始构建这种能力，但研究人员通常只会将此类 p 值视为“极小”；如果需要，您可以说“<1e-308” .我见过的唯一一个人们担心如此小的 p 值的精确值的研究领域是生物信息学，其中 p 值本身被视为度量标准，而不是统计清晰度的定性指标差异。

这是一个小示例，它测试样本量逐渐增大的非重叠集的 p 值，显示 p 值减小然后下溢到零（请参见位于下方 y 轴上的点图的右边缘）：

w <- function(n=20) {
    wilcox.test(1:n,1e6+1:n)$p.value
}
nvec <- seq(20,1000,by=10)
pvec <- sapply(nvec,w)

破解 log-p 值

深入研究 stats:::wilcox.test.default 中的代码，我们可以找到根据测试统计和分组样本大小计算 p 值的位置，并使用 log.p=TRUE 重新计算它们.下面的代码跳过了一些细节，例如考虑关系和允许不同的备选假设（即假设是双侧测试）。

这为您提供了 p 值的自然对数；您可以通过乘以 log10(exp(1)) ...

转换回 log10

wilcox_log_p <- function(x,y,exact=FALSE,correct=TRUE,...) {
    ## assume two-sided
    w <- wilcox.test(x,y,...)
    n.x <- length(x)
    n.y <- length(y)
    STATISTIC <- w$statistic
    if (exact) {
        if (STATISTIC > (n.x * n.y/2)) {
            return(pwilcox(STATISTIC - 1, n.x, n.y, 
                   lower.tail = FALSE, log.p=TRUE))
        }
        return(pwilcox(STATISTIC, n.x, n.y, log.p=TRUE))
    } else {
        NTIES <- 0 ## assume no ties!
        z <- STATISTIC - n.x * n.y/2
        SIGMA <- sqrt((n.x * n.y/12) * ((n.x + n.y + 1) - 
                 sum(NTIES^3 - NTIES)/((n.x + n.y) * (n.x + n.y - 
                  1))))
            if (correct) {
                CORRECTION <- sign(z) * 0.5
            }
            z <- (z - CORRECTION)/SIGMA
            PVAL <-  log(2) + min(pnorm(z, log.p=TRUE), 
                             pnorm(z, lower.tail = FALSE, log.p=TRUE))
        return(PVAL)
    }
}

w <- function(n=20) {
    wilcox.test(1:n,1e6+1:n, exact=FALSE)$p.value
}
w2 <- function(n=20) {
    wilcox_log_p(1:n,1e6+1:n)
}
nvec <- seq(20,1100,by=10)
pvec <- sapply(nvec,w)
pvec2 <- sapply(nvec,w2)

dd <- data.frame(n=rep(nvec,2),p=c(log(pvec),pvec2),
                 method=rep(c("default","log_p"),each=length(nvec)))
library(ggplot2); theme_set(theme_bw())
ggplot(dd, aes(n,p,colour=method)) + geom_point() + geom_line()
    scale_x_log10()

当 Wilcoxon 测试 returns 一些 0 p 值时该怎么办？

What to do when Wilcoxon test returns some 0 p-values?

r

statistical-test

破解 log-p 值