GLMMをRで実行する方法

 

GLMMをSASで実行する方法をすでにアップしましたが,次はRで実行する方法についてまとめます。

 

RでGLMMができる関数

RではGLMMを実行するためのプロシージャはいくつかあります。代表的なのは,glmmMLパッケージのglmmML関数と,lme4パッケージのglmer関数でしょうか。

glmmML関数は,ガウス-エルミート求積法による積分計算を行うので推定精度は高いようです。ただし,SASの時と同様に,変量効果は1種類しか指定できません。

一方,lme4パッケージのglmer関数は,ラプラス近似を用いますが,変量効果は複数推定することができます。

どちらがオススメか,というのは難しいところですが,変量効果が一つだけならglmmMLが,それ以上ならlme4がいいように思います。lme4はバージョンによっても結果が変わってきていて,現状の最新版のver 1.1-7はSASと一致する結果が得られましたが,ver 1.0-5では標準誤差がやや小さめに推定されていました。この辺りの不安定な感じが不安なら,glmmMLかSASを使うといいでしょう。

他にも,glmmパッケージなども出ているようですが,まだ触ってないので現状ではよくわかりません。

 

今回使うサンプルデータ

データはSASの時と同じ,一人から4回反復測定したデータを用います。目的変数yは2値なので二項分布に従うデータです。説明変数xは1~5の値を取る連続量です。tは試行を意味しています。個人を識別する変数はIDで,10人のデータです。以下のデータをcsvファイルに保存し,binomial.csvとしてRの実行ファイルと同じディレクトリに保存してください。

ID,t,y,x
1,1,0,1
1,2,0,1
1,3,0,1
1,4,0,1
2,1,0,1
2,2,1,1
2,3,0,1
2,4,0,1
3,1,1,2
3,2,0,2
3,3,1,2
3,4,0,2
4,1,1,2
4,2,1,2
4,3,0,2
4,4,0,2
5,1,0,3
5,2,0,3
5,3,0,3
5,4,0,3
6,1,1,3
6,2,1,3
6,3,1,3
6,4,1,3
7,1,1,4
7,2,1,4
7,3,1,4
7,4,1,4
8,1,0,4
8,2,1,4
8,3,1,4
8,4,0,4
9,1,1,5
9,2,1,5
9,3,1,5
9,4,1,5
10,1,1,5
10,2,1,5
10,3,1,5
10,4,1,5

 

glmmMLによるGLMM

glmmMLは変量効果を一つだけ指定できます。変量効果はcluster=で変数を指定します。今回はIDがそれです。よって,以下の様なコードになります。

dat

library(glmmML)

fit.ml summary(fit.ml)

y~xは,lmと同じ式を書きます。familyは,binomialとpoissonが指定できます。clusterは先ほど書いたように変量効果を推定する変数を指定します。method="ghq"はガウス-エルミート求積法を指定しています。method = "laplace"にすればラプラス近似による積分計算となります。

これを実行すると,

 

              coef se(coef)      z Pr(>|z|)
(Intercept) -3.734   1.9299 -1.935   0.0530
x            1.485   0.6917  2.147   0.0318

Scale parameter in mixing distribution:  1.503 gaussian 
Std. Error:                              0.9771 

        LR p-value for H_0: sigma = 0:  0.05682

 

という結果が得られます。SASの結果とほとんど同じですが,自由度を考慮に入れたt検定を使っていないので,p値は小さめに推定されています。

 

lme4のglmer関数

glmer関数は,複数の変量効果が指定できるところが便利です。変量効果の指定方法がglmmMLと違っていて,()の中に変量効果を指定します。

上と同じ分析をglmerで角と次のようになります。

dat

library(lme4)

fit.glmer summary(fit.glmer)

glmerでは,モデル指定のところに変量効果も同時に書きます。y ~ xまでは同じですが,+ (1|ID)が違います。最初の1は切片を意味していて,SASでいうところのinterceptと書いたものと同じです。"|"のあとにIDとあるのは,IDの変動についての変量効果を推定しろ,という意味です。

これを実行すると,

 

Formula: y ~ x + (1 | ID)
   Data: dat

     AIC      BIC   logLik deviance df.resid 
    43.5     48.5    -18.7     37.5       37 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-1.4344 -0.4761  0.1457  0.3948  2.2026 

Random effects:
 Groups Name        Variance Std.Dev.
 ID     (Intercept) 2.128    1.459   
Number of obs: 40, groups:  ID, 10

Fixed effects:
            Estimate Std. Error z value Pr(>|z|)  
(Intercept)  -3.7343     1.9185  -1.946   0.0516 .
x             1.4819     0.6879   2.154   0.0312 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

 

これもほぼSASと同じ結果です。laplaceによる方法なので少し違っていますが,SASでlaplaceを選んだ場合と同じ結果となっています。glmmMLと同様に,自由度を考慮に入れていないのでp値は小さめに推定されます。

二つ以上の変量効果を推定する場合は,次のように書きます。たとえばすべての参加者に同じ刺激セットを用いてデータを取ったとすると,変数tによる刺激間変動の推定が可能です。これをglmerで書くと,

fit.glmer2 summary(fit.glmer2)

となります。要は,変量効果の種類は()を増やしていけばいい,というだけです。この文法はなかなかシンプルでわかりやすいので個人的には好きです。

ただ,どうもglmerではSASのようにうまく刺激間変動を推定できませんでした。おそらくかなり小さい値だからだと思います。

 

これまでは最尤推定によるGLMMの方法を書きましたが,また別の記事でベイズ推定によるGLMMの話も書きたいと思います。

 

This entry was posted in 心理統計学. Bookmark the permalink.