用r语言绘制roc曲线（基于R语言的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之间，代表分类模型发现真正病人与非病人的总能力。约登指数越大，表示分类模型性能越好。约登指数的示意图如下：

youden index

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

本文采用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曲线上显示最佳阈值点

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

ROC曲线

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

# 获取最佳阈值 roc_result <- coords(rocobj, "best") # 计算在最佳阈值下混淆矩阵各项的值 TP <- dim(aSAH[as.numeric(aSAH$outcome)==2 & aSAH$s100b > result$threshold, ])[1] FP <- dim(aSAH[as.numeric(aSAH$outcome)==1 & aSAH$s100b > result$threshold, ])[1] TN <- dim(aSAH[as.numeric(aSAH$outcome)==1 & aSAH$s100b <= result$threshold, ])[1] FN <- dim(aSAH[as.numeric(aSAH$outcome)==2 & aSAH$s100b <= result$threshold, ])[1]

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

# 根据混淆矩阵计算特异度、敏感度以及准确度指标 TPR <- TP / (TP FN) TNR <- TN / (TN FP) ACC <- (TP TN) / (TP TN FP FN)

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