使用R包BACON提取P值时pnorm()返回0的问题及解决求助

问题根源 #

当检验统计量的绝对值超过约37.5时，pnorm()计算出的单侧概率会小于R双精度浮点数的最小可表示正数（~2.2e-308），直接返回0，导致结果丢失极小P值的有效信息。

解决方案 #

以下是三种可行的解决方法，可根据需求选择：

方法1：修改P值提取函数，避免下溢

重写原函数，通过对数P值计算再转换，对溢出为0的结果替换为R可表示的最小正数值：

get_bacon_pvalues <- function(object, corrected = TRUE) {
    if (!corrected | any(is.na(bias(object))) | any(is.na(inflation(object)))) {
        # 先计算对数形式的单侧P值
        log_p <- pnorm(-abs(object@teststatistics), log.p = TRUE)
        pvalues <- 2 * exp(log_p)
    } else {
        teststatistics <- t(t(object@teststatistics) - bias(object))
        teststatistics <- t(t(teststatistics)/inflation(object))
        log_p <- pnorm(-abs(teststatistics), mean = 0, sd = 1, log.p = TRUE)
        pvalues <- 2 * exp(log_p)
    }
    # 将溢出为0的P值替换为R能表示的最小正双精度数
    pvalues[pvalues == 0] <- .Machine$double.xmin
    return(pvalues)
}

方法2：直接使用对数P值分析

对数P值不会出现下溢问题，且适用于绝大多数后续分析（如多重检验校正、曼哈顿图绘制）：

get_bacon_logp <- function(object, corrected = TRUE) {
    if (!corrected | any(is.na(bias(object))) | any(is.na(inflation(object)))) {
        # 双侧检验的对数P值 = log(2) + 单侧对数P值
        log_pvalues <- log(2) + pnorm(-abs(object@teststatistics), log.p = TRUE)
    } else {
        teststatistics <- t(t(object@teststatistics) - bias(object))
        teststatistics <- t(t(teststatistics)/inflation(object))
        log_pvalues <- log(2) + pnorm(-abs(teststatistics), mean = 0, sd = 1, log.p = TRUE)
    }
    return(log_pvalues)
}

方法3：用高精度计算保留极小P值

如果需要极高精度的原始P值，使用Rmpfr包进行多精度浮点数运算：

# 先安装并加载包
install.packages("Rmpfr")
library(Rmpfr)

get_bacon_pvalues_mpfr <- function(object, corrected = TRUE, precBits = 1024) {
    if (!corrected | any(is.na(bias(object))) | any(is.na(inflation(object)))) {
        ts <- as(object@teststatistics, "mpfr")
        log_p <- pnorm(-abs(ts), log.p = TRUE, precBits = precBits)
        pvalues <- 2 * exp(log_p)
    } else {
        ts <- as(object@teststatistics, "mpfr")
        ts <- t(t(ts) - bias(object))
        ts <- t(t(ts)/inflation(object))
        log_p <- pnorm(-abs(ts), mean = 0, sd = 1, log.p = TRUE, precBits = precBits)
        pvalues <- 2 * exp(log_p)
    }
    return(pvalues)
}

注意事项 #

• 方法1和方法2是常规分析的优先选择，计算效率高且满足绝大多数需求；
• 方法3适合需要精确极小P值的场景，但计算速度较慢，结果为mpfr类型，需注意后续处理的兼容性。

内容的提问来源于stack exchange，提问作者chris333999