生態学のデータ解析 - lmer 紹介

• 線形混合モデルと一般化線形混合モデル (GLMM 参照) の推定計算する関数 lmer() を紹介します (R の関数)

• 複雑なモデル (例: 複数の種類の random effects をあつかう) を lmer() で推定させると， あまり正しくない推定結果が得られる場合が「よく」あります (2011-12-27 ……けっきょくこの欠点は何年たっても改善されなかった)
• このように複雑なモデル WinBUGS や JAGS， あるいは MCMCglmm() を使って推定したほうがよいと思います

[もくじ]

• lmer() とは何か
• 他の推定関数との比較
  • 実験: 乱数つかったデータ生成と GLMM 推定
  • 結果 (1): 推定値 b の比較
  • 結果 (2): deviance 比較 (glmmML() vs lmer())
  • とりあえずの，まとめ
  • Appendix: そもそも method = "Laplace" って?

lmer() とは何か 

• こまごまとしたことなどは lmer()雑 に
• 線形混合モデルと一般化線形混合モデルの推定計算をする関数
  • R の help: lmer()
• R の package lme4 に含まれている
  • R-2.4.0 より前の version の R では library(Matrix) に含まれている
• lmer() による最尤推定したいときの method 指定方法
  • family=gaussian: method = "ML"
  • family=binomial or poisson: method = "Laplace" (これは Laplace 近似による近似尤度計算)
• R に混合モデルを計算できる関数はいろいろある けれど，lmer() (method = "Laplace" 指定時の) 長所としては
  • 複雑な random effects -- 階層性のある random effects や「ランダム傾き」をあつかえる
  • Penalized quasi-likelihood (PQL) だけでなく (近似的な) 最尤推定も可能
  • モデル選択できるかもしれない (下の実験結果をみてください)
  • さらに mcmcsamp() 関数を使えば lmer, glmer オブジェクトを階層ベイズモデルとみなして MCMC 計算ができる → lmer MCMC 計算
• いっぽうで短所は
  • 近似的最尤推定の method が method = "Laplace" しか準備されておらず，より高精度ということになっている method = "AGQ" はまだ使えない (ver.0.995-8 時点において)
  • まだまだ未完成状態 -- "... not yet implemented" (未実装です) エラー多い
  • 「たいへん」(?) なデータで計算させると R ごと落ちる場合がある
• 以前は lme() あるいは lme4() だった関数が lmer() になっている?

他の推定関数との比較 

• ここでは「(よく使われている) glmmPQL() in library(MASS) より lmer() を method = "Laplace" で使ったほうがマシじゃないでしょうか」という趣旨で同じデータに対する推定結果比較をやってみました
• 例題としてとりあげるのは 自由集会2006 の久保 GLMM 紹介でつかった「架空植物の結実率は葉数 x にどう依存しているか」です
  • 自由集会ペイジにある PDF ファイルを参照してください
  • あるいは生態学会誌 (2006) 掲載の GLMM 解説記事 の説明でつかっている架空例とも同じ

実験: 乱数つかったデータ生成と GLMM 推定 

• 自由集会2006 の久保 GLMM 紹介と同じ状況
• (結実率) ~ a + b × (葉数) + (「個体差」) という混合ロジスティックモデルで表現できる結実率，パラメーター b の「真の値」は 1 なのでこれを正確に推定したい
• 「個体差」は平均ゼロで標準偏差 3 の正規乱数 (毎回生成，これはかなり「個体差」大きい)
• 個体内の 8 個の胚珠の結実は二項乱数で与える
• 乱数で生成されたデータから以下の四とおりの方法でパラメーター推定した
  • glm() -- random effects なし
  • glmmPQL() -- PQL 最大化によるパラメーター推定
  • glmmML() -- 最尤法
  • lmer() -- method = "Laplace" を指定してLaplace 積分による近似最尤法
• データ生成→GLMM を 100 回くりかえした

結果 (1): 推定値 b の比較 

• b = 1 と推定するのが正しい
• glm() の推定結果は b がかなり小さい方向にずれている
• glmmPQL() もやはり小さいほうにずれている
• glmmML() と lmer() の推定はほとんど偏っていない -- しかしばらつきは大きい

結果 (2): deviance 比較 (glmmML() vs lmer()) 

• 100 回乱数生成した各データセットに対して glmmML() と lmer() が計算した deviance (= -2 × 最大化対数尤度) を比較してみた
• 両者はだいたい同じ deviance (最大化対数尤度) を計算している
• また両者で b の推定値を比較すると

• 同じデータセットに対してだいたい同じ推定値が得られている
• これに対して fit<-glmmPQL(...) して -2*fit$logLik なる量と glmmML() の deviance を比較すると

• おそらく fit$logLik の符合が逆なのだろう ……
  • すると glmmML() の deviance とよい相関があるように見える (lmer() の場合よりはばらついているが)
  • しかしながら，その場合でも両者をつなぐ「傾き」が 1 ではないので AIC など尤度にもとづくモデル選択法は使えない

とりあえずの，まとめ 

• glmmML() でも推定できるようなこういう単純な問題に関しては lmer() は十分よい推定をしているのかもしれない (method = "Laplace" を指定した場合)
• さらに複雑な random effects のもとでの推定について実験が必要である
• いっぽうで glmmPQL() はこれほど単純な問題に関しても偏った推定結果をだしている
• AIC など使ったモデル選択について
  • 結果 (2) をみると Laplace 近似でだいたい正しい尤度が計算できているように見える (こういう単純な random effects のもとでは)
  • 調べたいモデルをすべて lmer() で計算するなら AIC によるモデル選択が可能かもしれない
  • ただしこのときに glmmML() など他の方法で得られた推定結果とその AIC をまぜてモデル選択してはいけない

Appendix: そもそも method = "Laplace" って? 

• GLMM 最尤推定推定計算の難点は「すべての random effects に関する積分」にあります
• random effects が単純な場合 (階層構造ナシの random intercept だけ) glmmML() などは「きあい」でこれを数値積分して尤度を計算しています (積分が一個だけなので)
• しかし random effects が複雑になると「きあい」だけでは数値積分できません (多重積分になるので)
• method = "Laplace" は「Laplace 近似」という (ラプラス変換をつかった) 数値積分技法を使って面倒な場合においてもできるだけ正確な尤度を何とか計算しようとする試みです (尤度の近似的計算)
  • したがって「近似が高精度」な計算に成功した場合，近似的尤度はホントの尤度にかなり近い値になります
  • 上の「結果 (2): deviance 比較」も参照してください
• これに対して quasi-likelihood は「(計算がラクな) 尤度みたいな量」を最大化します
  • こちらはホントの尤度との関係が明確ではありません