モデル評価基準 - ROCに関して - R{ROCR}
モデル評価をROCを用いて行うと仮定した場合の、考察とメモ。
RのROCRパッケージを使用。
混同行列（Confusion Matrix）の復習。
こんなデータがあったとする。
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="")
TPRを縦軸に、FPRを横軸にプロットしたものがROCとなる。
要するに、縦軸が「実際に1のものを1と予測できた割合」、横軸が「実際は0のものを1と予測してしまった割合」。
ある予測モデルのROCとは、その予測モデルによる予測値をさまざまな箇所で切った場合、どの程度0/1を正確に判断できるかを把握するためのグラフとなる。
上のROC TableのTPRとFPRを上から順にプロットして、点を繋いだたもの。
各点において、混同行列を作成できる。
誤分類率が最小になる5行目と、TPR-FPRが最大となる11行目の混同行列をROCプロットにメモしておくとこんな感じ。
閾値をどこにするか、すなわち、どこで切って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")
黒がpred1で、青がpred2である。
AUCの値は両方とも同じであるが、曲線の形状は異なる。
pred1はROCの下側が外に膨らんでおり、pred2はROCの上側が外に膨らんでいる。
このデータを別の書き方をしてみる。
データは、pred1と2の予測値の水準を同じとした擬似的なものなので、この予測値の大きい順に並び替えてみる。
predがpred1と2の値。obs of pred1がpred1の予測値に対応したobserved。obs of pred2がpred2に対応したobservedである。
少数のみ1と予測したくて、その1と予測したオブザベーションが実際に1である場合を多くしたい場合は、pred1型のROCとなるモデルを採用するのが良い。
沢山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
}
シミュレーションにより計算した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))
として、上のシミュレーションコードを実行。
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。
0% 25% 50% 75% 100% 0.1930 0.4345 0.5270 0.6275 0.8130
1が5個しかないので、縦軸のTPRは５段分割しかされない。
データは合計で205あるが、AUCも大きくばらつく。1/0を100づつでシミュレーションすると、こんなに大きくAUCがばらつくことはない。
