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の開発元だ。なんか関係あるのかな。
1. イントロダクション
一般化線型混合モデル(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.ってのを先に読んでおけばよかったのかな…]
2. 縦断二値アウトカムの回帰モデル
被験者\( 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(0, \Sigma_V)\)に従うとする。観察個数は\( n= \sum^N n_i\)である。
\(\Sigma_V = TT^\top \)とコレスキー分解すると $$ logit ( P(Y_{ij} = 1 | \theta_i) ) = x^\top_{ij} \beta^{ss} + z^\top_{ij} T \theta_i $$ とも書ける。\(\theta_i\)は標準化されたランダム効果で\(N(0, I)\)に従う。
\(\beta^{ss}\)はランダム効果をコントールしたときの共変量の効果を表しており、”subject-specific”に解釈される。その意味はランダム効果が何かによって変わってくる。
こんどはmarginalモデルについて考えよう。$$ logit( P(Y_{ij} = 1) ) = x^\top_{ij} \beta^{pa} $$ marginalモデルでは平均構造と連関構造が分離される。上のモデルは平均構造のモデルであり、このモデルを縦断データについて使う場合には、連関についてのなんらかの特徴付けも必要になるだろう。GEEにおけるworking correlationとか。
\(\beta^{pa}\)は”population-averaged”に解釈される。ここからpopulation-averaged確率が得られる。$$ \pi_{ij}^{pa} = \Psi(\lambda^{pa}) = logit^{-1}(\lambda^{pa})$$ \( \lambda^{pa} \)はpopulation-averagedロジット、\( \Psi\)は標準ロジスティック分布のCDF (いわゆるexpit)である。
3. Subject-specificな回帰係数の周辺化
さて、混合ロジスティック回帰モデルで\(\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 \)だからこうなるのである。
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) $$ これをガウス・エルミート求積法で近似すると、積分の各次元における求積点を\(Q\)として$$ \hat{\pi}^{pa}_{ij} \approx \sum_{q=1}^{Q^r} \Psi(x^\top_{ij} \hat{\beta}^{ss} + z^\top_{ij} \hat{T} B_q) A(B_q) $$ と書ける。単変量の標準正規分布の場合の最適点\(B_q\)と重み\(A(B_q)\)についてはStroud & Sechrest(1966)をみよ。多変量密度の場合は\(B_q\)はベクトルとなり\(A(B_q) = \prod_{h=1}^r A(B_{qh})\)となる。[恥ずかしながら、このあたり、ぜんぜんわからん…]
まあとにかく\( \hat{\pi}^{pa}_{ij} \)が出せたとする。縦に積み上げ長さ\(n\)のベクトルにして\( \hat{\pi}^{pa} \)と書く。その個々の要素のロジットからなるベクトルを \( \lambda^{pa} \) とする。共変量も\(n \times p\)行列にして\(X\)とすると \( \lambda^{pa} = X \beta^{pa} \)である。ここから $$ \beta^{pa} = (X^\top X)^{-1} X^\top \lambda^{pa} $$である。この式に \( \hat{\pi}^{pa} \)をロジットに変換した \( \hat{\lambda}^{pa} \)を放り込めば\( \hat{\beta}^{pa} \)が手に入る、という寸法である。
[つまるところ、個々の\(i,j\)についてGLMMで求めた\(\hat{\pi}^{ss}_{ij}\)を、\(\theta\)で積分して\(\hat{\pi}^{pa}_{ij}\)としたとき、うまいこと\(logit(\hat{\pi}^{pa}_{ij}) = x^\top_{ij} \beta\)となるような\(\beta\)を考えなさい、それが\( \hat{\beta}^{pa} \)だ、ってことね]
[\( \hat{\beta}^{pa} \)のSEについて。途中からややこしくなるのでパス]
[ロジットリンクに限らず一般化すると… 云々。読み飛ばした]
[シミュレーション。ランダム切片プロビットモデルについて、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という言葉には要注意だ。その人がどういう意図で使っているのかに気をつけないといけない。