在R中创建多个ROC曲线？

Question

我有150列得分和1列标签（1/0）。

我的目标是创建150个AUC得分。

以下是一个手动示例：

auc(roc(df$label, df$col1)),
auc(roc(df$label, df$col2)),

...

我可以使用Map / sapply / lapply，但是否有其他方法或函数？

- steves

---

apply家族函数有什么问题吗？在我看来，purrr包做类似的事情但更加直观。但是，它与apply家族函数相似。你的最终目标是什么，是一个包含所有分数的数据框吗？ - Amar

在 dplyr 中，你可以这样做：1. 使用 gather() 收集所有的 col1, col2, ... 分数 2. 对分数进行 group_by() 操作 3. 使用 summarise() 获得 AUC。 - Marius

@Marius 我知道这种方法，但是它是否比运行sapply更快或者一样快呢？ - steves

我没有测试过，所以很难说。我猜应该差不多，因为大部分运行时间都在 roc()/auc() 中，而不受你用来循环/迭代的方法影响，但也有可能会有一些差异。 - Marius

我看到只有一颗CPU在工作，我想让所有的16颗CPU都参与计算，我该如何告诉它们使用所有CPU？@ Marius - steves

Answer 1

This is a bit of an XY question. What you actually want to achieve is speed up your calculation. gfgm's answer answers it with parallelization, but that's only one way to go.

If, as I assume, you are using library(pROC)'s roc/auc functions, you can gain even more speed by selecting the appropriate algorithm for your dataset. pROC comes with essentially two algorithms that scale very differently depending on the characteristics of your data set. You can benchmark which one is the fastest by passing algorithm=0 to roc:

# generate some toy data
label <- rbinom(600000, 1, 0.5)
score <- rpois(600000, 10)

library(pROC)
roc(label, score, algorithm=0)
Starting benchmark of algorithms 2 and 3, 10 iterations...
  expr        min         lq       mean     median        uq      max neval
2    2 4805.58762 5827.75410 5910.40251 6036.52975 6085.8416 6620.733    10
3    3   98.46237   99.05378   99.52434   99.12077  100.0773  101.363    10
Selecting algorithm 3.

在这里，我们选择算法3，在阈值数量较少时效果非常好。但是，如果600000个数据点需要5分钟才能计算出来，我强烈怀疑您的数据非常连续（没有相同值的测量），并且您的阈值数量与数据点数量相当（600000）。在这种情况下，您可以直接跳转到算法2，随着ROC曲线中阈值数量的增加，它的可扩展性更好。

然后您可以运行：

auc(roc(df$label, df$col1, algorithm=2)),
auc(roc(df$label, df$col2, algorithm=2)),

在我的机器上，每次调用roc现在大约需要5秒钟，与阈值的数量相对独立。这样你应该在不到15分钟内完成。除非你有50个或更多核心，否则这比仅仅并行化要快。但当然你也可以两者都做...

- Calimo

Answer 2

如果您想并行计算，可以按如下方式进行：

# generate some toy data
label <- rbinom(1000, 1, .5)
scores <- matrix(runif(1000*150), ncol = 150)
df <- data.frame(label, scores)

library(pROC)
library(parallel)

auc(roc(df$label, df$X1))
#> Area under the curve: 0.5103

auc_res <- mclapply(df[,2:ncol(df)], function(row){auc(roc(df$label, row))})
head(auc_res)
#> $X1
#> Area under the curve: 0.5103
#> 
#> $X2
#> Area under the curve: 0.5235
#> 
#> $X3
#> Area under the curve: 0.5181
#> 
#> $X4
#> Area under the curve: 0.5119
#> 
#> $X5
#> Area under the curve: 0.5083
#> 
#> $X6
#> Area under the curve: 0.5159

由于大部分的计算时间似乎都是针对auc(roc(...))的调用，如果你有一台多核机器，那么这应该可以加速处理。

- gfgm

---

我完全按照你的方式做了，但是使用了parSapply。为什么你使用了mclapply？@gfgm - steves

说实话，因为我更好地记得语法，并且它可以节省一行代码在集群上加载环境。但实际上应该给出parSapply()/parLapply()的示例，因为mcapply不兼容Windows（这是我没有考虑到的）@steves - gfgm

这里讨论了Mclapply和Parlapply在R中的区别，请参考以下链接：https://dev59.com/CmQm5IYBdhLWcg3w-i9L?utm_medium=organic&utm_source=google_rich_qa&utm_campaign=google_rich_qa - gfgm

非常感谢，问题在于如果我使用parSapply，它会占用更少的内存，而mclapply则会占用更多的内存，有时会崩溃。@gfgm - steves

Answer 3

在cutpointr包中有一个函数可以实现这个功能。它还会计算截断点和其他指标，但是您可以将其丢弃。默认情况下，它将尝试除响应列以外的所有列作为预测变量。此外，您可以通过省略direction或手动设置来选择ROC曲线的方向（较大的值是否意味着正类或相反）。

dat <- iris[1:100, ]
library(tidyverse)
library(cutpointr)
mc <- multi_cutpointr(data = dat, class = "Species", pos_class = "versicolor", 
                silent = FALSE)
mc %>% select(variable, direction, AUC)

# A tibble: 4 x 3
  variable     direction   AUC
  <chr>        <chr>     <dbl>
1 Sepal.Length >=        0.933
2 Sepal.Width  <=        0.925
3 Petal.Length >=        1.00 
4 Petal.Width  >=        1.00  

顺便提一下，运行时间在这里不应该是一个问题，因为使用cutpointr或ROCR计算ROC曲线（甚至包括切点）对于一个变量和一百万个观测值只需要不到一秒钟的时间，所以您的任务大约需要一到两分钟才能完成。

如果内存是限制因素，那么并行化可能会使问题更加严重。如果上述解决方案占用太多内存，因为它在删除这些列之前返回所有变量的ROC曲线，您可以尝试在调用map时立即选择感兴趣的列。

# 600.000 observations for 150 variables and a binary outcome

predictors <- matrix(data = rnorm(150 * 6e5), ncol = 150)
dat <- as.data.frame(cbind(y = sample(0:1, size = 6e5, replace = T), predictors))

library(cutpointr)
library(tidyverse)

vars <- colnames(dat)[colnames(dat) != "y"]
result <- map_df(vars, function(coln) {
    cutpointr_(dat, x = coln, class = "y", silent = TRUE, pos_class = 1) %>%
        select(direction, AUC) %>%
        mutate(variable = coln)
})

result

# A tibble: 150 x 3
   direction   AUC variable
   <chr>     <dbl> <chr>   
 1 >=        0.500 V2      
 2 <=        0.501 V3      
 3 >=        0.501 V4      
 4 >=        0.501 V5      
 5 <=        0.501 V6      
 6 <=        0.500 V7      
 7 <=        0.500 V8      
 8 >=        0.502 V9      
 9 >=        0.501 V10     
10 <=        0.500 V11     
# ... with 140 more rows 

- thie1e