ロジスティック回帰をやるときに完全分離状態っていうのに陥るときがあって、例えば以下のようなデータセット
Y <- c(1,1,1,1,0,0,0,0) X <- c(8,7,3,1,0,-2,-1,-9) dat <- data.frame(Y,X) dat
> dat Y X 1 1 8 2 1 7 3 1 3 4 1 1 5 0 0 6 0 -2 7 0 -1 8 0 -9
これは目的変数Yが、Xの値が0より大きければ1、0より小さければ0というような例。
具体例では、例えばXが製品の質を表している数値で、この数値が0より大きければ合格、0以下なら欠陥品として処分。みたいな感じ。
で、このような例ではXに対し0てYの生起を「確率的」に評価するのが困難になる。だって手元のデータでは「X>0なら100%の確率で合格」「X≦0なら100%で不合格」のようになってしまって、Xの値がこれくらい大きくなると、Yが1になる確率はこのくらい大きくなる、みたいな表現ができないため。
で、これを実際にRで実行するとこうなる
glm(data = dat , Y ~ X , family=binomial(link = "logit"))
> glm(data = dat , Y ~ X , family=binomial(link = "logit"))
Call: glm(formula = Y ~ X, family = binomial(link = "logit"), data = dat)
Coefficients:
(Intercept) X
-21.71 43.55
Degrees of Freedom: 7 Total (i.e. Null); 6 Residual
Null Deviance: 11.09
Residual Deviance: 1.402e-09 AIC: 4
警告メッセージ:
1: glm.fit: アルゴリズムは収束しませんでした
2: glm.fit: 数値的に 0 か 1 である確率が生じました
数値的に0か1である確率が生じましたってのは、X>0なら確率1でYが起こるっていうことかな。
まあとりあえずこんな感じでエラーが起こるわけ。
他にも、準完全分離っていう状態では、以下のようなデータ
Y <- c(1,1,1,1,0,1,0,0) X <- c(8,7,3,1,0,0,-1,-9) dat2 <- data.frame(Y,X) dat2
> dat2 Y X 1 1 8 2 1 7 3 1 3 4 1 1 5 0 0 6 1 0 7 0 -1 8 0 -9
これは、X>0なら100%でYが1、X<0なら100%でYが0、そして境界線のX=0ではYが1になったりならなかったり。こういう状態は準完全分離と呼ばれて、やはりこれもYを確率的に評価するのが厳しくなる。
glm(data = dat2 ,Y ~ X , family = binomial(link="logit"))
> glm(data = dat2 ,Y ~ X , family = binomial(link="logit")) Call: glm(formula = Y ~ X, family = binomial(link = "logit"), data = dat2) Coefficients: (Intercept) X -3.333e-15 1.989e+01 Degrees of Freedom: 7 Total (i.e. Null); 6 Residual Null Deviance: 10.59 Residual Deviance: 2.773 AIC: 6.773 警告メッセージ: glm.fit: 数値的に 0 か 1 である確率が生じました
やはり警告が出る。まあ一応推定自体はできるけど。
でまあ、これを解決するためにFirthの方法だったりの各種補正法があるらしいんだけどそれは置いといて、以下のようなデータを見てみる。
Y <- c(1,1,1,1,0,0,0,0) X <- c(1,1,1,1,0,0,0,0) dat3 <- data.frame(Y,X) dat3
> dat3 Y X 1 1 1 2 1 1 3 1 1 4 1 1 5 0 0 6 0 0 7 0 0 8 0 0
Xがカテゴリカルで、Yもまたカテゴリカルで、かつ完全分離が起きているような状態。これはどう見ても完全分離状態だが、
Rで以下を実行すると、
glm(data = dat3 , Y ~ X , family = binomial(link="logit"))
> glm(data = dat3 , Y ~ X , family = binomial(link="logit"))
Call: glm(formula = Y ~ X, family = binomial(link = "logit"), data = dat3)
Coefficients:
(Intercept) X
-24.57 49.13
Degrees of Freedom: 7 Total (i.e. Null); 6 Residual
Null Deviance: 11.09
Residual Deviance: 3.429e-10 AIC: 4
警告メッセージが出てこない。……なんで？？？(マジで)
まあ、多分内部のアルゴリズムが上手く収束しているってことなんだろうけど、そんなことある？
警告メッセージが表示されないってことはRのglmで目的変数も説明変数も両方ともカテゴリカルであるようなロジスティック回帰をして変な推定が返ってきたとして、なんとそれに気づけないってことになる。
いやしかし、待ってほしい。以下のデータを見てみよう。
set.seed(1) X1 <- c(rep(1,10) , rep(0,10)) X2 <- c(rep(1,15) , rep(0,5)) X3 <- c(rep(1,3) , rep(0,17)) X4 <- c(rep(1,10) , rep(0,10)) X5 <- c(rep(1,18) , rep(0,2)) X6 <- c(rep(1,16) , rep(0,4)) #X1,X4はそこそこ出てきて、そこそこの効果 #X2は頻繁に出てきて、強めの効果 #X3はごくまれに出てきて、非常に強い効果 #X5,X6は頻繁に出てきて、弱い効果 X1 <- sample(X1 , replace = F) X2 <- sample(X2 , replace = F) X3 <- sample(X3 , replace = F) X4 <- sample(X4 , replace = F) X5 <- sample(X5 , replace = F) X6 <- sample(X6 , replace = F) logitP <- -2*X1 + 4*X2 + 15*X3 - 2.5*X4 + 0.2*X5 -1.2*X6 expitP <- exp(logitP)/(1+exp(logitP)) expitP Y <- rbinom(n=20,size=1, prob=expitP) data2 <- data.frame(Y ,X1,X2,X3,X4,X5,X6) data2
> data2 Y X1 X2 X3 X4 X5 X6 1 1 1 1 0 0 1 1 2 0 1 1 0 1 1 0 3 0 1 1 0 1 1 1 4 0 1 0 0 1 1 1 5 0 0 0 0 0 1 1 6 1 0 1 0 1 1 1 7 1 0 1 1 1 1 1 8 1 0 0 0 0 1 0 9 1 0 1 0 0 1 1 10 0 1 1 0 1 1 1 11 1 0 1 1 0 1 1 12 1 1 1 0 0 1 1 13 0 1 1 0 0 1 1 14 1 0 1 0 0 1 0 15 0 1 1 0 1 1 1 16 1 0 0 1 0 0 1 17 0 0 0 0 1 1 0 18 0 1 1 0 1 1 1 19 1 0 1 0 0 0 1 20 0 1 1 0 1 1 1
たとえばX3(あとX5も)に注目してみると、これは準完全分離が起きている。
これに対してロジスティック回帰を行ってもアルゴリズムが上手く推定できないと思われる。
glm(data = data2 , Y~X1+X2+X3+X4+X5+X6 , family="binomial")
> glm(data = data2 , Y~X1+X2+X3+X4+X5+X6 , family="binomial")
Call: glm(formula = Y ~ X1 + X2 + X3 + X4 + X5 + X6, family = "binomial",
data = data2)
Coefficients:
(Intercept) X1 X2 X3 X4 X5 X6
60.084 -81.116 101.752 2.379 -60.798 -39.962 -40.065
Degrees of Freedom: 19 Total (i.e. Null); 13 Residual
Null Deviance: 27.73
Residual Deviance: 3.819 AIC: 17.82
警告メッセージ:
glm.fit: 数値的に 0 か 1 である確率が生じました
あれ？警告メッセージが出てきた。
いやまあ、助かる。助かるけど、なんでさっきは警告が出てこなかったのに今度は警告が出てくるんだ？？？
マジで謎なんで、誰か詳しい人教えてください。
おまけ Firthの補正
せっかくなんでやっとこう。
Firthの補正は罰則項にフィッシャー情報量を利用した罰則付き最尤法になってるけど、まあ理論的なとこは各自で調べてください。
切片なしでまずやる。
#Firthの方法を用いた罰則つき最尤法を行う
install.packages("logistf")
library("logistf")
s <- logistf(data = data2 , Y ~0+X1+X2+X3+X4+X5+X6 , firth = T)
s
#切片なしなら一応上手く推定できた。
> s
logistf(formula = Y ~ 0 + X1 + X2 + X3 + X4 + X5 + X6, data = data2,
firth = T)
Model fitted by Penalized ML
Confidence intervals and p-values by Profile Likelihood
Coefficients:
X1 X2 X3 X4 X5 X6
-3.2825732 4.2851815 2.3987510 -3.0169298 0.9003723 -1.5259974
Likelihood ratio test=13.32006 on 6 df, p=0.0382258, n=20
まあよくわらんがとりあえずエラーもなく一応最尤推定値が出てきた。(まあ厳密にいえばこれは最尤推定値ではなく罰則付き最尤推定とでも呼ぶべきものだけど)
一応、識別率をチェックしてみる。
(predict(s)>0) == data2$Y #20個中1個だけ誤判別。
> (predict(s)>0) == data2$Y #20個中1個だけ誤判別。 [1] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE FALSE TRUE TRUE TRUE TRUE TRUE [19] TRUE TRUE
20個中1個だけ誤判別とのこと。
ま、切片なしなら上手くいった。しかし切片ありで推定すると……
logistf(data = data2 , Y ~X1+X2+X3+X4+X5+X6 , firth = T) #切片ありだと警告が出る。色々調整しなきゃらしいので調整する。
> logistf(data = data2 , Y ~X1+X2+X3+X4+X5+X6 , firth = T)
logistf(formula = Y ~ X1 + X2 + X3 + X4 + X5 + X6, data = data2,
firth = T)
Model fitted by Penalized ML
Confidence intervals and p-values by Profile Likelihood
Coefficients:
(Intercept) X1 X2 X3 X4 X5 X6
0.2255353 -2.5385945 3.1778146 1.7764449 -2.6034431 0.6743863 -1.2892589
Likelihood ratio test=12.38409 on 6 df, p=0.05392869, n=20
警告メッセージ:
1: logistf(data = data2, Y ~ X1 + X2 + X3 + X4 + X5 + X6, firth = T) で:
Maximum number of iterations exceeded. Try to increase the number of iterations or alter step size by passing 'logistf.control(maxit=..., maxstep=...)' to parameter control
2: logistf(data = data2, Y ~ X1 + X2 + X3 + X4 + X5 + X6, firth = T) で:
Nonconverged PL confidence limits: maximum number of iterations for variables: X6 exceeded. Try to increase the number of iterations by passing 'logistpl.control(maxit=...)' to parameter plcontrol
3: logistf(data = data2, Y ~ X1 + X2 + X3 + X4 + X5 + X6, firth = T) で:
Maximum number of iterations for variables: (Intercept), X1, X2, X3, X5, X6 exceeded. P-value may be incorrect. Try to increase the number of iterations by passing 'logistf.control(maxit=...)' to parameter control
なんかうるさいのがいっぱい出てきた。
どうやらイテレーションを増やせ、maxstepを増やせ、とのこと。PLってのはプロファイル尤度。信頼区間がどうのこうの言われてしまった。(ちなみにプロファイル尤度っていうのはパラメータが多次元の場合に、注目しているパラメータについて固定し、それ以外の余分なパラメータについて最大化した尤度のこと。特定のパラメータの方向から尤度関数の"輪郭(profile)"を見たときに、その輪郭っていうのはその他の余り物のパラメータについて最大化した値になってることから名づけられてる。わかりやすくて良いネーミングだと思う。逸脱度統計量(deviance statistic)ってのがカイ二乗分布に漸近的に従うからこれを利用していろいろできる。ネットに上がってるスライドだとhttps://www.researchgate.net/publication/341702223_purofairuyouduxinlaiqujian_Profile_likelihood_confidence_interval_yitengyaoer_muyaotanhuahui_2019nian3yue14riとかみたらわかりやすいかもしれん。尤度比検定とか勉強したことある人なら「あ～あれね」ってなるかも)
logistf()関数には引数としてcontrol と plcontrol があるので、ここをいじってまあとりあえず適当にイテレーションを増やしてみる。
#切片ありだと警告が出る。色々調整しなきゃらしいので調整する。 cont <- logistf.control(maxit=1000,maxstep = 100) plcont <- logistpl.control(maxit=1000) ss <- logistf(data=data2 , Y ~X1+X2+X3+X4+X5+X6 , firth = T , control = cont,plcontrol = plcont) ss
> ss
logistf(formula = Y ~ X1 + X2 + X3 + X4 + X5 + X6, data = data2,
control = cont, plcontrol = plcont, firth = T)
Model fitted by Penalized ML
Confidence intervals and p-values by Profile Likelihood
Coefficients:
(Intercept) X1 X2 X3 X4 X5 X6
0.2235775 -2.5389918 3.1783897 1.7748356 -2.6033374 0.6761207 -1.2892522
Likelihood ratio test=12.3841 on 6 df, p=0.05392857, n=20
警告が消えた。ただまあ、さっきの場合の推定値とほとんど変わらない。
せっかくなので信頼区間も見てみる
ss$ci.lower ss$ci.upper
> ss$ci.lower (Intercept) X1 X2 X3 X4 X5 X6 -8.5842939 -10.0082939 -0.2387445 -5.1621690 -8.8773211 -8.7860762 -6.6975334 > ss$ci.upper (Intercept) X1 X2 X3 X4 X5 X6 10.11215546 0.25278118 14.07936182 14.74288884 -0.05546432 10.02277278 2.32220824
X5は有意に負の効果を持つと推定されている。ちなみに真の構造ではX5はYの生起に対して非常に弱い正の影響を及ぼす。
……逆じゃねーか！（おわり）
追記：見間違えてました。X4が有意に負の効果を持つと推定されてます！真の構造でも確かに負です！よかった～
