R t.test 中的 P 值触底

Question

我有大约 8000 个基因的三个向量，每个向量都有相关的突变频率，例如：

alpha = c(0.84, 0.87, 0.91...)
beta = c(0.97, 0.94, 0.99...)
kappa = c(0.72, 0.68, 0.75...)

我正在使用 R 中的 t.test 函数来计算这些向量之间的 P 值。由于有多少个基因以及每个向量之间的差异有多一致，我的 P 值总是非常低。例如，当我跑步时：

ab = t.test(x = alpha, y = beta)
ak = t.test(x = alpha, y = kappa)
bk = t.test(x = beta, y = kappa)

所有测试的摘要输出相同的值 2.2e-16，这是我不喜欢的。当我直接提取P值时：

t.test(alpha, kappa)$p.value

R 只是输出 0。我假设这是由于一些内部浮点/双精度大小限制，但它看起来不太好。我们正在发布这些数据，无论我们采用相同的汇总 P 值还是 0，看起来都不好。 R 中有没有办法解决这个问题，以便我可以计算任意低的 P 值？如果另一个工具效果更好，我也会对此感兴趣。谢谢！

Answer 1

如果您的 t 统计量确实如此之大，以至于您的 p 值下溢为零（即 p 值小于大约 1e-308 (!!)），您可以从 this answer 调整机制以获得p 值...

x <- rep(1:5, 100)
y <- rep(10001:10005, 100)
tt <- t.test(x,y)

结果：

    Welch Two Sample t-test

data:  x and y
t = -111692, df = 998, p-value < 2.2e-16
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -10000.176  -9999.824
sample estimates:
mean of x mean of y 
        3     10003 

tt$p.value  ## 0

pvalue.extreme <- function(t, df) {
   log.pvalue <- log(2) + pt(abs(t), df = df, 
            lower.tail = FALSE, log.p = TRUE)
   log10.pvalue <- log.pvalue/log(10) ## from natural log to log10
   mantissa <- 10^(log10.pvalue %% 1)
   exponent <- log10.pvalue %/% 1
   ## or return(c(mantissa,exponent))
   return(sprintf("p value is %1.2f times 10^(%d)",mantissa,exponent))
}

pvalue.extreme(tt$statistic, tt$parameter)
[1] "p value is 1.11 times 10^(-3543)"