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

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

  

[もくじ]

• glmmML() とは何か
• glmmML() の AIC
• ver.0.65 以降いろいろな 改良 が進行しつつある
• 最尤推定計算が収束しない場合
• posterior.modes で個体差・場所差の事後分布を調べる
• start.sigma 変更わざ

glmmML() とは何か 

• 一般化線形混合モデルの推定計算をする関数
  • random effects 切片 (intercept) をいれたモデルの推定計算ができる．
  • random effects は nest できない; 一個だけしか指定できない
  • 二項分布 (binomial)・ポアソン分布 (poisson) があつかえる
• R の package glmmML に含まれている
  • R の help → library(help = glmmML)
• 作者: Göran Broström (url)
• MCMCglmm() 関数よりもおてがる

glmmML() の AIC 

• glmmML() の算出する AIC は glm() の AIC とは計算方法が異なります
  • glmmML() では AIC = (residual deviance) + 2×(パラメーター数)
  • glm() では AIC = (deviance) + 2×(パラメーター数)
    • もし両者の AIC を比較したいのであれば，とりあえず glm() の結果出力にふくまれている residual deviance を使って，glmmML() と同じ方法で AIC を評価するのが「お手軽」な方法でしょう

ver.0.65 以降いろいろな 改良 が進行しつつある 

• family=binomial において formula = cbind(生起した回数,生起しなかった回数)~... 書法が使えるようになった (ふつうの glm() と同じように書ける)
• stepAIC() (in library(MASS)) を使った自動モデル選択 (パラメーター選択 + 変数選択) ができる -- またできなくなった?
• 「切片」の最大事後確率 (MAP) 推定値も計算される

最尤推定計算が収束しない場合 

• 個体差・場所差がうまく推定できない場合
  • 例 1. 「ホントに」個体差・場所差がない
  • 例 2. 個体差・場所差を推定するにはサンプル数が少なすぎる場合
  • 対策: glmmML() による推定はあきらめて，glm() で推定する
• 個体差などのばらつきパラメーター sigma の初期値設定に問題がある場合
  • 対策: 下記の「start.sigma 変更わざ」を参照
  • しかし，最近の version ではこの問題はあまり発生しなくなってきた

posterior.modes で個体差・場所差の事後分布を調べる 

• fit <- glmmML(...) というふうにあてはめ結果を fit に格納したときに，fit$posterior.modes とすると各 cluster の「ランダム効果の値の分布の最頻値 (mode)」を取り出すことができる
  • つまり個体差・場所差などの個々の事後分布の「山のてっぺん」となる値がわかる
  • なぜかこの部分だけベイズ的に posterior とか記述されている

start.sigma 変更わざ 

glmmML() で「あてはまりが悪そうなモデル」を計算させると，以下のようなエラーがでて計算してくれなくなることがある．

[glmmml] fail = 1
以下にエラーglmmML.fit(X, Y, weights, start.coef, start.sigma, fix.sigma,  : 

こういったときは下のように start.sigma 値を設定してやるとうまくいく場合がある．

glmmML(あれこれ, start.sigma = 2)

• この start.sigma とは sigma の「最尤推定のための試行錯誤計算をスタートする時点での仮の値」である．じゃあ， sigma って何なのかと言えば，GLMM 参照 において「個体差」など random effects を表現させている正規分布の「幅」のことで，値の大きさが「個体差」の大きさをあらわしている．
• このようにユーザーが「脱エラー」目的で start.sigma を指定して計算させる場合，start.sigma の default 値は 0.5 なので help(glmmML) 参照)，0.5 より大きい値を設定してやる必要がある
• 最尤推定の収束計算がうまくいかない場合，start.sigma が小さすぎることが原因になっていることがある．
• データ (応答変数? 説明変数?) に欠測があると，パラメーター推定に失敗する