仕事の話なので抽象化して書くけれど、被験者内１要因の実験計画、被験者x要因内でさらに反復測定(反復回数は一様でない)、目的変数は二値。検定したいんだけどやり方がよくわからん、どうすればいい？ …という主旨のお問い合わせを、先日受けた。うーん、それは確かに、ちょっと困るかも。少なくとも市場調査のルーチンワークからは外れている。
それはもうGLMMなんじゃないっすかね、と説明しかけて、いや待てよ、こういうときにはGEEってのもあるよな、というのが頭をよぎり、どんよりした気分になった。GEE(一般化推定方程式)、それは過去なんどか勉強しようとしては挫折した、私にとっての鬼門のひとつなのである。
Hubbard, A.E., et al. (2010) To GEE or Not to GEE: Comparing Population Avarage and Mixed Models for Estimating the Associations Between Neighborhood Risk Factors and Health. Epidemiology, 21(4), 467-474.
というわけで、易しそうな文献で再チャレンジ。
いわく。
たとえば、コミュニティの特性と健康上のアウトカムとの連関を、住民の個人特性を調整した上で決めたいというようなとき、マルチレベル・アプローチがよく使われる。
連関の推定に使われる方法は主にふたつある。その１，混合モデルないしランダム効果モデルを最尤推定。その２，母集団平均モデル。[GEEの文脈でpopulation-levelっていうとき、populationってなんと訳せばいいんだろうか？ 母集団でいいの？ なんだかしっくりこないけど、以下では母集団と訳す]
これまでの研究では混合モデルが好まれてきた。おそらくその理由は、混合モデルはアウトカムがコミュニティ内で持つ共分散構造とコミュニティ間で持つ共分散構造を明示的に分割してモデリングするからだろう。しかし、分散を分割できるというのは、分析手法決定に際しての考慮事項のひとつにすぎない。
本論文は、コミュニティレベルの曝露と個人レベルのアウトカムとの間の連関を回帰で研究する際の手法を比較する。コミュニティ研究における因果推論に関するその他の諸問題については、Oakes (2004 Soc.Sci.Med.)を見て下さい。混合モデルと母集団平均モデルについての概観は、Heagerty & Zager (2000 Stat.Sci.), Gardiner, Luo, & Roman (2009 Stat.Med.)[←これ面白そう]、などを見て下さい。
この論文の目的は、それぞれのアプローチがどういう仮定の上に成り立っているかを概観することである。一般に、回帰モデルは真実の近似と考えるのがもっとも現実的だ。混合モデルはデータ生成分布について検証不能な仮定に依存している… [ひとしきり混合モデルのディスりが続くけど、ここは論文の予告編的なくだりなのでスキップ]
混合モデル
地区\(j\) における対象者\(i\)のアウトカムを\(Y_{ij}\)とする。共変量\(X_{ij}\)、固定係数\(\beta\)のもとでの平均的反応を\(\mu(X_{ij} | \beta)\)とする[ちょっと待って、この縦棒はなんなの…? \(\mu(X_{ij}, \beta)\)でよくない？]。近隣のランダム効果を\(\alpha_j\)とする。誤差項を\(U_{ij}(\alpha_j, X_{ij})\)、ただし\(E[U_{ij}(\alpha_j, X_{ij}) | X_{ij}] = 0\)とする。$$ E(Y_{ij} | X_{ij}) = g[ \mu(X_{ij}|\beta) + U_{ij}(\alpha_j, X_{ij})] $$ と書ける。ただし\(g\)はリンク関数で、線型なら\(g^{-1}(u) = u\), 対数なら\(g^{-1}(u) = \log(u)\), ロジスティックなら\(g^{-1}(u) = \log(u/(1-u))\)である。
具体例でいこう。地区\(j\)の対象者\(i\)が週2時間以上歩いてたら\(Y_{ij}=1\)、そうでないなら\(0\)とする。\(X_{ij}\)は地区の犯罪率とする(\(j\)のなかで等しい)。次の単純なランダム効果モデルを考える[第二項の\(X_{ij}\)は原文では\(X_{1j}\)。勝手に修正した]: $$ \log \frac{P(Y_{ij} = 1 | X_{ij}, \alpha_j)}{1-P(Y_{ij} = 1 | X_{ij}, \alpha_j)} = \beta_0 + \alpha_{0j} + (\beta_1 + \alpha_{1j}) X_{ij}, \ \ \alpha_j \sim MVN(0, \Sigma)$$ ただし\(\alpha_j = (\alpha_{0j}, \alpha_{1j}) \)。ついでに\(\beta = (\beta_0, \beta_1)\)と書くことにして $$ E(Y_{ij} | X_{ij}, \alpha_j) = P(Y_{ij} = 1 | X_{ij}, \alpha_j)$$ $$\mu(X_{ij} | \beta) = \beta_0 + \beta_1 X_{ij}$$ $$U_{ij}(\alpha_j, X_{ij}) = \alpha_{0j} + \alpha_{1j} X_{ij}$$ である。
このように、地区\(j\)のある対象者の「歩く時間」変数の平均が、\(X_{ij}, \beta, \alpha_j\)の関数としてモデル化される。\(\beta\)と誤差項の分布パラメータが最尤法ないし制約付き最尤法で導出される。[…中略…] \(\beta\)の推定値についての推論(SEとか)の精度は、誤差と\(\alpha_j\)の分布を正しく指定したかどうかに依存する。
線型モデルや対数線型モデルの場合、回帰係数は、\(X_{ij}, \alpha_j\)のもとでの\(Y_{ij}\)の条件付き平均の係数とも解釈できるし、\(Y_{ij}\)と\(X_{ij}\)のpopulation average association の係数としても解釈できる[\(E(Y_{ij} | X_{ij}, \alpha_j)\)のモデルの係数でもあるし\(E(Y_{ij}|X_{ij})\)のモデルの係数でもあるってことね]。ブートストラップやサンドイッチ推定量を使えば、\(E(Y_{ij}|X_{ij})\)のモデルが正しいことだけに依存する頑健推定が可能である。いっぽうロジスティック回帰の場合、ランダム効果モデルの係数と母集団平均モデルの係数は異なる。
観察データの尤度は $$ L(Y_{ij} | X_{ij}) = \int_{\alpha_j} f(Y_{ij} | X_{ij}, \alpha_j) h(\alpha_j | X_{ij}) d \alpha_j$$ つまり、\(X_{ij}, \alpha_j\)のもとでの\(Y_{ij}\)の確率密度を\(\alpha_j\)の確率密度を通じて平均したものである。しかし、\(X_{ij}\)のもとで\(Y_{ij}\)の同じ周辺分布を提供する\(f\)と\(h\)の組み合わせは無数にありうる。つまりモデルはnonparametrically nonidentifiableである。こういうモデルを使うときは、誤差とランダム効果の同時分布の誤指定のせいで、推定と推論が過度に誤りませんように、と祈るしかない。
母集団平均モデル
母集団平均モデル(周辺モデルともいう)は、共変量の変化のもとでの母集団平均の変化について述べる。係数推定には一般化推定方程式(GEE)が使われることが多いが、最尤法も使えるので、GEEとモデルを混同しないように。
GEEアプローチは分布の仮定を必要としない。必要とするのは、(1)関心あるパラメータ(ここでは\(E(Y_{ij}|X_{ij})\)のモデルのパラメータ)についての提案と、(2)もし真のパラメータをいれたら0になるような関数を見つけること、だけである。
線型モデルの場合、推定方程式アプローチは特定の混合モデルと実務的にみて同じ推定量を提供することが多い。たとえば、単純なランダム切片線型モデルは、すべての事例で地域内の分散共分散が等しく、地域をまたいだ事例間には相関がないと仮定する。推定値はGEEでの交換可能作業相関モデルと等しくなる。このように、特定の状況下では、GEEと混合モデルの推定値は同じである。
しかし\(\beta\)の推定値についての推論を行うときには話が異なる。推定方程式アプローチでは尤度を指定しないので、推定値についての最尤推論はできない。かわりに頑健推論かサンドイッチ推論を使う。これらの推定量は漸近的に線型である(推定量は漸近的に、iidな確率変数の合計として書ける。これをinfluence curveという)。influence curveを閉形式で書ければ、データの背後にある分布について仮定することなく、\(\hat{\beta}\)についての頑健推論が可能になる。ノンパラメトリック・ブートストラップという手もある。地区を復元抽出し、オリジナルのデータと同じサイズの擬似母集団をつくり、オリジナルのデータでやったのと同じ推定をやる。これを繰り返す。いずれにせよ、地区の数が十分多いことが必要。
パラメータの解釈: 混合モデル vs. 母集団平均モデル
ここからは、上で挙げた犯罪率と歩行の関係を例にして、パラメータ解釈上のちがいを示そう。モデルはもっと簡単に、単純なランダム切片モデルにして…[原文で\(\sigma\)とあるところを\(\sigma^2\)に修正] $$ logit ( P(Y_{ij} = 1 | X_{ij}, \alpha_j) ) = \beta_0 + \alpha_{0j} + \beta_1 X_{ij}, \ \ \alpha_{0j} \sim N(0, \sigma^2)$$ \(\beta_1\)は、地区の効果を固定したときの、犯罪率の1単位増大に伴う歩行の対数オッズ比である。
次に、母集団平均のロジスティック回帰モデルを考えよう(このモデルが上のモデルから含意されるわけではない点に注意): $$ logit( P(Y_{ij} = 1 | X_{ij})) = \beta^{*}_0 + \beta^{*}_1 X_{ij}$$ \(\beta^{*}_1\)は、犯罪率が1単位異なる地区で暮らす人の間での歩行の確率を比べている。
ロジスティック回帰の場合、\(\beta_0\)と\(\beta^{*}_0\)は理論的には異なる。なお線型モデル・対数線型モデルの場合は、数値的に同じ事である。なぜなら$$ E(Y_{ij}) = E_{\alpha_j}[E(Y_{ij} | X_{ij}, \alpha_j)] = \mu(X_{ij} | \beta) + E_{\alpha_j}[U_{ij}(\alpha_j, X_{ij})] = \mu(X_{ij} | \beta) $$ だからだ。
ロジスティック回帰に話を戻すと、真のモデルが単純なランダム切片モデルであるとき、母集団平均モデルでのオッズ比は混合モデルでのオッズ比に比べて1に近づく。もっと具体的にいうと、\(X_{ij}\)が犯罪率の高低を表す二値だとして、ランダム切片ロジスティック回帰モデルの傾きは、ある地区の犯罪率が高いときの歩行と低い時の歩行をを比べた対数オッズ比を地区内で平均したものである。いっぽう母集団平均ロジスティック回帰モデルの傾きは、犯罪率が高い地区と低い地区の間で歩行の確率を比べた対数オッズである。
典型的なクロス・セクショナルな地区レベルデータ(\(X_{ij}\)がある\(j\)のなかで共通)の場合、ランダム効果モデルにおける連関の推定値は、異なる地区を通じた比較から得るしかない(地区内の対数オッズは直接には推定できないから)。従って、傾きの推定値はランダム効果の分布についての仮定に依存する。たとえばランダム効果の正規性の仮定である。こうしたランダム効果モデルから得られる地区間変動の量の推定は、これらの仮定があまり破られていないという仮定に依存する。
projectionとしての回帰モデル
そもそも分析とは、(1)関心ある科学的問いを提案し、(2)その問いに関する関心あるパラメータを決め、(3)推定量を定め(ここで経験的には識別不能な仮定が必要になるかもしれない)、(4)推定値と適切な注意書きを報告する、というものである。
混合モデルのような潜在変数モデルは常に、データからは検証できない仮定を必要とする。尤度ベースの推論が信頼できる推論であるためには、モデルは近似であってはだめで、共変量と潜在変数のもとでの平均モデルの真の形式でなければならない。モデル選択を支える十分な理論がないとき、推定された回帰は、せいぜい条件付き平均の近似にすぎない。固定効果の部分が正しくモデル化されていないときはもちろん、潜在変数の部分が正しくモデル化されていないときも、固定効果の推定値をどう解釈したら良いのかわからなくなる。
いっぽう、真の母集団平均モデルの近似を推定することを正当化することは容易である。このアプローチなら、モデルの特定は不要である。[…中略…]
たとえば、単純な反復無し測定構造を考えよう。真のモデルを$$ Y = b_0 + b_1 X + b_2 X^2 + e, \ \ e \sim N(0, \sigma^2)$$とする。ここから得られたデータに線型モデルをあてはめてしまったとする。さらに、関心あるパラメータは、母集団全体に$$ Y = c_0 + c_1 X + \varepsilon$$をあてはめた場合の\(c_0, c_1\)であるとする。つまり、(1)真の平均モデル、(2)線型モデルの真実へのprojection, (3)真のモデルから得られた現実のデータ、(4)データにあてはめた線型モデル、は異なる。
仮に、関心あるパラメータが\(Y\)の真の条件つき平均であると考える人がいたら、その人はほぼ常に誤る。いっぽう、関心あるパラメータはなんらかの小さなモデル集合への真実のprojectionだと考える人であれば、合理的な推定値を得られる事が多い。
真のモデルが次の混合モデルである場合を考えよう。$$ logit(P(Y_{ij} = 1 | X_{ij}, b_{0i})) = b_0 + b_1 X_{ij} + b_{0i} \exp(b_2 X_{ij}), \ \ b_{0i} \sim N(0,1)$$ 次の5つは異なる。(1)アウトカムの真の周辺確率\(P(Y_{ij}=1 | X_{ij})。(2)単純なロジット線型モデルと独立作業相関によるモデルをあてはめて得た近似。(3)真のモデルから生成されたデータを用いた、GEEによる近似。(4)同じデータを用いた、単純なランダム切片を持つロジスティック回帰モデル。(5)ランダム効果モデルを、推定されたランダム効果の分布を通じて周辺化して得た母集団平均推定値。ご覧下さい、GEEのほうがいいでしょう。[と図を示している。いや、全部けっこうずれてますやんか…]
誤指定されたランダム効果モデルから得られる母集団平均の推定値は、もちろんバイアスを受けている。これは、パラメータを推定するためにはランダム効果の部分を正しく指定する必要はないわけで、その意味で不必要なバイアスである。もしランダム効果の部分を誤指定したら、地区特定的な効果(たとえばランダム効果モデルにおける係数)も、そこから含意される母集団平均モデルも、推定の対象の量と比べて予測不能なバイアスを持ちうる。条件つきで指定されたモデルにおける回帰係数(ランダム効果モデルにおける固定効果)は、それに対応する母集団平均モデルに比べて、ランダム効果についての仮定に対してはるかに敏感である。
[…後略]
考察
[くたびれたのでメモ省略]
—-
易しい解説かと思ったら、これが意外にそうでもなかった…
ときどき、数学が得意な人があえて数式を使わずに説明しようとしてかえって曖昧模糊となってしまうことがあるけど、これもそれだと思う。我々のように数学が苦手な馬鹿はですね、単にフォーマルな記述を減らせば喜ぶわけじゃなくて、フォーマルな記述とインフォーマルな説明との間を丁寧に対応づけて欲しいと思っているのに、そのことがいまいちわかって頂けないのである。いやー、涙ちょちょぎれますね。
うううむ。勉強不足のせいだろうけど、謎が深まるばかりだ。著者らのGEE推しに反する感想で申し訳ないんだけど…
知りたい事柄が、混合モデルの固定効果パラメータとしても母集団平均モデルのパラメータとしても表現できるとき、母集団平均モデルのほうが(余計な仮定がないぶんだけ)有利だ、というのはわかる。でも、知りたい事柄をどちらのモデルでも表現できるというのは果たしてどんな状況なのか、うまく想像できないんですよね。
冒頭に挙げた、1要因反復測定デザインの実験で目的変数が二値という事例では、要因の水準\(i\), 被験者\(j\)の\(k\)回目の測定を\(Y_{ijk}\)として、混合モデルならば $$ g^{-1}(P(Y_{ijk} = 1 | \alpha_j)) = \beta_0 + \alpha_j + A_i$$ 母集団平均モデルなら $$ g^{-1} ( P(Y_{ijk} = 1)) = \beta^{*}_0 + A^{*}_i $$ だと思う(\(A_i, A^{*}_i\)は要因の主効果で\(\sum_{i} A_i = \sum_{i} A^{*}_i = 0\)というつもり)。実験者が関心を持つのは、任意の\(j\)さんに対する介入がもたらす平均的な因果効果を表すはずの\(A_i\)であって、\(A^{*}_i\)じゃないと思うんですよね… こういうときに母集団平均モデルを使うのはおかしい、なんらか\(\alpha_j\)の分布について仮定した上で混合モデルを推定せざるをえない、という理解で正しいだろうか？
本文中の地区犯罪率\(X_{ij}\)と歩行\(Y_{ij}\)の例でいうと、人が近所を散歩するかどうかには犯罪率の他にも、坂道の多さとか遊歩道が整備されているかとか、そういう地理的な特性も影響していることがわかっており、しかしそれらは測定できていない、とする。仮に、犯罪率と歩行確率との関係をそうした未知の特性を調整した上で検討したいならば(そりゃ調整したいですよね。だってそれらは犯罪率と歩行確率の両方の共変量かもしれないから)、その場合は、万難を乗り越えて\(E(Y_{ij} | \alpha_j, X_{ij})\)のモデルを推定すべきであって、\(E(Y_{ij} | X_{ij})\)のモデルをGEEで推定するのはおかしい… という理解で合っているだろうか？
こうしてみると、知りたい事柄を混合モデルの固定効果パラメータとしても母集団平均モデルのパラメータとしても表現できるような状況ってのがいったいどんな状況なのか、だんだんわかんなくなってきた… もうちょっときちんと勉強しないと、よくわからんなあ。