ロジスティック回帰分析の適合度検定の手法の一つであるホスマー・レメショウ検定（Hosmer-Lemeshow test）*1について、Rにおける実装について考えてみる。
要旨としては
- RでHosmer-Lemeshow
testを実行する
{ResourceSelection}のhoslem.test()では0.3.5を含む以前のバージョンではロジスティック回帰分析のモデルによっては正しいp値が計算されていない可能性がある。 {performance}のhosmer_test()は{ResourceSelection}のhoslem.test()で正しいp値が計算されないロジスティック回帰分析のモデルの場合は計算が完了せずエラーでストップする。
Hosmer-Lemeshow testについては邦文だと次の総説が詳しい。
Wikipediaも参照のこと。
Hosmer-Lemeshow testは、データをn個のグループに分けて（多くの場合、十分位で）観測結果と回帰式からの予測値（予測確率・期待値・リスク値）との残差を用いて統計量の計算を行い、その統計量が自由度n-2のカイ二乗分布に従うというアイディアに依る。
Hosmer氏が次のように述べている通り*2、Hosmer-Lemeshow testには次のようなメリットがある。
The advantage of the Hosmer-Lemeshow type tests is that they are based on groupings of theestimated probabilities that are intuitively appealing and easily understood by subject matterscientists.
自画自賛である。
ただし、Hosmer氏が言う通り、
The disadvantage is that the value of the statistic depends on the choice of cutpointsthat define the groups.
というデメリットもある。グループ分けがメリット・デメリットになっているというわけで、Hosmer-Lemeshow testではグループ分けにかかわる処理というのは重要ということがわかる。
まずはn個のグループに分け統計量を計算し、その統計量の自由度n-2のカイ二乗分布に従うという点について考えてみる。
現在、Rで比較的広く用いられているHosmer-Lemeshow
testの実装は{ResourceSelection}パッケージのhoslem.test()だと思う（Rによる多変量解析入門(Rによる多変量解析入門 データ分析の実践と理論 | Ohmsha)などでも紹介されている）。あとは、比較的広く用いられているパッケージとしては回帰分析等のモデルのパフォーマンスを調べるための関数群の{performance}パッケージのperformance_hosmer()がある。
これらは2つの関数はともに、ロジスティック回帰分析の回帰式の期待値の分位数でn個のグループをつくり、統計量を計算している。ここで既定値では十分位を取るようになっている。
とりあえずmtcarsデータでglm()で適当にロジスティック回帰分析を行ってみる。
model <- glm(vs ~ wt + mpg, data = mtcars, family = binomial())
このロジスティック回帰分析の適合度検定を{ResourceSelection}パッケージのhoslem.test()で計算してみる。既定値でgは10、つまり十分位でのグループ分けになっているが一応明示しておく。
library(ResourceSelection)
## ResourceSelection 0.3-5 2019-07-22
HL <- hoslem.test(model$y, fitted.values(model), g = 10) HL
##
## Hosmer and Lemeshow goodness of fit (GOF) test
##
## data: model$y, fitted.values(model)
## X-squared = 5.1369, df = 8, p-value = 0.7429
既定値だと10個のグループに分けられているので、自由度8のカイ二乗分布でp値が計算されている。
このとき、観測値と期待値のグループ分けされた結果は次のようにして取り出せる。
HL$observed
##
## cutyhat y0 y1
## [0.0174,0.051] 4 0
## (0.051,0.0828] 3 0
## (0.0828,0.158] 3 0
## (0.158,0.243] 2 1
## (0.243,0.369] 2 1
## (0.369,0.52] 2 1
## (0.52,0.615] 1 2
## (0.615,0.868] 0 3
## (0.868,0.982] 1 2
## (0.982,0.998] 0 4
HL$expected
##
## cutyhat yhat0 yhat1
## [0.0174,0.051] 3.87984951 0.12015049
## (0.051,0.0828] 2.77779066 0.22220934
## (0.0828,0.158] 2.67732615 0.32267385
## (0.158,0.243] 2.38429554 0.61570446
## (0.243,0.369] 2.06107259 0.93892741
## (0.369,0.52] 1.67780793 1.32219207
## (0.52,0.615] 1.35146089 1.64853911
## (0.615,0.868] 0.90896607 2.09103393
## (0.868,0.982] 0.24958571 2.75041429
## (0.982,0.998] 0.03184495 3.96815505
たしかに、10個のグループに分けられていることが確認できる。
次にTitanicデータでロジスティック回帰分析を考えてみる。
data("Titanic") df <- data.frame(Titanic) df <- data.frame(Class = rep(df$Class, df$Freq), Sex = rep(df$Sex, df$Freq), Age = rep(df$Age, df$Freq), Survived = rep(df$Survived, df$Freq)) model <- glm(Survived ~ . ,data = df, family = binomial())
このHosmer-Lemeshow testは次のようになる。
HL <- hoslem.test(model$y, fitted.values(model), g = 10) HL
##
## Hosmer and Lemeshow goodness of fit (GOF) test
##
## data: model$y, fitted.values(model)
## X-squared = 16.733, df = 8, p-value = 0.03301
このときの観測値と期待値のグループ分けは次のようになっている。
HL$observed
##
## cutyhat y0 y1
## [0.104,0.225] 1211 281
## (0.225,0.407] 153 70
## (0.407,0.566] 89 87
## (0.566,0.736] 13 85
## (0.736,0.957] 24 188
HL$expected
##
## cutyhat yhat0 yhat1
## [0.104,0.225] 1216.20514 275.79486
## (0.225,0.407] 139.71270 83.28730
## (0.407,0.566] 77.99548 98.00452
## (0.566,0.736] 26.21904 71.78096
## (0.736,0.957] 29.86764 182.13236
10個のグループ分けになるようにしているのに、5個のグループ分けになっている。
これは{ResourcesSelection}のhoslme.test()のコード中でグループ分けする範囲（期待値のどこからどこまでを1つにするのか）を決めている次の箇所による。
yhat <- y qq <- unique(quantile(yhat, probs = seq(0, 1, 1/g))) cutyhat <- cut(yhat, breaks = qq, include.lowest = TRUE)
yはモデルの期待値が入力され、quantile(yhat, probs = seq(0, 1, 1/g))で十分位の値が計算されることとなる。
ここでunique()をとっており、被っていると消えてしまう。
Titanicでのロジスティック回帰分析の期待値は次のように分布している。
yhat <- fitted.values(model) hist(yhat)
ということで分布が著しく偏っている。
十分位を見てみると次のようになる。
quantile(yhat, seq(0.1, 1, 0.1))
## 10% 20% 30% 40% 50% 60% 70% 80%
## 0.1039594 0.1039594 0.2254997 0.2254997 0.2254997 0.2254997 0.4070382 0.5661291
## 90% 100%
## 0.7360897 0.9571141
十分位を取るときに同じ値が複数でてくる。そのため、これをunique()した結果からグループ分けすると10個未満のグループになってしまう。
少し前にもどって結果をTitanicのロジスティック回帰分析の結果は本来は5個のグループに分けられて計算された統計量は、自由度3のカイ二乗分布でp値が計算されなければならないので、この場合の本当のp値は次のようになる。
1 - pchisq(HL$statistic, 3)
## X-squared
## 0.0008018872
先に上げた結果とぜんぜん違う（もともと適合していないという結果ではある）。
なんでunique()しているのかわからない（十分位を計算できなかったときにエラーになるのを避けたかったんだろうと推察できるが）（英語版Wikipediaでもそういう実装のコードが参照されていた、コードの発想はメーリングリストによるとDocumentsに記載されているので同源なのかもしれない）。
上記の問題についてはGithub上でISSUEで報告しており、Github上で公開されているバージョンのパッケージでは指定されたグループ分け数ではなく実際に分けられたグループ数でカイ二乗分布の自由度を決める処理になっている（CRANには現在のところまだ上がっていない）。
ところで、もう一つの{performance}パッケージのhosmer_test()はこのように一意に分位数を決定できない場合にはエラーが出るようになっている。
library(performance) performance_hosmer(model, n_bins = 10)
## Error in cut.default(yhat, breaks = stats::quantile(yhat, probs = seq(0, : 'breaks' が一意的ではありません
これは、{ResourcesSelection}のhoslem.test()で分位で分割する範囲をを計算している箇所に該当する次のコードの所でエラーが出ている。
cutyhat <- cut(yhat, breaks = stats::quantile(yhat, probs = seq(0, 1, 1/n_bins)), include.lowest = TRUE)
{ResourcesSelection}のhoslem.test()とは異なり、unique()していないため一意に分割できないためにエラーが出ている。つまり指定された分位でちゃんと計算できないときにはエラーが出てストップするという硬派な処理になっている。
ところで、このエラー文だけではなんでエラーになっているのかわかりにくいので、「指定されたn_binsでは分割できないから計算できません。」みたいなエラー文にしてほしいな～とISSUEでお願いしてみているが、どうなるかはわからない。
少し長くなってきたので、グループ分けについてもう少し考えてみたが次の稿に分けることにする。
*1:David W. Hosmer, Stanley Lemesbow (1980). Goodness of fit tests for the multiple logistic regression model, Communications in Statistics - Theory and Methods, 9:10, 1043-1069, doi:10.1080/03610928008827941
*2:HOSMER, D.W., HOSMER, T., LE CESSIE, S. and LEMESHOW, S. (1997), A COMPARISON OF GOODNESS-OF-FIT TESTS FOR THE LOGISTIC REGRESSION MODEL. Statist. Med., 16: 965-980. doi:10.1002/(SICI)1097-0258(19970515)16:9<965::AID-SIM509>3.0.CO;2-O
