大样本下R语言计算Cliffs delta遇整数溢出的解决方法咨询

问题背景 #

需要在R中计算Cliffs Delta，核心逻辑为：将A组n个观测值与B组m个观测值逐一比较，统计A组元素大于B组的次数、小于B组的次数，最终通过(大于次数 - 小于次数)/(n*m)得到结果。

但当A组有66208个观测值、B组有228691个观测值时，66208*228691的结果（约151亿）超过R中32位整数的最大值（约21亿），直接触发整数溢出：

set.seed(123)
a <- runif(228691)
b <- runif(66208)

x <- length(a) * length(b)
# 警告信息：In length(a) * length(b) : NAs produced by integer overflow

使用effsize包的cliff.delta()函数计算时，同样会抛出溢出警告，但仍能输出估计值：

library(effsize)
x <- cliff.delta(a, b)
# 警告信息：
# 1: In n1 * n2 : NAs produced by integer overflow
# 2: In n1 * n2 : NAs produced by integer overflow

x$estimate
# [1] -0.0005022877

自行实现计算逻辑时，无论是计算分母n*m还是分子的总比较次数，都会遇到相同的溢出问题。

原因分析 #

R默认整数类型为32位有符号整数，最大值为2^31 - 1（约2147483647）。本次数据量的乘积228691*66208=15141173728远大于该阈值，直接计算会触发溢出并返回NA。

解决方案 #

1. 复用effsize包的高效计算逻辑（无溢出风险） #

提取effsize::cliff.delta()的核心计算代码可见，它并未直接计算n1*n2，而是通过排序+二分查找分步计算均值，完全避开超大整数运算，因此不会触发溢出：

set.seed(123)
vector1 <- runif(228691)
vector2 <- runif(66208)

.bsearch.partition <- function(x, a, b = 1, e = length(a)) {
  n <- length(x)
  low <- rep(NA, n)
  L <- rep(b, n)
  H <- rep(e, n)
  
  repeat {
    M <- as.integer((L + H) / 2)
    left <- x <= a[M]
    H[left] <- M[left]
    L[!left] <- M[!left] + 1
    if (all(H <= L)) {
      break
    }
  }
  
  H <- L
  repeat {
    below <- a[H] == x
    below[is.na(below)] <- FALSE
    if (!any(below)) 
      break
    H[below] <- H[below] + 1
  }
  
  repeat {
    L.clean <- L
    L.clean[L.clean < 1] <- NA
    above <- a[L.clean] >= x
    above[is.na(above)] <- FALSE
    if (!any(above)) 
      break
    L[above] <- L[above] - 1
  }
  
  if (any(L == H)){
    H[H == L] <- L[H == L] + 1
  }
  H[H > length(a) + 1] <- length(a) + 1
  cbind(below = L, above = H)
}

treatment <- sort(vector1)
control <- sort(vector2)

n1 <- length(treatment)
n2 <- length(control)

partitions <- .bsearch.partition(treatment, control)
partitions[, 2] <- n2 - partitions[, 2] + 1L
partitions[partitions[, 1] > n2, 1] <- n2

d_i. <- mean(partitions %*% c(1L, -1L) / n2)
d <- mean(d_i.)

d
# [1] -0.0005022877

该代码核心思路：

• 对两组数据排序，用二分查找快速定位每个元素在另一组中的位置
• 对每个元素计算(小于等于当前元素的数量 - 大于等于当前元素的数量)/n2后取均值
• 最终结果等价于(总大于次数 - 总小于次数)/(n1*n2)，全程无需计算超大整数乘积

2. 手动处理大整数运算 #

若需自行实现完整逻辑，可通过以下方式避免溢出：

• 转换为64位整数：使用bit64包的integer64类型存储大整数：

  library(bit64)
  n1 <- as.integer64(length(a))
  n2 <- as.integer64(length(b))
  total <- n1 * n2  # 不会溢出
  

• 转换为浮点数计算：双精度浮点数可精确表示2^53以内的整数（约9e15），远大于本次乘积，直接转换为numeric类型即可：

  total <- as.numeric(length(a)) * as.numeric(length(b))
  

关于effsize包结果的可信性 #

effsize包的溢出警告仅出现在计算置信区间等辅助统计量的步骤中，核心的Cliffs Delta估计值是通过上述高效逻辑计算的，与手动提取代码的结果完全一致，因此该估计值可信。

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