R语言roc曲线切换阈值设定 r语言roc曲线代码

前言：以前使用Matlab绘制ROC曲线常常是工具箱有就画，没有就不画，而且在想画的时候工具箱恰恰就没有，很纳闷。然后无意间发现了一篇用R语言绘制ROC曲线的文章，赶紧学了并分享出来，以备不时之需。

先通过一个例子来讲解一下参数的作用，使用的数据是大名鼎鼎的Iris数据集，R语言自带。

1.数据处理

第一步当然得处理一下数据。默认的Iris数据集有三类鸢尾花，我目前的理解是只有二分类才画的出ROC曲线，所以才去一定的手段处理一下数据：

• 输入

# 数据准备
iris2 <- iris
iris2$label[iris2$Species == 'setosa'] <- 1
iris2$label[iris2$Species == 'versicolor'] <- 2
iris2 <- iris2[-which(iris2$Species == 'virginica'), ]  # 剔除类型为virginica的数据
iris2$Species <- NULL   # 去除Species列
head(iris2,10)  # 显示前10个数据

• 输出

> head(iris2,10)
   Sepal.Length Sepal.Width Petal.Length Petal.Width label
1           5.1         3.5          1.4         0.2     1
2           4.9         3.0          1.4         0.2     1
3           4.7         3.2          1.3         0.2     1
4           4.6         3.1          1.5         0.2     1
5           5.0         3.6          1.4         0.2     1
6           5.4         3.9          1.7         0.4     1
7           4.6         3.4          1.4         0.3     1
8           5.0         3.4          1.5         0.2     1
9           4.4         2.9          1.4         0.2     1
10          4.9         3.1          1.5         0.1     1

2.参数设置

以label与Sepal.Length之间的ROC曲线为例：

auc1 <- roc(label~Sepal.Length, data=iris2, smooth=FALSE)

接着通过实际图形看看使用plot绘制ROC曲线的一些参数：

• 原始图形

plot(auc1)

• print.auc：在图中显示AUC的值，AUC的大小等于ROC曲线下方的面积大小；

plot(auc1, print.auc=TRUE)

• print.thres：在图中显示ROC曲线的阈值(threshold)，大概为ROC曲线中最尖的那个点；

plot(auc1, print.thres=TRUE)

• print.thres.col：设置阈值数据的颜色。

plot(auc1, print.thres=TRUE, print.thres.col="blue")

• col：设置ROC曲线的颜色。

plot(auc1,col="blue")

• identity.col：设置对角线的颜色。

plot(auc1, identity.col="blue")

• identity.lty：设置对角线的类型，取数字。

plot(auc1,identity.lty=2)

• identity.lwd：设置对角线的线宽，默认宽度为1。

plot(auc1, identity.lwd=2)

• 综合以上参数，得到如下图形

auc1 <- roc(label~Sepal.Length, data=iris2, smooth=FALSE)
plot(auc1, print.auc=TRUE, print.thres=TRUE, main="多组ROC曲线比较",
     col="blue", print.thres.col="blue", identity.col="blue",
     identity.lty=2, identity.lwd=1)

3.多变量ROC曲线比较

# ROC曲线的绘制
auc1 <- roc(label~Sepal.Length, data=iris2, smooth=FALSE)
plot(auc1, print.auc=TRUE, print.thres=TRUE, main="多组ROC曲线比较",
     col="blue", print.thres.col="blue", identity.col="blue",
     identity.lty=2, identity.lwd=1)
auc2 <- roc(label~Sepal.Width, data=iris2,smooth=FALSE)
auc3 <- roc(label~Petal.Length, data=iris2,smooth=FALSE)
auc4 <- roc(label~Petal.Width, data=iris2, smooth=FALSE)
lines(auc2, col="red")
lines(auc3,col="green")
lines(auc4,col="yellow")

4.ROC检验——DeLong’s Test

• 输入

roc.test(auc1,auc2)

• 输出

Bootstrap test for two correlated ROC curves

data:  auc1 and auc2
D = 0.19387, boot.n = 2000, boot.stratified = 1, p-value =
0.8463
alternative hypothesis: true difference in AUC is not equal to 0
sample estimates:
AUC of roc1 AUC of roc2 
     0.9326      0.9248

5.四组间整体检验

• 输入

compareROCdep(iris2[,-5], iris2$label, method="auc")  ##四组间整体检验

• 输出

> compareROCdep(iris2[,-5], iris2$label, method="auc")  ##四组间整体检验
In the considered database there are 50 controls and 50 cases.

Method considered: AUC comparison (DeLong, DeLong and Clarke-Pearson, 1988)

Progress bar: Estimation of statistic value in each variable (k = 4)
  |=======================================================================================================================================| 100%
Error in solve.default(M) : 
  Lapack例行程序dgesv: 系统正好是奇异的: U[3,3] = 0

6.多元统计课程案例实践

先上代码：

#----------------------------绘制ROC曲线图并展示预测效果评判指标-------------------------------------#
plot.myROC <- function(class.true,class.pred){   # class.pred必须为numeric或者order
  # 绘制ROC曲线
  library(pROC,quietly = T)
  modelroc <- roc(class.true,class.pred,levels=base::levels(as.factor(class.true)),direction=c("<"))  
  plot(modelroc, print.auc=TRUE, auc.polygon=TRUE, grid=c(0.1, 0.2),
       grid.col=c("green", "red"), max.auc.polygon=TRUE,
       auc.polygon.col="skyblue", print.thres=TRUE)
  
  # 展示混淆矩阵与准确率
  tab <- table(class.true,class.pred)   # 混淆矩阵
  cat("该方法的混淆矩阵为：",tab,"(从左至右依次为TP、FP、FN、TN)")
  cat("\nFPR = ",tab[2,1]/(tab[2,1]+tab[2,2]),"\t","TPR = ",tab[1,1]/(tab[1,1]+tab[2,2]))
  cat("\n该方法的准确率为：",sum(diag(prop.table(tab)))*100,"%")
}

# Step1：准备数据
path <- 'E:/桌面文档/学习文件/大三/多元统计/课程设计/guocaiyang.csv'
data <- read.table(path,header = F,sep=',')
names(data) <- c('性别','年级','时长','使用频率','使用数量','是否为重度抑郁情绪','专业类别')

x <- data[,-7]  # 输入变量
y <- data[,7]  # 输出变量

# Step2：主成分分析
PCA <- princomp(x,cor=T)   # 从R出发
PCA$loadings    # 载荷矩阵
summary(PCA)

PCA2 <- princomp(x)  # 从Sigma出发
PCA2$loadings

summary(PCA2)

# 绘制碎石图
screeplot(PCA2,type='lines')

# Step3：判别分析
data.PCA <- data.frame(y=data[,7],PCA2$scores)
class.true <- data.PCA$y   # 真实类别

library(MASS,quietly = T)
#### Fisher's线性判别分析      
ld1 <- lda(y~Comp.1+Comp.2+Comp.3+Comp.4,data=data.PCA)   
ld1

Z1 <- predict(ld1) #预测所属类别
plot.myROC(class.true,as.numeric(Z1$class))   

#### 二次判别
qd <- qda(y~Comp.1+Comp.2+Comp.3+Comp.4,data=data.PCA)
qd

Z2 <- predict(qd) # 预测
plot.myROC(class.true,as.numeric(Z2$class))   

#### 贝叶斯判别(不等的初始概率)
ld2 <- lda(y~Comp.1+Comp.2+Comp.3+Comp.4,data=data.PCA,prior=c(0.82,0.18))   
ld2

Z3 <- predict(ld2) # 预测
plot.myROC(class.true,as.numeric(Z3$class))   

# Step4：聚类分析
#### k-means聚类
km <- kmeans(scale(PCA2$scores[,1:4]), 2)
km

plot.myROC(class.true,as.numeric(km$cluster)) 

#### 系统聚类法(完全连接法)
hc <- hclust(dist(scale(PCA2$scores[,1:4])),'complete')  # 完全连接法
plot(hc,hang=-1,sub='');rect.hclust(hc,2)  # 聚类树状图

classpredict <- cutree(hc,2)  # 分2类
plot.myROC(class.true,as.numeric(classpredict))

#### 系统聚类法(简单连接法)
hc2 <- hclust(dist(scale(PCA2$scores[,1:4])),'single') 
plot(hc2,hang=-1,sub='');rect.hclust(hc2,2)  # 聚类树状图

classpredict2 <- cutree(hc2,2)  # 分2类
plot.myROC(class.true,as.numeric(classpredict2))

#### 系统聚类法(ward.D)
hc3 <- hclust(dist(scale(PCA2$scores[,1:4])),'ward.D')  
plot(hc3,hang=-1,sub='');rect.hclust(hc3,2)  # 聚类树状图

classpredict3 <- cutree(hc3,2)  # 分2类
plot.myROC(class.true,as.numeric(classpredict3))

#### 系统聚类法(平均距离法)
hc4 <- hclust(dist(scale(PCA2$scores[,1:4])),'average')  
plot(hc4,hang=-1,sub='');rect.hclust(hc4,2)  # 聚类树状图

classpredict4 <- cutree(hc4,2)  # 分2类
plot.myROC(class.true,as.numeric(classpredict4))

#### 系统聚类法(centroid，重心法)
hc4 <- hclust(dist(scale(PCA2$scores[,1:4])),'centroid') 
plot(hc4,hang=-1,sub='');rect.hclust(hc4,2)  # 聚类树状图

classpredict5 <- cutree(hc4,2)  # 分2类
plot.myROC(class.true,as.numeric(classpredict5))

得到的ROC曲线和AUC值为：