elsur.jpn.org >

« 読了:「中国 人口減少の真実」 | メイン | 覚え書き: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の因子型として評価される。縦棒は、効果とグルーピング因子との交互作用を表すと思えばよい。
 いくつか例を挙げると、

FEexprから$\mathbf{X}$, その右から$\mathbf{Z}$と$\Lambda_\theta$がつくられる。[以下、$\mathbf{Z}$と$\Lambda_\theta$の作り方について詳しい説明があったけど、半分くらいで力尽きた。メモ省略]

 モジュール2, 目的関数モジュールは...[罰則付き最小二乗法によるあてはめの説明とかが12頁にわたって説明されている。読んだところで理解できそうにないのでまるごとパス]

 モジュール3. 非線形最適化モジュールは... [ここもどうせわからんので読んでないけど、混合モデルに限らない汎用的な最適化関数らしい。ふーん]

 モジュール4. 出力モジュールは... [単に出力の話かと思ったら、意外に仕事が多くて驚いた。まあこれは困ったときに読めばいいかな]

 ... 途中から面倒になってぱらぱらめくっただけなんだけど、まあいいや。
 lme4パッケージの全貌についてのユーザ向け紹介を期待していたんだけど、このビニエットはlmer()の仕組みについて詳しーく説明するものであった。しょうがないのでメモしておくけど、lme4の主力関数は次の3つである。ときに、末尾のerってなんの略なんすかね?

 nlmeパッケージでいうところのcorrelation = corAR1()のように、誤差項に自己相関を持たせる方法が知りたかったんだけどな... おそらく、lme4でもできなかないけど、めんどくさいのであろう。

論文:データ解析(2018-) - 読了:Bates, et al. (2015?) lme4::lmer()でたのしい線形混合モデリング

rebuilt: 2020年4月20日 18:52
validate this page