当 Wilcoxon 检验返回一些 0 p 值时该怎么办？

Question

我正在一个大型列表中使用 R 执行 Wilcoxon 测试（包含 86 个数据帧，行数和值可变）。我不明白为什么在 p-values

我在这里报告我的脚本：

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 说：

请注意，在大多数平台上，正值小于 '.Machine$double.xmin' 可能会出现。在典型的 R 平台上 最小的正双精度约为'5e-324'。

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

虽然如果您非常需要这种能力，可能可以从头开始构建，但研究人员通常只是将这种 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()