Pavlou, M., Ambler, G., Seaman, S., Omar, R. (2015) A note on obtaining correct marginal predictions from a random intercepts model for binary outcomes. BMC Medical Research Methodology, 15:59.
たまたまどこかで題名をみかけて、そうそうこれ前から疑問に思ってた話だよなと気づき、移動中にざざーっと目を通した。
いわく。
階層データの分析にはランダム効果モデル(混合モデル)とGEEがよく用いられる。前者は特定のクラスタのもとでの条件付き推論を、後者はpopulationを通じて平均した周辺推論を提供する。線形モデルの場合、populationレベルの回帰係数は両者で一致する。しかしアウトカムが二値の場合は必ずしも一致しない。
Bouwmeester et al.(2013)は、ランダム切片付きのロジスティック混合モデルで条件付きリスクを予測した。彼らがいうには、分析データに含まれているクラスタに属する新メンバーでイベントが生じる確率は、周辺予測よりも条件付き予測のほうがうまくいった。しかしよく見ると、彼らのいう周辺予測というのはランダム切片をゼロにしただけである。これは実は正しくないやり方で、級内相関(ICC)が小さいときに損失が生じる。
この論文では、なぜ損失が生じるのか、正しい周辺予測を得るにはどうしたらよいか、を説明しよう。
クラスタ\(i (=1, \ldots, K)\)のメンバー\(j (=1, \ldots, N_i)\)の二値アウトカムを\(Y_{ij}\), 共変量の\(p\)次元ベクトルを\(X_{ij}\)とする。
ランダム切片モデルは $$logit(P(Y_{ij}=1|X_{ij}=x_{ij}, u_i) = \alpha^{RE} + u_i + \sum_m^p x_{ijm} \beta^{RE}_m $$ \(u_i\)はクラスタのランダム切片で、ふつう\(u_i \sim N(0, \sigma^2_u)\)と仮定する。[原文では\(\alpha, \beta\)のREは下添字だが、インデクスの添字と変数名の添字が混在してややこしいので上添字にする。次の式の添字Mも同様。\(X_{ij}, x_{ij}, \beta^{RE}\)はベクトルなのでほんとは太字だけど、面倒なので太くせずに書く]
一方周辺モデルは $$ logit(P(Y_{ij} = 1 | X_{ij} = x_{ij})) = \alpha^M + \sum_m^p \beta^M_m $$ このモデルはクラスタを無視してML推定できる。推定量を\( \hat{\alpha}^{ML}, \hat{\beta}^{ML} \)とする。[原文では添字はMLでなくてstandardなんだけど、長くて書きにくいので変更した]
しかし理想的には、クラスタ内相関を説明する必要があるので、上の式は適切な’working’相関行列を使ってGEE推定するべきである。推定量を\( \hat{\alpha}^{GEE}, \hat{\beta}^{GEE} \)とする。
\( \hat{\alpha}^{ML}, \hat{\beta}^{ML} \)と\( \hat{\alpha}^{GEE}, \hat{\beta}^{GEE} \)はいずれも\( \alpha^M, \beta^M \)の一致推定量である(欠損がMCARであれば)。またworking相関行列として独立な行列を使えばML推定量とGEE推定量は一致する。
さてここで重要なのは、ランダム切片モデルの固定効果は、周辺モデルの固定効果より絶対値が大きい、という点である。つまり $$ |\alpha^{RE}| \geq |\alpha^M|$$であり、すべての要素\(m\)について $$|\beta^{RE}_m| \geq |\beta^M_m|$$である。その差はランダム効果の分散が大きいほど大きくなる。Zegar et al. (1988 Biometrics)によれば、$$ \alpha^M \approx \frac{\alpha^{RE}}{\sqrt{c^2 \sigma^2_u + 1}}, \ \beta^M_m \approx \frac{\beta^{RE}_m}{\sqrt{c^2 \sigma^2_u + 1}} $$ $$ c = \frac{16\sqrt{3}}{15 \pi} $$という性質がある。[へえええ！ なお原文では\(\alpha\)に下添字\(m\)がついてんだけど間違いじゃないかと思う]
いっぽうBauwmeesterさんたちの周辺予測は $$ P(Y_{ij} | X_{ij} = x_{ij}) = logit^{-1} (\hat{\alpha}^{RE} + \sum_m^p x_{ijm} \hat{\beta}_m^{RE} ) $$ となっておるわけです。しかしほんとは、\(u\)の密度関数を\(f(u)\)として $$ P(Y_{ij} | X_{ij} = x_{ij}) = \int logit^{-1} (\hat{\alpha}^{RE} + u + \sum_m^p x_{ijm} \hat{\beta}_m^{RE} )f(u) du $$ としないといけない。[あー、なるほど、納得… でもどうせ\(f(u)\)は\(N(0,\sigma^2_u) \)だから、\(\sigma^2_u\)が十分に小さかったら気にしなくていいような気もするけど…]
以下のシミュレーションでは次の5つを比較する。
- No clustering. クラスタを無視して周辺モデルをML推定。
- RE-zero. ランダム切片モデル、ランダム効果を0とする。
- Re-integ. ランダム切片モデル、ランダム効果の推定された分布を積分する。
- Re-approx. ランダム切片モデル、推定された条件付き係数をZegerらのの近似でリスケールする。[うん、なんだかそれでよさそうに聞こえますね]
- GEE. exchangebleな相関行列を使ったGEEであてはめた周辺モデル。
ではシミュレーションしてみまーす。
母集団はクラスタ数100、それぞれのサイズは\( N_i \sim Poisson(\exp(\lambda_i)), \lambda_i \sim N(5.7, 0.3^2) \)として生成した。全部で32502人、クラスタ別人数の中央値は312となりました。
予測子は6つ。うち3つは連続量で正規分布, 平均0, SDは順に1, 0.3, 0.2。残りの3つは二値で、生起率は0.2, 0.3, 0.4。回帰係数はすべて1。[あ、説明変数は互いに独立なのね…]
さて懸案のクラスタリング。\(u\)の分散\(\sigma^2_u\)を動かす。\(ICC = \frac{\sigma^2_u}{\pi^2 / 3 + \sigma^2_u} \)となる。
というわけで、\(\eta_{ij} = u_i + x_{ij}\)ができます[\(x_{ij}\)は6次元ベクトルね]。\(P(Y_{ij} = 1) = logit^{-1} (\eta_{ij}) \)として\(Y_{ij}\)を生成しました。
こっから2段階抽出します。まず20クラスタ抽出する。それらからランダムに1000人選ぶ。これを訓練データセットとする。残りをテストデータセットとする[残りって? 選ばれなかった80クラスタのことだよね？ 選ばれた20クラスタのなかの選ばれなかった人は入らないんだよね？]。これを100回繰り返しました。
[ここの段落は全訳。Calibration in the large/Calibration slopeって医学系の言い回しなんですかね？]
周辺リスクを得るための異なる諸手法の予測性能をcalibrationの観点から検討する。テストデータに2つの異なる回帰モデルを適合させることによって、Calibration in the largeとcalibration slopeを算出した。Calibration slopeとは、推定された線形予測子(予測確率の逆ロジット関数。リスクスコア、予測されたログオッズともいう)を二値アウトカムに回帰するロジスティック回帰モデルの傾きである。Calibration in the largeとは、線形予測子セットをオフセット項として持つ切片のみのロジスティック回帰モデルにおける切片項である。完全にcalibrateされたモデルは、calibration in the largeの値が0となり(予測確率の平均がアウトカムのprevalenceと等しいことを意味する)、calibration slopeの値が1になるはずである。calibration slopeが1より小さいということは、予測の範囲が広すぎ、なんらかの縮小が必要だということを示す。
[ううう。慣れない話でよくわからん。えーと、モデルを学習して\( \hat{u}_i, \hat{\beta} \) が得られましたよね？ で、テストデータの個体について、どうやるのは別にして、ログオッズ\( \hat{\eta}_{ij} \)が得られますよね？ たとえばRE-zero方式なら \(\hat{\eta}_{ij} = 0 + (\hat{\beta}^{RE})^{\top} x_{ij}\)っていう風に。で、\(P(Y_{ij} = 1) = logit^{-1}(a + \hat{\eta}_{ij}) \)というモデルと\(P(Y_{ij} = 1) = logit^{-1}(b \hat{\eta}_{ij}) \)というモデルを別々にあてはめて、\(\hat{a}\)が0に近いか、\(\hat{b}\)が1に近いかを調べる、ってこと？ なぜ\(P(Y_{ij} = 1) = logit^{-1}(a + b \hat{\eta}_{ij}) \)をあてはめないんだろうか？ それにさ、素人考えだけど、単純に\(\hat{Y}_{ij}\)と\(Y_{ij}\)のクロス表をみればよくない？ なにか話のポイントを理解し損ねてるんだろうか…]
結果。
RE-zero はcalibrationのinterceptが大きくslopeが小さい(ICCが大きいと特に)。これはoverfittingのせいではなくて、アプローチが間違ってるせいであろう。
いっぽうRE-integ, RE-approxはNo ClusteringやGEEと比べても成績が良い。[←へー！]
X2以降を0に固定し、X1を横軸、予測された周辺確率を縦にとった曲線をみるとシグモイド曲線になるんだけど、RE-zeroだけ勾配がきつく、あとはほぼ同じ。[なるほど… つまりRE-zeroはいつも確率を過大推定するわけじゃなくて、予測子の効果を過大に捉えちゃうってことね]
考察。
学習時に未知だったクラスタに属する個人については周辺予測に関心がもたれるし、学習時に既知だったクラスタに属する新しい個人については条件付き予測に関心がもたれるだろう。条件付き予測についての分析結果はAppendixをみよ。
ランダム切片モデルは予測子がランダム効果から独立だと仮定している。これが破られていると係数にバイアスがかかるんだけど、予測という観点からはむしろ利点がある(クラスタ間のprevalenceのちがいを部分的に説明できるから。
同様に、クラスタサイズもランダム効果から独立だと仮定されているが、これが破られているとき(情報的クラスタサイズ)はクラスタサイズを予測子にいれるとよい。
結論。
現実問題として、病院と患者のデータではICCはふつう0.15以下で、そういうときにはまちがった周辺予測法でもあんまし気にしなくてよい。
云々。
… 細かいところがちょっと腑に落ちないんだけど(おそらく私の知識不足のせいで)、でも勉強になりましたです。
階層回帰モデルで、学習時にはなかったクラスタに属する個体について予測するときどうすればいいの？ というのは前から疑問だったので、その点について頭がクリアになったのが収穫であった。要するにあれですか、線形モデルならランダム切片にゼロを突っ込んでそのまま予測していいけど、一般化線形モデルの場合は回帰係数に修正が必要だってことでしょうか。
正直なところ、学習時になかったクラスタについて予測したいときはGEEを使うべきなのだろうか… でもあれ難しいので関わり合いを持ちたくない… というジレンマがあったんだけど、少なくともこの論文のセッティングについていえば、わざわざGEEに乗り換える必要はなく、階層ロジスティック回帰モデルで推定した係数をきちんと修正すればよい模様である。ありがたい。
なお、シミュレーションはlme4::glmer()でやったそうだ。Re-integ方式は自力でコードを書いたそうで、Appendixに載っている由。