基于R语言的ROC曲线绘制及最佳阈值点(Cutoff)选择

https://zhuanlan.zhihu.com/p/159201640

ROC曲线

在介绍ROC曲线之前，我们首先需要介绍混淆矩阵(Confusion Matrix)。在统计分类模型的评估过程中分别统计分类模型归错类，归对类的观测值个数，然后把结果放在一个表里展示出来的表格就是混淆矩阵。混淆矩阵的示意图如下：

 



在混淆矩阵中：

TP代表的是真实值是positive，模型分类为positive的样本数量。

FP代表的是真实值是negative，模型分类为positive的样本数量。

TN代表的是真实值是negative，模型分类为negative的样本数量。

FN代表的是真实值是positive，模型分类为negative的样本数量。

根据混淆矩阵我们可以计算出模型分类的真阳性率(TPR, True Positive Rate)，真阴性率(TNR, True Negative Rate)以及假阴性率(FNR, False Negative Rate)：

TPR = TP / (TP + FN)

TNR = TN / (TN + FP)

FNR = FN / (TN + FP) = 1- TNR

其中TPR又称敏感度(sensitivity)，TNR又称特异度(specificity)。

而我们知道，通过给定分类模型不同的阈值，其分类结果是变化的。即混淆矩阵中的TP，FP，TN，FN的取值根据不同的阈值会发生变化，那么相应的TPR和FNR也将随之发生变化。随着阈值的变化，以FNR(1-specificity)和TPR(sensitivity)分别为横轴和纵轴绘制的曲线即为ROC曲线。

 

最佳阈值点选择

在实际使用中，ROC曲线上的最佳阈值点所对应的混淆矩阵将是我们计算敏感度(sensitivity)、特异度(specificity)以及准确度等指标的依据。那么ROC曲线上的哪一个点对应的阈值是最佳阈值点呢？通常情况下我们会通过约登指数(Youden index)进行选择。约登指数也称正确指数，是指敏感度和特异度之和减去1：

Youden index = Sensitivity + Specificity − 1

约登指数指数范围取值介于0-1之间，代表分类模型发现真正病人与非病人的总能力。约登指数越大，表示分类模型性能越好。约登指数的示意图如下：

 



 

 

图中的C表示最佳阈值点，红色线段J的长度表示约登指数的取值。

 

R语言实践

本文采用R语言中pROC包用于ROC曲线对象的生成与绘制，如果你在之间没有安装过pROC包，需要在RStudio的Console中运行下面的命令进行安装：

 

install.packages("pROC")

为了简化步骤，本文没有利用线性回归方法构造分类模型然后再完成ROC曲线的生成与绘制，而是采用pROC包中内置的“aSAH”数据集直接进行ROC曲线的生成与绘制。

首先我们需要导入pROC包和aSAH数据集，代码如下：

 

library(pROC)

# 导入数据

data(aSAH)

然后我们以aSAH中的s100b作为预测值(s100b这里通常是使用分类模型的预测值)，以outcome作为真实值构造roc曲线对象并进行绘制，代码如下：

 

# 生成roc曲线对象

rocobj <- roc(aSAH$outcome, aSAH$s100b)

# 绘制roc曲线

plot(rocobj,

     legacy.axes = TRUE,

     main="ROC曲线最佳阈值点",

     thresholds="best", # 基于youden指数选择roc曲线最佳阈值点

     print.thres="best") # 在roc曲线上显示最佳阈值点

上述代码的绘制结果如下图：

 



通过上图中的信息我们可以了解到，最佳阈值是0.205，相应的特异度和敏感度分别为0.806和0.634。至此，ROC曲线的生成与绘制基本已经结束，但是为了获取更多指标我们还需要进行进一步的处理。

我们需要根据最佳阈值计算相应混淆矩阵各项的取值，代码如下：

 

# 获取最佳阈值

roc_result <- coords(rocobj, "best")

# 计算在最佳阈值下混淆矩阵各项的值

TP <- dim(aSAH[as.numeric(aSAH$outcome)==2 & aSAH$s100b > roc_result$threshold, ])[1]

FP <- dim(aSAH[as.numeric(aSAH$outcome)==1 & aSAH$s100b > roc_result$threshold, ])[1]

TN <- dim(aSAH[as.numeric(aSAH$outcome)==1 & aSAH$s100b <= roc_result$threshold, ])[1]

FN <- dim(aSAH[as.numeric(aSAH$outcome)==2 & aSAH$s100b <= roc_result$threshold, ])[1]

最后，在获取了混淆矩阵各项的取值之后，我们便可以计算我们需要的指标，如特异度、敏感度以及准确度。

 

# 根据混淆矩阵计算特异度、敏感度以及准确度指标

TPR <- TP / (TP + FN)

TNR <- TN / (TN + FP)

ACC <- (TP + TN) / (TP + TN + FP + FN)

通过混淆矩阵计算出的特异度、敏感度以及准确度分别为0.6341463、0.8055556以及0.7433628。当然，更复杂的指标也是能够计算的，本文在这里就不继续展开了。