作者:白介素2
相关阅读:
R语言-multiROC package
R语言生存分析04-Cox比例风险模型诊断
R语言生存分析03-Cox比例风险模型
R语言生存分析-02-ggforest
R语言生存分析-01
生存曲线
如果没有时间精力学习代码,推荐了解:零代码数据挖掘课程
pROC package
以下是本包中常用的一些缩写
ROC: receiver operating characteristic,ROC曲线
AUC: area under the ROC curve,曲线下面积
pAUC: partial area under the ROC curve 部分曲线下面积
CI: confidence interval 可信区间
SP: specificity 特异度
SE: sensitivity 灵敏度
require(pROC)
data(aSAH)
if(!require(DT)) install.packages(DT)
DT::datatable(aSAH)
aSAH[1:5,1:5]
roc函数建立roc曲线
- 支持在管道中运行
- 参数分别为data, event, predict marker
library(dplyr)
aSAH %>%
filter(gender == "Female") %>%
roc(outcome, s100b)
Call:
roc.data.frame(data = ., response = outcome, predictor = s100b)
Data: s100b in 50 controls (outcome Good) < 21 cases (outcome Poor).
Area under the curve: 0.72
coords函数中筛选有效的的坐标
- transpose参数指返回值的格式,FALSE 为row
- 这样筛选出了敏感度和特异度>0.6的坐标
library(dplyr)
aSAH %>%
filter(gender == "Female") %>%
roc(outcome, s100b) %>%
coords(transpose=FALSE) %>%
filter(sensitivity > 0.6,
specificity > 0.6)
threshold specificity sensitivity
1 0.155 0.68 0.6666667
2 0.165 0.74 0.6666667
3 0.175 0.76 0.6666667
4 0.185 0.78 0.6666667
5 0.215 0.80 0.6666667
6 0.245 0.82 0.6666667
7 0.255 0.82 0.6190476
建立roc 对象的方法
# Build a ROC object and compute the AUC
roc(aSAH$outcome, aSAH$s100b)
roc(outcome ~ s100b, aSAH)
建立光滑曲线
# Smooth ROC curve
roc(outcome ~ s100b, aSAH, smooth=TRUE)
Call:
roc.formula(formula = outcome ~ s100b, data = aSAH, smooth = TRUE)
Data: s100b in 72 controls (outcome Good) < 41 cases (outcome Poor).
Smoothing: binormal
Area under the curve: 0.74
可信区间与绘图
# more options, CI and plotting
roc1 <- roc(aSAH$outcome,
aSAH$s100b, percent=TRUE,
# arguments for auc
partial.auc=c(100, 90), partial.auc.correct=TRUE,
partial.auc.focus="sens",
# arguments for ci
ci=TRUE, boot.n=100, ci.alpha=0.9, stratified=FALSE,
# arguments for plot
plot=TRUE, auc.polygon=TRUE, max.auc.polygon=TRUE, grid=TRUE,
print.auc=TRUE, show.thres=TRUE)
## 在原有图形上继续绘制
roc2 <- roc(aSAH$outcome, aSAH$wfns,
plot=TRUE, add=TRUE, percent=roc1$percent)
找出感兴趣的坐标
## Coordinates of the curve ##
coords(roc1, "best", ret=c("threshold", "specificity", "1-npv"),transpose = FALSE
)
coords(roc2, "local maximas", ret=c("threshold", "sens", "spec", "ppv", "npv"),transpose = FALSE)
threshold sensitivity specificity ppv npv
local.maximas -Inf 100.00000 0.00000 36.28319 NaN
local.maximas.1 1.5 95.12195 51.38889 52.70270 94.87179
local.maximas.2 2.5 65.85366 79.16667 64.28571 80.28169
local.maximas.3 3.5 63.41463 83.33333 68.42105 80.00000
local.maximas.4 4.5 43.90244 94.44444 81.81818 74.72527
local.maximas.5 Inf 0.00000 100.00000 NaN 63.71681
计算AUC可信区间
# CI of the AUC
ci(roc2)
95% CI: 74.85%-89.88% (DeLong)
plot在原有图形上增加
- add=TRUE参数
roc1 <- roc(aSAH$outcome,
aSAH$s100b, percent=TRUE,
# arguments for auc
partial.auc=c(100, 90), partial.auc.correct=TRUE,
partial.auc.focus="sens",
# arguments for ci
ci=TRUE, boot.n=100, ci.alpha=0.9, stratified=FALSE,
# arguments for plot
plot=TRUE, auc.polygon=TRUE, max.auc.polygon=TRUE, grid=TRUE,
print.auc=TRUE, show.thres=TRUE)
plot(roc2, add=TRUE)
比较AUC
- 看是否有统计学意义
# Test on the whole AUC
roc.test(roc1, roc2, reuse.auc=FALSE)
DeLong's test for two correlated ROC curves
data: roc1 and roc2
Z = -2.209, p-value = 0.02718
alternative hypothesis: true difference in AUC is not equal to 0
sample estimates:
AUC of roc1 AUC of roc2
73.13686 82.36789
绘制ROC曲线-基于ggplot2
- 创建roc对象
- ggroc绘图
# Create a basic roc object
data(aSAH)
rocobj <- roc(aSAH$outcome, aSAH$s100b)
rocobj2 <- roc(aSAH$outcome, aSAH$wfns)
绘图
- 基础绘图
library(ggplot2)
g <- ggroc(rocobj)
g
- 美化参数设置
ggroc(rocobj, alpha = 0.5, colour = "red", linetype = 2, size = 2)
支持gglot2语法的美化
# You can then your own theme, etc.
g + theme_minimal() + ggtitle("My ROC curve") +
geom_segment(aes(x = 1, xend = 0, y = 0, yend = 1), color="grey", linetype="dashed")
修改横纵坐标
# And change axis labels to FPR/FPR
gl <- ggroc(rocobj, legacy.axes = TRUE)
gl
gl + xlab("FPR") + ylab("TPR") +
geom_segment(aes(x = 0, xend = 1, y = 0, yend = 1), color="darkgrey", linetype="dashed")
绘制多条曲线
- ggroc以list格式包裹roc对象
# Multiple curves:
g2 <- ggroc(list(s100b=rocobj, wfns=rocobj2, ndka=roc(aSAH$outcome, aSAH$ndka)))
g2
- 也可先构建好公式,再绘制
# This is equivalent to using roc.formula:
roc.list <- roc(outcome ~ s100b + ndka + wfns, data = aSAH)
g.list <- ggroc(roc.list)
g.list
美化修改
- size设置线条粗细
- alpha设置透明度
# with additional aesthetics:
g3 <- ggroc(roc.list, size = 1.2,alpha=.6)
g3+ggsci::scale_color_lancet()
改变参数
- aes即按什么属性进行区分
g4 <- ggroc(roc.list, aes="linetype", color="red")
g4
按多种属性区分ROC曲线
# changing multiple aesthetics:
g5 <- ggroc(roc.list, aes=c("linetype", "color"))
g5
分面绘制ROC曲线
# OR faceting
g.list + facet_grid(.~name) + theme(legend.position="none")
所有曲线有相同颜色
- group参数
# To have all the curves of the same color, use aes="group":
g.group <- ggroc(roc.list, aes="group",color="red")
g.group
g.group + facet_grid(.~name)
我是白介素2,本期内容就到这里,下期再见。