読者です 読者をやめる 読者になる 読者になる

東京に棲む日々

データ分析、統計、ITを勉強中。未だ世に出ず。

モデル評価基準 - ROCに関して - R{ROCR}

モデル評価をROCを用いて行うと仮定した場合の、考察とメモ。

RのROCRパッケージを使用。


混同行列(Confusion Matrix)の復習。

f:id:High_School_Student:20150524081619j:plain

 

こんなデータがあったとする。

f:id:High_School_Student:20150524081841j:plain

pred1、pred2はそれぞれ、モデル1、モデル2による予測値。
observedは実測値。n=19で、Prob(observed=1)=5/19=0.263。


後で、比較を行いやすくするためモデル1、2は同水準の予測値を出力したと仮定してこのデータは作成している。

 


Rへデータの読み込み。
d_roc1 <- read.delim("clipboard", header=TRUE)
d_roc1

   pred1 pred2 observed
1   0.95  0.90        1
2   0.90  0.80        1
3   0.80  0.75        1
4   0.60  0.70        1
5   0.50  0.60        1
6   0.85  0.95        0
7   0.75  0.85        0
8   0.70  0.65        0
9   0.65  0.55        0
10  0.55  0.50        0
11  0.45  0.45        0
12  0.40  0.40        0
13  0.35  0.35        0
14  0.30  0.30        0
15  0.25  0.25        0
16  0.20  0.20        0
17  0.15  0.15        0
18  0.10  0.10        0
19  0.05  0.05        0

 

observation=1の割合。

base_ratio <- sum(d_roc1$observed==1)/length(d_roc1$observed)
0.2631579となる。

 

 

ROCRパッケージのロード。
library(ROCR)

 

pred1に関して、ROC分析を行う。

prediction(予測値,観測値)関数で、予測値(pred1)の大きい順に並び替え、TP,FP,FN,TNを計算。
pred_p1 <- prediction(predictions=d_roc1$pred1, labels=d_roc1$observed)

An object of class "prediction"
Slot "predictions":
1
 [1] 0.95 0.90 0.80 0.60 0.50 0.85 0.75 0.70 0.65 0.55 0.45 0.40 0.35 0.30 0.25 0.20
[17] 0.15 0.10 0.05


Slot "labels":
1
 [1] 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0
Levels: 0 < 1


Slot "cutoffs":
1
 [1]  Inf 0.95 0.90 0.85 0.80 0.75 0.70 0.65 0.60 0.55 0.50 0.45 0.40 0.35 0.30 0.25
[17] 0.20 0.15 0.10 0.05


Slot "fp":
1
 [1]  0  0  0  1  1  2  3  4  4  5  5  6  7  8  9 10 11 12 13 14


Slot "tp":
1
 [1] 0 1 2 2 3 3 3 3 4 4 5 5 5 5 5 5 5 5 5 5


Slot "tn":
1
 [1] 14 14 14 13 13 12 11 10 10  9  9  8  7  6  5  4  3  2  1  0


Slot "fn":
1
 [1] 5 4 3 3 2 2 2 2 1 1 0 0 0 0 0 0 0 0 0 0


Slot "n.pos":
1
[1] 5


Slot "n.neg":
1
[1] 14


Slot "n.pos.pred":
1
 [1]  0  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19


Slot "n.neg.pred":
1
 [1] 19 18 17 16 15 14 13 12 11 10  9  8  7  6  5  4  3  2  1  0

 

performance(predictionオブジェクト,測定基準1,測定基準2)で 感度(真陽性率)、1-特異度(偽陽性率)を計算。
perf_p1 <- performance(pred_p1, measure="tpr", x.measure="fpr") 

# tpr: True Positive Rate, fpr: False Positive Rate

An object of class "performance"
Slot "x.name":
[1] "False positive rate"

Slot "y.name":
[1] "True positive rate"

Slot "alpha.name":
[1] "Cutoff"

Slot "x.values":
1
 [1] 0.00000000 0.00000000 0.00000000 0.07142857 0.07142857 0.14285714 0.21428571
 [8] 0.28571429 0.28571429 0.35714286 0.35714286 0.42857143 0.50000000 0.57142857
[15] 0.64285714 0.71428571 0.78571429 0.85714286 0.92857143 1.00000000


Slot "y.values":
1
 [1] 0.0 0.2 0.4 0.4 0.6 0.6 0.6 0.6 0.8 0.8 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0


Slot "alpha.values":
1
 [1]  Inf 0.95 0.90 0.85 0.80 0.75 0.70 0.65 0.60 0.55 0.50 0.45 0.40 0.35 0.30 0.25
[17] 0.20 0.15 0.10 0.05

y.valuesがTPRで、x.valuesがFPRとなる。

 


ROC Tableの作成。MissClassRは後分類率。
ROCtable_p1 <- data.frame(Cutoff=unlist(pred_p1@cutoffs),
   TP=unlist(pred_p1@tp), FP=unlist(pred_p1@fp),
   FN=unlist(pred_p1@fn), TN=unlist(pred_p1@tn),
   TPR=unlist(perf_p1@y.values), FPR=unlist(perf_p1@x.values),
   diff_TPR_FPR=unlist(perf_p1@y.values)-unlist(perf_p1@x.values),
   MissClasR=(unlist(pred_p1@fp)+unlist(pred_p1@fn))/(unlist(pred_p1@tp)+unlist(pred_p1@fp)+unlist(pred_p1@fn)+unlist(pred_p1@tn))
)

   Cutoff TP FP FN TN TPR        FPR diff_TPR_FPR MissClasR
1     Inf  0  0  5 14 0.0 0.00000000   0.00000000 0.2631579
2    0.95  1  0  4 14 0.2 0.00000000   0.20000000 0.2105263
3    0.90  2  0  3 14 0.4 0.00000000   0.40000000 0.1578947
4    0.85  2  1  3 13 0.4 0.07142857   0.32857143 0.2105263
5    0.80  3  1  2 13 0.6 0.07142857   0.52857143 0.1578947
6    0.75  3  2  2 12 0.6 0.14285714   0.45714286 0.2105263
7    0.70  3  3  2 11 0.6 0.21428571   0.38571429 0.2631579
8    0.65  3  4  2 10 0.6 0.28571429   0.31428571 0.3157895
9    0.60  4  4  1 10 0.8 0.28571429   0.51428571 0.2631579
10   0.55  4  5  1  9 0.8 0.35714286   0.44285714 0.3157895
11   0.50  5  5  0  9 1.0 0.35714286   0.64285714 0.2631579
12   0.45  5  6  0  8 1.0 0.42857143   0.57142857 0.3157895
13   0.40  5  7  0  7 1.0 0.50000000   0.50000000 0.3684211
14   0.35  5  8  0  6 1.0 0.57142857   0.42857143 0.4210526
15   0.30  5  9  0  5 1.0 0.64285714   0.35714286 0.4736842
16   0.25  5 10  0  4 1.0 0.71428571   0.28571429 0.5263158
17   0.20  5 11  0  3 1.0 0.78571429   0.21428571 0.5789474
18   0.15  5 12  0  2 1.0 0.85714286   0.14285714 0.6315789
19   0.10  5 13  0  1 1.0 0.92857143   0.07142857 0.6842105
20   0.05  5 14  0  0 1.0 1.00000000   0.00000000 0.7368421

 
例えば、TPRは以下でも計算可能。
ROCtable_p1$TP / (ROCtable_p1$TP+ROCtable_p1$FN)

 [1] 0.0 0.2 0.4 0.4 0.6 0.6 0.6 0.6 0.8 0.8 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0

 
TPR-FPRの最大値。
max(ROCtable_p1$diff_TPR_FPR)
0.6428571となる。

 

TPR-FPRの最大値とは、真陽性率を高めつつ偽陽性率を小さく押さえる最適な判断基準として、0/1の予測を分ける閾値(Cutt-off)として良く使われる。


この行を出力してみる。
for(i in 1:nrow(ROCtable_p1)){
     mx <- max(ROCtable_p1$diff_TPR_FPR) # max value of TPR-FPR
     if(ROCtable_p1$diff_TPR_FPR[i]==mx){ print(ROCtable_p1[i,]) }
}

   Cutoff TP FP FN TN TPR       FPR diff_TPR_FPR MissClasR
11    0.5  5  5  0  9   1 0.3571429    0.6428571 0.2631579

 11行目が最大値となる。


上で、base_ratio=0.2631579と計算した。これは、すべて0と予測したときの誤分類率となる。
誤分類率がこの値より低い行を出力してみる。
for(i in 1:nrow(ROCtable_p1)){
     if(ROCtable_p1$MissClasR[i]<=base_ratio){ print(ROCtable_p1[i,]) }
}

  Cutoff TP FP FN TN TPR FPR diff_TPR_FPR MissClasR
1    Inf  0  0  5 14   0   0            0 0.2631579
  Cutoff TP FP FN TN TPR FPR diff_TPR_FPR MissClasR
2   0.95  1  0  4 14 0.2   0          0.2 0.2105263
  Cutoff TP FP FN TN TPR FPR diff_TPR_FPR MissClasR
3    0.9  2  0  3 14 0.4   0          0.4 0.1578947
  Cutoff TP FP FN TN TPR        FPR diff_TPR_FPR MissClasR
4   0.85  2  1  3 13 0.4 0.07142857    0.3285714 0.2105263
  Cutoff TP FP FN TN TPR        FPR diff_TPR_FPR MissClasR
5    0.8  3  1  2 13 0.6 0.07142857    0.5285714 0.1578947
  Cutoff TP FP FN TN TPR       FPR diff_TPR_FPR MissClasR
6   0.75  3  2  2 12 0.6 0.1428571    0.4571429 0.2105263
  Cutoff TP FP FN TN TPR       FPR diff_TPR_FPR MissClasR
7    0.7  3  3  2 11 0.6 0.2142857    0.3857143 0.2631579
  Cutoff TP FP FN TN TPR       FPR diff_TPR_FPR MissClasR
9    0.6  4  4  1 10 0.8 0.2857143    0.5142857 0.2631579
   Cutoff TP FP FN TN TPR       FPR diff_TPR_FPR MissClasR
11    0.5  5  5  0  9   1 0.3571429    0.6428571 0.2631579

5行目が最小値。
min(ROCtable_p1$MissClasR)
0.1578947となる。

 

 

AUC(ROCの下の面積)の計算。モデルの精度を数値変換したもので、良く使われるモデル評価値の一つ。
auc.temp <- performance(pred_p1,"auc")
auc_p1 <- unlist(auc.temp@y.values)
0.8571429となる。

 


ROCをプロットする。
plot(perf_p1)
45度線も重ねる。
abline(a=0, b=1, col="yellow")
グリッド線も重ねる。
abline( h=(1:10)/10, v=(1:10)/10, lty=3 )

各行でのTPR、FPRの点も重ねておく。
par(new=TRUE)
plot(ROCtable_p1$FPR, ROCtable_p1$TPR, xlab="", ylab="")

f:id:High_School_Student:20150524083834j:plain

TPRを縦軸に、FPRを横軸にプロットしたものがROCとなる。

要するに、縦軸が「実際に1のものを1と予測できた割合」、横軸が「実際は0のものを1と予測してしまった割合」

 


ある予測モデルのROCとは、その予測モデルによる予測値をさまざまな箇所で切った場合、どの程度0/1を正確に判断できるかを把握するためのグラフとなる。

上のROC TableのTPRとFPRを上から順にプロットして、点を繋いだたもの。

各点において、混同行列を作成できる。

誤分類率が最小になる5行目と、TPR-FPRが最大となる11行目の混同行列をROCプロットにメモしておくとこんな感じ。

f:id:High_School_Student:20150524084028j:plain

閾値をどこにするか、すなわち、どこで切って0/1を判断するかは、分析の目的による。

 


pred2に関してもROC分析を行い、この点に関して考えてみる。


上の各手続きを関数化しておく。

 

# ROCテーブルを返す関数
fn_ROC_Table <- function(pred, obs){

     pred1 <- prediction(predictions=pred, labels=obs)
     perf1 <- performance(pred1, measure="tpr", x.measure="fpr")

     # ROC Tableの作成
     # MissClassRは後分類率
     ROCtable <- data.frame(
     Cutoff=unlist(pred1@cutoffs),
     TP=unlist(pred1@tp), FP=unlist(pred1@fp),
     FN=unlist(pred1@fn), TN=unlist(pred1@tn),
     TPR=unlist(perf1@y.values), FPR=unlist(perf1@x.values),
     diff_TPR_FPR=unlist(perf1@y.values)-unlist(perf1@x.values),
     MissClasR=(unlist(pred1@fp)+unlist(pred1@fn))/
(unlist(pred1@tp)+unlist(pred1@fp)+unlist(pred1@fn)+unlist(pred1@tn))
     )

     return(ROCtable)
}

 

fn_ROC_Table(d_roc1$pred1, d_roc1$observed)でROCtable_p1と同じとなる。

 

pred2に関して。
fn_ROC_Table(d_roc1$pred2, d_roc1$observed)

   Cutoff TP FP FN TN TPR        FPR diff_TPR_FPR MissClasR
1     Inf  0  0  5 14 0.0 0.00000000   0.00000000 0.2631579
2    0.95  0  1  5 13 0.0 0.07142857  -0.07142857 0.3157895
3    0.90  1  1  4 13 0.2 0.07142857   0.12857143 0.2631579
4    0.85  1  2  4 12 0.2 0.14285714   0.05714286 0.3157895
5    0.80  2  2  3 12 0.4 0.14285714   0.25714286 0.2631579
6    0.75  3  2  2 12 0.6 0.14285714   0.45714286 0.2105263
7    0.70  4  2  1 12 0.8 0.14285714   0.65714286 0.1578947
8    0.65  4  3  1 11 0.8 0.21428571   0.58571429 0.2105263
9    0.60  5  3  0 11 1.0 0.21428571   0.78571429 0.1578947
10   0.55  5  4  0 10 1.0 0.28571429   0.71428571 0.2105263
11   0.50  5  5  0  9 1.0 0.35714286   0.64285714 0.2631579
12   0.45  5  6  0  8 1.0 0.42857143   0.57142857 0.3157895
13   0.40  5  7  0  7 1.0 0.50000000   0.50000000 0.3684211
14   0.35  5  8  0  6 1.0 0.57142857   0.42857143 0.4210526
15   0.30  5  9  0  5 1.0 0.64285714   0.35714286 0.4736842
16   0.25  5 10  0  4 1.0 0.71428571   0.28571429 0.5263158
17   0.20  5 11  0  3 1.0 0.78571429   0.21428571 0.5789474
18   0.15  5 12  0  2 1.0 0.85714286   0.14285714 0.6315789
19   0.10  5 13  0  1 1.0 0.92857143   0.07142857 0.6842105
20   0.05  5 14  0  0 1.0 1.00000000   0.00000000 0.7368421 

 

 

# AUCを返す関数
fn_AUC <- function(pred, obs){

     pred1 <- prediction(predictions=pred, labels=obs)
     auc.temp <- performance(pred1, "auc")
     auc <- unlist(auc.temp@y.values)
     return(auc)
}

 

fn_AUC(d_roc1$pred1, d_roc1$observed)でauc_p1と同じとなる。

 

pred2に関して。
fn_AUC(d_roc1$pred2, d_roc1$observed)
0.8571429と、AUCはpred1と同じである。

 


# ROC曲線を描く関数
fn_plot_ROC <- function(pred, obs, line45=TRUE, gridL=TRUE, color="black"){
     pred1 <- prediction(predictions=pred, labels=obs)
     perf1 <- performance(pred1, measure="tpr", x.measure="fpr")

     plot(perf1, col=color)

     # 45 degree line
     if(line45==TRUE){ abline(a=0, b=1, col="yellow") }

     # grid
     if(gridL==TRUE){ abline( h=(1:10)/10, v=(1:10)/10, lty=3 ) }
}


fn_plot_ROC(d_roc1$pred1, d_roc1$observed)でplot(perf_p1)と同じになる。

 

pred1とpred2の曲線を重ねてみる。
fn_plot_ROC(d_roc1$pred1, d_roc1$observed)
par(new=TRUE)
fn_plot_ROC(d_roc1$pred2, d_roc1$observed, color="blue")

f:id:High_School_Student:20150524084753j:plain

がpred1で、がpred2である。

AUCの値は両方とも同じであるが、曲線の形状は異なる。


pred1はROCの下側が外に膨らんでおり、pred2はROCの上側が外に膨らんでいる。

 

このデータを別の書き方をしてみる。

データは、pred1と2の予測値の水準を同じとした擬似的なものなので、この予測値の大きい順に並び替えてみる。

f:id:High_School_Student:20150524084908j:plain

predがpred1と2の値。obs of pred1がpred1の予測値に対応したobserved。obs of pred2がpred2に対応したobservedである。


少数のみ1と予測したくて、その1と予測したオブザベーションが実際に1である場合を多くしたい場合は、pred1型のROCとなるモデルを採用するのが良い。
沢山1と予測してしまってはまずい場合。


なるべく1である場合をもらざず、がさっと取ってきたいときはは、pred2型のROCとなるモデルを採用するのが良い。
1を逃してしまうとまずい場合。

 

 

 

 

では、AUCがどの程度の精度を持つか、予測値をランダム出力し、シミュレーションしてみる。一様分布により、0から1の間の予測値を出力する。


observation(1/0)の指定。この割合をいろいろと変えてシミュレーションしてみる。
obs1 <- c(rep(1,10), rep(0,15))

 [1] 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

length1 <- length(obs1)

# AUCを格納するベクトル
auc1 <- c()

# シミュレーション回数(100回)
n_sim <- 100

for(i in 1:n_sim){

     pred1 <- runif(length1)      # 一様分布

     # plot ROC
     fn_plot_ROC(pred1, obs1)
     if(i>=n_sim){ break }
     par(new=TRUE)

     # AUC
     auc1 <- append(auc1, fn_AUC(pred1, obs1))

     i <- i+1
}

f:id:High_School_Student:20150524085525j:plain

 

シミュレーションにより計算したAUC(100個)のQuantileを求める。
quantile(auc1)

       0%       25%       50%       75%      100% 
0.1933333 0.4633333 0.5266667 0.5900000 0.7866667 

 

 

データ数が25くらいだと、全く意味のないモデルでもAUCがかなり高くなることがある。

 

数を増やしてみる。合計10,000行のデータ。
obs1 <- c(rep(1,1000), rep(0,9000))

として、上のシミュレーションコードを実行。

f:id:High_School_Student:20150524085815j:plain

       0%       25%       50%       75%      100% 
0.4732207 0.4911747 0.4994109 0.5058587 0.5229493 

これくらいデータがあると、ランダム予測AUCは0.5という理論が成り立つ。

 


ちょっと1/0のバランスが悪い場合。
obs1 <- c(rep(1,5), rep(0,200))
5/205=0.024。

f:id:High_School_Student:20150524090214j:plain

    0%    25%    50%    75%   100% 
0.1930 0.4345 0.5270 0.6275 0.8130 

1が5個しかないので、縦軸のTPRは5段分割しかされない。

データは合計で205あるが、AUCも大きくばらつく。1/0を100づつでシミュレーションすると、こんなに大きくAUCがばらつくことはない。