読了:Hadeker et al. (2018) ロジスティック混合回帰モデルで得た回帰係数を集団レベルの係数に変換する方法

Hedeker, D., du Tout, S., Demirtas, H., Gibbons, R.D. (2018) A note on marginalization of regression parameters from mixed models of binary outcomes. Biometrics, 74(1), 354-361.

 ロジスティック回帰混合モデルからランダム切片を取っ払ったとき(marginalizeしたとき)、固定効果の係数をどう修正すれば良いかという解説。仕事の都合で必要になりそうな話題なので目を通した。
 第二著者の所属はScientific Software International。ここはたしかHLMの開発元だ。なんか関係あるのかな。

 いわく。
 一般化線型混合モデル(GLMM)の回帰係数はsubject-specific(SS)な係数で、ランダム効果をコントロールしたときの共変量の効果としてconditionalに解釈できる。いっぽう、GEEとかmarginalized multilevel modelでは、係数はpopulation-averaged (PA)でありmarginalに解釈できる。[← marginalized multilevel modelというのは、この論文とは別の、いったんmultilevel modelを組んでそれをmarginalizeするというアプローチのことらしい。Heagerty & Zeger(2000 Stat.Sci.)というのがreferされている]
 SSモデルは個人レベルの推論に焦点をあてており、PAモデルはpopulationレベルの推論に焦点を当てている。本論文では、縦断データ、二値アウトカム、ロジットリンクの場合に、SSモデルからPAモデルを数値的に求める方法を考えます。

 PAモデルの場合、GEEは欠損がMCARのときにのみ妥当である。MARだったらweighted GEEという手があるが、縦断アウトカムと欠損の同時モデルをつくり、欠損モデルを正しく指定しなければならない。いっぽうSSモデルなら、MARであってもfull-likelihood mixed modelが組めるので [混合モデルを完全情報最尤法で推定するっていうこと?]、いったんSSモデルを組んでおいてそれをPAモデルに変換するという手がある。もっとも分散共分散構造を正しく指定しないといけませんけどね。[←ど素人なのでこのへんの感覚がうまくつかめない。分散共分散構造を正しく指定するのって、欠損モデルを正しく指定するのと同じくらいの無理難題ではなかろうか]
 また、混合モデルは時間間隔がirregularなデータをより柔軟に扱える。分散共分散構造は連続時間変数に依存するからである[←ここはほんとに知識不足でわからないんだけど、縦断のGEEってのは時間を離散的に捉えてんですか?]。
 というわけで、いったんmultilevelモデルを組んでそれをmarginalizeするというのは、MARの下でとか時間間隔不均一のときとかにPA推定値を得るための良い代替案である。ソフトとかコードについて知りたければGriswold et al.(2013 Stat.)を見よ。[ああ、このGriswold et al.ってのを先に読んでおけばよかったのかな…]

 被験者\( i (=1,\ldots,N) \), 時点 \( j(=1,\ldots,n_j) \)の二値アウトカムを\( Y_{ij} \)とする。次のモデルを考える[原文ではlogと分数で書いてあるけど、メモではlogit関数を使って略記する。\(x, \beta, z, v\)はベクトルだけど面倒なのでそのまま書く]。$$ logit ( P(Y_{ij} = 1 | v_i) ) = x^\top_{ij} \beta^{ss} + z^\top_{ij} v_i $$ \(v_i\)は長さ\(r\)のランダム効果で、\(N(, \Sigma_V)\)に従うとする。\(\Sigma_V = TT^\top \)とコレスキー分解すると $$ logit ( P(Y_{ij} = 1 | v_i) ) = x^\top_{ij} \beta^{ss} + z^\top_{ij} T \theta_i $$ \(\theta_i\)は標準化されたランダム効果で\(N(0, I)\)に従う。

 こんどはmarginalモデルについて考えよう。$$ logit( P(Y_{ij} = 1) ) = x^\top_{ij} \beta^{pa} $$ このモデルを使う際には、縦断データにおける連関についてのなんらかの特徴付けも必要になるだろう。GEEにおけるworking correlationとか。
 このモデルからpopulation-averaged確率が得られる。$$ \pi_{ij}^{pa} = \Psi(\lambda^{pa}) = logit^{-1}(\lambda^{pa})$$ \( \lambda^{pa} \)はpopulation-averagedロジット、\( \Psi\)は標準ロジスティック分布のCDF (いわゆるexpit)である。

 さて、混合ロジスティック回帰モデルで\(\hat{\beta}^{ss}, \hat{T}\)を得たとして、そこから\(\hat{\beta}^{pa}\)を得るためにはどうしたらよいか。

 もしこれがランダム切片プロビットモデルだったら閉形式で書けてしまう。ランダム切片の分散を\( \sigma^2_v \)として $$ \hat{\beta}^{pa} = \hat{\beta}^{ss} / \sqrt{ 1 + \hat{\sigma}^2_v } $$ である。[ふーん]
 しかるにランダム切片ロジスティックモデルの場合は、$$ \hat{\beta}^{pa} \approx \hat{\beta}^{ss} / \sqrt{ \frac{\hat{\sigma}^2_v + \pi^2/3}{\pi^2/3} } $$である。標準ロジスティック分布の分散が\( \pi^2/3 \)だからこうなるのである。文句のあるやつぁ Agresti “Categorical Data Analysis” (2002)をみやがれ。[はい、そうします… 会社の本棚だから出社しないと読めないけど…]
 Zegar et al.(1988)いわく、\( \pi^2/3 \)を15/16倍したほうがよい。これはロジスティック関数の累積正規分布による近似から得られる値である。[よくわからん…]
 右辺の \( \sqrt{ \frac{\hat{\sigma}^2_v + \pi^2/3}{\pi^2/3} } \)は、conditionalなSDとmarginalなSDとの比を表している。つまり、SS係数をPA係数に変換する”marginalization”ファクターと捉えることができる。これはスカラーで時間不変である。

 ここまでは、ランダム効果が切片のみだった場合の話である。ランダム効果が複数あると、話は俄然ややこしくなる。$$ \hat{\pi}^{pa}_{ij} = \int_\theta \Psi(x^\top_{ij} \hat{\beta}^{ss} + z^\top_{ij} \hat{T}\theta_i) dF(\theta) $$ これをガウス・エルミート求積法で近似すると…[メモ省略]。
 またとにかく\( \hat{\pi}^{pa} \)が出せたとする。ときに、\( \pi^{pa} \)の個々の要素のロジットからなるベクトルを \( \lambda^{pa} \) とすると、\( \lambda^{pa} = X \beta^{pa} \)である。ここから $$ \beta^{pa} = (X^\top X)^{-1} X^\top \lambda^{pa} $$である。この式に \( \hat{\pi}^{pa} \)をロジットに変換した \( \hat{\lambda}^{pa} \)を放り込めば\( \hat{\beta}^{pa} \)が手に入る、という寸法である。
 
 [\( \hat{\beta}^{pa} \)の標準誤差はどうなるかという話。読み飛ばした]
 [ロジットリンクに限らず一般化すると… 云々。読み飛ばした]
 [シミュレーション。ランダム切片プロビットモデルについて、SSモデルを数値積分でmarginalizeした結果は閉経式でmarginalizeした結果とほぼ一致。ランダム切片と時間トレンドがはいったロジスティック混合モデルについて、SSモデルを数値積分でmarginalizeした結果、Griswold & Zegar (2004)という別のやり方で得た結果、GEEの3つを比較するとほぼ同じ、とかなんとか。読み飛ばした]
 [事例。力尽きたので丸ごと読み飛ばした]

 考察。
 本論文で提案した数値積分によるmarginalizeは簡単。GEEと変わらないし、MAR欠損の下ではGEEより良い。
 弱点は、混合モデルにおいて共分散構造を正しく指定しなければならないという点。
 本論文の提案はmarginalized multilevel modelの代替案である。marginalized multilevel modeを実装しているソフトは少ないけど、混合モデルなら簡単である。
 なお、提案手法は我々のつくったSupermixというソフトに載ってまっせ。
 云々。

 …途中から「こりゃあソフトがないと辛いわ」という気がしてきて、斜め読みになってしまったんだけど、ま、勉強になりました。
 自分なりにまとめておくと… 回帰モデルを組みました。データにクラスタ構造があるので階層回帰モデルにしました。説明変数の効果は偏回帰係数に表れておりますと。ところがそれは、ランダム係数の推定ができている個々のクラスタに属する個体についての推定結果なのであって、クラスタを超えた集団レベルについての推定結果ではないかもしれない。モデルに交互作用項や二次の項が入っていたり、リンク関数が入っていたりすると、個体レベルの効果と集団レベルの効果はずれる。なぜなら後者は、モデルの右辺のランダム係数入りの式を、(リンク関数で変換した後に)ランダム係数で積分したらどうなりますかという話だからである。
 オーケー、じゃ係数を修正しましょうと。その場合、ランダム係数が切片のみでロジットリンクないしプロビットリンクだったら手計算でもなんとかなる。そうでないときはかなーりややこしいと覚悟する必要がある。ってことであろうか。

 ところで… この前にmarginsパッケージのvignetteを読んでいて、marginalとconditionalという言葉の意味について頭が混乱してしまったんだけど、ここでの用語法は明確であった。
 ここでは、marginal-conditionalは対義語で、ある共変量の係数が、ランダム効果をコントロールしたもとでの係数だったらconditional, でなかったらmarginalである。他の共変量の係数をコントロールしているかという意味ではない(その意味ではいずれもconditionalだということになると思う)。回帰平面の傾きを指してmarginalと呼んでいるわけでもない (この論文の用語法では、おそらく Y=(ランダム切片)+bX+cX^2 というモデル のbはconditional, Y=(全員共通の切片)+bX+cX^2というモデルのbはmarginalだが、しかしどちらのモデルでもbは回帰平面の傾きではない)。嗚呼、腑に落ちる… ありがとう著者の先生たち…
 わたしゃ十分な統計学教育を受けていないので、何事に関しても自信が持てないんだけど、どうやらmarginal/conditionalという言葉には要注意だ。その人がどういう意図で使っているのかに気をつけないといけない。