基于基础R函数编写2×2列联表卡方被积函数遇integrate错误求助

背景与问题描述 #

手动实现给定边际的2×2列联表卡方统计量概率密度函数（PDF），不调用R原生函数以深入理解原理。列联表结构如下：

a   b     ab
c   d     cd
ac  bd   abcd

变量定义：

• a、b、c、d：列联表单元格值
• ab = a+b：基因集总基因数
• ac = a+c：目标组基因数
• abcd：检测到的总基因数

编写了被积函数integrand()，要求对参数a向量化以适配integrate()，但调用时频繁触发roundoff error is detected in the extrapolation table错误，排查两天未定位问题。已排除小单元格卡方效能、积分限设置的影响，未使用Yates连续性校正。

原实现代码 #

integrand <- function(a, ...)
{
  args <- list(...)
  sapply(a, \(x, ...)
  {
    args <- list(...)
    observed <- c(x, args[[1]] - x, args[[2]] - x, args[[3]] - args[[2]] - args[[1]] + x)
    expected <- c(args[[1]] * args[[2]], args[[1]] * (args[[3]] - args[[2]]), (args[[3]] - args[[1]]) * args[[2]], (args[[3]] - args[[1]]) * (args[[3]] - args[[2]])) / args[[3]]
    s <- sum((observed - expected)^2 / expected)
    (s^-0.5) * exp(-s / 2) / (2^0.5 * gamma(0.5))  #  df fixed to 1 appropriately for 2x2 contingency table.
  }
  , args[[1]], args[[2]], args[[3]])
}

ab <- 11
ac <- 10
abcd <- 25

# plot the pdf between 0 and 10
curve(integrand(x, ab, ac, abcd), 0, 10)
#  integrate between 1 and 5
integrate(integrand, lower = 1, upper = 5, subdivisions = 1000, ab = ab, ac = ac, abcd = abcd)

错误原因分析 #

1. 参数传递混乱：嵌套使用list(...)和sapply，导致内层args覆盖外层参数，传递逻辑混乱，integrate()无法正确处理向量输入。
2. 无合法范围限制：当a超出合理区间（a < max(0, ab+ac-abcd)或a > min(ab, ac)）时，observed单元格会出现负数，计算卡方时产生NaN/无穷大，导致函数不连续。
3. 除以0问题：当观测值等于期望值时，s=0，s^-0.5会触发无穷大，破坏函数连续性。
4. 统计逻辑混淆：固定边际的2×2列联表中，a服从超几何分布，直接用卡方(1)分布PDF映射到a的PDF是近似处理，但原实现未处理这种映射的定义域问题。

修正后的代码 #

integrand <- function(a, ab, ac, abcd) {
  # 计算a的合法取值范围，超出范围返回0
  a_min <- max(0, ab + ac - abcd)
  a_max <- min(ab, ac)
  res <- numeric(length(a))
  
  # 仅处理合法范围内的a值
  valid_idx <- a >= a_min & a <= a_max
  a_valid <- a[valid_idx]
  
  if (length(a_valid) > 0) {
    # 向量化计算观测值矩阵
    observed <- cbind(a_valid, 
                      ab - a_valid, 
                      ac - a_valid, 
                      abcd - ab - ac + a_valid)
    # 向量化计算期望值矩阵
    expected <- cbind(ab * ac / abcd,
                      ab * (abcd - ac) / abcd,
                      (abcd - ab) * ac / abcd,
                      (abcd - ab) * (abcd - ac) / abcd)
    # 计算卡方统计量
    s <- rowSums((observed - expected)^2 / expected)
    # 替换s=0的情况，避免除以0
    s[s == 0] <- .Machine$double.eps
    # 计算卡方(1)分布的PDF值
    res[valid_idx] <- (s^-0.5) * exp(-s / 2) / (sqrt(2) * gamma(0.5))
  }
  res
}

ab <- 11
ac <- 10
abcd <- 25

# 绘制合法区间内的PDF
curve(integrand(x, ab, ac, abcd), 0, 10, ylab = "PDF")
# 在合法区间内积分
integrate(integrand, lower = 1, upper = 5, ab = ab, ac = ac, abcd = abcd)

修正说明 #

• 向量化实现：用矩阵运算直接处理向量输入，替代sapply，适配integrate()的要求，同时提升效率。
• 合法范围限制：确保所有单元格观测值非负，避免NaN/无穷大，保证函数连续性。
• 除以0处理：用R内置的极小值替代s=0，避免计算错误。
• 简化参数传递：显式传递参数，避免...导致的参数混乱。

注意：固定边际的2×2列联表中，a的真实分布是超几何分布，上述实现是基于卡方统计量近似构建的PDF，仅为教学演示用途，实际场景建议直接使用超几何分布相关函数。

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