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は１～５の値を取る連続量です。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 <- read.csv("binomial.csv")

library(glmmML)

fit.ml <- glmmML (y ~x, data =dat, family =binomial, cluster = ID, method = "ghq")

summary(fit.ml)

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

これを実行すると，

 

 

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

 

lme4のglmer関数

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

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

dat <- read.csv("binomial.csv")

library(lme4)

fit.glmer <- glmer(y ~ x + (1|ID), data = dat, family = binomial)

summary(fit.glmer)

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

これを実行すると，