1からみなおす線形モデル (3) - 一般化線形混合モデル

はじめに

ということで，引き続き線形モデルの拡張について説明していきます．前回と前々回からさらに発展して，今回は一般化線形混合モデルについて説明していきます．

一般化線形混合モデルとは

前回で説明した一般化線形モデルは，以下のような数式でした．

$log\frac{p}{1-p}=\beta_{\fs{-1}0}+\beta_{\fs{-1}1}x_{\fs{-1}1}+\beta_{\fs{-1}2}x_{\fs{-1}2}$

これに対して，一般化線形混合モデルは次のように，変数がさらに1つ増えます．

$log\frac{p}{1-p}=\beta_{\fs{-1}0}+\beta_{\fs{-1}1}x_{\fs{-1}1}+\beta_{\fs{-1}2}x_{\fs{-1}2}+r}$

この増えた $r$ は，個体差のようなものを意味するパラメタです．一般化線形モデルでは，推定したパラメタ（上の数式でいう $\beta_{\fs{-1}0}$ とか $\beta_{\fs{-1}1}$ とか $\beta_{\fs{-1}2}$ ですね）はどのサンプルに対しても同じ値を取るものでした．それに対してここで加えた $r$ は，サンプルごとに異なる値が推定されます．

そのため，前者の全サンプルに同じ値が当てはめられるものを固定効果，サンプルごとに異なる値が当てはめられるものをランダム効果という呼び方で区別します．一般化線形混合モデルとは，これまでの固定効果に加えてランダム効果がモデルに含まれている，つまり複数の効果がモデルに含まれているという意味で混合なのです．

ちなみに，上の説明ではランダム効果を個体差といいましたが，それ以外の場合もあります．たとえば1学年4クラスの生徒について，勉強時間と成績の関係を線形モデルで表す際に，クラス間の差をモデルに組み込むような場合です．この場合，生徒が所属しているクラス単位で異なる値が $r$ に当てはめられます．

ランダム効果の定式化

一般化線形混合モデルでは，ランダム効果の各変数について，平均ゼロで標準偏差が $s$ の正規分布に従うという仮定を立てます．そして，先ほどの $\beta_{\fs{-1}1}x_{\fs{-1}1}$ とか $\beta_{\fs{-1}2}x_{\fs{-1}2}}$ と合わせて，この $s$ を推定することで，よりデータに合致したパラメタ推定を行うことができるようになるわけです．

Rによるパラメタ推定

今回は，ネタ本の久保先生のサポートサイトからダウンロードしてきた，サンプルデータを使って説明をします（データのリンクはこちら）．

まずは分析に先立って，一般化線形混合モデルのパッケージであるglmmMLをCRANからダウンロードして，読み込んでおきます．その次にダウンロードしたデータを開いてみると，N, y, x, dの4変数があります．これは順に，Nが植物の種子数，yが種子のうち生存していた数，xが葉の数，idが個体ごとの識別番号になります．

そしてglmmML関数で，glm関数と同じように推定を行っていきます．この際に，コマンドの最後にcluster=idとあるのが，個体差の指定になるわけです．その推定結果を見ると，定数項もxの係数も有為であることがみてとれます．また先ほど説明したランダム効果の正規分布分散項のsが2.408になっています．

> install.packages("glmmML")
Installing package(s) into '/Library/Frameworks/R.framework/Versions/2.15/Resources/library'
(as 'lib' is unspecified)
trying URL 'http://cran.rstudio.com/bin/macosx/leopard/contrib/2.15/glmmML_0.82-1.tgz'
Content type 'application/x-gzip' length 508486 bytes (496 Kb)
opened URL
==================================================
downloaded 496 Kb

tar: Failed to set default locale

The downloaded binary packages are in
	/var/folders/69/4pf7w3g930l6n3d7tz04yxr00000gn/T//RtmpqOTnnX/downloaded_packages
> library(glmmML)
> data <- read.csv("~/Downloads/data.csv")
> head(data)
  N y x id
1 8 0 2  1
2 8 1 2  2
3 8 2 2  3
4 8 4 2  4
5 8 1 2  5
6 8 0 2  6
> fit <- glmmML(cbind(y, N-y)~x, data=data, family=binomial, cluster=id)
> summary(fit)

Call:  glmmML(formula = cbind(y, N - y) ~ x, family = binomial, data = data,      cluster = id) 


              coef se(coef)      z Pr(>|z|)
(Intercept) -4.190   0.8777 -4.774 1.81e-06
x            1.005   0.2075  4.843 1.28e-06

Scale parameter in mixing distribution:  2.408 gaussian 
Std. Error:                              0.2202 

        LR p-value for H_0: sigma = 0:  2.136e-55 

Residual deviance: 269.4 on 97 degrees of freedom 	AIC: 275.4