2

Rでロジスティック回帰を行った際のresidual devianceを取得したいと考えています。
glm.fit関数を用いたところ、以下のような完全分離の警告が出ました。

glm.fit: fitted probabilities numerically 0 or 1 occurred 

そこで、こちらを参考にlogistf関数を利用してみました。

cont <- logistf.control(maxit=1000,maxstep = 100)
plcont <- logistpl.control(maxit=1000)
res <- logistf(y ~ x,firth = TRUE, control = cont,plcontrol = plcont)

summary(res)の中身は以下のとおりです。

> summary(res)
logistf(formula = dat[, j] ~ Dat_res[, i], control = cont, plcontrol = plcont, firth = TRUE)

Model fitted by Penalized ML
Coefficients:
                   coef   se(coef)  lower 0.95 upper 0.95    Chisq        p
(Intercept)   0.1375199 0.11492692 -0.08696857  0.3642894 1.439322 0.230249
Dat_res[, i] -0.7204668 0.05044381 -0.82565105 -0.6269034      Inf 0.000000
             method
(Intercept)       2
Dat_res[, i]      2

Method: 1-Wald, 2-Profile penalized log-likelihood, 3-None

Likelihood ratio test=895.6586 on 1 df, p=0, n=1000
Wald test = 203.9916 on 1 df, p = 0

この結果から、こちらのサイトでも説明されているdeviance残差(residual deviance)を計算しようと考えていますが、導出する方法がよく分かっていません。
理解不足で申し訳ありませんが、ご教授よろしくお願いします。

1 件の回答 1

1

一般化線形モデル (GLM) の deviance (逸脱度) の定義は「-2 x モデルの対数尤度」です (別の定義もありますがそのスライドや質問の文脈からこれを意図していると判断します), この場合, 通常 residual deviance は推定したモデルと飽和モデル (フルモデル) の逸脱度の差のことを指します.

この場合, logisf のオブジェクトの loglik には罰則推定モデルとフルモデルの対数尤度が含まれているので, -2 * (res$loglik[1] - res$loglik[2]) で計算できることになります.

そしてこれはデフォルト設定での尤度比検定の検定統計量と等価なので, summary() で表示したときの Likelihood ratio test=... の値と同じになるはずです.

おせっかいかも知れませんが, ネット上の記述は多くの場合断片的なので, 体系的に説明している教科書を1冊は読んでおいたほうがいいと思います. 例えば日本語で, かつ逸脱度について言及がある本の1つに, 久保拓哉『データ解析のための統計モデリング入門』 (2012, 岩波書店) があります.

ただし, 免責事項というか, 以下の点に注意してください.

  1. 私は logistf パッケージを今回初めて知ったので, 適切に動作するかは確認していません (説明書通りの方法を提示しただけです)
  2. そちらが何を目的にしているのかはこちらにはわからないため, このアプローチが適切なのかまでは保証できません (i.e., 罰則付き推定と罰則なし推定は別物ですし, 完全分離がどのような理由で発生しているかは私には知りえません)

回答から時間が経ってしまいましたが, 誤解を招きそうなので以下追記

質問文とコメントを読み返して気づきましたが, 私もあなたも residual deviance と deviance residuals を混同して書いています. 私の方は「deviance残差」と表記してる箇所もありますが, 意図していたのは全て residual deviance です. モデルに対して定義されるのがこの residual deviance で, -2*(対数尤度 - フルモデルの対数尤度) と定義され, これは上記通りの方法で数値を得られます.

しかし, 参考にされたスライドでは, 上記の residual deviance を単に "deviance" と呼んでいます. deviance を -2*(対数尤度) とするか -2*(対数尤度 - フルモデルの対数尤度) と定義するかは文献によってまちまちなので紛らわしいです. R の glm のヘルプや上記の参考書では, 参考になさったスライドとは異なり, 前者を deviance とし, 後者を residual deviance と呼んでいます. (以降はこの区別のため, 元スライドと違う呼び方になってることに注意してください).

一方で観測ごとに計算する deviance residuals / deviance 残差 は, 上記の (residual) deviance の平方根に符号をつけたものを指します. deviance = -2*(対数尤度) と residula deviance = -2*(対数尤度 - フルモデルの対数尤度) のどちらを使うかもまた, 文献によってまちまちなので紛らわしいです. (参照したスライドも式を明示していないので私も読み違えていたようです) つまり,

  1. deviance の定義には複数あり, しかもどれが主流か決まってない
  2. (モデルに対する) residual deviance と (観測点に対する) deviance residuals 名前が似ているが全く別物である

という二重の紛らわしい箇所があることに注意してください (後者は読み違えないように "unit --" とか "casewise ---" とかを頭につけて区別することがあります).

そちらはタイトルでは「deviance残差」と書いて本文では「residual deviance」「deviance残差(residual deviance)」などと書いています. 私は, 表記にゆらぎがあるものの, モデルに対する residual deviance を意図していると判断してして回答しましたが, 実際に欲しいのが residual deviance ではなくこちらの (casewise な) deviance residuals であるなら, glm では residual(<glm オブジェクト>) で得られます. ただしスライドの定義とは異なり, R では -2*(対数尤度) を deviance と定義しており, この関数で計算されるものもこれに準じています (GLM で使われる residuals の定義は他にもいくつもありますが, どれを使うべきかは目的によると思います). 一方で logistf には residuals を計算する機能はありません. よって観測点ごとの対数尤度か事後対数確率から計算することになります. 例えば以下に計算例が挙げられています. これは尤度を基にした場合で, かつ deviance は -2*(対数尤度) と定義しているようです.

https://github.com/georgheinze/logistf/issues/37

なお, コメントに書いたように, 線型回帰モデルの残差とGLMのdeviance 残差は必ずしも同じふるまいをするわけではないため, 常に残差と同じように扱えると思い込まないでほしい, という点は変わりません.

7
  • コメントが遅くなりました、すみません。 なるほど、確かに尤度比検定の検定統計量として求められる値と一致しますね。 ただ、その値は、こちらのサイトの34枚目で述べられているdevianceに相当すると思いますが、その場合、このdeviance残差に相当する値(それぞれのデータにおける期待値と実測値の差の平方根に相当)はどのようにして求められるのでしょうか。 線形回帰ですと、実測値と予測値の差というものが感覚的にも分かりやすいですが、一般化線形モデルにおいてはイメージがしづらく困っておりました。 体系的な学習についてもご教授いただきありがとうございます。 そちらの書籍について、一通り学習をさせていただきました。 また、免責事項として述べて頂いた完全分離の理由ですが、人工的なデータ生成の際のデータの偏りによるものです。そのため、logistf関数の適用以外にもそもそものデータ生成の方法についても見直してみようかと考えております。
    – katht
    2021年9月20日 0:39
  • こちらも返事が遅れました: まず, ルートを求めたいのであれば R には sqrt 関数があります. 次に34ページの記述について, そのスライドの30ページに書かれているのが (数式ではなく数値で書かれていますが), 一般的な deviance 残差の定義だと思います. 34ページは「各々差を取って, 2倍して」という記述までがこれと対応していると思いますが, 「平方根を取ったもの」というのは deviance 残差の定義とは異なります. 同ページ上部の式とも一致しないので, スライドを見ただけではどういう意図なのかなんとも言えません. (もしかすると後述する unit deviance residuals と混同しているのかもしれません) 2021年9月24日 12:06
  • 最後に, GLMのdeviance残差がイメージしづらいという気持ちはわかりますが, 例えばロジスティック回帰だとyがゼロか1しかないため, 通常の回帰モデルのようにデータの点ごとの残差を計算してもうまくいきません. GLM でも線形回帰で計算できるようなデータの点個別の残差に対応するものとして unit deviance residuals というのが定義されており, R でも glm 関数ならば residuals 関数で求められますが, GLMは非線形モデルなので単純な差とみなすことはできず, 見方には注意が必要です (定義も複数存在します). 一方で, 線形回帰モデルの残差平方和は対数尤度と一致するので, モデル全体の尤度で評価するのはこの点では普遍的な方法だと思います. 2021年9月24日 12:06
  • ご連絡が大変遅くなり,申し訳ありません. 仰る通り,求めたいものがcasewiseなdeviance residualsのことでした. 私自身,似た用語が多く混乱していましたが,その背景には,整理して頂いたように二重の紛らわしい点があったためでした. 線形回帰モデルの残差平方和が対数尤度に一致するため,モデル全体の尤度もしくはそこから計算される(casewiseでない)residual devianceを用いることが望ましいかと思いますが,別の手法を用いるために,個々のデータ点における数値を得る必要があったため,このcasewiseなdeviance residualsを計算する方法が分かり大変助かりました. 拙い文章にもかかわらず,何度も重ねて丁寧にご説明いただきありがとうございます.また問題が起きる可能性もありますが,一度試してみます.
    – katht
    2021年10月20日 7:12
  • こちらこそ, 前提の認識共有に時間がかかってしまいすみません. 関心があるのが casewise な残差となると, (deviance 残差に限らず) GLM での残差診断全般の解説資料を紹介したほうがよい気がします (そちらが参考になさったスライドも, 私が紹介した教科書も残差診断にはあまり詳しくないです). GLM の残差診断の解説は, 英語なら有用なものが結構ありますが, 日本語だと教科書であってもかなり少ないです. よって, 英語になってしまうのですが, 数種類の残差の定義とその違いを紹介したものとして, こんなページがRコードつきで良いと思います: stats.stackexchange.com/questions/1432/… 線型回帰モデルの残差の延長として使いたいならば, deviance 残差ではなく pearson 残差のほうが良いかもしれませんが, 使用するデータによっては0と1ばかりになり使いにくいかもしれません. (このページ内でも Harrel という統計学者がコメントしているように, 一般にロジスティック回帰での残差は得られる情報が少ないという問題があります) 2021年10月20日 8:20

この質問に回答するには、ログインする必要があります。

求めていた回答ではありませんか? のタグが付いた他の質問を参照する。