So-net無料ブログ作成
検索選択

ROC曲線をggplot2で描画 [R]

ROC曲線をRの美麗な可視化パッケージであるggplot2で表示するプログラムを作ってみました。

ROCを描画するパッケージとしてEpiというものがあるのだが、ちょっと出力される画像が好みではないのでggplot2でやってみました。三重大学の奥村先生のサイトを参考にしました。

ROC <- function(score, actual) {
  ###ROC-AUC計算
  o = order(score, decreasing = T) 
  fp = tp = fp_prev = tp_prev = 0
  nF = sum(actual == F)
  nT = sum(actual == T)
  score_prev = -Inf
  ber_min = Inf
  area = 0
  rx = ry = numeric(length(o))
  n = 0
  
  #cat("i", "\t", "n", "\t", "score[j]", "\t", "tp", "\t", "fp", "\t", "rx[n]", "\t", "ry[n]", "\t", "area", "\t", "AUC", "\n")
  
  for (i in seq_along(o)) {
    j = o[i]
    if (score[j] != score_prev) {
      area = area + (fp - fp_prev) * (tp + tp_prev) / 2
      n = n + 1
      rx[n] = fp/nF
      ry[n] = tp/nT
      ber = (fp/nF + 1 - tp/nT)/2
      AUC = area/(nF*nT) 
      
      #cat(i, "\t", n, "\t", score[j], "\t", tp, "\t", fp, "\t", rx[n], "\t", ry[n], "\t", area, "\t", AUC, "\n")
      
      if (ber < ber_min) {
        ber_min = ber
        th = score_prev
        rx_best = fp/nF
        ry_best = tp/nT
      }
      score_prev = score[j]
      fp_prev = fp
      tp_prev = tp
    }
    if (actual[j] == T) {
      tp = tp + 1
    } else {
      fp = fp + 1
    }
  }
  area = area + (fp - fp_prev) * (tp + tp_prev) / 2
  AUC = area/(nF*nT)
  n = n + 1
  rx[n] = fp/nF  # = 1
  ry[n] = tp/nT  # = 1
  
  ###ggplot2で可視化
  AUC <- round(AUC, digits = 4)
  tmp <- data.frame(x=rx, y=ry)
  label <- paste("AUC = ", AUC)
  gg <- ggplot(tmp) + geom_path(aes(x,y, color = "red"), size = 1.0) + 
    geom_path(data=data.frame(x = c(0,1), y = c(0,1)),aes(x,y), colour = "gray", size = 1.0) +
    scale_x_continuous("False_Positive Rate (1 - Specificity)", labels=percent, limits = c(0, 1)) +
    scale_y_continuous("True_Positive Rate (Sensivity or Recall)", labels=percent,limits = c(0, 1)) +
    theme(legend.position = "none") +
    ggtitle("ROC Curve") +
    geom_text(data = NULL, x = 0.85, y = 0.15, label = label, colour = "black")
  print(gg)
  
  ###ROC-AUC計算
  cat("\n")
  cat("AUC =", area/(nF*nT), "th =", th, "\n")
  cat("BER =", (rx_best + (1-ry_best))/2,
      "OR =", (ry_best/(1-ry_best))/(rx_best/(1-rx_best)), "\n")
  print(table(score >= th, actual, dnn=c("Predicted","Actual")))
}

####実行例
library(ggplot2)
library(scales) #ggplot2のscaleを%表示にするのに必要

s = c(16,15,14,13,12,11,10, 9, 8, 8, 8, 8, 7, 6, 5) #score
a = c( T, T, F, T, T, T, F, T, T, T, T, F, F, T, F) #actual
ROC(s, a)


【実行結果】
ROC_Curve.png
やっぱりggplot2はキレイでいいですね!

◆参考文献
・マシンラーニング (Rで学ぶデータサイエンス 6)

マシンラーニング (Rで学ぶデータサイエンス 6)

マシンラーニング (Rで学ぶデータサイエンス 6)



ROCやAUCの基本的な考え方はこの本の1章にあります。


nice!(0)  コメント(0)  トラックバック(1) 

nice! 0

コメント 0

コメントを書く

お名前:
URL:
コメント:
画像認証:
下の画像に表示されている文字を入力してください。

トラックバック 1

この記事のトラックバックURL:
メッセージを送る

この広告は前回の更新から一定期間経過したブログに表示されています。更新すると自動で解除されます。

×

この広告は1年以上新しい記事の更新がないブログに表示されております。