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

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

• ROCRパッケージのインストール（Ubuntu)

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

1	$ sudo apt-get install r-cran-rocr

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

• ROCRを使うための準備

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

01 02 03 04 05 06 07 08 09 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40

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)で起動できます。

1 2	$ R > library(ROCR)

• データの読み込み

先ほどのroc_data.txtをrocdataという変数に読み込みます。変数名は何でもいいのですが、ここではそうします。
read.tableという関数で表を読み込めるので、それを使います。

1	rocdata <- read.table("roc_data.txt")

• ROCR ステップ1: prediction

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

1	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)で定義)を求めます。以下のようにタイプします。

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

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

• ROCR ステップ3: グラフの描画

それでは、ROC曲線を描きます。非常に簡単です。

1	plot(perf)

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



• グラフの保存

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

1 2 3

png("roc-curve.png") plot(perf) dev.off()

• AUCの算出

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

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

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

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

1 2	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

1 2 3 4 5 6 7

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を表示させると、以下のように表示されます。

01 02 03 04 05 06 07 08 09 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43

> 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がもっとも高いところを見つけるには、

1	max(table$Accuracy)

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

1 2	> max(table$Accuracy) [1] 0.775

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

• おまけ：感度特異度曲線

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

1 2 3 4

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

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

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