« 読了：「中国 人口減少の真実」 | メイン | 覚え書き：Rのformulaのなかで使える表現 »
2020年2月29日 (土)
仕事の都合でときどき混合効果モデルを組みたくなることがある。本腰を入れる場合はMplusを使うけど、ほんのちょっとした用事の場合など、ちょっぴり億劫に感じることもある。RStudioのなかでちゃちゃっと済ませちゃいたいな、なあんて... ううむ、Mplus信者としては敬虔さが足りない、すみませんすみません。
Rでなにかのモデルを組むとき、「どのパッケージを使えばいいんだ...」と困ることがあるけど、混合効果モデルはその典型例だと思う。古株としてはnlme (たしか標準パッケージだ)、lme4があり、さらにはglmmMLだの、glmmだの(あ、2018年で開発が止まったようだ)、かの複雑怪奇なMCMCglmmとか... 参っちゃいますね。
Bates, D., Machler, M., Bolker, B.M., Walker, S.C. Fitting Linear Mixed-Effects Models Using lme4.
というわけで、緊急事態に対応するため、どれでもいいからちょっと土地勘をつけておきたいと思い、lme4パッケージのビニエットに目を通した。発行年が書いてないんだけど、これの古い版は2015年にJ. Statistical Softwareに載っている模様。
ほんとはnlmeパッケージについて調べたかったのだけど、適切な文献が見当たらなかったのである。開発者による解説書、ピネイロ&ベイツ「S-PLUSによる混合効果モデル解析」を読めってことでしょうか。うーん、それは勘弁してほしい。あれは滅茶苦茶むずかしい。
ところが、今回はじめて知ったんだけど、lme4パッケージの筆頭開発者はDouglas Batesさん。ピネイロ&ベイツ本の第二著者、かつてピネイロさんとともにnlmeパッケージを開発した人であった。どうやら、Batesさんは2007年にnlme開発チームを離れ、ライバルとなる新パッケージlme4を作り始めた、ということらしい。へええ。
いわく。
lme4はnlmeと比べ...[さまざまな美点が挙げてある]。いっぽうnlmeは、残差に自己相関のような構造があるモデルとか、ランダム効果の共分散行列に構造があるモデルとかを組むときのUIが優れている。
まずは線形混合モデルから。
線形モデルは、ベクトル値ランダム反応変数$\mathcal{Y}$の分布
$\mathcal{Y} \sim N(\mathbf{X} \beta + \mathbf{o}, \sigma^2 \mathbf{W}^{-1})$
として記述できる。$\mathbf{W}$は対角行列で既知のウェイト、$\mathbf{o}$は既知のオフセット項。パラメータは$\beta$と$\sigma^2$である。
同様に、線形混合モデルは$\mathcal{Y}$とランダム効果ベクトル$\mathcal{B}$の分布として記述できる。$\mathcal{B}=\mathbf{b}$の下での条件つき分布は
$(\mathcal{Y} | \mathcal{B} = \mathbf{b}) = N(\mathbf{X} \beta + \mathbf{Zb} + \mathbf{o}, \sigma^2 \mathbf{W}^{-1})$
$\mathcal{B}$の分布は
$\mathcal{B} \sim N(\mathbf{0}, \mathbf{\Sigma})$
ただし$\mathbf{\Sigma}$は半正定値[すべての固有値が非負ってことね]。これは、分散成分パラメータ$\theta$に依存する相対共分散ファクター$\Lambda_\theta$を考えて
$\Sigma_\theta = \sigma^2 \Lambda_\theta \Lambda_\theta^\top$
と書くと便利である。[こういわれると抽象的で困っちゃいますが... ]
さて、lme4パッケージがご提供する線形混合モデルの関数は lmer() であります。
たとえば...[対象者の睡眠を9日間制限し、毎日反応時間課題をやったというデータの紹介...]。この場合、
lmer(Reaction ~ Days + (Days | Subject), data = sleepstudy)
と書く。
実は、lmer() は以下の4つのモジュールの組み合わせである。
その1, Formula.
parsedFormula <- lFormula(formula = Reaction ~ Days + (Days | Subject), data = sleepstudy)
その2, 目的関数.
devianceFunction <- do.call(mkLmerDevfun, parsedFormula)
その3, 最適化.
optimizerOutput <- optimizeLmer(devianceFunction)
その4, 出力. (めんどくさいのでコードは省略するが、mkMerMod()という関数)
モジュール1, Formulaについて。
formula引数は
resp ~ FEexpr + (REexpr1 | factor1) + (REexpr2 | factor2) + ...
という形式をとる。FEexprは固定効果モデル行列$X$を決める表現。問題はその右側、ランダム効果の項である。
個々のREexprは線形モデルのformulaとして評価される。つまりlm()なんかのformulaと同じである。1と書いたら切片である。factorはRの因子型として評価される。縦棒は、効果とグルーピング因子との交互作用を表すと思えばよい。
いくつか例を挙げると、
(1 | g)は「gの個々の水準がランダム切片を持つ」という意味になる。ランダム切片の平均とSDがパラメータになる。1 + (1 | g)と書いてもよい。- ランダム切片の平均が既知の変数oであるならば
0 + offset(o) + (1|g)と書く。先頭の0は-1と書いてもよい。 - g1の下にg2があって、切片がg1の水準間でもg2の水準間でも動くなら
(1|g1) + (1|g1:g2)と書く。これは(1 | g1/g2)とも書ける。[←へー] - g1とg2がネストじゃなくてクロスしているなら
1 + (1|g1) + (1|g2)と書く。IRTなんかがそうだ(被験者と項目がクロスする)。先頭の1は略記できる。 - 切片と傾きの両方が動くなら
1 + x + (1+x|g)と書く。これはx + (x | g)と略記できる[←えっそうなの？ 知らなかった]。lme4では同じランダム効果項に出てくる係数はすべて相関することに注意。つまり切片と傾きは相関する。 - 切片と傾きを無相関にしないなら、
1 + x + (1|g) + (0+x|g)と書けばよい。これはx + (x || g)と略記できる[←へー！]。なお、こういう制約は予測子が比率尺度の場合に限るのが望ましい。この制約の下では予測子をシフトするだけで尤度が変わるから。[←ここ、ちょっとピンとこなかった... いずれ必要になったらゆっくり考えよう]
FEexprから$\mathbf{X}$, その右から$\mathbf{Z}$と$\Lambda_\theta$がつくられる。[以下、$\mathbf{Z}$と$\Lambda_\theta$の作り方について詳しい説明があったけど、半分くらいで力尽きた。メモ省略]
モジュール2, 目的関数モジュールは...[罰則付き最小二乗法によるあてはめの説明とかが12頁にわたって説明されている。読んだところで理解できそうにないのでまるごとパス]
モジュール3. 非線形最適化モジュールは... [ここもどうせわからんので読んでないけど、混合モデルに限らない汎用的な最適化関数らしい。ふーん]
モジュール4. 出力モジュールは... [単に出力の話かと思ったら、意外に仕事が多くて驚いた。まあこれは困ったときに読めばいいかな]
... 途中から面倒になってぱらぱらめくっただけなんだけど、まあいいや。
lme4パッケージの全貌についてのユーザ向け紹介を期待していたんだけど、このビニエットはlmer()の仕組みについて詳しーく説明するものであった。しょうがないのでメモしておくけど、lme4の主力関数は次の3つである。ときに、末尾のerってなんの略なんすかね？
- lmer(): 線形混合モデル
- glmer(): 一般化線形混合モデル
- nlmer(): 非線形混合モデル
nlmeパッケージでいうところのcorrelation = corAR1()のように、誤差項に自己相関を持たせる方法が知りたかったんだけどな... おそらく、lme4でもできなかないけど、めんどくさいのであろう。
論文：データ解析(2018-) - 読了：Bates, et al. (2015?) lme4::lmer()でたのしい線形混合モデリング