临床预测模型之二分类资料ROC曲线绘制
模型 分类 预测 绘制 资料 曲线 之二 临床
2023-06-13 09:15:06 时间
ROC曲线是评价模型的重要工具，曲线下面积AUC可能是大家最常见的模型评价指标之一。
如果你还不太了解关于ROC曲线中的各种指标，请看下面这张图，有你需要的一切（建议保存）：
混淆矩阵
混淆矩阵计算
R语言中有非常多的方法可以实现ROC曲线，但是基本上都是至少需要2列数据，一列是真实结果，另一列是预测值，有了这两列数据，就可以轻松使用各种方法画出ROC曲线并计算AUC。
这篇文章带大家介绍最常见的并且好用的二分类变量的ROC曲线画法。
- 方法1
- 方法2
- 方法3
方法1
使用pROC包，不过使用这个包需要注意，一定要指定direction，否则可能会得出错误的结果。
这个R包计算AUC是基于中位数的，哪一组的中位数大就计算哪一组的AUC，在计算时千万要注意！
使用pROC包的aSAH数据，其中outcome列是结果变量，1代表Good，2代表Poor。
library(pROC)
## Type 'citation("pROC")' for a citation.
##
## 载入程辑包：'pROC'
## The following objects are masked from 'package:stats':
##
## cov, smooth, var
data(aSAH)
dim(aSAH)
## [1] 113 7
str(aSAH)
## 'data.frame': 113 obs. of 7 variables:
## $ gos6 : Ord.factor w/ 5 levels "1"<"2"<"3"<"4"<..: 5 5 5 5 1 1 4 1 5 4 ...
## $ outcome: Factor w/ 2 levels "Good","Poor": 1 1 1 1 2 2 1 2 1 1 ...
## $ gender : Factor w/ 2 levels "Male","Female": 2 2 2 2 2 1 1 1 2 2 ...
## $ age : int 42 37 42 27 42 48 57 41 49 75 ...
## $ wfns : Ord.factor w/ 5 levels "1"<"2"<"3"<"4"<..: 1 1 1 1 3 2 5 4 1 2 ...
## $ s100b : num 0.13 0.14 0.1 0.04 0.13 0.1 0.47 0.16 0.18 0.1 ...
## $ ndka : num 3.01 8.54 8.09 10.42 17.4 ...
计算AUC及可信区间：
res <- roc(aSAH$outcome,aSAH$s100b,ci=T,auc=T)
## Setting levels: control = Good, case = Poor
## Setting direction: controls < cases
res
##
## Call:
## roc.default(response = aSAH$outcome, predictor = aSAH$s100b, auc = T, ci = T)
##
## Data: aSAH$s100b in 72 controls (aSAH$outcome Good) < 41 cases (aSAH$outcome Poor).
## Area under the curve: 0.7314
## 95% CI: 0.6301-0.8326 (DeLong)
plot(res,legacy.axes = TRUE)
plot of chunk unnamed-chunk-3
可以显示最佳截点，比如AUC最大的点：
plot(res,
legacy.axes = TRUE,
thresholds="best", # AUC最大的点
print.thres="best")
plot of chunk unnamed-chunk-4
可以显示AUC的可信区间：
rocobj <- plot.roc(aSAH$outcome, aSAH$s100b,
main="Confidence intervals",
percent=TRUE,ci=TRUE,
print.auc=TRUE
)
## Setting levels: control = Good, case = Poor
## Setting direction: controls < cases
ciobj <- ci.se(rocobj,
specificities=seq(0, 100, 5)
)
plot(ciobj, type="shape", col="#1c61b6AA")
plot(ci(rocobj, of="thresholds", thresholds="best"))
plot of chunk unnamed-chunk-5
多条ROC曲线画在一起：
rocobj1 <- plot.roc(aSAH$outcome, aSAH$s100,percent=TRUE, col="#1c61b6")
## Setting levels: control = Good, case = Poor
## Setting direction: controls < cases
rocobj2 <- lines.roc(aSAH$outcome, aSAH$ndka, percent=TRUE, col="#008600")
## Setting levels: control = Good, case = Poor
## Setting direction: controls < cases
legend("bottomright", legend=c("S100B", "NDKA"), col=c("#1c61b6", "#008600"), lwd=2)
plot of chunk unnamed-chunk-6
两条ROC曲线的比较，可以添加P值：
rocobj1 <- plot.roc(aSAH$outcome, aSAH$s100,percent=TRUE, col="#1c61b6")
## Setting levels: control = Good, case = Poor
## Setting direction: controls < cases
rocobj2 <- lines.roc(aSAH$outcome, aSAH$ndka, percent=TRUE, col="#008600")
## Setting levels: control = Good, case = Poor
## Setting direction: controls < cases
legend("bottomright", legend=c("S100B", "NDKA"), col=c("#1c61b6", "#008600"), lwd=2)
testobj <- roc.test(rocobj1, rocobj2)
text(50, 50, labels=paste("p-value =", format.pval(testobj$p.value)), adj=c(0, .5))
plot of chunk unnamed-chunk-7
方法2
使用ROCR，如果你只是为了画一条ROC曲线，这是我最推荐的方法了，美观又简单！
library(ROCR)
使用非常简单，3句代码，其中第2句是关键，可以更改各种参数，然后就可以画出各种不同的图形：
pred <- prediction(aSAH$s100b,aSAH$outcome)
perf <- performance(pred, "tpr","fpr")
auc <- round(performance(pred, "auc")@y.values[[1]],digits = 4)
plot(perf,lwd=2,col="tomato")
abline(0,1,lty=2)
legend("bottomright", legend="AUC of s100b: 0.7314", col="tomato", lwd=2,bty = "n")
plot of chunk unnamed-chunk-9
添加箱线图：
perf <- performance(pred, "tpr", "fpr")
perf
## A performance instance
## 'False positive rate' vs. 'True positive rate' (alpha: 'Cutoff')
## with 51 data points
plot(perf,
avg="threshold",
spread.estimate="boxplot")
plot of chunk unnamed-chunk-10
还可以绘制PR曲线，召回率recall为横坐标，精确率precision 为纵坐标：
perf <- performance(pred, "prec", "rec")
plot(perf,
avg= "threshold",
colorize=TRUE,
lwd= 3,
main= "Precision-Recall plot")
plot(perf,
lty=3,
col="grey78",
add=TRUE)
plot of chunk unnamed-chunk-11
还可以把特异度为横坐标，灵敏度为纵坐标：
perf <- performance(pred, "sens", "spec")
plot(perf,
avg= "threshold",
colorize=TRUE,
lwd= 3,
main="Sensitivity/Specificity plots")
plot(perf,
lty=3,
col="grey78",
add=TRUE)
plot of chunk unnamed-chunk-12
这个包还可以计算非常多其他的指标，各种图都能画，大家可以自己探索。
方法3
使用tidymodels。这个包很有来头，它是R中专门做机器学习的，我很快就会详细介绍它，它也是目前R语言机器学习领域两大当红辣子鸡之一！另一个是mlr3。
suppressPackageStartupMessages(library(tidymodels))
它很优雅，如果你要计算AUC，那么就是roc_auc()函数：
aSAH %>% roc_auc(outcome, s100b,event_level="second")
## # A tibble: 1 x 3
## .metric .estimator .estimate
## <chr> <chr> <dbl>
## 1 roc_auc binary 0.731
如果你是要画ROC曲线，那么就是roc_curve()函数：
aSAH %>% roc_curve(outcome, s100b,event_level="second") %>%
ggplot(aes(x = 1 - specificity, y = sensitivity)) +
geom_path(size=1.2,color="firebrick") +
geom_abline(lty = 3) +
coord_equal() +
theme_bw()
plot of chunk unnamed-chunk-15
还有太多方法可以画ROC了，不过pROC和ROCR基本上技能解决99%的问题了。
最后，给大家看看cran中比较常见的画ROC曲线的包，大家有兴趣可以自己探索：
library(pkgsearch)
rocPkg <- pkg_search(query="ROC",size=200)
rocPkgShort <- rocPkg %>%
filter(maintainer_name != "ORPHANED") %>%
select(score, package, downloads_last_month) %>%
arrange(desc(downloads_last_month))
head(rocPkgShort,20)
## # A data frame: 20 x 3
## score package downloads_last_month
## * <dbl> <chr> <int>
## 1 4711. caTools 172954
## 2 10868. pROC 123923
## 3 932. ROCR 82311
## 4 1551. plotROC 9719
## 5 2471. PRROC 7209
## 6 282. riskRegression 6622
## 7 173. InformationValue 4710
## 8 1950. survivalROC 4096
## 9 421. PresenceAbsence 3534
## 10 2068. cvAUC 3144
## 11 1594. timeROC 2570
## 12 779. risksetROC 2253
## 13 316. mlr3viz 2243
## 14 109. RcmdrPlugin.EZR 1953
## 15 198. logcondens 1947
## 16 427. ROCit 1744
## 17 1469. precrec 1497
## 18 183. WVPlots 1491
## 19 205. hmeasure 1013
## 20 67.4 wrProteo 965
方法太多了，千万不要迷失哦~
以上就是今天的内容，希望对你有帮助哦！
相关文章
- 无人机集群自组织搜索仿真模型设计与实现
- CSS小技能:常用样式属性、选择器分类、盒子模型
- 用GFlowNets统一生成模型，Bengio等人数页论文给讲通了
- 扩散模型爆火，这是首篇综述与Github论文分类汇总
- 预测三分类变量模型的ROC介绍
- 分类模型的f1值大概是多少_准确率召回率f1分数
- 机器学习模型评估的方法总结（回归、分类模型的评估）
- python分类模型_nlp模型评估指标
- 测试数据增强_预测模型最佳cutoff值
- 应用实践：Paddle分类模型大集成者[PaddleHub、Finetune、prompt]
- 使用阈值调优改进分类模型性能
- 【视频】R语言逻辑回归（Logistic回归）模型分类预测病人冠心病风险|数据分享|附代码数据
- Paddlenlp之UIE分类模型【以情感倾向分析新闻分类为例】含智能标注方案）
- 【视频】R语言逻辑回归（Logistic回归）模型分类预测病人冠心病风险|数据分享|附代码数据
- Meta发布首个「非参数化」掩码语言模型NPM：吊打500倍参数量的GPT-3
- 语言模型自己学会用搜索引擎了？Meta AI提出API调用自监督学习方法Toolformer
- 【数据挖掘】拉普拉斯修正 ( 判别模型 | 概率模型 | 贝叶斯分类 | 拉普拉斯修正 | 朴素贝叶斯分类应用场景 | 朴素贝叶斯优缺点 )
- 【组合数学】生成函数 ( 正整数拆分 | 正整数拆分基本模型 | 有限制条件的无序拆分 )
- 【计算理论】计算复杂性 ( 阶段总结 | 计算理论内容概览 | 计算问题的有效性 | 语言与算法模型 | 可计算性与可判定性 | 可判定性与有效性 | 语言分类 ) ★
- 机器学习模型以出色的精度进行有机反应机理分类
- 分群思维（三）基于帕累托模型的渠道分类
- 领英创始人里德·霍夫曼退出OpenAI；微软推出多模态模型 Kosmos-1；谷歌推出PaLM-E
- 吴恩达的2022年终盘点：生成式AI、ViT、大模型
- 分割一切后，Segment Anything又能分辨类别了：Meta/UTAustin提出全新开放类分割模型
- K-3D是基于GNU/Linux和Win32的一个三维建模、动画和绘制系统，是一款免费、开放原始码的 3D 模型和动画制作与渲染 (rendering) 工具，它强大的功能可以满足专业人士的需求。它可以创建和编辑 3D 几何图形，提供极具弹性的面向对象的插件增强功能及以节点作基础的可视化管线架构，所有参数和选项的调整，都会立即显现结果，而且可以无限次数地复原与取消复原。此外，它使用与 RenderMan 相符的渲染引擎 (render engine)，可创作出电影质量的 3D 动画。