R语言——ROC 终结者- - 科研猫

# ROC 曲线——科研猫

数据解释

第一列为二分类结局变量，第2~n列为模型预测概率。


代码如下：

options(stringsAsFactors = F)

options(warn = -1)

# 安装包

#install.packages("pROC")

#install.packages("wesanderson")

#install.packages("openxlsx")

# 加载包

library(pROC)

library(wesanderson)

library(openxlsx)

# 数据要求，第一列为 y，后面为 x_1~x_n

data = read.xlsx("data.xlsx",sheet=1,rowNames = F,colNames = T,startRow = 1,

                detectDates = F,na.strings = "#NA")

# check groups (should be <= 5)

num.value = ncol(data) - 1

if(num.value > 5){

  print("Error: Too many groups of observations!")

  q(save="no")

}

name.value = colnames(data)[2:ncol(data)]

# define favoritable colors

col = wes_palette("Darjeeling1",num.value,type=c("discrete"))

# ROC curve

pdf("1.ROC_Combined.pdf",height = 4,width = 4)

for (i in 1:num.value){

  if(i == 1){

    roc.data = roc(data[,1],data[,i+1],

                  percent=T,plot=T, grid=T,lty=i,

                  print.auc=F,col=col[i])

    text(30,50,"AUC",font = 2,col="darkgray")

    text(30,50-10*i,

        paste(name.value[i],":",sprintf("%0.4f",as.numeric(roc.data$auc))),

        col=col[i])

  }else{

    roc.data = roc(data[,1],data[,i+1],

                  percent=T,plot=T, grid=T, add=T,lty=i,

                  print.auc=F,col=col[i])

    text(30,50-10*i,

        paste(name.value[i],":",sprintf("%0.4f",as.numeric(roc.data$auc))),

        col=col[i])

  }

}

dev.off()

# smooth ROC curve

pdf("2.ROC_with_smooth_curve_Combined.pdf",height = 4,width = 4)

for (i in 1:num.value){

  if(i == 1){

    roc.data = roc(data[,1],data[,i+1],

                  percent=T,plot=T, grid=T,lty=i,smooth = T,

                  print.auc=F,col=col[i])

    text(30,50,"AUC",font = 2,col="darkgray")

    text(30,50-10*i,

        paste(name.value[i],":",sprintf("%0.4f",as.numeric(roc.data$auc))),

        col=col[i])

  }else{

    roc.data = roc(data[,1],data[,i+1],

                  percent=T,plot=T, grid=T, add=T,lty=i,smooth = T,

                  print.auc=F,col=col[i])

    text(30,50-10*i,

        paste(name.value[i],":",sprintf("%0.4f",as.numeric(roc.data$auc))),

        col=col[i])

  }

}

dev.off()

## ROC with Confidence intervals ##

transparentColor <- function(col,alpha=200){

  # alpha is an integer >= 1 and <= 255


  col.rgb = as.numeric(col2rgb(col))

  col.rgb.alpha = rgb(col.rgb[1],col.rgb[2],col.rgb[3],alpha=alpha,maxColorValue = 255)

  return(col.rgb.alpha)

}

## CI of the curve

# curve shape

pdf("3.ROC_with_confidence_level_ribbon.pdf",height = 3,width = 3*num.value)

par(mfrow=c(1,num.value))

for (i in 1:num.value){

  roc.data = roc(data[,1],data[,i+1],

                percent=T,plot=T, grid=T,

                print.auc=F,col=transparentColor("white",255))

  sens.ci <- ci.se(roc.data, boot.n=100,conf.level = 0.95,specificities=seq(0, 100, 5))

  plot(sens.ci,type="shape",col=transparentColor(col[i],100))


  text(30,40,"AUC",font = 2,col="darkgray")

  text(30,30,

      paste(name.value[i],":",sprintf("%0.4f",as.numeric(roc.data$auc))),

      col=col[i])

}

dev.off()

#curve bar

pdf("4.ROC_with_confidence_level_bars_Combined.pdf",height = 3,width = 3*num.value)

par(mfrow=c(1,num.value))

for (i in 1:num.value){

  roc.data = roc(data[,1],data[,i+1],

                percent=T,plot=T, grid=T,

                print.auc=F,col=col[i])

  sens.ci <- ci.se(roc.data, boot.n=100,conf.level = 0.95,specificities=seq(0, 100, 5))

  plot(sens.ci,type="bars",col=col[i])


  text(30,40,"AUC",font = 2,col="darkgray")

  text(30,30,

      paste(name.value[i],":",sprintf("%0.4f",as.numeric(roc.data$auc))),

      col=col[i])

}

dev.off()

## Performance ##

getROCPerformance <- function(rocdata){

  performance = coords(roc.data, "best", best.method = "youden",

                      ret=c("threshold", "sensitivity","specificity",

                            "npv","ppv","tpr","fpr",

                            "tnr","fnr","fdr","accuracy",

                            "precision","youden"),

                      transpose = F)

  res = t(as.data.frame(performance))

  auc = as.numeric(roc.data$auc)

  res = rbind(auc,res)

  return(res)

}

perf.all = NULL

roc.data.all = NULL

for (i in 1:num.value){

  roc.data = roc(data[,1],data[,i+1],

                percent=T,plot=F)

  roc.data.all = c(roc.data.all,roc.data)

  perf = getROCPerformance(roc.data)

  colnames(perf) = name.value[i]

  perf.all = cbind(perf.all,perf)

}

perf.all = as.data.frame(perf.all)

perf.all$Index = c("AUC","Best Cut-off Value","Sensitivity","Specificity",

                  "Negative Predictive Value","Positive Predictive Value",

                  "True Positive Rate","False Positive Rate",

                  "True Negatice Rate","False Negative Rate",

                  "False Discovery Rate","Accuracy","Precision","Youden Index")

perf.all = perf.all[,c(ncol(perf.all),1:(ncol(perf.all)-1))]

write.table(perf.all,"5.Model.Performance.csv",row.names = F,col.names = T,quote=F,sep=",")

## Delong Comparisons ##

pair = t(combn(2:ncol(data),2))

compair.result = NULL

for (i in 1:nrow(pair)){

  name1 = colnames(data)[pair[i,1]]

  compair.result$Model1.ROC = c(compair.result$Model1.ROC,name1)

  roc1 = roc(data[,1],data[,pair[i,1]],

            percent=T,plot=F)


  name2 = colnames(data)[pair[i,2]]

  compair.result$Model2.ROC = c(compair.result$Model2.ROC,name2)

  roc2 = roc(data[,1],data[,pair[i,2]],

            percent=T,plot=F)


  p = roc.test(roc1, roc2, reuse.auc=FALSE, boot.n = 1000, boot.stratified = F)$p.value

  compair.result$Pvalue.Delong.test = c(compair.result$Pvalue.Delong.test,p)

}

compair.result = as.data.frame(compair.result)

write.table(compair.result,"6.Delong.Comparision.csv",row.names = F,col.names = T,quote=F,sep=",")

print('finish')

最后编辑于 ：
©著作权归作者所有,转载或内容合作请联系作者
  • 人面猴
    序言：七十年代末，一起剥皮案震惊了整个滨河市，随后出现的几起案子，更是在滨河造成了极大的恐慌，老刑警刘岩，带你破解...
    沈念sama阅读 215,874评论 6赞 498
  • 死咒
    序言：滨河连续发生了三起死亡事件，死亡现场离奇诡异，居然都是意外死亡，警方通过查阅死者的电脑和手机，发现死者居然都...
    沈念sama阅读 92,102评论 3赞 391
  • 救了他两次的神仙让他今天三更去死
    文/潘晓璐 我一进店门，熙熙楼的掌柜王于贵愁眉苦脸地迎上来，“玉大人，你说我怎么就摊上这事。” “怎么了？”我有些...
    开封第一讲书人阅读 161,676评论 0赞 351
  • 道士缉凶录：失踪的卖姜人
    文/不坏的土叔 我叫张陵，是天一观的道长。 经常有香客问我，道长，这世上最难降的妖魔是什么？ 我笑而不...
    开封第一讲书人阅读 57,911评论 1赞 290
  • 港岛之恋（遗憾婚礼）
    正文 为了忘掉前任，我火速办了婚礼，结果婚礼上，老公的妹妹穿的比我还像新娘。我一直安慰自己，他们只是感情好，可当我...
    茶点故事阅读 66,937评论 6赞 388
  • 恶毒庶女顶嫁案：这布局不是一般人想出来的
    文/花漫 我一把揭开白布。 她就那样静静地躺着，像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上，一...
    开封第一讲书人阅读 50,935评论 1赞 295
  • 城市分裂传说
    那天，我揣着相机与录音，去河边找鬼。 笑死，一个胖子当着我的面吹牛，可吹牛的内容都是我干的。 我是一名探鬼主播，决...
    沈念sama阅读 39,860评论 3赞 416
  • 双鸳鸯连环套：你想象不到人心有多黑
    文/苍兰香墨 我猛地睁开眼，长吁一口气：“原来是场噩梦啊……” “哼！你这毒妇竟也来了？” 一声冷哼从身侧响起，我...
    开封第一讲书人阅读 38,660评论 0赞 271
  • 万荣杀人案实录
    序言：老挝万荣一对情侣失踪，失踪者是张志新（化名）和其女友刘颖，没想到半个月后，有当地人在树林里发现了一具尸体，经...
    沈念sama阅读 45,113评论 1赞 308
  • 护林员之死
    正文 独居荒郊野岭守林人离奇死亡，尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 37,363评论 2赞 331
  • 白月光启示录
    正文 我和宋清朗相恋三年，在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 39,506评论 1赞 346
  • 活死人
    序言：一个原本活蹦乱跳的男人离奇死亡，死状恐怖，灵堂内的尸体忽然破棺而出，到底是诈尸还是另有隐情，我是刑警宁泽，带...
    沈念sama阅读 35,238评论 5赞 341
  • 日本核电站爆炸内幕
    正文 年R本政府宣布，位于F岛的核电站，受9级特大地震影响，放射性物质发生泄漏。R本人自食恶果不足惜，却给世界环境...
    茶点故事阅读 40,861评论 3赞 325
  • 男人毒药：我在死后第九天来索命
    文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹，春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 31,486评论 0赞 21
  • 一桩弑父案，背后竟有这般阴谋
    文/苍兰香墨 我抬头看了看天上的太阳。三九已至，却和暖如春，着一层夹袄步出监牢的瞬间，已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 32,674评论 1赞 268
  • 情欲美人皮
    我被黑心中介骗来泰国打工， 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留，地道东北人。 一个月前我还...
    沈念sama阅读 47,513评论 2赞 368
  • 代替公主和亲
    正文 我出身青楼，却偏偏与公主长得像，于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子，可洞房花烛夜当晚...
    茶点故事阅读 44,426评论 2赞 352