2024-06-16

pROC带你轻松搞定ROC曲线绘制

小师妹生信果 2024-02-26 19:01:06

转自公众号：R语言学徒
http://mp.weixin.qq.com/s?__biz=Mzg5MDk3Mzg4OA==&mid=2247487165&idx=1&sn=1223985c1bcd8abb6c13c9fa1cb6510e

为了能够计算和评估医学诊断试验和预测模型性能，我们常常用到一个强大的R包工具：pROC包，它是一个R语言中的统计学分析包，它提供了计算、绘制和分析接受者操作特征曲线（ROC曲线）的函数和方法，并计算某个阈值下的真阳性率和假阳性率。

了解了pROC的功能的强大之后，下面就由小果带大家来学习一下，如何使用使用pROC来绘制ROC曲线吧！

一、数据的准备与pROC包的加载

在这里，我们将采用aSAH数据集作为示例，它是脑动脉瘤研究中的一个经典数据集，其中包含了200个亚实质性蛛网膜下腔出血（aSAH）病例的相关信息。该数据集包括了11个变量，示例包括年龄、颅内压、术前状态等。这个数据集通常用于分析不同变量与患者病情恢复的关系，并且也常用于介绍如何在R中做ROC分析。

# 加载MASS包            library(MASS)            

            # 加载pROC包            library(pROC)            

            # 加载aSAH数据集            data(aSAH)

需要注意的是，aSAH数据集是MASS包中的一个示例数据集，在使用前，需要加载安装并MASS包。以下为aSAH的数据的展示。

二、使用roc()进行ROC分析

在aSAH数据集中，我们用age这个数值型预测变量来预测outcome二分类因变量的取值，并且指定因变量的两个水平分别为”Good”和”Poor”。roc1是二分类预测变量age在预测outcome二分类因变量时的ROC曲线和AUC值。我们可以通过这个分析来评估age对outcome的预测能力。

#使用roc()来建立ROC对象            roc1 <- roc(aSAH$outcome, aSAH$age,                        levels=c("Good", "Poor"))            

            # 绘制ROC曲线并进行个性化设置            plot(roc1,                  # 打印AUC值                 print.auc=TRUE,                  # 设置打印AUC值的x轴位置                 print.auc.x=0.4,                  # 设置打印AUC值的y轴位置                 print.auc.y=0.5,                 # 不生成ROC曲线下方多边形                 auc.polygon=FALSE,                  # 设置两个轴间隔的大小                 grid=c(0.4, 0.3),                  # 设置两个轴间隔的颜色                 grid.col=c("#FF8C00", "#FF8C00"),                  # 不显示最佳阈值的位置                 print.thres=FALSE,                  # 图形标题                 main="ROC curve for aSAH dataset",                 # ROC曲线颜色                 col="#FF2E63",                  # 绘制旧式的轴刻度标签                 legacy.axes=TRUE,                  # 打印AUC值时使用橙色文本                 print.auc.col="#FF8C00")

接着，我们使用plot()，它是一个R语言基础内置函数，它的作用是用于将数据可视化，并生成各种类型的图形或图表。在画出roc1对象所代表的ROC曲线图之后，我们可以得到如图所示的ROC曲线图。

三、曲线的比较与p值的添加

为了检验两条曲线的p值，我们再次引入另一条曲线，参数更改为aSAH$s100b，并使用 add=TRUE将它与上一条曲线放置在同一图上，代码以及图例如下：

#使用roc()来建立ROC对象            roc2 <- roc(aSAH$outcome, aSAH$s100b, levels = c("Good", "Poor"))            

            # 绘制 ROC 曲线            plot(roc2, col = "blue", add=TRUE, # 增加曲线，                 print.auc = TRUE,                 main = "ROC curve for S100B",                 xlab = "1 - Specificity", ylab = "Sensitivity",                 print.auc.x=0.6,print.auc.y=0.6)

紧接着，我们使用 roc.test() 函数对前面得到的两条 ROC 曲线 roc1 和 roc2 进行比较，并得出它们是否具有显著差异。函数返回的结果存储在变量 testp 中。然后我们使用 text() 函数在图表上添加了一个注释，显示了检验结果中的 p 值。

testp <- roc.test(roc1,roc2)   # 检验两条曲线            

            text(0.4, 0.2,  # 设置添加p值的位置(x,y)                 labels=paste("P value =", format.pval(testp$p.value)),                  adj=c(0, 0.5)) # 在图上添加P值