Josephy, H., Loeys, Tom, Rosseel, Y. (2016) A Review of R-packages for Random-Intercept Probit Regression in Small Clusters. Frontiers in Applied Mathematics and Statistics. 13.
題名の通り、一般化線型混合モデル、アウトカムは二値、リンクはプロビット、ランダム切片付き、クラスタサイズはめっちゃ小さい、という際のRに使えるパッケージを比較しましたという論文。
正直なところ、そういう局面になったらそのとき悩めばいいわけで、別に読まなくてもいいんだけど、なんだかなあこんなんで論文一本書けちゃうんだなあ(すいません)…などと呟きながら眺めていて、つい最後まで読んでしまった。だって気楽じゃないですかこういう話題。いささか心がなごむのであります。
イントロは飛ばして…
例として、Vandeweghe et al.という人のRCT研究を取り上げる。
お題は、未就学児が野菜好きになるのをどう促進するか。次の4つの手続きを比較した。(1)チコリを食べるのを推奨する[encourgementの訳。具体的になにをしたのかわからない。映画「鬼畜」で緒形拳が息子の口に毒入りのあんパンを押しつけて「喰え！利一、喰えーっ」っていう場面があったけど、ああいう感じの手続きでしょうか]、(2)チコリを食べたら報酬、(3)チコリの繰り返し曝露 [いやいや曝露って…]、そしてコントロール。プリテスト、ポストテスト、フォローアップの3時点で、野菜が好きかどうか(二値)を測定。個々のこどもがクラスタで、各クラスタが値が3つ持つ。ほんとは学校というクラスタもあるんだけど、3レベルになっちゃうので、ある学校のデータについて分析する。話を簡単にするため、推奨群とコントロール群にのみ注目。標本サイズは両群あわせて37。
モデルは次のとおり。測定時点を\(i (=0,1,2)\), 個体番号を\(j\)として $$ P(y_{ij} = 1 | x_{ij}, b_j) = \Phi(\beta_0 + \beta_1 x_{ij} + b_j) $$ \(y_{ij}\)は野菜が好きか。\(x_{ij}\)は推奨したか。関心があるのは\( H_0: \beta_1 = 0 \)の検定である。
[本来、推奨有無は時点とは独立な共変量で、かつ時点と交互作用するはずだから \(\beta_{1,i} x_{j}\)だと思うんだけど、\(\beta_{1,0}=0, \beta_{1,1}=\beta_{1,2}\)と決め打ちしてパラメータを減らし、データ\(x_{ij}\)の側は推奨群の\(x_{1j}, x_{2j}\)だけ1にしてあとはみな0にするのだと思う]
お待たせしました、選手入場です。
選手1番、GLMM。$$ E(Y_{ij} | x_{ij}, b_j) = g^{-1}(\beta_0 + \beta_1 x_{ij} + b_j), \ \ b_j \sim N(0, \tau) $$ \(g^{-1}\)は標準正規累積分布とする。書き換えると$$ P(Y_{ij} = 1| x_{ij}, b_j) = \Psi(\beta_0 + \beta_1 x_{ij} + b_j)$$ 尤度関数はランダム効果を積分して $$ l(\beta, \tau | y_{ij}) = \prod_j \int \prod_i f(y_ij|\beta, b_j)\phi(b_j|\tau) db_j$$ \(\phi\)はランダム切片の密度、ここでは正規である。
これを最大化するのは解析的には無理で、以下の方法が考えられる。
1の1、ラプラス近似。二次テーラー近似でもって周辺尤度を閉形式で書く。Rではlme4::glmer()がこれにあたる。
1の2、罰則付き擬似尤度(PQL)。反応関数をテイラー展開してLMMで近似してしまう。MASS::glmmPQL()がこれにあたる。
1の3、Adaptive Gaussian Quadrature(AGQ)という手法を使って積分記号をsummationで近似する。lme4::glmer()でできる。
選手2、ベイジアン。
2の1、MCMCglmmパッケージ。
2の2、INLA。[説明が書いてあるけど、毎度の事ながらさっぱり理解できないのでメモ省略…] R-inlaパッケージ。
選手3、SEM。測定モデルを $$ y_j = v + \Lambda \eta_j + K x_j + \eta_j $$ 構造モデルを $$\eta_j = \zeta_j $$ とする。\( v = (v \ v)^\top \)は切片ベクトル[転置記号を勝手に付け加えた]、\(\eta_j\)は因子で負荷は\(\Lambda = (1 1)^\top\)。\(x_j = (x_{1j} \ x_{2j})^\top\), \(K\)は2×2の対角行列で要素は\(k, k\)。ここでは目的変数が二値なので、いったん連続潜在変数\(y^{*}\)を置く。\(\eta_j\)の分散を1に固定するか(thetaパラメータ化)、\(\zeta\)の分散と\(\eta\)の分散の和を1に固定する(deltaパラメータ化)。でもって、\(y^{*}\)が閾値\(\kappa\)を超えたら\(y=1\)とする。
推定方法は大きく分けてMLとWLSがある。ここでは小サンプルでもうまくいくといわれているDWLSを使う。ソフトはlavaanパッケージ。チェックのためにMplusも試す。
実データで実験。クラスタサイズが2(フォローアップを捨てる)のときと3のときを比べた。2のときは手法間で推定値のばらつきが大きくなった。時間はラプラス近似が最速、MCMCは他の10倍くらいかかった。
よろしい、ではシミュレーションだ。
[なんだかどうでもよい話のような気がしてきたので(すいません)、一気に考察までスキップ]
考察。
クラスタ数2 (つまりペアード・データ)の場合、SEMとAGQの成績がよかった[←へー]。ラプラス近似はあんまりよくなかった(正規分布か大クラスタ用の手法なのでこれは予期した結果)。PQLもあまりよくない。ベイジアンはMCMCだろうがINLAだろうが最悪[←うわー]。
云々、云々。
… 申し訳ないくらいに適当にしか読んでないんだけど… 著者らも考察のところで書いているように、たとえばベイジアンの推定精度がひどいというのは、事前分布の与え方のせいかもしれないし、サンプラーの問題かもしれないし、MCMCglmmの実装のせいかもしれない。正直こういうシミュレーションって、はあそうでしたか、としかいえないなあ…という感想である。
こういう論文のありがたみはむしろsupplementだという気がしますね。付録にコードをつけてくれているそうだ。R-INLAのサンプルコードもついているはずで(あれ世紀末的にわかりにくいのである)、ありがたいことである。