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

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

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

如果你还不太了解关于ROC曲线中的各种指标，请看下面这张图（请至公号查看），有你需要的一切（建议保存）

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

文章目录

• • 方法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, vardata(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)

可以显示最佳截点，比如AUC最大的点：

plot(res,legacy.axes = TRUE,thresholds="best", # AUC最大的点print.thres="best") 

可以显示AUC的可信区间：

rocobj <- (aSAH$outcome, aSAH$s100b,main="Confidence intervals", percent=TRUE,ci=TRUE, print.auc=TRUE) 
## Setting levels: control = Good, case = Poor
## Setting direction: controls < casesciobj <- ci.se(rocobj,specificities=seq(0, 100, 5))plot(ciobj, type="shape", col="#1c61b6AA")
plot(ci(rocobj, of="thresholds", thresholds="best")) 

多条ROC曲线画在一起：

rocobj1 <- (aSAH$outcome, aSAH$s100,percent=TRUE, col="#1c61b6")
## Setting levels: control = Good, case = Poor
## Setting direction: controls < cases
rocobj2 <- (aSAH$outcome, aSAH$ndka, percent=TRUE, col="#008600")
## Setting levels: control = Good, case = Poor
## Setting direction: controls < caseslegend("bottomright", legend=c("S100B", "NDKA"), col=c("#1c61b6", "#008600"), lwd=2)

两条ROC曲线的比较，可以添加P值：

rocobj1 <- (aSAH$outcome, aSAH$s100,percent=TRUE, col="#1c61b6")
## Setting levels: control = Good, case = Poor
## Setting direction: controls < cases
rocobj2 <- (aSAH$outcome, aSAH$ndka, percent=TRUE, col="#008600")
## Setting levels: control = Good, case = Poor
## Setting direction: controls < caseslegend("bottomright", legend=c("S100B", "NDKA"), col=c("#1c61b6", "#008600"), lwd=2)testobj <- st(rocobj1, rocobj2)text(50, 50, labels=paste("p-value =", format.pval(testobj$p.value)), adj=c(0, .5))

方法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")

添加箱线图：

perf <- performance(pred, "tpr", "fpr")
perf
## A performance instance
##   'False positive rate' vs. 'True positive rate' (alpha: 'Cutoff')
##   with 51 data pointsplot(perf,avg="threshold",spread.estimate="boxplot")

还可以绘制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)

还可以把特异度为横坐标，灵敏度为纵坐标：

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)

这个包还可以计算非常多其他的指标，各种图都能画，大家可以自己探索。

方法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()

还有太多方法可以画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