我们在做数据分析的时候经常需要用到ROC曲线，我们一般都在下面的两种情况用到ROC曲线：
ROC曲线有一些特定的值都是需要在做ROC的时候给出的。例如：灵敏度；特意度；阈值；曲线下面积以及约登指数等等。
R语言可以实现绘制ROC曲线的有很多，例如：pROC;RROC;cutpointr都可以对一个变量进行ROC计算以及相关数据的统计。这次我们就通过pROC包的相关函数来批量计算ROC曲线的结果。
pROC包当中提供了一个函数叫roc，我们可以通过这个函数来计算两组之间的ROC结果ci.auc函数了计算ROC曲线的曲线下面积以及95%CIcoords来计算ROC曲线的相关信息。这个函数可以接受的参数包括：“threshold”, “specificity”, “sensitivity”, “accuracy”, “tn” (true negative count), “tp” (true positive count), “fn” (false negative count), “fp” (false positive count), “npv” (negative predictive value), “ppv” (positive predictive value), “precision”, “recall”. “1-specificity”, “1-sensitivity”, “1-accuracy”, “1-npv” and “1-ppv”`。我们可以选择不同的参数来返回不同的结果灵敏度+特意度-1。我们可以自己计算。通过上面的需求，我们可以创建一个函数来返回相关的ROC分析的结果
### 加载使用的包
library(tidyverse)## ── Attaching packages ──────────────────────── tidyverse 1.2.1 ──## ✔ ggplot2 3.2.1 ✔ purrr 0.3.2
## ✔ tibble 2.1.3 ✔ dplyr 0.8.3
## ✔ tidyr 0.8.3 ✔ stringr 1.4.0
## ✔ readr 1.3.1 ✔ forcats 0.4.0## ── Conflicts ─────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()library(pROC)## Type 'citation("pROC")' for a citation.##
## Attaching package: 'pROC'## The following objects are masked from 'package:stats':
##
## cov, smooth, varROCStatFunc <- function(dat, group, var,retype = c("threshold", "specificity", "sensitivity"),
auc = T,youden = T, digit = 3){
subgroup <- levels(as.factor(dat[[group]]))
subgroup1 <- paste0(subgroup[2], " vs ", subgroup[1])
rocmodel <- roc(dat[[group]], dat[[var]])
other <- coords(rocmodel, "b", ret = retype)
other <- round(other, digit)
if(auc == T){
auc <- round(ci.auc(rocmodel),digit)
auc <- paste0(auc[2],"(",auc[1],"-",auc[3],")")
if(youden == T){
abc <- coords(rocmodel, "b", ret = c("specificity", "sensitivity"))
youdenres <- abc[1] + abc[2] - 1
youdenres <- round(youdenres, digit)
result <- c(group, subgroup1, auc, other, youdenres)
names(result) <- c("group", "subgroup","auc(95%CI)", retype, "youden")
}else{
result <- c(group, subgroup1, auc, other)
names(result) <- c("group", "subgroup", "auc(95%CI)", retype)
}
}else{
if(youden == T){
abc <- coords(rocmodel, "b", ret = c("specificity", "sensitivity"))
youdenres <- abc[1] + abc[2] - 1
youdenres <- round(youdenres, digit)
result <- c(group, subgroup1, other, youdenres)
names(result) <- c("group","subgroup", retype, "youden")
}else{
result <- c(group, subgroup1,other)
names(result) <- c("group", "subgroup",retype)
}
}
return(result)
}我们构建的函数一共包括7个参数。分别是：
由于不是自己写的函数，我们在使用pROC包的时候会返回很多信息。我们可以去掉这些结果
quiteROCFunc <- quietly(ROCStatFunc)我们使用示例数据来查看结果
data("aSAH")
head(aSAH)## gos6 outcome gender age wfns s100b ndka
## 29 5 Good Female 42 1 0.13 3.01
## 30 5 Good Female 37 1 0.14 8.54
## 31 5 Good Female 42 1 0.10 8.09
## 32 5 Good Female 27 1 0.04 10.42
## 33 1 Poor Female 42 3 0.13 17.40
## 34 1 Poor Male 48 2 0.10 12.75### 计算outcome为结局变量的age的相关信息
quiteROCFunc(aSAH, group = "s100b", var = "age")$result## group subgroup auc(95%CI) threshold specificity
## "s100b" "0.04 vs 0.03" "NA(NA-NA)" "39" "1"
## sensitivity youden
## "0.6" "0.6"# 批量计算变量的ROC结果
## 定义group
multigroup <- c("age", "s100b", "ndka")
rocRes <- lapply(multigroup, function(x) quiteROCFunc(aSAH, "outcome", x)$result)
rocResDat <- do.call(rbind, rocRes)
rocResDat## group subgroup auc(95%CI) threshold specificity
## [1,] "outcome" "Poor vs Good" "0.615(0.508-0.722)" "50.5" "0.569"
## [2,] "outcome" "Poor vs Good" "0.731(0.63-0.833)" "0.205" "0.806"
## [3,] "outcome" "Poor vs Good" "0.612(0.501-0.723)" "11.08" "0.514"
## sensitivity youden
## [1,] "0.634" "0.204"
## [2,] "0.634" "0.44"
## [3,] "0.707" "0.221"PS: 1. 这个函数依赖的是tidyververse以及pROC包，运行之前需要提前加载 2. 函数对于分组的话必须是两组。 3. 结果当中的subgroup是case vs control。