著者別アーカイブ： kiyotaka

研究で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)

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



• グラフの保存

グラフを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という名前で保存されます。

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