为什么我无法获得小于 2.2e-16 的 p 值？

Question

我在 R 中的 t 检验和卡方中发现了这个问题，但我认为这个问题通常适用于其他测试。如果我这样做：

a <- 1:10
b <- 100:110
t.test(a,b) 

我得到：t = -64.6472，df = 18.998，p值< 2.2e-16。我从评论中知道 2.2e-16 是 .Machine$double.eps 的值 - 满足 1 + x != 的最小浮点数1，但当然 R 可以表示比它小得多的数字。我还从 R FAQ 中知道 R 必须将浮点数舍入为 53 位二进制数字精度：R 常见问题解答。

有几个问题：(1) 我是否正确地将其读取为 精度 的 53 个二进制数字，或者是 R 中的值？ .Machine$double.eps 计算不准确？ (2) 为什么在进行此类计算时，即使精度有所损失，R 也不提供显示较小 p 值的方法？ (3) 有没有办法显示更小的 p 值，即使我失去了一些精度？对于单个测试，2 位小数有效数字就可以了，对于我要 Bonferroni 纠正的值，我需要更多。当我说“失去一些精度”时，我想< 53 个二进制数字，但是 (4) 我完全错了，任何 p 值 < .Machine$double.eps 非常不准确？ (5) R 是否只是诚实的，而其他统计数据包则不然？

在我的领域，非常小的 p 值是常态，一些例子：http://www. ncbi.nlm.nih.gov/pubmed/20154341，http://www.plos Genetics.org/article/info%3Adoi %2F10.1371%2Fjournal.pgen.1002215 这就是为什么我想代表这么小的p 值。

感谢您的帮助，抱歉问了这么曲折的问题。

Answer 1

尝试类似 t.test(a,b)$p.value 的操作，看看是否能够满足您所需的准确性。我相信它与结果的打印有关，而不是与实际存储的计算机值有关，后者应该具有必要的精度。

Answer 2

一些 R 包解决了这个问题。最好的方法是通过 pspearman 包。

source("http://www.bioconductor.org/biocLite.R")
biocLite("pspearman")
library("pspearman")
a=c(1:110,110)
b=1:111
out <- spearman.test(a, b, alternative = "greater", approximation="t-distribution")
out$p.value

[1] 3.819961e-294

Answer 3

最近有同样的问题。统计学家研究员建议：

A <- cor.test(…)
p <- 2* pt(A$statistic,  df = A$parameter, lower.tail=FALSE)

Answer 4

这是一个流行的问题，但令人惊讶的是没有提到使用对数表示作为解决方案的答案。

在一些研究领域，特别是生物信息学（尤其是基因组学，但在其他组学领域越来越多）精确的 log10（p 值）被用来将证据与无效证据进行比较。通过将 log.p=TRUE 传递给适当的分位数分布函数，可以在 R 中获得 p 值的对数以进行常见测试。

t 测试

a = 1:10
b = 110:120

log10_t_test = function(...) {
    model = t.test(...)
    # note: you need to modify below if passing `alternative` arg
    log_e_p = log(2) + unname(pt(abs(model$statistic), df=model$parameter, lower.tail=FALSE, log.p=TRUE))
    model$log10_pvalue = log_e_p / log(10)
    model
}

model = log10_t_test(a, b)
model$log10_pvalue  # gives  -24.6699

，您可以根据 log10(p) 的朴素计算进行评估：

t(sapply(seq(2, 7), function(order_of_magnitude) {
    n = 10 ** order_of_magnitude
    a = rnorm(n, mean=0)
    b = rnorm(n, mean=0.05)
    model = log10_t_test(a, b)
    c(
        proper_log10p=model$log10_pvalue,
        naive_log10p=log10(model$p.value)
    )
}))

proper_log10p	naive_log10p
-0.2239816	-0.2239816
-1.4719561	-1.4719561
-0.7009232	-0.7009232
-30.7052283	-30.7052283
-250.8593000	-250.8593000
-2737.2759952	-Inf

Correlation

log10_cor_test = function(x, y, ..., method='pearson') {
    model = cor.test(x, y, ..., method=method)
    if (method == 'spearman') {
        r = model$estimate
        n = length(x)  # note: this assumes no missing values
        statistic = r * sqrt((n - 2) / (1 - r**2))
        df = n - 2
    } else {
        statistic = model$statistic
        df = model$parameter
    }
    log_e_p = log(2) + unname(pt(abs(statistic), df=df, lower.tail=FALSE, log.p=TRUE))
    model$log10_pvalue = log_e_p / log(10)
    model
}

Pearson

方便地，这也使用 t 统计量，并且可以直接从 cor.test 结果中提取统计量以及自由度参数。

比较：

t(sapply(seq(2, 7), function(order_of_magnitude) {
    n = 10 ** order_of_magnitude
    a = seq(1, n)
    b = a + rnorm(n) * 10**7 # add strong noise
    model = log10_cor_test(a, b)
    c(
        proper_log10p=model$log10_pvalue,
        naive_log10p=log10(model$p.value)
    )
}))

proper_log10p	naive_log10p
-2.435782e+00	-2.4357824
-6.848772e-01	-0.6848772
-1.073320e-01	-0.1073320
-7.630908e-01	-0.7630908
-1.815829e+02	-181.5829491
-1.735492e+05	-Inf

Spearman

这需要更多的手动工作，因为我们需要手动计算自由度 (n-2) 和统计数据。

如果您对 t 分布近似感到满意，则可以使用以下公式计算检验统计量：r * sqrt((n - 2) / (1 - r**2)) 公式并重复使用相同的 pt 函数。

比较：

t(sapply(seq(2, 7), function(order_of_magnitude) {
    n = 10 ** order_of_magnitude
    a = seq(1, n)
    b = a + rnorm(n) * 10**7 # add strong noise
    model = log10_cor_test(a, b, method='spearman')
    c(
        proper_log10p=model$log10_pvalue,
        naive_log10p=log10(model$p.value)
    )
}))

proper_log10p	naive_log10p
-6.525561e-02	-0.06535951
-3.388555e-01	-0.33891807
-7.684660e-03	-0.00768466
-1.337620e-02	-0.01337620
-1.874304e+02	-187.43039010
-1.682590e+05	-Inf

Answer 5

在这里交换答案和评论时，我对一些事情感到困惑。

首先，当我尝试OP的原始示例时，我没有得到像这里讨论的那样小的p值（几个不同的2.13.x版本和R-devel）：

a <- 1:10
b <- 10:20
t.test(a,b)
## data:  a and b 
## t = -6.862, df = 18.998, p-value = 1.513e-06

第二，当我使组之间的差异变得更大时，我实际上得到了@eWizardII建议的结果：

a <- 1:10
b <- 110:120
(t1 <- t.test(a,b))
# data:  a and b 
# t = -79.0935, df = 18.998, p-value < 2.2e-16
#
> t1$p.value
[1] 2.138461e-25

t.test中打印输出的行为是由它的调用驱动的stats:::print.htest （也被其他统计测试函数调用，例如 chisq.test，如 OP 所示），它又调用 format.pval，它呈现的 p 值小于其 eps 值（即 .Machine$double.eps）默认）为 <每股收益。我很惊讶地发现自己不同意这些普遍精明的评论者......

最后，尽管担心非常小的p值的精确值似乎很愚蠢，但OP是正确的，这些值是在生物信息学文献中经常被用作证据强度的指数——例如，人们可能会测试 100,000 个候选基因并查看结果p值的分布（搜索“火山图”作为一个例子这种 程序）。

Answer 6

两个问题：

1) 1e-16 和 1e-32 的 p 值之间的统计含义可能存在什么差异？如果您确实可以证明它的合理性，那么使用记录的值就是正确的方法。

2）当您对 R 的数值准确性感兴趣时，为什么使用维基百科？

R-FAQ 表示“其他[即非整数]数字必须四舍五入到（通常）53 位二进制数字的精度。” 16 位数字已经是极限了。这是在控制台时获得精度限制的方法：

> .Machine$double.eps
[1] 2.220446e-16

当在 [0,1] 范围内解释时，该数字实际上为零

Answer 7

您链接到的维基百科页面适用于 R 不使用的 Decimal64 类型 - 它使用标准发行的双精度数。

首先，来自 .Machine 帮助页面的一些定义。

double.eps：最小的正浮点数“x”，使得
'1 + x！= 1'。 ...通常为“2.220446e-16”。

double.xmin：最小的非零标准化浮点数
...通常为“2.225074e-308”。

因此，您可以表示小于 2.2e-16 的数字，但它们的精度会降低，并且会导致计算问题。尝试一些数字接近最小可表示值的示例。

2e-350 - 1e-350
sqrt(1e-350)

您在评论中提到您想要进行 bonferroni 修正。我建议您使用 p.adjust(your_p_value, method = "bonferroni") 而不是为此滚动自己的代码。 pairwise.t.test 使用这个。