【问题标题】:Calculate AUC in R?计算R中的AUC?
【发布时间】:2011-06-21 15:29:37
【问题描述】:

给定一个分数向量和一个实际类标签向量,如何计算 R 语言或简单英语中二元分类器的单数 AUC 度量?

"AUC: a Better Measure..."的第9页似乎需要知道类标签,这里是an example in MATLAB我不明白的地方

R(Actual == 1))

因为 R(不要与 R 语言混淆)被定义为向量但用作函数?

【问题讨论】:

标签: r machine-learning data-mining auc


【解决方案1】:

使用包pROC,您可以使用函数auc(),就像帮助页面中的这个例子:

> data(aSAH)
> 
> # Syntax (response, predictor):
> auc(aSAH$outcome, aSAH$s100b)
Area under the curve: 0.7314

【讨论】:

    【解决方案2】:

    The ROCR package 将计算 AUC 以及其他统计信息:

    auc.tmp <- performance(pred,"auc"); auc <- as.numeric(auc.tmp@y.values)
    

    【讨论】:

    • 我使用 ROCR 来绘制性能,但我看不到它如何计算“单数 AUC 指标”(来自原始问题)。
    • auc.tmp &lt;- performance(pred,"auc"); auc &lt;- as.numeric(auc.tmp@y.values)
    【解决方案3】:

    正如其他人所提到的,您可以使用 ROCR 包计算 AUC。使用 ROCR 包,您还可以绘制 ROC 曲线、提升曲线和其他模型选择措施。

    您可以直接计算 AUC,而无需使用任何软件包,因为 AUC 等于真阳性得分高于真阴性的概率。

    例如,如果pos.scores 是一个包含正例得分的向量,neg.scores 是一个包含负例的向量,则 AUC 近似为:

    > mean(sample(pos.scores,1000,replace=T) > sample(neg.scores,1000,replace=T))
    [1] 0.7261
    

    将给出 AUC 的近似值。您还可以通过自举估计 AUC 的方差:

    > aucs = replicate(1000,mean(sample(pos.scores,1000,replace=T) > sample(neg.scores,1000,replace=T)))
    

    【讨论】:

    • 对于我的测试数据集,您的复制值与@jonw 的(为 0.8504,您的为 0.850591)非常相似,但我不需要安装 pROC。谢谢
    • @Andrew @eric 这是一个糟糕的答案。您确实估计 AUC 的方差 - 您只估计重采样过程的方差。为了说服自己,请尝试更改sample 中的样本量...除以 10,您的方差乘以 10。乘以 10,您的方差除以 10。这当然不是计算AUC 的方差。
    • 此外,答案应注意估计与重复次数一样好。去无穷大,你会得到实际的 AUC。
    • 同意@Calimo,这不是引导程序。要引导,您必须用替换 M 次重新采样 N 个数据点,其中 N 是原始数据集的总大小,M 可以是任何值(通常是几百或更多)。 N 不是任意的。如果 N 未设置为完整的数据集大小,您将获得有偏差的统计信息。
    • 我对显示的基本 R 方法有点不清楚。可以纯粹从混淆矩阵中计算出来吗?在给定混淆矩阵的上下文中,pos.scoresneg.scores 会是什么?
    【解决方案4】:

    无需任何额外的包:

    true_Y = c(1,1,1,1,2,1,2,1,2,2)
    probs = c(1,0.999,0.999,0.973,0.568,0.421,0.382,0.377,0.146,0.11)
    
    getROC_AUC = function(probs, true_Y){
        probsSort = sort(probs, decreasing = TRUE, index.return = TRUE)
        val = unlist(probsSort$x)
        idx = unlist(probsSort$ix)  
    
        roc_y = true_Y[idx];
        stack_x = cumsum(roc_y == 2)/sum(roc_y == 2)
        stack_y = cumsum(roc_y == 1)/sum(roc_y == 1)    
    
        auc = sum((stack_x[2:length(roc_y)]-stack_x[1:length(roc_y)-1])*stack_y[2:length(roc_y)])
        return(list(stack_x=stack_x, stack_y=stack_y, auc=auc))
    }
    
    aList = getROC_AUC(probs, true_Y) 
    
    stack_x = unlist(aList$stack_x)
    stack_y = unlist(aList$stack_y)
    auc = unlist(aList$auc)
    
    plot(stack_x, stack_y, type = "l", col = "blue", xlab = "False Positive Rate", ylab = "True Positive Rate", main = "ROC")
    axis(1, seq(0.0,1.0,0.1))
    axis(2, seq(0.0,1.0,0.1))
    abline(h=seq(0.0,1.0,0.1), v=seq(0.0,1.0,0.1), col="gray", lty=3)
    legend(0.7, 0.3, sprintf("%3.3f",auc), lty=c(1,1), lwd=c(2.5,2.5), col="blue", title = "AUC")
    

    【讨论】:

    • 如果您复制粘贴此代码并收到Error in plot.window(...) : need finite 'xlim' values,可能是因为您的标签是0-1,而@AGS 使用的标签是1-2。
    • 如果两个观测值具有相同的概率并且观测值的顺序不是随机的,则它不会给出真实的 AUC。否则代码又好又快。
    • 不知道为什么此解决方案不适用于我的数据,我的概率未标准化为在 [0,1] 内
    【解决方案5】:

    我发现这里的一些解决方案很慢和/或令人困惑(其中一些不能正确处理关系),所以我在我的 R 包mltools 中编写了我自己的基于data.table 的函数auc_roc()

    library(data.table)
    library(mltools)
    
    preds <- c(.1, .3, .3, .9)
    actuals <- c(0, 0, 1, 1)
    
    auc_roc(preds, actuals)  # 0.875
    
    auc_roc(preds, actuals, returnDT=TRUE)
       Pred CountFalse CountTrue CumulativeFPR CumulativeTPR AdditionalArea CumulativeArea
    1:  0.9          0         1           0.0           0.5          0.000          0.000
    2:  0.3          1         1           0.5           1.0          0.375          0.375
    3:  0.1          1         0           1.0           1.0          0.500          0.875
    

    【讨论】:

    • 这个解决方案比 pROC 包中的 auc() 方法快得多!如果必须计算多类或多输出回归问题的 auc 分数,则 pROC 包中的 auc() 方法非常慢。
    【解决方案6】:

    您可以在Miron Kursa 的这篇博文中了解有关 AUROC 的更多信息:

    https://mbq.me/blog/augh-roc/

    他为 AUROC 提供了一个快速函数:

    # By Miron Kursa https://mbq.me
    auroc <- function(score, bool) {
      n1 <- sum(!bool)
      n2 <- sum(bool)
      U  <- sum(rank(score)[!bool]) - n1 * (n1 + 1) / 2
      return(1 - U / n1 / n2)
    }
    

    让我们测试一下:

    set.seed(42)
    score <- rnorm(1e3)
    bool  <- sample(c(TRUE, FALSE), 1e3, replace = TRUE)
    
    pROC::auc(bool, score)
    mltools::auc_roc(score, bool)
    ROCR::performance(ROCR::prediction(score, bool), "auc")@y.values[[1]]
    auroc(score, bool)
    
    0.51371668847094
    0.51371668847094
    0.51371668847094
    0.51371668847094
    

    auroc()pROC::auc()computeAUC() 快 100 倍。

    auroc()mltools::auc_roc()ROCR::performance() 快 10 倍。

    print(microbenchmark(
      pROC::auc(bool, score),
      computeAUC(score[bool], score[!bool]),
      mltools::auc_roc(score, bool),
      ROCR::performance(ROCR::prediction(score, bool), "auc")@y.values,
      auroc(score, bool)
    ))
    
    Unit: microseconds
                                                                 expr       min
                                               pROC::auc(bool, score) 21000.146
                                computeAUC(score[bool], score[!bool]) 11878.605
                                        mltools::auc_roc(score, bool)  5750.651
     ROCR::performance(ROCR::prediction(score, bool), "auc")@y.values  2899.573
                                                   auroc(score, bool)   236.531
             lq       mean     median        uq        max neval  cld
     22005.3350 23738.3447 22206.5730 22710.853  32628.347   100    d
     12323.0305 16173.0645 12378.5540 12624.981 233701.511   100   c 
      6186.0245  6495.5158  6325.3955  6573.993  14698.244   100  b  
      3019.6310  3300.1961  3068.0240  3237.534  11995.667   100 ab  
       245.4755   253.1109   251.8505   257.578    300.506   100 a   
    

    【讨论】:

    • 对于更大的样本量,bigstatsr::AUC() 甚至更快(在 C++ 中实现)。免责声明:我是作者。
    【解决方案7】:

    将来自ISL 9.6.3 ROC Curves 的代码与@J 组合在一起。 Won.'s answer to this question和其他几个地方,下面绘制ROC曲线并在图的右下角打印AUC。

    probs 下面是二进制分类预测概率的数字向量,test$label 包含测试数据的真实标签。

    require(ROCR)
    require(pROC)
    
    rocplot <- function(pred, truth, ...) {
      predob = prediction(pred, truth)
      perf = performance(predob, "tpr", "fpr")
      plot(perf, ...)
      area <- auc(truth, pred)
      area <- format(round(area, 4), nsmall = 4)
      text(x=0.8, y=0.1, labels = paste("AUC =", area))
    
      # the reference x=y line
      segments(x0=0, y0=0, x1=1, y1=1, col="gray", lty=2)
    }
    
    rocplot(probs, test$label, col="blue")
    

    这给出了这样的情节:

    【讨论】:

      【解决方案8】:

      我通常使用 DiagnosisMed 包中的函数ROC。我喜欢它产生的图表。 AUC 连同它的置信区间一起返回,并且在图表上也提到了。

      ROC(classLabels,scores,Full=TRUE)
      

      【讨论】:

      【解决方案9】:

      按照 erik 的回答,您还应该能够通过比较 pos.scores 和 neg.scores 中所有可能的值对来直接计算 ROC:

      score.pairs <- merge(pos.scores, neg.scores)
      names(score.pairs) <- c("pos.score", "neg.score")
      sum(score.pairs$pos.score > score.pairs$neg.score) / nrow(score.pairs)
      

      肯定比示例方法或 pROC::auc 效率低,但比前者更稳定,并且比后者需要更少的安装。

      相关:当我尝试这个时,它给出了与 pROC 值相似的结果,但不完全相同(相差 0.02 左右);结果更接近 N 非常高的示例方法。如果有人知道为什么我会感兴趣。

      【讨论】:

      • 不准确的一个来源是处理关系。从技术上讲,您应该采用阳性病例分数严格大于阴性分数的概率 + 1/2 * 它们相等的概率。如果所有分数都是唯一的,这将不是问题。
      【解决方案10】:

      当前投票最多的答案不正确,因为它忽略了关系。当正负分相等时,AUC 应为 0.5。下面是更正的例子。

      computeAUC <- function(pos.scores, neg.scores, n_sample=100000) {
        # Args:
        #   pos.scores: scores of positive observations
        #   neg.scores: scores of negative observations
        #   n_samples : number of samples to approximate AUC
      
        pos.sample <- sample(pos.scores, n_sample, replace=T)
        neg.sample <- sample(neg.scores, n_sample, replace=T)
        mean(1.0*(pos.sample > neg.sample) + 0.5*(pos.sample==neg.sample))
      }
      

      【讨论】:

        【解决方案11】:

        使用Metrics 包计算 AUC 非常简单直接:

        library(Metrics)
        
        actual <- c(0, 0, 1, 1)
        predicted <- c(.1, .3, .3, .9)
        
        auc(actual, predicted)
        
        0.875
        

        【讨论】:

          猜你喜欢
          • 2013-05-08
          • 2011-12-10
          • 2013-12-29
          • 2018-03-28
          • 2021-08-07
          • 2021-11-24
          • 2016-08-30
          • 1970-01-01
          • 2021-04-26
          相关资源
          最近更新 更多