研究でROC解析を行う必要があり、Rでどうやったらできるのか調べてみました。
そうしたところ、ROCRというパッケージが公開されており、比較的簡単にROC解析を行い、グラフを作成できることがわかりました。
- ROCRパッケージのインストール（Ubuntu)
既にRはインストールされているとします。Ubuntuの場合、ROCRパッケージはapt経由で簡単に入手できます。
$ sudo apt-get install r-cran-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はRを立ち上げた後に、library(ROCR)で起動できます。
$ R > library(ROCR)
先ほどのroc_data.txtをrocdataという変数に読み込みます。変数名は何でもいいのですが、ここではそうします。
read.tableという関数で表を読み込めるので、それを使います。
rocdata <- read.table("roc_data.txt")
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です。ここでは、感度、すなわち真陽性率 (TP/(TP+FN)で定義)と、1-特異度、すなわち偽陽性率(FP/(FP+TN)で定義)を求めます。以下のようにタイプします。
perf <- performance(pred, "tpr", "fpr")
ここでtprはtrue positive rateを、fprはfalse positive rateを意味します。
それでは、ROC曲線を描きます。非常に簡単です。
plot(perf)
そうすると、次のようなグラフが現れるはずです。
グラフをPNG形式で保存するには、次のように行うことで、roc-curve.pngという名前でワーキングディレクトリに保存されます。
png("roc-curve.png")
plot(perf)
dev.off()
先ほどの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という名前で保存されます。
比較的簡単に求められるので便利です。