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

尝试这样的事情，看看这是否能为您提供所需的准确性。我相信它与结果的打印有关，而不是实际存储的计算机值，而实际存储的计算机值应该具有必要的精度。t.test(a,b)$p.value

两个问题：

1） 1e-16 和 1e-32 的 p 值之间在统计意义上可能存在什么差异？如果您真的可以证明它的合理性，那么使用记录的值是要走的路。

2）当你对R的数值精度感兴趣时，你为什么要使用维基百科？

R-FAQ说：“其他[表示非整数]数字必须四舍五入到（通常）53位二进制数字的精度。16 位数字大约是限制。这是在控制台上获得准确性限制的方法：

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

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

您链接到的维基百科页面是 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

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

首先，当我尝试 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

打印输出 in 的行为由其调用（也由其他统计测试函数调用，如 OP 所指出的）驱动，而 又调用 ，它表示的 p 值小于其值 （默认为）为 。我很惊讶地发现自己不同意这些普遍精明的评论者......t.teststats:::print.htestchisq.testformat.pvaleps.Machine$double.eps< eps

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

某些 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

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

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

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

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

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)
    )
}))

相关

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
}

皮尔森

方便的是，这也使用了 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)
    )
}))

斯皮尔曼

这需要更多的手动工作，因为我们需要手动计算自由度 （） 和统计量。n-2

如果您对 t 分布近似感到满意，则可以使用 公式 计算检验统计量 formula，并重用相同的函数。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)
    )
}))