阿越1229

V1

2022/09/08阅读:21主题:自定义主题1

临床预测模型之二分类资料的ROC曲线绘制

本文首发于公众号:医学和生信笔记

医学和生信笔记,专注R语言在临床医学中的使用,R语言数据分析和可视化。主要分享R语言做医学统计学、meta分析、网络药理学、临床预测模型、机器学习、生物信息学等。

ROC曲线是评价模型的重要工具,曲线下面积AUC可能是大家最常见的模型评价指标之一。

如果你还不太了解关于ROC曲线中的各种指标,请看下面这张图,有你需要的一切(建议保存):

混淆矩阵
混淆矩阵
混淆矩阵计算
混淆矩阵计算

R语言中有非常多的方法可以实现ROC曲线,但是基本上都是至少需要2列数据,一列是真实结果,另一列是预测值,有了这两列数据,就可以轻松使用各种方法画出ROC曲线并计算AUC。

这篇文章带大家介绍最常见的并且好用的二分类变量的ROC曲线画法。

方法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
plot of chunk unnamed-chunk-3

可以显示最佳截点,比如AUC最大的点:

plot(res,
     legacy.axes = TRUE,
     thresholds="best"# AUC最大的点
     print.thres="best"
plot of chunk unnamed-chunk-4
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(01005)
               )

plot(ciobj, type="shape", col="#1c61b6AA")
plot(ci(rocobj, of="thresholds", thresholds="best")) 
plot of chunk unnamed-chunk-5
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
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(5050, labels=paste("p-value =", format.pval(testobj$p.value)), adj=c(0.5))
plot of chunk unnamed-chunk-7
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
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
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
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
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
plot of chunk unnamed-chunk-15

还有太多方法可以画ROC了,不过pROCROCR基本上技能解决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

太多了,千万不要迷失哦~

本文首发于公众号:医学和生信笔记

医学和生信笔记,专注R语言在临床医学中的使用,R语言数据分析和可视化。主要分享R语言做医学统计学、meta分析、网络药理学、临床预测模型、机器学习、生物信息学等。

分类:

其他

标签:

医学

作者介绍

阿越1229
V1

黄金矿工。