加速分位数计算

Question

我正在使用Hmisc包计算两个连续变量的分位数，并在交叉表中比较结果。以下是我的代码。

我的问题是，如果观察数量增加，分位数的计算需要相当长的时间。

有没有可能通过使用data.table、ddply或其他任何包来加速这个过程？

谢谢。

library(Hmisc)

# Set seed
set.seed(123)

# Generate some data
a <- sample(1:25, 1e7, replace=TRUE)
b <- sample(1:25, 1e7, replace=TRUE)
c <- data.frame(a,b)

# Calculate quantiles
c$a.quantile <- cut2(a, g=5)
c$b.quantile <- cut2(b, g=5)

# Output some descriptives
summaryM(a.quantile ~ b.quantile, data=c, overall=TRUE)

# Time spent for calculation:
#       User      System verstrichen 
#      25.13        3.47       28.73 

- majom

Answer 1

如 jlhoward 和 Ricardo Saporta所述，data.table 在这种情况下似乎并没有太大的加速效果。显然，cut2 函数是瓶颈。我使用另一个函数来计算分位数（参见R 中是否有更好的创建分位数“哑变量”/因子的方法？），并且成功将计算时间减少了一半：

qcut <- function(x, n) {
  if(n<=2)
    { 
    stop("The sample must be split in at least 3 parts.")
  }
  else{
    break.values <- quantile(x, seq(0, 1, length = n + 1), type = 7)
    break.labels <- c(
      paste0(">=",break.values[1], " & <=", break.values[2]),
      sapply(break.values[3:(n)], function(x){paste0(">",break.values[which(break.values == x)-1], " & <=", x)}),
      paste0(">",break.values[(n)], " & <=", break.values[(n+1)]))
    cut(x, break.values, labels = break.labels,include.lowest = TRUE)
  }
}

c$a.quantile.2 <- qcut(c$a, 5)
c$b.quantile.2 <- qcut(c$b, 5)
summaryM(a.quantile.2 ~ b.quantile.2, data=c, overall=TRUE)

# Time spent for calculation:
#       User      System verstrichen 
#      10.22        1.47       11.70 

使用data.table可以再节省一秒的计算时间，但我更喜欢Hmisc包提供的摘要总结。

- majom

Answer 2

您可以使用`data.table`内置的变量`.N`来快速制表。

library(data.table)
library(Hmisc)

DT <- data.table(a,b)
DT[, paste0(c("a", "b"), ".quantile") := lapply(.SD, cut2, g=5), .SDcols=c("a", "b")]

DT[, .N, keyby=list(b.quantile, a.quantile)][, setNames(as.list(N), as.character(b.quantile)), by=a.quantile]

你可以将最后一行分解为两个步骤，以了解正在发生的情况。第二个"[ "只是将数据重新格式化为清晰的格式。

DT.tabulated <- DT[, .N, keyby=list(b.quantile, a.quantile)]
DT.tabulated

DT.tabulated[, setNames(as.list(N), as.character(b.quantile)), by=a.quantile]

- Ricardo Saporta

---

感谢您的回复。很高兴看到使用data.table可以实现这一点。然而，正如jlhoward所提到的，加速效果相对较小。我想知道为什么会这样？与quantile()相比，cut2()似乎比较慢。 - majom

是的，一般来说，在将函数应用于数据框时，有几种情况会看到通过将相同的函数应用于数据表而获得显着增加。在这里，cut2将成为瓶颈，无论是df还是dt。如果你只想要分位数，那么绝对可以使用它，而不是cut2。 - Ricardo Saporta

我发布了一个答案，使用另一个函数来计算分位数，这个函数是 Hadley Wickham 一段时间前发布的，确实可以显著加快处理速度。 - majom

Answer 3

数据表似乎在这里并没有改善情况：

library(Hmisc)
set.seed(123)
a <- sample(1:25, 1e7, replace=TRUE)
b <- sample(1:25, 1e7, replace=TRUE)

library(data.table)
# original approach
system.time({
  c <- data.frame(a,b)
  c$a.quantile <- cut2(a, g=5)
  c$b.quantile <- cut2(b, g=5)
  smry.1 <-summaryM(a.quantile ~ b.quantile, data=c, overall=TRUE)
})
   user  system elapsed 
  72.79    6.22   79.02 

# original data.table approach
system.time({
  DT <- data.table(a,b)
  DT[, paste0(c("a", "b"), ".quantile") := lapply(.SD, cut2, g=5), .SDcols=c("a", "b")]
  smry.2 <- DT[, .N, keyby=list(b.quantile, a.quantile)][, setNames(as.list(N), as.character(b.quantile)), by=a.quantile]
})
   user  system elapsed 
  66.86    5.11   71.98 

# different data.table approach (simpler, and uses table(...))
system.time({
  dt     <- data.table(a,b)
  smry.3 <- table(dt[,lapply(dt,cut2,g=5)])
})
   user  system elapsed 
  67.24    5.02   72.26 

- jlhoward