当Z得分很大时(pvalue远低于零),如何计算R中z分数的p值？

Question

在遗传学中,非常小的p值是常见的(例如10 ^ -400),当我在R中得分很大时,我正在寻找一种获得非常小的p值(双尾)的方法,例如:

z=40
pvalue = 2*pnorm(abs(z), lower.tail = F)

这给了我零而不是非常小的值,这是非常重要的.

Answer 1

无法处理小于约10 ^( - 308)(.Machine$double.xmin)的p值并不是R的错误,而是使用双精度(64位)浮点数存储数字信息的任何计算系统的一般限制.

通过计算对数刻度来解决问题并不难,但是你不能将结果存储为R中的数值; 相反,您需要将结果存储(或打印)为尾数加指数.

pvalue.extreme <- function(z) {
   log.pvalue <- log(2) + pnorm(abs(z), 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(5)
## [1] "p value is 5.73 times 10^(-7)"
2*pnorm(5,lower.tail=FALSE)
## [1] 5.733031e-07

更极端:

pvalue.extreme(40)
## [1] "p value is 7.31 times 10^(-350)"

在R(Brobdingnag,Rmpfr,...)中有各种各样的包处理极大/小数字并具有扩展的精度.例如,

2*Rmpfr::pnorm(mpfr(40, precBits=100), lower.tail=FALSE, log.p = FALSE)
## 1 'mpfr' number of precision  100   bits 
## [1] 7.3117870818300594074979715966414e-350

但是,在使用任意精度系统时,您将在计算效率和方便性方面付出巨大代价.