RでのROC解析:ROCRパッケージを使ったROC曲線とAUCの求め方

研究でROC解析を行う必要があり、Rでどうやったらできるのか調べてみました。
そうしたところ、ROCRというパッケージが公開されており、比較的簡単にROC解析を行い、グラフを作成できることがわかりました。

  • ROCRパッケージのインストール(Ubuntu)
  • 既にRはインストールされているとします。Ubuntuの場合、ROCRパッケージはapt経由で簡単に入手できます。

    $ sudo apt-get install r-cran-rocr
    

    これでインストール完了です。

  • ROCRを使うための準備
  • ROC解析に必要なものは、何らかの指標と、それが属するグループの一覧です。具体例を挙げると、以下のようになります。
    第1列に指標、第2列に属するグループ(0か1)が記載されています。
    これをroc_data.txtとという名前で保存することとします。保存したディレクトリをRのワーキングディレクトリとします。

    0.9706	1
    0.9572	1
    0.4854	1
    0.8003	1
    0.1419	1
    0.4218	1
    0.9157	1
    0.7922	1
    0.9595	1
    0.6557	1
    0.0357	1
    0.8491	1
    0.934	1
    0.6787	1
    0.7577	1
    0.7431	1
    0.3922	1
    0.6555	1
    0.1712	1
    0.706	1
    0.4797	0
    0.4551	0
    0.0374	0
    0.081	0
    0.2984	0
    0.7597	0
    0.1404	0
    0.3853	0
    0.0238	0
    0.5513	0
    0.0551	0
    0.306	0
    0.4991	0
    0.6909	0
    0.7593	0
    0.3472	0
    0.0614	0
    0.0507	0
    0.0575	0
    0.6407	0
    
  • ROCRの起動
  • ROCRはRを立ち上げた後に、library(ROCR)で起動できます。

    $ R
    > library(ROCR)
    
  • データの読み込み
  • 先ほどのroc_data.txtをrocdataという変数に読み込みます。変数名は何でもいいのですが、ここではそうします。
    read.tableという関数で表を読み込めるので、それを使います。

    rocdata <- read.table("roc_data.txt")
    
  • ROCR ステップ1: prediction
  • ROCRのステップ1はpredictionで、値と属するグループを指定します。
    今、rocdataは20行2列の行列になっています。1列目はrocdata[,1]で、2列目はrocdata[,2]であらわすことができますので、以下のように記載します。

    pred <- prediction(rocdata[,1], rocdata[,2])
    

    ここで、変数predはpredictionの頭文字4文字です。もちろん、別の名前でもかまいません。
    私はよくカンマを忘れるので、カンマも忘れないようにしましょう。

    ここで、何が行われているかというと、データを大きい順にソートし、真陽性(TP)、偽陽性(FP)、偽陰性(FN)、真陰性(TN) の数を算出します。

  • ROCR ステップ2: performance
  • ROCRのステップ2はperformanceです。ここでは、感度、すなわち真陽性率 (TP/(TP+FN)で定義)と、1-特異度、すなわち偽陽性率(FP/(FP+TN)で定義)を求めます。以下のようにタイプします。

    perf <- performance(pred, "tpr", "fpr")
    

    ここでtprはtrue positive rateを、fprはfalse positive rateを意味します。

  • ROCR ステップ3: グラフの描画
  • それでは、ROC曲線を描きます。非常に簡単です。

    plot(perf)
    

    そうすると、次のようなグラフが現れるはずです。

    roc-curve

  • グラフの保存
  • グラフをPNG形式で保存するには、次のように行うことで、roc-curve.pngという名前でワーキングディレクトリに保存されます。

    png("roc-curve.png")
    plot(perf)
    dev.off()
    
  • AUCの算出
  • 先ほどのperformanceの際に”auc”と指定するとAUCも計算されます。ただ、1クッション入れる必要があります。具体的な方法はこちらのサイトに記載されていましたが、それを転載します。

    auc.tmp <- performance(pred,"auc")
    auc <- as.numeric(auc.tmp@y.values)
    

    performance(pred,”auc”)の結果をauc.tmpという変数に代入し、
    auc.tmpの中からy.valuesの値を取り出して、その値を変数aucに代入します。

    最後に変数aucを表示させてみます。

    auc
    [1] 0.8
    

    これでAUCが0.8だということがわかります。

  • 正診率の算出
  • ここからさらに一歩踏み込んで、正診率を求めたいと思います。

    正診率、感度、特異度は以下で定義されます。

    正診率=(TP+TN)/総数
    感度(真陽性率)=TP/(TP+FN)
    特異度=TN/(FP+TN)

    今、総数は、rocdataの行数を求めればよいですから、nrow(rocdata)で求められます

    カットオフ値を少しずつずらした時に、TP, FP, FN, TNは変わっていきますので、その一覧を表に出力しましょう。

    表の列は以下のようにしたいと思います

    Cutoff TP FP FN TN Sensitivity Specificity Accuracy

    table <- data.frame(Cutoff=unlist(pred@cutoffs),
      TP=unlist(pred@tp), FP=unlist(pred@fp),
      FN=unlist(pred@fn), TN=unlist(pred@tn),
      Sensitivity=unlist(pred@tp)/(unlist(pred@tp)+unlist(pred@fn)),
      Specificity=unlist(pred@tn)/(unlist(pred@fp)+unlist(pred@tn)),
      Accuracy=((unlist(pred@tp)+unlist(pred@tn))/nrow(rocdata))
      )
    

    これで、tableを表示させると、以下のように表示されます。

    > table
       Cutoff TP FP FN TN Sensitivity Specificity Accuracy
    1     Inf  0  0 20 20        0.00        1.00    0.500
    2  0.9706  1  0 19 20        0.05        1.00    0.525
    3  0.9595  2  0 18 20        0.10        1.00    0.550
    4  0.9572  3  0 17 20        0.15        1.00    0.575
    5  0.9340  4  0 16 20        0.20        1.00    0.600
    6  0.9157  5  0 15 20        0.25        1.00    0.625
    7  0.8491  6  0 14 20        0.30        1.00    0.650
    8  0.8003  7  0 13 20        0.35        1.00    0.675
    9  0.7922  8  0 12 20        0.40        1.00    0.700
    10 0.7597  8  1 12 19        0.40        0.95    0.675
    11 0.7593  8  2 12 18        0.40        0.90    0.650
    12 0.7577  9  2 11 18        0.45        0.90    0.675
    13 0.7431 10  2 10 18        0.50        0.90    0.700
    14 0.7060 11  2  9 18        0.55        0.90    0.725
    15 0.6909 11  3  9 17        0.55        0.85    0.700
    16 0.6787 12  3  8 17        0.60        0.85    0.725
    17 0.6557 13  3  7 17        0.65        0.85    0.750
    18 0.6555 14  3  6 17        0.70        0.85    0.775
    19 0.6407 14  4  6 16        0.70        0.80    0.750
    20 0.5513 14  5  6 15        0.70        0.75    0.725
    21 0.4991 14  6  6 14        0.70        0.70    0.700
    22 0.4854 15  6  5 14        0.75        0.70    0.725
    23 0.4797 15  7  5 13        0.75        0.65    0.700
    24 0.4551 15  8  5 12        0.75        0.60    0.675
    25 0.4218 16  8  4 12        0.80        0.60    0.700
    26 0.3922 17  8  3 12        0.85        0.60    0.725
    27 0.3853 17  9  3 11        0.85        0.55    0.700
    28 0.3472 17 10  3 10        0.85        0.50    0.675
    29 0.3060 17 11  3  9        0.85        0.45    0.650
    30 0.2984 17 12  3  8        0.85        0.40    0.625
    31 0.1712 18 12  2  8        0.90        0.40    0.650
    32 0.1419 19 12  1  8        0.95        0.40    0.675
    33 0.1404 19 13  1  7        0.95        0.35    0.650
    34 0.0810 19 14  1  6        0.95        0.30    0.625
    35 0.0614 19 15  1  5        0.95        0.25    0.600
    36 0.0575 19 16  1  4        0.95        0.20    0.575
    37 0.0551 19 17  1  3        0.95        0.15    0.550
    38 0.0507 19 18  1  2        0.95        0.10    0.525
    39 0.0374 19 19  1  1        0.95        0.05    0.500
    40 0.0357 20 19  0  1        1.00        0.05    0.525
    41 0.0238 20 20  0  0        1.00        0.00    0.500
    

    Accuracyがもっとも高いところを見つけるには、

    max(table$Accuracy)
    

    とします。そうすると、今は、

    > max(table$Accuracy)
    [1] 0.775
    

    となりますので、該当するところをみると、感度70%、特異度85%、正診率77.5%達成できるということがわかりました。

  • おまけ:感度特異度曲線
  • 感度と特異度の曲線も簡単に書けます。predictionまで行った後に、次のようにします。

    perf <- performance(pred, "sens", "spec")
    png("sens-spec-curve.png")
    plot(perf)
    dev.off()
    

    これで、下図のような感度、特異度の曲線がsens-spec-curve.pngという名前で保存されます。

    sens-spec-curve

比較的簡単に求められるので便利です。