3.8 判別の評価

モデルの評価基準として逸脱度を紹介したが,実際目的変数に対して推定が上手くいっているのかを確認することも非常に重要である.

ロジスティック回帰モデルの当てはめによって\(i\)番目のデータに対して\(\hat y_i = 1\)となる確率の推定値\(\hat \pi_i\)が得られ,実際Rでは推定結果のオブジェクトに対してfitted関数を当てはめることで推定値が得られる.

fitted(result) %>% head
##         1         2         3         4         5         6 
## 0.9997697 0.9975913 0.9999961 0.9806641 0.9995125 0.6299903

これは1番目のデータから順番にこれらの確率で\(\hat y_i\)は1になるだろうという推定でもある.

次に分析者に要求されるのは,これらの得られた確率の推定値\(\hat \pi_i\)の値を用いてどのように\(\hat y_i\)を計算するかであり.これらは最終的には分析者が決定することになりる.

例えば

  • \(\hat \pi_i \geq 0.5\)のとき\(\hat y_i = 1\),それ以外は\(\hat y_i = 0\)

などが考えられる.この\(0.5\)という数値は\(0 \sim 1\)を自由に動かしても良いわけなので,ここでは便宜上\(T\)としよう.もし\(T = 0\)とすると全てを1と予測し,\(T=1\)とすると全てを0と予測することになる.

さて問題はどのように\(T\)を決定するかになるが,これは分析状況によって異なるだろう. 推定および予測には必ず誤判別の可能性がある.誤判別とは本来の値とは異なる予測をすることで,今回であれば,

  • 本当は良性腫瘍なのに悪性と判別
  • 本当は悪性のはずなのに良性と判別

することを意味する.\(T\)の値が低ければ悪性腫瘍を見逃すことは減るが,一方で本当は良性の人を悪性と判断する可能性が高まる.逆に\(T\)の値を高くすると本当に悪性の人だけを悪性と判別できるようになるはずだが,本当は悪性の人を良性と判断してしまうケースが増えるだろう.

これは一般的にどちらが良いということはなく,分析目的によってどちらを許容するかは変わる.例えば,乳がんは生命に直結する事態であり,実は良性だが悪性と判断してしまうことより,実は悪性だが良性だと判断してしまう方が深刻な事態に繋がりかねない.一方,大した疾患ではなく軽症であれば治るようなものであればコスト削減の意味で逆の判断となる可能性もある.

それでは,これらがトレードオフであるという状況を確認してみよう.

# T = 0.5の場合
data.frame(
  Diagnosis = df$Diagnosis,
  hat_pi = fitted(result),
  color_value = ifelse(df$Diagnosis == 1, "orange", "skyblue")
) %>%
  arrange(hat_pi) %>% 
ggplot() + 
  geom_jitter(
    height = 0.08, width = 0,
    aes(x=hat_pi, y=Diagnosis, color = color_value)) +
  geom_vline(xintercept = c(0.5)) + 
  ggtitle("T = 0.5")

# T = 0.7の場合
data.frame(
  Diagnosis = df$Diagnosis,
  hat_pi = fitted(result),
  color_value = ifelse(df$Diagnosis == 1, "orange", "skyblue")
) %>%
  arrange(hat_pi) %>% 
ggplot() + 
  geom_jitter(
    height = 0.08, width = 0,
    aes(x=hat_pi, y=Diagnosis, color = color_value)) +
  geom_vline(xintercept = c(0.7)) + 
  ggtitle("T = 0.7")

このように考えると判別結果を考察するときに次のような4つの指標を見るべきだろうということがわかる.

実際の値\推定値 0 1
0 FN:真陰性(true negative) FP:偽陽性(false positive)
1 TN:偽陰性(false negative) TP:真陽性(true positive)

この表に基づいて数を集計したものを混合行列と呼ぶ.いま便宜上\(T=0.5\)として集計した混合行列を作成してみよう.

# 実際の値
y <- df$Diagnosis
hat_y <- ifelse(fitted(result) >= 0.5, 1, 0)
table(y,hat_y)
##    hat_y
## y     0   1
##   0 346  11
##   1  19 193

これらの指標を用いて得られる次の二つの比率もよく分析に用いられる.

  • 正解率(ACC: accuracy rate)
  • 誤判別立(ERR: error rate)

\[ \begin{align} {\rm ACC} &= \frac{\rm TP + TN}{\rm TP + FN + FP + TN}, \\ {\rm ERR} &= 1 - {\rm ACC} \end{align} \]

  • 真陽性率(TRP, true positive rate):陽性と判断したものの中で本当に陽性であるものの割合
  • 偽陰性率(FPR, false positive rate):陰性と判断したものの中で本当に陰性であるものの割合

混合行列とそれに関する指標の計算に関しては便利なパッケージも提供されている. caretというパッケージのconfusionMatrix関数を紹介しよう.

install.packages("caret")を実行しcaretパッケージをインストールしておこう. また,この関数の引数に与えるベクトルはファクター型にしておく必要があることに注意されたい.

library(caret)
confusionMatrix(data=as.factor(y), reference=as.factor(hat_y), positive = "1")
## Confusion Matrix and Statistics
## 
##           Reference
## Prediction   0   1
##          0 346  11
##          1  19 193
##                                           
##                Accuracy : 0.9473          
##                  95% CI : (0.9256, 0.9641)
##     No Information Rate : 0.6415          
##     P-Value [Acc > NIR] : <2e-16          
##                                           
##                   Kappa : 0.8864          
##                                           
##  Mcnemar's Test P-Value : 0.2012          
##                                           
##             Sensitivity : 0.9461          
##             Specificity : 0.9479          
##          Pos Pred Value : 0.9104          
##          Neg Pred Value : 0.9692          
##              Prevalence : 0.3585          
##          Detection Rate : 0.3392          
##    Detection Prevalence : 0.3726          
##       Balanced Accuracy : 0.9470          
##                                           
##        'Positive' Class : 1               
## 

3.8.1 ROC曲線

\(T\)の値を\(0 \sim 1\)まで少しずつ変化させ,TPRとFPRの値をプロットしたものをROC曲線(receiver operating characteristic curve)と言う.またROC曲線における曲線下側面積をAUC(area under the curve)と呼ぶ.

ROCの計算と作図およびAUCの計算は次のような手順で計算することができる.

y <- df$Diagnosis # 実際の値
hat_pi <- fitted(result) # piの推定値
t_values <- seq(0,1,0.001)

# TPR,FPRを格納するベクトルを生成
tp_count <- c()
fp_count <- c()


for(i in 1:length(t_values)){
  # Tの値はhat_piの値を昇順に並べて小さい方から採用する
  t_value <- t_values[i] # Tの値を更新
  hat_y <- ifelse(hat_pi >= t_value, 1 , 0) # yの推定値を更新
  cm <- table(y,hat_y) # 混合行列
  
  # False Positiveの集計 : 実際の値が0で推定値が1のもの
  fp_count[i] <- sum({y == 0} & {hat_y == 1})
  
  # True Positiveの集計 : 実際の値が1で推定値が1のもの
  tp_count[i] <- sum({y == 1} & {hat_y == 1})
}
# TPR, FPRの計算
true_count <- sum(df$Diagnosis == 1) # y=1のカウント
false_count <- sum(df$Diagnosis == 0) # y=0のカウント

tpr <- tp_count / true_count
fpr <- fp_count / false_count
plot(fpr, tpr, type="l")

# AUCの計算
1 - sum( tpr * c(fpr[1], diff(fpr)))
## [1] 0.9866815

また,ROC,AUCについてはpROCという便利なパッケージも提供されており分析の際はこちらを使うのが良いだろう.

install.packages("pROC")を実行しておこう.

library(pROC)
plot(roc(y, hat_pi))
## Setting levels: control = 0, case = 1
## Setting direction: controls < cases

roc(y, hat_pi)
## Setting levels: control = 0, case = 1
## Setting direction: controls < cases
## 
## Call:
## roc.default(response = y, predictor = hat_pi)
## 
## Data: hat_pi in 357 controls (y 0) < 212 cases (y 1).
## Area under the curve: 0.9867

モデルの評価についてはAUCの値を利用することもある. AUCの最大面積は1であり.すなわちAUCの値は1に近いほど良いと考えられる. AUCの意味でモデル選択をする場合は,得られたモデルの中で最もAUCが大きいものを選ぶことになる.