備忘ログ

チラシの裏的備忘録&メモ

ロジスティック回帰分析の適合度検定であるHosmer-Lemeshow test(ホスマー・レメショウ検定)のRでの実装について考えてみる-その2

ロジスティック回帰分析の適合度検定であるHosmer-Lemeshow test(ホスマー・レメショウ検定)のRでの実装について考えてみる-その2。

要旨としては

  • 単純にデータを期待確率の小さい順に並べて上から等分に指定されたサブグループ数で分けた場合でのHosmer-Lemeshow testについて考えてみた。

Hosmer-Lemeshow testはグループ分けが重要である。

この点については次の解説が詳しい。

回帰分析を用いた疫学データの解析 基本編(2)

このグループ分けについて、Rで多く使用されている{ResourcesSelection}パッケージのhoslem.test()の実装では回帰式の予測確率の分位数を用いてグループ分けをする手法が用いられているが、意図したグループ数で分けられない可能性があること(と現在CRANで公開されている仕様では誤ったp値を算出している可能性があること)をこの記事で示した。

ところで、総説にも記載されているとおり、グループ分けの手法は十分位に限られない。(ただし、予測確率の十分位が手法として優れている可能性が示唆されてはいる。)

この解説の

代表的なものはモデルが予測する目的変数の値を小さいほうから順に並べ各グループに含まれるデータポイント数が大体同じになるようグループ分けしていく方法である。

という方法についてついて。

つまり、同じ予測確率を示すデータは同じグループになるようにしつつ、各グループの数が同じになるように隣接するグループと合体させて一つのグループにしていくということである。

具体例は解説にある。

たとえば、予測確率とそれぞれのデータのn数が次のようにあった場合、

expect n
0.01 21
0.11 15
0.13 5
0.22 18
0.33 21
0.45 20
0.56 20
0.59 20
0.68 20
0.78 20
0.89 20

n数がおおよそ同じように10このグループに分けるときに、隣接するグループで統合してグループ分けし次のようなサブグループを作る作業をすると良いということになる。例の場合で言えば、期待確率0.11と0.13のグループを合体させるとよさそう。

つまり次のようなサブグループとすると良さそう。

subgroup expect n
1 0.01 21
2 0.11 - 0.13 20
3 0.22 18
4 0.33 21
5 0.45 20
6 0.56 20
7 0.59 20
8 0.68 20
9 0.78 20
10 0.89 20

この塩梅で、まずは最もシンプルな場合を考えてみる。

mtcarsのデータで次のようにロジスティック回帰分析を行ってみる。

model <- glm(vs ~ wt + mpg, data = mtcars, family = "binomial")

この時の予測確率は一様に分散しており、すべて異なる値となっている。

# 値がすべて異なっているということは、値から重複する値を削除したものと、元の値の長さは同じになるはず
length(unique(fitted.values(model))) == length(fitted.values(model))
## [1] TRUE

というわけで、この場合は予測確率を昇順に並べて可能な限り均等に10分割してもグループ分けとしてはそれほどおかしな結果にはならないはずである。

ざっと書いてみる。

hosmer_test <- function(model, g = 10){
  
  df <- g - 2
  out <- cbind.data.frame(y = model$y,
                          fitted.values = model$fitted.values)
  out <- out[order(out$fitted.values, out$y, decreasing = FALSE),]
  step_by <- nrow(out) / g
  subgroup_step <- round(seq(step_by, nrow(out), step_by))
  n_step <- c(subgroup_step[1], subgroup_step[-1] - subgroup_step[1:(length(subgroup_step) - 1)])
  out$subgroup <- unlist(mapply(function(x, y) rep(x, y),
                                1:length(n_step),
                                n_step,
                                SIMPLIFY = FALSE))
  
  obs <- stats::xtabs(cbind(y0_obs = 1 - y, y1_obs = y) ~ subgroup, data = out)
  expect <- stats::xtabs(cbind(y0_expect = 1 - fitted.values, y1_expect = fitted.values) ~ subgroup, data = out)
  chisq <- sum((obs - expect)^2 / expect)
  p.value <- 1 - stats::pchisq(chisq, df)
  
  names(chisq) <- "X-squared"
  names(df) <- "df"
  
  HL <- list(method = "Hosmer and Lemeshow goodness of fit test",
             parameter = df,
             p.value = p.value,
             statistic = chisq,
             data.name = paste(deparse(model$call)),
             observed = obs,
             expected = expect)
  class(HL) <- "htest"
  HL
}

処理としては、予測確率を昇順に並べて小さい順に等間隔になるように1~10の番号を降っていくという処理をしそれをサブグループとして、統計量を統計量とするものである。

データがグループ数の整数倍ではない場合、どこかでズレが生じるので適当に四捨五入している。

このときのHosmer-Lemeshow testの結果は次の通り。

HL <- hosmer_test(model, 10)
HL
## 
##  Hosmer and Lemeshow goodness of fit test
## 
## data:  glm(formula = vs ~ wt + mpg, family = "binomial", data = mtcars)
## X-squared = 7.5898, df = 8, p-value = 0.4745

この時の観測値と予測確率のグループ分けは次のとおりとなっている。

HL$observed
##         
## subgroup y0_obs y1_obs
##       1       3      0
##       2       3      0
##       3       4      0
##       4       2      1
##       5       2      1
##       6       2      1
##       7       1      2
##       8       0      4
##       9       1      2
##       10      0      3
HL$expected
##         
## subgroup  y0_expect  y1_expect
##       1  2.92882139 0.07117861
##       2  2.81097948 0.18902052
##       3  3.59516545 0.40483455
##       4  2.38429554 0.61570446
##       5  2.06107259 0.93892741
##       6  1.67780793 1.32219207
##       7  1.35146089 1.64853911
##       8  1.01755249 2.98244751
##       9  0.15475854 2.84524146
##       10 0.01808572 2.98191428

{ResourcesSelection}で十分位で計算した時の結果は次の通りとなる。

HL <- ResourceSelection::hoslem.test(model$y, fitted.values(model), 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
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の整数倍ではないためどこかでズレが生じているためである。

つまり、データ数が10の整数倍であれば同じ結果になる。たとえばmtcarsデータの上から20個のデータで見てみると同じになるはずである。

model <- glm(vs ~ wt + mpg, data = mtcars[1:20,], family = "binomial")
# このモデルも20個の予測確率をもち、それらは同じ値がない。
length(fitted.values(model))
## [1] 20
length(unique(fitted.values(model))) == length(fitted.values(model))
## [1] TRUE
# 予測確率の小さい方から順に並べてサブグループのサンプル数が同じになるようにしたHosmer-Lemeshow testの場合
HL <- hosmer_test(model, g = 10)
HL
## 
##  Hosmer and Lemeshow goodness of fit test
## 
## data:  glm(formula = vs ~ wt + mpg, family = "binomial", data = mtcars[1:20,     ])
## X-squared = 2.4293, df = 8, p-value = 0.9649
HL$observed
##         
## subgroup y0_obs y1_obs
##       1       2      0
##       2       2      0
##       3       2      0
##       4       1      1
##       5       1      1
##       6       1      1
##       7       1      1
##       8       0      2
##       9       0      2
##       10      0      2
HL$expected
##         
## subgroup    y0_expect    y1_expect
##       1  1.9814201264 0.0185798736
##       2  1.9098763622 0.0901236378
##       3  1.6152333372 0.3847666628
##       4  1.5059246274 0.4940753726
##       5  1.2942785467 0.7057214533
##       6  0.8876493896 1.1123506104
##       7  0.5171035709 1.4828964291
##       8  0.2618654900 1.7381345100
##       9  0.0263991603 1.9736008397
##       10 0.0002493892 1.9997506108
# ResourcesSelection}のhoslem.test()で十分位で計算した時
HL <- ResourceSelection::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 = 2.4293, df = 8, p-value = 0.9649
HL$observed
##                  
## cutyhat           y0 y1
##   [0.00844,0.025]  2  0
##   (0.025,0.162]    2  0
##   (0.162,0.228]    2  0
##   (0.228,0.284]    1  1
##   (0.284,0.449]    1  1
##   (0.449,0.641]    1  1
##   (0.641,0.809]    1  1
##   (0.809,0.929]    0  2
##   (0.929,0.999]    0  2
##   (0.999,1]        0  2
HL$expected
##                  
## cutyhat                  yhat0        yhat1
##   [0.00844,0.025] 1.9814201264 0.0185798736
##   (0.025,0.162]   1.9098763622 0.0901236378
##   (0.162,0.228]   1.6152333372 0.3847666628
##   (0.228,0.284]   1.5059246274 0.4940753726
##   (0.284,0.449]   1.2942785467 0.7057214533
##   (0.449,0.641]   0.8876493896 1.1123506104
##   (0.641,0.809]   0.5171035709 1.4828964291
##   (0.809,0.929]   0.2618654900 1.7381345100
##   (0.929,0.999]   0.0263991603 1.9736008397
##   (0.999,1]       0.0002493892 1.9997506108

すべてのサグブループが均等に2ずつになり、統計量も小さい方から均等に取り出した場合と、十分位で計算した場合で同じになる。

とりあえずまずは、最もシンプルな10分割について考えてみた。もちろん、この手法では予測確率が一様に分散していなければ、異なるサブグループに同じ予測確率が割りふりされることになるので不適切である。

もう少し複雑なパターンについては次項で考える。