R语言机器学习 | 9 ROC曲线

Posted 2021-04-18 PsychRun

ROC曲线(receiver operating characteristic curve)全称为接受者操作特征曲线，是衡量模型表现的有力工具。在实验心理学里的信号检测论部分通常也会提到ROC曲线，实质和这边的是一样的东西，但是习惯用语上会有一些差异，信号检测论和ROC曲线的解释在这里不多介绍。

ROC曲线：曲线下覆盖的面积越大越好，曲线越接近左上角越好

信号检测论的基本参数

在介绍ROC曲线时，会先介绍一些容易混淆的概念：

# 真阳性（True Positive , TP）被模型预测为正的正样本；

# 假阴性（False Negative , FN）被模型预测为负的正样本；

# 假阳性（False Positive , FP）被模型预测为正的负样本；

# 真阴性（True Negative , TN）被模型预测为负的负样本。

# 真阳性率（TPR）也叫 灵敏度（sensitivity）也叫 召回率（recall） 也叫击中率；TPR = TP /（TP + FN） 

# 真阴性率（TNR) 也叫 特异度（specificity）；TNR = TN /（TN + FP）

# 假阳性率（1-TNR） 也叫 虚报率 ；1-TNR = FP/(TN+FP)

通常：ROC曲线 y轴为TPR（击中率）；x轴为1-TNR（虚报率）。有时也有ROC曲线的X轴会是TNR（特异度），要注意区别。

## 其他指标

# 准确率(Precision) =  TP/（TP+FP）

# F1 = 2 * P * R / (P + R)  同时描绘了准确率和灵敏度

# AUC：ROC曲线下面积

在机器学习建模后，ROC曲线和上面一系列指标通常都会报告，当然除了以上的一些，还有一些别的模型指标，详见：  

https://blog.csdn.net/sinat_26917383/article/details/51114244

目前有的软件在做机器学习时可以直接输出ROC和一系列模型参数，如JASP可以做不少算法，非常方便。那么在R中，可以利用一些已有的package如pROC等，也可以较为简便的实现ROC的绘制，接下来用各个之前提过的算法为例，尝试画一下这些算法的ROC曲线。画ROC的关键是得到“预测概率”和“标签”，不同算法的概率值获得的方法不同，详细见下面代码部分。

（1）逻辑回归

library(pROC) #输入一列真实的分类，一列模型的预测概率
## 实战1：逻辑回归ROClibrary(rio);library(tidyverse)library(car)library(mlogit)Cured = c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1)Intervention = c("NO", "NO", "YES", "NO", "NO", "NO", "YES","NO", "NO", "NO", "NO", "NO", "YES","YES","YES","YES", "NO","YES", "YES", "YES", "YES", "NO" )Duration = c(7, 7,6, 5, 7, 4, 7, 6, 3, 7, 6, 3, 8, 8, 7, 9, 9, 7, 8, 5, 10, 9)A = data.frame(Cured,Intervention,Duration)M1 = glm(data = A, Cured~Intervention,family = binomial()) #二项分布（二元）M2 = glm(data = A, Cured~Intervention+Duration,family = binomial()) # JASP中的factor和cov类似## ROC曲线的绘制Yhat1 = predict(M1,A,type = 'response') #预测概率值 注意逻辑回归里面type选responseYhat2 = predict(M2,A,type = 'response') #预测概率值A$cured;Yhat1;Yhat2modelroc = roc(A$Cured,Yhat1)modelroc2 = roc(A$Cured,Yhat2)
plot(modelroc,print.auc=T,auc.polygon=T,max.auc.polygon=T,print.thres=T,auc.polygon.col='skyblue' ,main="ROC", percent=TRUE)lines(modelroc2,col=2)auc(A$Cured,Yhat1);auc(A$Cured,Yhat2) #计算AUC值
#ggplot画风的ROC (g=ggroc(list(modelroc,modelroc2),legacy.axes = T,lwd=1.5))g+theme_minimal() + ggtitle("My ROC curve") +  geom_segment(aes(x = 0, xend = 1, y = 0, yend = 1), color="grey", linetype="dashed")

注意这里横坐标是特异度

与上图一样的ROC，但是ggplot画风，且注意这里横坐标是虚报率

（2）决策树

## 实战2：决策树ROCrm(list=ls())library(rpart.plot);set.seed(100)index <- sample(nrow(iris),0.7*nrow(iris)) #在x里面抽size个样本；这里直接是抽index，以便接下来使用。train <- iris[index,]test <- iris[-index,]tree = rpart(Species~Sepal.Width+Sepal.Width,data = train,method = 'class')Yhat= predict(tree,test,type = 'prob') #预测概率值 modelroc1 = roc(test$Species,Yhat[,1]) #注意一次只能一列modelroc2 = roc(test$Species,Yhat[,2])modelroc3 = roc(test$Species,Yhat[,3])(g=ggroc(list(modelroc1,modelroc2,modelroc3),legacy.axes = T,lwd=1.5,alpha=0.5,size=1.2))g=g+theme_minimal() + ggtitle("My ROC curve") +  geom_segment(aes(x = 0, xend = 1, y = 0, yend = 1), color="grey", linetype="dashed")# OR facetingg + facet_grid(.~name) + theme(legend.position="none")

facet+ggplot画风

（3）逻辑回归与决策树ROC的比较

## 实战3：决策树与逻辑回归的对比（二分类）df = import('tree.csv') # 为了之后与逻辑回归比较，这边用了一个关于肿瘤的二分类的数据，如果没有可以用别的二分类数据集代替df[df$class == 'benign',10] = 0df[df$class == 'malignant',10] = 1df$class=as.numeric(df$class)# 逻辑回归glm1 = glm(formula = class~.,family = binomial(),data = df)glm_AIC = step(glm1,k=2) #利用AIC准则逐步选出最佳的方程Yhat_glm = predict(glm_AIC,df,type = 'response')roc_glm = roc(df$class,Yhat_glm)#决策树tree2 = rpart(class~.,data = df,method = 'class')rpart.plot(tree2,cex = 0.45,type = 4)Yhat= predict(tree2,df,type = 'prob') #预测值roc_tree2=roc(df$class,Yhat[,2])
## ROC对比1plot(roc_tree2,print.auc=F,print.thres=F,col=2, percent=TRUE)lines(roc_glm,print.auc=F,print.thres=F)## ROC对比2(g=ggroc(list(roc_glm,roc_tree2),legacy.axes = T,lwd=1.5,alpha=0.5,size=1.2))(g=g+theme_minimal() + ggtitle("My ROC curve") +  geom_segment(aes(x = 0, xend = 1, y = 0, yend = 1), color="grey", linetype="dashed"))

（4）随机森林

## 实战4：随机森林rm(list=ls())library(randomForest)index <- sample(nrow(iris),0.7*nrow(iris)) #在x里面抽size个样本；这里直接是抽index，以便接下来使用。train <- iris[index,]test <- iris[-index,]forest = randomForest(Species~Sepal.Width+Sepal.Width,data = train) #默认500棵树Yhat=predict(forest,test,type = 'prob')#预测值roc_for1=roc(test$Species,Yhat[,1])roc_for2=roc(test$Species,Yhat[,2])roc_for3=roc(test$Species,Yhat[,3])
(g=ggroc(list(roc_for1,roc_for2,roc_for3),legacy.axes = T,lwd=1.5,alpha=0.5,size=1.2))(g=g+theme_minimal() + ggtitle("My ROC curve") +  geom_segment(aes(x = 0, xend = 1, y = 0, yend = 1), color="grey", linetype="dashed"))# OR facetingg + facet_grid(.~name) + theme(legend.position="none")

（5）SVM

## 实战5：SVMlibrary(e1071)SVM1 = svm(data=train,Species~.,kernal='sigmoid',probability = TRUE) #prob=T 来计算Yhat_svm = predict(SVM1,test,probability = TRUE) #prob=T 使导出预测值Yhat_svm = attributes(Yhat_svm)$probabilities #把预测值单独放出来
roc_for1=roc(test$Species,Yhat_svm[,1])roc_for2=roc(test$Species,Yhat_svm[,2])roc_for3=roc(test$Species,Yhat_svm[,3])(g=ggroc(list(roc_for1,roc_for2,roc_for3),legacy.axes = T,lwd=1.5,alpha=0.5,size=1.2))(g=g+theme_minimal() + ggtitle("My ROC curve") +  geom_segment(aes(x = 0, xend = 1, y = 0, yend = 1), color="grey", linetype="dashed"))

（6）NaiveBayes

## 实战6：朴素贝叶斯library(e1071) rm(list=ls())set.seed(100)index <- sample(nrow(iris),0.7*nrow(iris)) train <- iris[index,]test <- iris[-index,]nb1 <- naiveBayes(Species ~Sepal.Width, data =train )#计算NBYhat <- predict(nb1, test,type = 'raw')roc_for1=roc(test$Species,Yhat[,1])roc_for2=roc(test$Species,Yhat[,2])roc_for3=roc(test$Species,Yhat[,3])
(g=ggroc(list(roc_for1,roc_for2,roc_for3),legacy.axes = T,lwd=1.5,alpha=0.5,size=1.2))(g=g+theme_minimal() + ggtitle("My ROC curve") +  geom_segment(aes(x = 0, xend = 1, y = 0, yend = 1), color="grey", linetype="dashed"))# OR facetingg + facet_grid(.~name) + theme(legend.position="none")

以上是初步的对各算法的ROC进行绘制，如果要更好看的或有任何特殊要求的，可以自行DIV处理！