虽然比较基础,但还是要讲的,毕竟名字叫:零基础说科研,还是要照顾一下基础差的粉丝,也算不忘初心把
ROC曲线也叫受试者工作曲线,原来用在军事雷达中,后面广泛应用于医学统计中。ROC曲线是根据一系列不同的二分类方式(分界值或决定阈),以真阳性率(灵敏度)为纵坐标,假阳性率(1-特异度)为横坐标绘制的曲线。
ROC曲线主要应用于二分类结局,比如是否死亡,疾病诊断,肿瘤复发等等,可以用于自变量为连续变量的截点判定。
可以做ROC曲线的软件很多,如SPSS,R语言,Stata, SAS等等,既往我们已经介绍了使用SPSS制作ROC曲线,今天我们来视频介绍怎么使用R语言绘制多条ROC曲线,继续使用我们的肺炎数据(公众号回复:肺炎,可以获得数据)
手把手教你R语言基于pROC包绘制多条ROC曲线
代码
library(pROC)
#导入数据:公众号回复:肺炎数据,可以获得这个数据
bc<-read.csv("E:/r/test/feiyan.csv",sep=',',header=TRUE)
names(bc)
##生成roc数据
roc1<- roc(bc$fy, bc$TNF)
plot(roc1)
plot(smooth(roc1))
plot(smooth(roc1),col="red")
########
plot(roc1,col="red",legacy.axes=T)##更改X轴格式
plot(roc1,col="red",legacy.axes=T,print.auc=TRUE)##显示AUC面积
plot(roc1,col="red",legacy.axes=T,print.auc=TRUE,print.thres=TRUE)##添加截点和95%CI
###另一种画法
plot(roc1, print.auc=TRUE, auc.polygon=TRUE, grid=c(0.1, 0.2),
grid.col=c("green", "red"), max.auc.polygon=TRUE,
auc.polygon.col="lightblue", print.thres=TRUE)
####
plot(roc1, print.thres="best", print.thres.best.method="youden")
plot(roc1, print.thres="best", print.thres.best.method="closest.topleft")
#####多个ROC曲线
roc2<- roc(bc$fy, bc$IL6)
roc3<- roc(bc$fy, bc$sTREM1)
roc4<- roc(bc$fy, bc$PCT)
####
plot(roc1,col="red",legacy.axes=T)
plot(roc2, add=TRUE, col="blue")
plot(roc3, add=TRUE, col="green")
plot(roc4, add=TRUE, col="black")
legend("bottomright", legend=c("roc1","roc2","roc3","roc4"),col=c("red","blue","green","black"),lty=1)
#####
round(auc(roc1),3)##AUC
round(ci(roc1),3)##95%CI
round(auc(roc2),3)##AUC
round(ci(roc2),3)##95%CI
round(auc(roc3),3)##AUC
round(ci(roc3),3)##95%CI
round(auc(roc4),3)##AUC
round(ci(roc4),3)##95%CI
legend("bottomright", legend=c("roc1-auc0.637","roc2-auc0.657","roc3-auc0.774","roc4-auc0.872"),
col=c("red","blue","green","black"),lty=1)
###########可信区间
#######
rocobj <- plot.roc(bc$fy, bc$TNF)
ci.auc(rocobj)
ci.sp.obj <- ci.sp(rocobj, sensitivities=seq(0, 1, .01), boot.n=100)
plot(rocobj) # restart a new plot
plot(ci.sp.obj, type="shape", col="blue")
legend("bottomright",legend="AUC=0.707(0.6313-0.7829 )", lty = 1, col = 'black')
function (n, roc, se)
{
controls <- sample(roc$controls, replace = TRUE)
cases <- sample(roc$cases, replace = TRUE)
thresholds <- roc.utils.thresholds(c(cases, controls), roc$direction)
perfs <- roc$fun.sesp(thresholds = thresholds, controls = controls,
cases = cases, direction = roc$direction)
roc$sensitivities <- perfs$se * ifelse(roc$percent, 100,
1)
roc$specificities <- perfs$sp * ifelse(roc$percent, 100,
1)
roc$thresholds <- thresholds
return(sapply(se, function(x) coords.roc(roc, x, input = "sensitivity",
ret = "specificity", transpose = TRUE)))
}