ロジスティック回帰分析の適合度検定であるHosmer-Lemeshow
test(ホスマー・レメショウ検定)のRでの実装について考えてみる-その3。
要旨としては
- Hosmer-Lemeshow testについて、十分位でサブグループを分割する方法がうまく行かない場合のサブグループの分割方法について、期待確率を小さい順に並べたときに同じ期待確率の値のデータは同じサブグループになるように調整しつつ、各サブグループのデータポイント数がおおよそ等しくなる分け方について考えてみた。
- 1つ目はデータポイント数が最小のサブグループは隣接するサブグループと統合する、というのをサブグループ数が任意のグループ数になるまで繰り返す処理の場合。
- 2つ目は隣接するサブグループと統合しサブグループ数が任意の数になるパターンを総当りで計算し、その中でサブグループ数の分散が最も小さい組み合わせを取り出す方法の場合。
Hosmer-Lemeshow testはグループ分けが重要である。
この回帰分析を用いた疫学データの解析 基本編（2）の解説の
代表的なものはモデルが予測する目的変数の値を小さいほうから順に並べ各グループに含まれるデータポイント数が大体同じになるようグループ分けしていく方法である。
という方法についてついて、前回は期待確率の小さい順に並べて、単純に任意のサブグループ数で分割するという単純な礼につついて考えた。
ただ、期待確率が一様に分散していない場合は上記の単純に分割するだけでは、同じ期待確率のものが別のサブグループに含まれることになってしまうので不適切になってしまうことが容易に想像される。
たとえば、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())
このモデルの各期待確率の数は次のようにして取り出すことができる。
expect_logit <- sapply(sort(unique(fitted.values(model))),
function(x) sum(fitted.values(model) == x))
names(expect_logit) <- sort(unique(fitted.values(model)))
expect_logit
## 0.103959413496167 0.198719327281383 0.225499724406818 0.251158568967817
## 462 168 862 48
## 0.407038204021279 0.417565455179058 0.56612912085157 0.664924908154581
## 175 11 165 5
## 0.736089652810118 0.766053806997806 0.790446281655777 0.885323441946764
## 93 23 31 144
## 0.889661176752348 0.957114111822833
## 13 1
このように、期待確率は偏って分布している。
ヒストグラムにすると次の通り。
hist(fitted.values(model))
このような場合のサブグループの分割方法について考えてみる。
これを10個のサブグループにするときに、各グループに含まれるデータポイントの数は
sum(expect_logit) / 10
## [1] 220.1
を目指したいが、隣接する期待確率のグループを足して10個のサブグループを作ろうとしても、そのような組み合わせはない。
ということで、もう少し条件を考えてグループ分けできるようにしてみる。
ここからは各サブグループに含まれる個数だけで考えてみる。各サブグループに含むべきデータポイント数さえわかれば、小さい順に並べてその数だけ取り出せばいい。
データポイントの数が最小のサブグループのデータポイントの数が最大になる（ようにしたい）
サブグループに含まれているデータポイント数の最小値が最大になる組み合わせで和を取り、長さが指定されたサブグループ数のベクトルを作るように考える。
vec_g <- function(x, g){
if(g <= 3) stop("If g is less than 3, it cannot be calculated well.")
ans <- x
ans_len <- length(ans)
while (ans_len > g){
ans_min <- which(ans == min(ans))
ans_list <- lapply(ans_min, function(ans_min){
if(ans_min == 1){
ans_ <- c(ans[ans_min] + ans[ans_min + 1], ans[(ans_min + 2):ans_len])
}else if(ans_min == ans_len){
ans_ <- c(ans[1:(ans_min - 2)], ans[ans_min - 1] + ans[ans_min])
}else{
if(ans_min == 2){
ans_pre <- c(ans[ans_min - 1] + ans[ans_min], ans[(ans_min + 1):ans_len])
}else{
ans_pre <- c(ans[1:(ans_min - 2)], ans[ans_min - 1] + ans[ans_min], ans[(ans_min + 1):ans_len])
}
if((ans_min + 2) == ans_len){
ans_pos <- c(ans[1:(ans_min - 1)], ans[ans_min] + ans[ans_min + 1], ans[ans_len])
}else if((ans_min + 1) == ans_len){
ans_pos <- c(ans[1:(ans_min - 1)], ans[ans_min] + ans[ans_min + 1])
}else{
ans_pos <- c(ans[1:(ans_min - 1)], ans[ans_min] + ans[ans_min + 1], ans[(ans_min + 2):ans_len])
}
ans_pre_var <- stats::var(ans_pre)
ans_pos_var <- stats::var(ans_pos)
if(ans_pre_var < ans_pos_var || is.na(ans_pos_var)) ans <- ans_pre
else ans_ <- ans_pos
}
})
ans <- unlist(unname(ans_list[which.min(lapply(ans_list, stats::var))]))
ans_len <- length(ans)
}
ans
}
これは、まずベクトル内の最小値を求め、その後、前後の値との和を取り、分散が小さい方を採用するというのを、ベクトルの長さが指定された長さになるまで繰り返すという処理になっている。
ところで分散についてはRで計算される分散は不偏分散になるが、結果に影響はない。
試しに、Titanicデータのモデルの場合で数を確認してみる。
vec_g(unname(expect_logit), 10)
## [1] 462 168 862 48 175 176 98 23 31 158
結構いい感じにサブグループを分割できるようになっていると思う。
ただし、この方法ではいい感じに計算できないパターンがあり、例えば次のようなモデルを考えてみる。
model_1 <- glm(vs ~ wt + mpg, data = mtcars, family = binomial())
このときの期待確率は次のようになっている。
explect_1 <- sapply(sort(unique(fitted.values(model_1))), function(x) sum(fitted.values(model_1) == x))
explect_1
## [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
length(explect_1)
## [1] 32
と、1が32個のベクトルになっている。
vec_g(explect_1, 10)
## [1] 4 4 4 4 4 4 2 2 2 2
この1が32個の場合は、4が2つでのこり8個が3であるのが最も分散が小さく、各サブグループに含まれるデータポイントの数が平均に近づくはず。
次に別の方法を考えてみる。
すべての組み合わせを計算してデータポイントの数の分散が最小のものを選択する
すべてのデータポイントの数が平均に近い組み合わというのは、データポイントの数の分散が最も小さい組み合わとみなせると思う。
隣接するグループ同士で合体させて、任意の長さのベクトルになる組み合わせを総当りで計算し、分散が最も小さい組み合わせを取り出す方法を考えてみる。
vec_n <- function(x, g){
if(g <= 3) stop("If g is less than 3, it cannot be calculated well.")
vec_add <- function(x){
lapply(1:(length(x) - 1), function(a){
if(a == 1) c(x[1] + x[2], x[3:length(x)])
else if(a + 2 == length(x)) c(x[1:(a - 1)], x[a] + x[a + 1], x[a + 2])
else if(a + 1 == length(x)) c(x[1:(a - 1)], x[a] + x[a + 1])
else if(a == 2) c(x[1], x[2] + x[3], x[4:length(x)])
else c(x[1:(a - 1)], x[a] + x[a + 1], x[(a + 2):length(x)])
})
}
ans <- list(x)
while(length(ans[[1]]) > g){
ans <- lapply(ans, vec_add)
ans <- unique(purrr::list_flatten(ans))
}
ans_n <- sapply(ans, stats::var)
if(length(ans_n) == 1) ans <- unlist(ans)
else if(length(min(ans_n)) == 1) ans <- unlist(ans[which.min(ans_n)])
else{
ans <- ans[ans_n == min(ans_n)]
ans_n <- sapply(ans, min)
ans <- unlist(ans[which.max(ans_n)])
}
ans
}
この場合のTitanicデータのモデルの個数を考えてみる。
vec_n(unname(expect_logit), 10)
## [1] 462 168 862 48 175 176 98 54 144 14
これは先程のサブグループに含まれるデータポイントの数の最小値が最大になるようにした手法と一致する場合もあるが、今回の場合は分散がこちらの方が小さくなっている。
vec_g(unname(expect_logit), 10) |> var()
## [1] 66692.77
vec_n(unname(expect_logit), 10) |> var()
## [1] 66403.21
ただし、この計算は指定された長さにいたるまで計算を続けるのでメモリや計算量が多くなる。
たとえば、mtcarのデータでのモデルの場合、長さ32のベクトルから10のベクトルに至るまでのすべての組み合わせを計算する（一部同じベクトルがでれば削除するが）のでかなり重たく、とりあえず自分の環境では計算が終わらなかった。
もっと効率の良い計算方法があると思うがこれが一番簡単にかけた。
それぞれの組み合わせでHosmer-Lemeshow testしてみる
これらの組み合わせでサブグループをHosmer-Lemeshow
testを計算するRの関数をgithubにあげている{infun}パッケージに入れてみたのでそちらで計算してみる。
github.com
simpleをTRUEにすると前項の単純に切り分ける手法になる。simpleをFALSEにしてforceをFALSEにするとサブグループに含まれるデータポイント数の最小値が最大になる組み合わせで統計量を計算する。simpleをFALSEにしてforceをTRUEにすると総当りでサブグループに含まれるデータポイント数の分散が最小になる組み合わで統計量を計算する。
となっている。
データポイント数の最小値が最大の組み合わせでサブグループを計算した場合は次のとおりになる。
HL_1 <- infun::hosmer_test(model, g = 10, force = FALSE)
HL_1
##
## Hosmer and Lemeshow goodness of fit test
##
## data: glm(formula = Survived ~ ., family = binomial(), data = df)
## X-squared = 82.909, df = 8, p-value = 1.266e-14
HL_1$observed
##
## subgroup y0_obs y1_obs
## 0.104 387 75
## 0.199 154 14
## 0.225 670 192
## 0.251 35 13
## 0.407 118 57
## 0.418 - 0.566 89 87
## 0.665 - 0.736 13 85
## 0.766 3 20
## 0.79 17 14
## 0.885 - 0.957 4 154
HL_1$expected
##
## subgroup y0_expect y1_expect
## 0.104 413.970751 48.029249
## 0.199 134.615153 33.384847
## 0.225 667.619238 194.380762
## 0.251 35.944389 12.055611
## 0.407 103.768314 71.231686
## 0.418 - 0.566 77.995475 98.004525
## 0.665 - 0.736 26.219038 71.780962
## 0.766 5.380762 17.619238
## 0.79 6.496165 24.503835
## 0.885 - 0.957 17.990715 140.009285
総当りでサブグループに含まれているデータポイント数の分散が最小になる組み合わせで統計量を計算すると次のとおりになる。
HL_2 <- infun::hosmer_test(model, g = 10, force = TRUE)
HL_2
##
## Hosmer and Lemeshow goodness of fit test
##
## data: glm(formula = Survived ~ ., family = binomial(), data = df)
## X-squared = 67.253, df = 8, p-value = 1.727e-11
HL_2$observed
##
## subgroup y0_obs y1_obs
## 0.104 387 75
## 0.199 154 14
## 0.225 670 192
## 0.251 35 13
## 0.407 118 57
## 0.418 - 0.566 89 87
## 0.665 - 0.736 13 85
## 0.766 - 0.79 20 34
## 0.885 4 140
## 0.89 - 0.957 0 14
HL_2$expected
##
## subgroup y0_expect y1_expect
## 0.104 413.970751 48.029249
## 0.199 134.615153 33.384847
## 0.225 667.619238 194.380762
## 0.251 35.944389 12.055611
## 0.407 103.768314 71.231686
## 0.418 - 0.566 77.995475 98.004525
## 0.665 - 0.736 26.219038 71.780962
## 0.766 - 0.79 11.876928 42.123072
## 0.885 16.513424 127.486576
## 0.89 - 0.957 1.477291 12.522709
結構統計量が異なる事がわかる。ただ、いずれの場合もそれぞれのサブグループに含まれるデータポイント数が均等とは言い難い。
そもそも、サブグループに含まれるデータポイント数がある程度10個のサブグループで均等にできる場合には十分位で充分に均等にデータが分割できる状態にあると思うので、そちらの方でサブグループに分けたほうがいいと思う。こういう手法じゃないと均等にサブグループのデータポイント数が均等にならないパターンを考えてみたけど思い至らなかった。
ところで、この場合、サブグループ数は期待確率のユニークな値の数を超えることはできない。
なので、たとえば、今回の場合でいうと15個以上のサブグループをつくって計算しようとするとエラーが出るようにしている。
infun::hosmer_test(model, 17)
## Error in infun::hosmer_test(model, 17): Cannot be calculated because the unique value of expected value is less than the number of groups specified.
また、サブグループ数と期待確率のユニークな値の数が同じ場合は、期待確率のグループのまま計算される。
HL_3 <- infun::hosmer_test(model, 14)
HL_3
##
## Hosmer and Lemeshow goodness of fit test
##
## data: glm(formula = Survived ~ ., family = binomial(), data = df)
## X-squared = 103.83, df = 12, p-value < 2.2e-16
HL_3$observed
##
## subgroup y0_obs y1_obs
## 0.104 387 75
## 0.199 154 14
## 0.225 670 192
## 0.251 35 13
## 0.407 118 57
## 0.418 0 11
## 0.566 89 76
## 0.665 0 5
## 0.736 13 80
## 0.766 3 20
## 0.79 17 14
## 0.885 4 140
## 0.89 0 13
## 0.957 0 1
HL_3$expected
##
## subgroup y0_expect y1_expect
## 0.104 413.97075096 48.02924904
## 0.199 134.61515302 33.38484698
## 0.225 667.61923756 194.38076244
## 0.251 35.94438869 12.05561131
## 0.407 103.76831430 71.23168570
## 0.418 6.40677999 4.59322001
## 0.566 71.58869506 93.41130494
## 0.665 1.67537546 3.32462454
## 0.736 24.54366229 68.45633771
## 0.766 5.38076244 17.61923756
## 0.79 6.49616527 24.50383473
## 0.885 16.51342436 127.48657564
## 0.89 1.43440470 11.56559530
## 0.957 0.04288589 0.95711411
これはつまりどういうことかと言うと、今回の場合は当てはまらないが、期待確率のユニークな値の数が10個を下回るときには10個のサブグループは作れないよということ。
もちろん、これを指定されたサブグループ数が大きすぎるから計算出来ないよ～とワーニングメッセージを出しながら計算可能な最大のサブグループ数で計算し直した結果を返す処理にしても良いと思うが、勝手に自由度が変わる処理は混乱を招きかねないかなと思っている。
最後の方は自作関数の紹介になってしまった。
蛇足
こちらで現在CRANで公開されている{ResourceSelection}のhoslme.test()では期待確率の出現頻度が偏っているときには、正しい答えを出していない可能性について示した（この問題についてGithub上で正しい答えを出していない可能性について、Issueで報告し修正してもらっている）。
indenkun.hatenablog.com
エラーが出ているわけでもないので一見正しい結果を返しているように見えるが、$observedや$expectedで中身を確認すると期待したサブグループ数で分割されていないのに、期待されたサブグループ数で分割された場合の自由度のχ二乗分布で検定されているという問題だった。
例えばHosmer-Lemeshow
testで計算された統計量が15のとき、自由度8（サブグループ数10）の場合と自由度6（サブグループ数8）の場合では次のようにp値が異なる。
pchisq(15, 8, lower.tail = FALSE)
## [1] 0.05914546
pchisq(15, 6, lower.tail = FALSE)
## [1] 0.02025672
つまり、サブグループ数10を指定し十分位を取るときに重複する値が有り、サブグループ数が実際には8になってしまった場合に、本来は自由度6でp値が計算されるべきなのに、自由度8でp値が計算されてしまった場合にはモデルは適合していると判断される（されていた）可能性があるということになる。
このことについて色々思うところはあるけれど、それなりに長く使われていて知名度のあるパッケージの関数でこういう単純な謎挙動があるのは、きっと誰か気付いていたんじゃないかなーとは思わななくもない（調べた範囲ではこのことに言及しているものは見当たらなかったけど）。
