研究でROC解析を行う必要があり、Rでどうやったらできるのか調べてみました。
そうしたところ、ROCRというパッケージが公開されており、比較的簡単にROC解析を行い、グラフを作成できることがわかりました。
$ 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という名前で保存されます。
比較的簡単に求められるので便利です。