elsur.jpn.org >

« 読了：「砦に拠る」 | メイン | 読了：Bail, et al. (2018) Twitterでリベラルなツイートを読むと保守派は少しは歩み寄るか？ →それどころかもっと保守的になる »

2018年8月12日 (日)

　GLMM(一般化線形混合モデル)のRパッケージMCMCglmmというのがあるんだけど、その使い方があまりにヤヤコシイ(少なくとも私にとっては)。必要になるたびに資料をめくり、用事が終わるたびにすぐに忘れてしまう、という繰り返しで、流石にうんざりしてきたので、まとめてメモをとることにする。
　参照した資料は以下の通り。

概念
　最初に基本的な概念の説明から。
　MCMCglmmでは、まず反応変数$y_i$の確率密度を
　$f_i (y_i | l_i)$
とする。そういわれるとわかりにくいが、たとえばポアソン分布、ログリンクであれば、ポアソン密度関数を$f_P$として
　$f_P(y_i | \exp(l_i))$
である。リンク関数で変換する前の変数$l$を潜在変数と呼ぶ。

　潜在変数のベクトルを$\mathbf{l}$とする。これは個体レベルの話じゃなくて、たとえ単純なガウシアン回帰でも、100人いたら長さ100である。さらに反応変数が複数あったり、(反応変数が単一でも)潜在変数が複数あったりする場合、それを全部縦に積む。だから$\mathbf{l}$はすごく長い。
　固定効果を持つ説明変数の計画行列を$\mathbf{X}$, ランダム効果を持つ説明変数の計画行列を$Z$として、線型モデル
　$\mathbf{l} = \mathbf{X \beta} + \mathbf{Z u} + \mathbf{e}$
を考える。潜在変数に残差項$\mathbf{e}$が入っている点にご注目。これは測定誤差ではない。ポアソン・モデルであれば、この項でover-dispersionを説明することになる。

　MCMCglmmは、係数はすべてMVNに従うと考え
　$\mathbf{u} \sim N(\mathbf{0},\mathbf{G})$
　$\mathbf{e} \sim N(\mathbf{0}, \mathbf{R})$
とする(これはデータ生成メカニズム)。$\mathbf{G}, \mathbf{R}$は巨大な正方行列である。なお、事前分布は逆ウィシャート分布である模様。
　さらに、MCMCglmmはベイジアン・アプローチに立つので、固定効果の係数$\mathbf{\beta}$でさえ確率変数であり、
　$\mathbf{\beta} \sim N(\beta_0, \mathbf{B})$
だと考える(これは事前分布)。

実行例
　MCMCglmmパッケージ所収のサンプルデータBTdataには、Blue tit(アオガラ)という鳥についてのデータがはいっている。828行, 行はひなの個体を表す。変数は：

　さらに、BTpedには鳥の親子関係がはいっている。1040行、行は個体を表す。変数は:

サンプルコード。わかりやすいようにすべての引数を書いてみる。

m <- MCMCglmm(

  fixed  = cbind(tarsus, back) ~ trait:sex + trait:hatchdate -1, 

  random = ~ us(trait):animal + us(trait):fosternest, 

  rcov   = ~ us(trait):units, 

  family = c("gaussian", "gaussian"), 

  mev    = NULL, 

  data   = BTdata, 

  start  = NULL, 

  prior  = list(

    R = list(V = diag(2)/3, n=2), 

    G = list(

      G1 = list(V = diag(2)/3, n=2), 

      G2 = list(V = diag(2)/3, n=2)

    ),

  ),

  tune        = NULL, 

  pedigree    = BTped, 

  nodes       = "ALL", 

  scale       = TRUE, 

  nitt        = 60000, 

  thin        = 25, 

  burnin      = 10000, 

  pr          = FALSE, 

  pl          = FALSE,

  verbose     = TRUE, 

  DIC         = TRUE, 

  singular.ok = FALSE, 

  saveX       = TRUE, 

  saveZ       = TRUE, 

  saveXL      = TRUE, 

  slice       = FALSE, 

  ginverse    = NULL, 

  trunc       = FALSE

)

　この実行例に沿って、主な引数の指定の仕方をメモしていく。

fixed引数
　fixed引数で固定効果を指定する。

fixed  = cbind(tarsus, back) ~ trait:sex + trait:hatchdate -1

　fixed引数はformulaをとる。formulaの左辺には目的変数を書く。この例のように目的変数が複数ある場合にはcbind()でつなぐ。

　ややこしいのはformulaの右辺である。
　上の例では目的変数が2個ある。MCMCglmmは目的変数の2 x n行列を、(値, 目的変数名, 行番号)の3列からなる2n x 3行列に展開する(reshape2パッケージでいうmelt(), dplyrパッケージでいうgather()ね)。で、目的変数名の列に trait, 行番号の列に unitsという名前を付けるのである。そうです、traitというのは予約語なのです！ 勘弁してくださいな、もっと予約語っぽい名前にしてくださいよ... _TRAIT_ とかさ...。
　上の例では、目的変数tarsusとbackそれぞれについて、固定効果を持つ説明変数としてsexとhatchdateを指定したい。つまり、素直に書けば
　tarsus ~ sex + hatchdate
　back ~ sex + hatchdate
である。もちろん係数はtarsusとbackで異なると考えたい。
　tarsusとbackを縦に積んだ、長さ2 x nの反応変数ベクトル y についてモデルを組むことを考えよう。
　y ~ sex + hatchdate
とすると、係数がtarsusとbackで共通になってしまう。だから、反応変数名とsexの交互作用項 trait:sexと、反応変数名とhatchdateの交互作用項 trait:hatchdate を投入して
　y ~ trait:sex + trait:hatchdate
と書かないといけない... という理屈である。

　最後に -1 として切片項を除外しているのは、そうしないとbackに対するパラメータ推定値がtarsusとの対比になってしまうからである。
　この理屈もしばらく腑におちなかったのだが、上で得たオブジェクトm1について
　View(as.matrix(m1$X))
を観察してようやく意味がわかった(Xは固定効果の計画行列だが、疎行列のdgCMatrixクラスになっているのでas.matrix()しないとわかりにくい)。計画行列は828x2行、8列になる。列の内訳は、まずsexのダミー行列、2(traitの水準数)x3(sexの水準数)=6列。そして、hatchdateの値を格納した2列である(うち828行は左列のみ、残り828行は右列のみ値を持ち、残りは0で埋めてある)。
　仮に
　fixed = cbind(tarsus, back) ~ trait:sex + trait:hatchdate
とすると、計画行列にはすべての行で1を持つ切片列が加わり、かわりにダミー行列の右端の一列(trait=back:sex=UNK)が削除される。そ、それはわかりにくいわ...

　いまこのメモを書いていてふと気になったんだけど、仮に説明変数がhatchdateだけだったらどうするんだろう。hatchdateはダミー行列ではないので、
　fixed = cbind(tarsus, back) ~ trait:hatchdate - 1,
だと切片が消えてしまう。
　fixed = cbind(tarsus, back) ~ trait:hatchdate + trait - 1,
とかかな？ こんど試してみよう。

random引数
　いよいよ本丸である。random引数でランダム効果を指定する。はっきりいってすごくわかりにくいので、深呼吸するように。

random = ~ us(trait):animal + us(trait):fosternest

　random引数はformulaをとる。その書式は
　random = ~ 分散関数(formula):リンク関数(ランダム項) + ...
である。おそらくは左辺はいらないのだと思う。分散関数、formula, リンク関数は省略できる。

　まずは、いくつかの例からみていこう。
　その1, ランダム項だけ指定する場合。
　random = ~ fosternest
もし目的変数がtarsusしかなかったらこれでよい。tarsusに対するfosternestのランダム効果の分散$\sigma^2_f$を推定することになる。$\mathbf{u}$は長さ104のベクトルとなり、$\mathbf{R}$は$104 \times 104$の対角行列で、対角要素はすべて$\sigma^2_f$となる。
　しかし、目的変数がふたつあるときは、これではいけない。なぜか。この指定だと、tarsusに対するfosternestのランダム効果の分散は$\sigma^2_f$だ、backに対するランダム効果の分散も$\sigma^2_f$だ、そして「tarsusに対するfosternestのランダム効果とbackに対するfosternestのランダム効果の共分散」も$\sigma^2_f$だ、と推定することになるのである。
　なにそれ！わけがわからないよ！いやー、そうなる理由がぜんぜん腑に落ちない...

　その2, formulaを付け加える。
　random = ~ trait:fosternest
　目的変数が複数あるときは、この式が基本である。
　いま、$\mathbf{u}$は長さ828x2の長ーいベクトルである。ランダム項としてfosternestを指定したので、計画行列$\mathbf{Z}$は828x2行、104x2列の行列となる。$\mathbf{G}$は828x2行、828x2列の壮大な正方行列である。さて、この$\mathbf{G}$から、ある特定の巣に関わる行と列だけを取り出すと、2x2の正方行列になる。行と列はtraitである。
　random = ~ trait:fosternestとは、この2x2の分散共分散行列を、非対角要素は0、対角要素は同一とします、という指定である。これ、おそらくidv(trait):fosternest の略記だと思う。

　その3, 分散関数を付け加える。
　random = ~ idh(trait):fosternest
とすると、分散共分散行列の対角要素だけが別々に推定される。つまり、tarsusに対するfosternestのランダム効果の分散$\sigma^2_{f:tarsus}$とbackに対するfosternestのランダム効果の分散$\sigma^2_{f:back}$の2つのパラメータが推定される。非対角要素は0。
　random = ~ us(trait):fosternest
とすると、非対角要素を含めた4つが別々に推定される。
　この事例では、巣が体長に与える効果と体色に与える効果のあいだに共分散があると考えているので、~ us(trait): fosternest となっているわけだ。

　さて、書式の詳細について。
　分散関数には以下がある。

　formulaは、上の例ではtraitになっているが、量的変数を指定しても良い。たとえば
　random = ~ us(1+age):individual
とすると、切片とageのスロープがindividualごとにランダムに決まることになる。
　random = ~ us(1+poly(age,2)):individual
もアリである。うーむ、ややこしい。説明の流れのせいでわかりにくくなっているが、formulaでは「ある個体における計画行列」を指定し、ランダム項には個体を表す変数を指定する、という理屈であろう。

　リンク関数には次の2種類がある。いずれも、複数のランダム項を+でつないで記述する。

　なお、上の例でus(trait):animalというのをいれているのは、pedigree引数を指定しているからなのだそうである。ここはよくわからない。

rcov引数
　rcov引数で$\mathbf{R}$を指定する。

rcov   = ~ us(trait):units

　書式はrandom引数と同じだが、項はひとつしかとれない。
　おそらく、多変量なら ~ us(traits):units とするのが基本であろう。単変量なら ~ units かな？

family引数
　family引数で分布を指定する。

family = c("gaussian", "gaussian")

　family引数は、反応変数の数だけの文字列ベクトルをとる。指定できる分布はは以下の通り。

mev引数
　各データ点における測定誤差分散を表すベクトルを指定できるらしい。
　この引数はメタ分析の際に使うことを想定しているらしいのだが... まてよ、これを使えば小地域推定のFay-Harriotモデルを推定できるのではなかろうか。

prior引数　もうひとつの山場、事前分布の指定である。

prior  = list(

    R = list(V = diag(2)/3, n=2), 

    G = list(

      G1 = list(V = diag(2)/3, n=2), 

      G2 = list(V = diag(2)/3, n=2)

    ),

  )

　prior引数はリストをとる。リストには3つの要素をいれることができる。

　さて、上述のR、ならびにGの各要素は、以下の要素を持つリストである。

　わかりにくいので例に基づいて考えよう。
　上の例では、R, G1, G2のいずれについても、Vは2x2の対角行列(対角要素は1/3)、nは2としている。Rに注目しよう。rcov引数により、体長の潜在変数における残差項の分散$\sigma^2_{e:tarsus}$、体色の潜在変数の残差項の分散$\sigma^2_{e:back}$, 体長の潜在変数における残差項と体色の潜在変数における残差項との共分散$\sigma_{e:tarsus,back}$の3つを推定することになっている。しかしここでは、Rの持つVは対角行列で、対角要素は1/3である。つまり、$\sigma_{e:tarsus,back}$の事前分布の期待値は0だと指定していることになる。
　対角要素を1/3としているのは、この事例では反応変数がすでに標準化されていて分散1であり、それを3つの要素に均等に分けている、という意味合いである。

　ところで、prior引数の要素R、ならびに要素Gの各要素は、デフォルトではlist(V=1, nu=0)なのだそうである。うーん、よくわからん。

starts引数
　初期値の指定。リストをとる。次の4つの要素をいれることができる。

　... と、この辺で力尽きたので、メモはここまでにしておく。さらに難解なpedigree引数、ginverse引数には、残念ながら踏み込めなかった。またの機会に。
　いやー、自分の能力を棚に上げていうけど、なんてわかりにくいんだろう。SASのproc mixedのほうがはるかにわかりやすい。参ったなあ。

雑記：データ解析 - 覚え書き：RのMCMCglmmパッケージの使用方法の難解さを目の前にして途方に暮れるの巻

rebuilt: 2020年11月16日 22:54
validate this page