【統計】ロジスティック回帰分析（ホスマー・レメショウ検定、ROC、AUC等）

R、R Studio Python SAS 統計

ロジスティック回帰の評価について。

【目次】

 

前回で理論部分を取り扱った。その続きでモデルの評価方法。

cochineal19.hatenablog.com

 

計算式等

モデル評価には、大きく分けて「Discrimination（判別能力）」と「Calibration（較正）」がある。
「Discrimination（判別能力）」では、モデルが正しくイベントの発生有無を判別できるかの予測性能を評価する。
「Calibration（較正）」では、サブグループ毎に実測値と予測値のズレ具合を評価し、部分的に予測精度の悪さがないかを評価する。

 

 Discrimination（判別能力）

① 真陽性率（感度）、真陰性率（特異度）等

実測と予測の関係を混同行列（Confusion Matrix）に表すと次のとおり。

  実測
(+) (-)
予測 (+) 真陽性
Ture Positive (TP)		 偽陽性
False Positive (FP)
(-) 偽陰性
False Negative (FN)		 真陰性
True Negative (TN)


この混同行列を参考として、各種評価指標を次のように計算する。 

項目 式 類語
真陽性率 (True Positive Rate) TP / (TP + FN) 感度 (Sensitivity)、検出率 (Recall)
真陰性率 (True Negative Rate) TN / (FP + TN) 特異度 (Specificity)
偽陽性率 (False Positive Rate) FP / (FP + TN) 1 - 特異度 (Specificity)
偽陰性率 (False Negative Rate) FN / (TP + FN) 1 - 感度 (Sensitivity)
陽性適中率 (Positive Predictive Value) TP / (TP + FP) 適合率 (Precision)
陰性適中率 (Negative Predictive Value) TN / (FN + TN) ー
偽発見率 (False Discovery Rate) FP / (TP + FP) 1 - 適合率 (Precision)
正診率 (Accuracy) (TP + TN) / N ー

どの指標を使うかは目的に応じる。

・感度（検出率）：実際に陽性だったもののうち、どれだけ陽性と言えたか。
・特異度：実際に陰性だったもののうち、どれだけ陰性と言えたか。
・適合率：陽性と予測したもののうち、本当に陽性だったのはどのくらいか。
・正診度：予測全体のうち、本当に陽性または陰性だったのはどのくらいか。

なお、例えば、めったに起きないイベントの有無を予測したい場合、無と答えていれば大体正解できてしまうため予測性能が高いモデルのように見えてしまう。
このような不均衡データに対応する指標として、マシューズ相関係数（MCC）やF1スコアなどがある。F1スコアは検出率と適合率の調和平均。

\quad MCC=\dfrac{TP\times TN-FP\times FN}{\sqrt{\left( TP+FP\right) \left( TP+FN\right) \left( TN+FP\right) \left( TN+FN\right) }}

\quad F1\ Score=2\times\dfrac{Recall \times Precision}{Recall + Precision}=2\times\dfrac{\{TP / \left(TP + FN\right)\}\times \{TP / \left(TP + FP\right)\} }{ \{TP / \left(TP + FN\right)\} + \{TP / \left(TP + FP\right)\} }

 

② ROC曲線（Receiver Operating Characteristic Curve）、
　 AUC（Area Under the Curve）

ROC曲線は、縦軸に真陽性率（感度）、横軸に偽陽性率（１－特異度）をプロットした曲線。横軸を真陰性率（特異度）とすることもある。
AUCは、ROC曲線下の面積。０～１の範囲を取り、一つの基準として " 0.7 以上で可、0.8 以上で良、0.9 以上で優" とされる。

f:id:cochineal19:20210130225421p:plain

 

Calibration（較正）

① ホスマー・レメショウ検定（Hosmer and Lemeshow goodness of fit (GOF) test）

ホスマー・レメショウ検定は、データセットをk個に分割してサブグループ（リスクグループ）を作成する。このサブグループ毎に "実測度数 O_{k} と期待度数 E_{k} のズレ" を算出し、足し合わせた \chi _{HL}^{2} 統計量を用いてモデルの適合度を評価する手法。
※ 死亡有無を予測するモデルを例にすると、死亡率０～９％、10％台、20％台… とリスクグループを分割することで各リスクグループでの予実の差（ズレ）を評価できる。これにより、特定のリスクグループに予測精度の悪さがないかを評価できる。特に目立ったズレがなければ全体的にモデルが適合していると言える。

\qquad \chi _{HL}^{2}=\sum \dfrac{O_{k}-E_{k}}{N_{k}\widehat{p} _{k}\left( 1-\widehat{p}_{k}\right) }=\sum \dfrac{O_{k}-N_{k}\widehat{p} _{k}}{N_{k}\widehat{p} _{k}\left( 1-\widehat{p}_{k}\right) }

　　※ \widehat{p} は予測値 \widehat{y} の平均値。
　　※ グループ数は明確な基準はない。５～１０くらい？

\chi _{HL}^{2} は自由度 k-2 の \chi ^{2} 分布に近似することができ、p>0.05 の時、帰無仮説を棄却せず、モデルが適合していると評価する。
　帰無仮説H_{0}: p=\widehat{p}（モデルがどのサブグループにも適合している）
　対立仮説H_{1}: p≠\widehat{p}（モデルがいずれかのサブグループに適合していない）

 

② Calibration plot

予測値 \widehat{p}_{k} と実測値 p_{k} をプロットしたグラフ。
上のホスマー・レメショウ検定の分割したデータセットを使って作ることができる。

 

---------
その他にも、McFaddenの疑似決定係数、Brier Score、Log Lossなどの指標がある。

 

計算例

ロジスティック回帰の理論部分で使ったサンプルデータを用いる。
（SASのサンプルデータ作成は 前回記事 ）
※機械学習では訓練データ：バリデーションデータ：テストデータなどに分割するが、ここではこれら手順は行わない（また別のタイミングで）。

ads <- data.frame(
  x1=c(rep(5,20),rep(5,22),rep(10,25),rep(10,30),rep(15,30),rep(15,30))
  ,x2=c(rep(0,20),rep(1,22),rep(0,25),rep(1,30),rep(0,30),rep(1,30))
  ,y=c(rep(0,19),rep(1,1),rep(0,16),rep(1,6),rep(0,16),rep(1,9),rep(0,14),rep(1,16),rep(0,11),rep(1,19),rep(0,2),rep(1,28))
)
fit1 <- glm(y ~ x1 + x2, family=binomial, data=ads)
ads$yhat <- fit1$fitted.values


------------------------------------

Discrimination（判別能力）

ROC曲線を描画し、AUCを算出する。 

Rでの実行：

library(tidyverse)
library(pROC)
(roc1 <- pROC::roc(response=fit1$y, predictor=fit1$fitted, ci=T))
roc1.res1 <- (pROC::coords(roc1, x="best", best.method="youden", ret=c("threshold","tpr","fpr","tnr","fnr","npv","ppv"), transpose=F))
roc1.res2 <- (pROC::coords(roc1, x="best", best.method="closest.topleft", ret=c("threshold","tpr","fpr","tnr","fnr","npv","ppv"), transpose=F))
AUC1 <- round(roc1$ci, 4)

ggroc(list(Model1=roc1), legacy.axes=T, size=1.5, color="darkgreen", alpha=.7) +
  geom_abline(intercept=0, slope=1, linetype="dashed") +
  annotate("text", x=1-roc1.res1$tnr, y=roc1.res1$tpr, hjust=0.5, vjust=0.5, size=10, label="×", color="blue") +
  annotate("text", x=1-roc1.res2$tnr, y=roc1.res2$tpr, hjust=0.5, vjust=0.5, size=10, label="〇", color="red") +
  annotate("text", x=0.75, y=0.1, size=6
           ,label=paste0("Area Under the Curve ( 95%CI ) :\n",AUC1[2]," ( ",AUC1[1],", ",AUC1[3]," )")) +
  labs(title="ROC曲線")

f:id:cochineal19:20210130230434p:plain

SASでの実行：

ods output LackFitChiSq=HOSLEM1;
proc logistic data=ads descending;
  model y(event='1') = x1 x2 / outroc=ROC1 lackfit;
  output out=OutDS;
run;

f:id:cochineal19:20210130230629p:plain

Rのグラフでは最適カットオフ値を載せている。
Youden Indexが青色の×印、closest.topleft が赤色の〇印。
Youden Indexは45度線から一番離れているROC曲線のポイント、closest.topleftは左上の角（AUC=1.0の位置）から一番近いROC曲線のポイントをそれぞれ最適カットオフ値とする。今回のデータでは、両方法とも同じ最適カットオフ値（threshold=0.4497）だった。

今回のデータでは、AUC=0.8085 [ 0.7435, 0.8735] で0.8を超える良好な結果。
また、最適カットオフ値での混同行列（Confusion Matrix）と各評価指標は次のとおり。

  実測
(+) (-)
予測 (+) ６３ ２７
(-) １６ ５１

・感度（検出率）=63/(63+13)=0.7975
・特異度 =51/(27+51)=0.6538
・適合率 =63/(63+27)=0.7000
・正診率 =(63+51)/(63+27+16+51)=0.7261 
・F1スコア=2×0.7975×0.7/(0.7975+0.7)=0.7456
・MCC =(63×51-27×16)/sqrt{(63+27)(63+16)(51+27)(51+16)}=0.4562

 

------------------------------------

Calibration（較正）

ホスマー・レメショウ検定をRとSASで行う。 

Rでの実行：（※グループ数はSASの実行結果に合わせて 6 にしている）

> library(tidyverse)
> library(ResourceSelection)
> (HL1 <- ResourceSelection::hoslem.test(fit1$y, fit1$fitted, g=6))
	Hosmer and Lemeshow goodness of fit (GOF) test
data:  fit1$y, fit1$fitted
X-squared = 2.1597, df = 4, p-value = 0.7064

> cbind(HL1$observed, HL1$expected)
               y0 y1     yhat0     yhat1
[0.0763,0.237] 35  7 35.261211  6.738789
(0.237,0.292]  16  9 17.705136  7.294864
(0.292,0.608]  14 16 11.772441 18.227559
(0.608,0.673]  11 19  9.821233 20.178767
(0.673,0.885]   2 28  3.439979 26.560021

SASでの実行は先ほどのコードから。

f:id:cochineal19:20210130231534p:plain

残念ながら結果は一致せず、どうやら bin の設定の仕方が異なるよう。
ということで、Rのコードを自作。

> #-----------------------------------------------------------
> #-- 自作 Hosmer and Lemeshow test, and Calibration Plot
> #-----------------------------------------------------------
> grp1 <- 6 #-- グループ数
> sep1 <- unique(c(0,quantile(fit1$fitted.values, probs = seq(0, 1, 1/grp1)))) #-区間
> out1 <- data.frame(y=fit1$y
+                    ,yhat=fit1$fitted.values
+                    ,ycut=cut(fit1$fitted.values, breaks=sep1, include.lowest=T))
> out2 <- out1 %>% 
+   dplyr::group_by(ycut) %>%
+   dplyr::summarise(Act=sum(y)/n()
+                    ,Pred=mean(yhat)
+                    ,N=n()
+                    ,O=sum(y)
+                    ,E=n()*mean(yhat)
+                    ,ChiSq=(sum(y) - n()*mean(yhat))^2 / (n()*mean(yhat)*(1-mean(yhat))))
> out3 <- data.frame(Chi_Square = sum(out2$ChiSq)
+                    ,P_value = 1 - pchisq(sum(out2$ChiSq), grp1 - 2))
> print(out2)
  ycut             Act   Pred     N     O     E ChiSq
                  
1 [0,0.0763]     0.05  0.0763    20     1  1.53 0.197
2 (0.0763,0.237] 0.273 0.237     22     6  5.21 0.156
3 (0.237,0.292]  0.36  0.292     25     9  7.29 0.563
4 (0.292,0.608]  0.533 0.608     30    16 18.2  0.694
5 (0.608,0.673]  0.633 0.673     30    19 20.2  0.210
6 (0.673,0.885]  0.933 0.885     30    28 26.6  0.681
> print(out3$Chi_Square)
[1] 2.500154
> print(out3$P_value)
[1] 0.6446081

自作プログラムではSASの結果と一致。
今回は p≒0.645 で帰無仮説は却下せず、モデルが適合していると評価する。

Rでキャリブレーションプロットも描いてみる。
先ほどのホスマー・レメショウ検定で使用したグループ毎の予測値と実測値をプロットする。45° 線付近にプロットがあり、大きなズレはなさそう。

ggplot(data=out2) +
  geom_abline(intercept=0, slope=1, linetype="dashed") +
  scale_x_continuous(limits=c(0,1), breaks=seq(0,1,0.25)) +
  scale_y_continuous(limits=c(0,1), breaks=seq(0,1,0.25)) +
  geom_point(mapping=aes(x=Pred, y=Act)) +
  labs(title="キャリブレーションプロット", x="予測値", y="実測値")

f:id:cochineal19:20210130232214p:plain

 

プログラムコード

■ Rのコード

#-----------------------------------------------------------------------------------
# Discrimination（判別能力）
# ROC、AUC、感度、特異度
#-----------------------------------------------------------------------------------
library(tidyverse)
library(pROC)
(roc1 <- pROC::roc(response=fit1$y, predictor=fit1$fitted, ci=T)) #-- glmの結果を投入
(roc1.res1 <- (pROC::coords(roc1, x="best", best.method="youden", ret=c("threshold","tpr","fpr","tnr","fnr","npv","ppv"), transpose=F)))
(roc1.res2 <- (pROC::coords(roc1, x="best", best.method="closest.topleft", ret=c("threshold","tpr","fpr","tnr","fnr","npv","ppv"), transpose=F)))
(AUC1 <- round(roc1$ci, 4))

ggroc(list(Model1=roc1), legacy.axes=T, size=1.5) +
  geom_abline(intercept=0, slope=1, linetype="dashed") +
  annotate("text", x=1-roc1.res1$tnr, y=roc1.res1$tpr, hjust=0.5, vjust=0.5, size=10, label="〇", color="blue") +
  annotate("text", x=1-roc1.res2$tnr, y=roc1.res2$tpr, hjust=0.5, vjust=0.5, size=10, label="×", color="blue") +
  annotate("text", x=1-roc1.res1$tnr, y=roc1.res1$tpr-0.05, hjust=0, vjust=0.5, size=5, label="youden",) +
  annotate("text", x=1-roc1.res2$tnr, y=roc1.res2$tpr-0.05, hjust=0, vjust=0.5, size=5, label="closest.topleft") +
  annotate("text", x=0.8, y=0.1, size=4
           ,label=paste0("Area Under the Curve ( 95%CI ) :\n",AUC1[2]," ( ",AUC1[1],", ",AUC1[3]," )")) +
  labs(title="ROC曲線")

#-----------------------------------------------------------------------------------
# Calibration（較正）
# Hosmer and Lemeshow goodness of fit (GOF) test
# Calibration plot
#-----------------------------------------------------------------------------------
grp1 <- 6 #-- グループ数
sep1 <- unique(c(0,quantile(fit1$fitted.values, probs = seq(0, 1, 1/grp1)))) #-- 区間
out1 <- data.frame(y=fit1$y                  #-- glmの結果を投入（実測値）
                   ,yhat=fit1$fitted.values  #-- glmの結果を投入（予測値）
                   ,ycut=cut(fit1$fitted.values, breaks=sep1, include.lowest=T) )
out2 <- out1 %>% 
  dplyr::group_by(ycut) %>%
  dplyr::summarise(Act=sum(y)/n()
                   ,Pred=mean(yhat)
                   ,N=n()
                   ,O=sum(y)
                   ,E=n()*mean(yhat)
                   ,ChiSq=(sum(y) - n()*mean(yhat))^2 / (n()*mean(yhat)*(1-mean(yhat))) )
out3 <- data.frame(Chi_Square = sum(out2$ChiSq)
                   ,P_value = 1 - pchisq(sum(out2$ChiSq), grp1 - 2))
print(out2)
print(out3$Chi_Square)
print(out3$P_value)

ggplot(data=out2) +
  geom_abline(intercept=0, slope=1, linetype="dashed") +
  scale_x_continuous(limits=c(0,1), breaks=seq(0,1,0.25)) +
  scale_y_continuous(limits=c(0,1), breaks=seq(0,1,0.25)) +
  geom_point(mapping=aes(x=Pred, y=Act)) +
  labs(title="キャリブレーションプロット", x="予測値", y="実測値")

 

■ SASのコード

ods output LackFitChiSq=HOSLEM1;
proc logistic data=ads descending;
  model y(event='1') = x1 x2 / outroc=ROC1 lackfit;
  output out=OutDS;
run;

 

■ Pythonのコード

整備中

 

参考

https://www.ism.ac.jp/~noma/2018-12-07%20JBS%20Seminar%20Kyoto.pdf

 

★サイトマップ

cochineal19.hatenablog.com

本ブログは個人メモです。 本ブログの内容によって生じた損害等の一切の責任を負いかねますのでご了承ください。