2

如何处理 R 中的 p 值?

我期待非常低的 p 值,例如:

1.00E-80

我需要-log10

-log10(1.00E-80)

-log10(0) 是 Inf,但 Inf 也具有舍入感。

但似乎在 1.00E-308 之后,R 产生 0。

1/10^308  
[1] 1e-308

 1/10^309 
[1] 0

功能与截止点1e-308相同的p值显示的准确性lm,还是只是设计成我们需要一个截止点,我需要考虑一个不同的截止点-例如1e-100(对于例如)用 <1e-100 替换 0。

4

2 回答 2

8

有多种可能的答案——哪一个最有用取决于上下文:

  • 在通常情况下,R 确实无法存储比 更接近零的浮点值.Machine$double.xmin,后者因平台而异,但通常(如您所发现的)大约为1e-308. 如果您真的需要处理这么小的数字并且找不到直接处理对数刻度的方法,则需要搜索 Stack Overflow 或 R wiki 以获取处理任意/扩展精度值的方法(但您可能应该尝试在对数尺度上工作——这将不那么麻烦)
  • 在许多情况下,R 实际上在内部计算(自然)对数标度上的 p 值,并且可以在请求时返回对数值,而不是在给出答案之前对它们求幂。例如,dnorm(-100,log=TRUE)给出 -5000.919。log10您可以通过除以 =-2171 直接转换为 log10 比例(无需取幂,然后使用 )log(10),因为dnorm(-100,log=TRUE)/log(10)它太小而无法用浮点表示。对于p***(累积分布函数)函数,使用log.p=TRUE而不是log=TRUE. (这一点在很大程度上取决于您的特定上下文。即使您没有使用内置的 R 函数,您也可以找到一种方法来提取对数刻度上的结果。)
  • <2.2e-16在某些情况下,即使已知更精确的值,R 也会将 p 值结果呈现为:(t1 <- t.test(rnorm(10,100),rnorm(10,80)))

印刷

....
t = 56.2902, df = 17.904, p-value < 2.2e-16

但您仍然可以从结果中提取精确的 p 值

> t1$p.value
[1] 1.856174e-18

(在许多情况下,此行为由format.pval()函数控制)

说明所有这些将如何工作lm

d <- data.frame(x=rep(1:5,each=10))
set.seed(101)
d$y <- rnorm(50,mean=d$x,sd=0.0001)
lm1 <- lm(y~x,data=d)

summary(lm1)将斜率的 p 值打印为<2.2e-16,但如果我们使用coef(summary(lm1))(不使用 p 值格式),我们可以看到该值为 9.690173e-203。

一个更极端的例子:

set.seed(101); d$y <- rnorm(50,mean=d$x,sd=1e-7)
lm2 <- lm(y~x,data=d)
coef(summary(lm2))

表明 p 值实际上已下溢为零。但是,我们仍然可以在对数尺度上得到答案:

tval <- coef(summary(lm2))["x","t value"]
2*pt(abs(tval),df=48,lower.tail=FALSE,log.p=TRUE)/log(10)

给出 -692.62(您可以使用前面的示例检查此方法,其中 p 值不会溢出,并看到您得到的答案与摘要中打印的相同)。

于 2012-07-04T12:38:55.903 回答
2

小数字通常很难处理。

R 中无限的限制是由使用双精度浮点引起的:

?double 所有 R 平台都需要使用符合 IEC 60559(也称为 IEEE 754)标准的值。这基本上以 53 位的精度工作,并代表从大约 2e-308 到 2e+308 的绝对值范围。

http://en.wikipedia.org/wiki/Double_precision_floating-point_format

您可能会发现Rmpfr包在这里很有帮助,因为它允许您创建多个精度数字。

install.packages("Rmpfr")
require(Rmpfr)

log(mpfr(1/10^309, precBits=500))
于 2012-07-04T13:03:42.193 回答