计算具有非常大的组的 Cliffs delta 会导致整数溢出

Question

在 R 中，我需要计算 Cliffs delta。公式如下：

其中 x_i 是 A 组中的观测值，x_j 是 B 组中的观测值，并且 [x_i > x_j] 为 1，如果 x_i > x_j 为是的。

这是 Cliff (1993) 的非数学解释：

...一组中的每个 n x 与以下 m 个中的每个进行比较 另一个，并计算该成员的次数 第一组高了多少倍，低了多少倍。

然后除以成对比较的数量，即每组中观察数量的乘积。

问题来了：在我的数据中，A 组有 66208 个观测值，B 组有 228691 个观测值。 66208*228691 导致 R 中整数溢出。

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

x <- length(a) * length(b)

Warning message:
In length(a) * length(b) : NAs produced by integer overflow

因此，我不确定在计算悬崖增量时是否可以信任

effsize

的结果，因为它会抛出类似的警告：

library(effsize)
x <- cliff.delta(a, b)

Warning messages:
1: In n1 * n2 : NAs produced by integer overflow
2: In n1 * n2 : NAs produced by integer overflow

> x$estimate
[1] -0.0005022877

我可能可以找到一种不使用

effsize

来实现计算的方法，但我没有看到任何解决需要除以 15141173728 (228691*66208) 的问题的方法，我不知道该怎么做，因为它会导致整数溢出。另外，分子也可能导致整数溢出。

有什么方法可以让我处理大量数据，以便我可以计算数据中的悬崖增量？

编辑：

我遵循了 Pbulls 的建议，找到了

effsize:cliff.delta

代码中计算 d 的部分。

effsize:cliff.delta

返回 10 个列表，如果我 d 唯一需要的东西，则不会有任何警告：

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

#Remove comments to replicate subsample sanity check
#vector1 <- vector1[1:1000]
#vector2 <- vector2[1:1000]

.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

Answer 1

在处理不可能的大乘法时，对数是你最好的朋友。

这个实现不是高效——你可能会很好地从

effsize

本身复制更好的成对比较方法——但它应该在数字上很好直到

log(n1) + log(n2) ~ 702

或直到添加溢出：

log.cliff.delta <- function(x, y) {
   ## bigger <- sum(outer(x, y, `>`))   ## Needs 113 Gb memory
   bigger <- sum(vapply(y, \(yo) sum(x > yo), numeric(1)))
   smaller <- sum(vapply(y, \(yo) sum(x < yo), numeric(1)))
   sign <- c(1, -1)[1+(bigger < smaller)]
   log_d <- log(abs(bigger - smaller)) - log(length(x)) - log(length(y))
   sign*exp(log_d)
}

## Subsample sanity check
effsize::cliff.delta(a[1:1000], b[1:1000])
#> delta estimate: -0.016404 (negligible)

log.cliff.delta(a[1:1000], b[1:1000])
#> -0.016404

## The whole nine yards
effsize::cliff.delta(a, b)
#> delta estimate: -0.0005022877 (negligible)
#> Warning messages:
#> 1: In n1 * n2 : NAs produced by integer overflow
#> 2: In n1 * n2 : NAs produced by integer overflow

log.cliff.delta(a, b)   ## is *not* fast
#> -0.0005022877

从中您可以看到，来自

effsize::cliff.delta

的警告并不像您想象的那么令人担忧：它实际上在计算 delta 估计本身时实现了一些巧妙的矩阵乘法和平均，

n1 * n2

分母不参与其中（而且它比我天真的尝试要快得多）。 但是，我不能代表置信极限，这些可能会受到未知程度的影响，但您没有为它们提供公式，并且在这些样本大小的情况下，任何置信区间都可能会如此狭窄，以至于毫无意义.