線形混合効果モデルメモ

線形混合効果モデルについて調べてみたメモ。

線形混合効果モデルは分野によって呼び方が異なるらしい。

R で線形混合効果モデルを使うには、lme4 パッケージが標準的みたいだけど、ここでは mlmRev パッケージを使う。
(mlmRev パッケージの中では結局 lme4 パッケージを使っている)

library(mlmRev)
## Loading required package: lme4
## Loading required package: Matrix
## Loading required package: Rcpp

データは mlmRev パッケージの Exam データ。
ロンドンの 69 の学校における 4509 人の学生の試験の点数のデータらしい。

head(Exam)
##   school normexam schgend schavg      vr     intake standLRT sex type
## 1      1   0.2613   mixed 0.1662 mid 50% bottom 25%   0.6191   F  Mxd
## 2      1   0.1341   mixed 0.1662 mid 50%    mid 50%   0.2058   F  Mxd
## 3      1  -1.7239   mixed 0.1662 mid 50%    top 25%  -1.3646   M  Mxd
## 4      1   0.9676   mixed 0.1662 mid 50%    mid 50%   0.2058   F  Mxd
## 5      1   0.5443   mixed 0.1662 mid 50%    mid 50%   0.3711   F  Mxd
## 6      1   1.7349   mixed 0.1662 mid 50% bottom 25%   2.1894   M  Mxd
##   student
## 1     143
## 2     145
## 3     142
## 4     141
## 5     138
## 6     155

normexam が試験の得点を標準化した数値で、standLRT が入学時に実施された読解力テストの得点を標準化したもの。
この2つの関係を見てみると、かなりきれいな相関が見える。

data <- Exam
plot(normexam ~ standLRT, data)

plot of chunk unnamed-chunk-3

これに通常の線形モデル

\[y_i = \beta_0 + \beta_1 x_i + \varepsilon_i\]

を適用すると、次のような結果になる。

linear.model <- lm(normexam ~ standLRT, data)
summary(linear.model)
## 
## Call:
## lm(formula = normexam ~ standLRT, data = data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.6562 -0.5185  0.0126  0.5440  2.9740 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -0.00119    0.01264   -0.09     0.92    
## standLRT     0.59506    0.01273   46.74   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.805 on 4057 degrees of freedom
## Multiple R-squared:  0.35,   Adjusted R-squared:  0.35 
## F-statistic: 2.19e+03 on 1 and 4057 DF,  p-value: <2e-16

これに対して、混合効果モデルは、学校ごとの成績の差異をモデルに組み込むことができる。
学校ごとに試験の得点と読解力テストの関係に違いがあるのか図示してみる。

library(lattice)
library(dplyr)
sample <- data %>% filter(as.integer(school) <= 15) # 15 校のみ抽出
xyplot(normexam ~ standLRT|school, sample)

plot of chunk unnamed-chunk-5

よくわからないので、級内相関係数 ICC(Intraclass correlation cofficient) を調べる。
まず、最小モデル

\[y_{ij} = \beta_{0j} + \varepsilon_{ij}\] \[\beta_{0j} = \gamma_{00} + u_{0j}\]

を考える。
最小モデルは、変量効果を持つ切片のみを仮定する最も単純な混合効果モデルである。

null.model <- lmer(normexam ~ 1 + (1 | school), data=data)
summary(null.model)
## Linear mixed model fit by REML ['lmerMod']
## Formula: normexam ~ 1 + (1 | school)
##    Data: data
## 
## REML criterion at convergence: 11015
## 
## Scaled residuals: 
##    Min     1Q Median     3Q    Max 
## -3.947 -0.649  0.012  0.699  3.657 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  school   (Intercept) 0.172    0.414   
##  Residual             0.848    0.921   
## Number of obs: 4059, groups: school, 65
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept)  -0.0133     0.0541   -0.25

この最小モデルにおいて、\(u_{0j}\) の分散は級間分散と呼ばれ、全体の分散に対する級間分散の占める割合を級内相関係数 ICC と呼ぶ。

\[{\rm ICC} = \frac{{\rm VAR}(u_{0j})}{{\rm VAR}(\varepsilon_{ij}) + {\rm VAR}(u_{0j})}\]

今回のデータでは、

\[{\rm ICC} = \frac{0.172}{0.172 + 0.848} = 0.1683\]

となり、全体の分散の 17% が学校間の違いにより説明される。
17% が大きいのかどうかはよくわからないが、ICC が大きい場合には混合効果モデルを適用したほうが良い。
(10%以上あれば混合効果モデルを適用した方がいいという説があるらしい。)

さて、本格的に混合効果モデルを作成していこう。
まずは、切片のみに変量効果を持ち、係数は固定効果とするモデル

\[y_{ij} = \beta_{0j} + \beta_{1} x + \varepsilon_{ij}\] \[\beta_{0j} = \gamma_{00} + u_{0j}\]

を考える。このモデルをランダム切片モデル(RI)と呼ぶことにする。

ri.model <- lmer(normexam ~ standLRT + (1 | school), data=data)
summary(ri.model)
## Linear mixed model fit by REML ['lmerMod']
## Formula: normexam ~ standLRT + (1 | school)
##    Data: data
## 
## REML criterion at convergence: 9369
## 
## Scaled residuals: 
##    Min     1Q Median     3Q    Max 
## -3.717 -0.630  0.029  0.685  3.267 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  school   (Intercept) 0.0938   0.306   
##  Residual             0.5659   0.752   
## Number of obs: 4059, groups: school, 65
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept)  0.00232    0.04035     0.1
## standLRT     0.56331    0.01247    45.2
## 
## Correlation of Fixed Effects:
##          (Intr)
## standLRT 0.008

次に、切片にも係数にも変量効果を持つモデル

\[y_{ij} = \beta_{0j} + \beta_{1j} x + \varepsilon_{ij}\] \[\beta_{0j} = \gamma_{00} + u_{0j}\] \[\beta_{1j} = \gamma_{10} + u_{1j}\]

を考える。このモデルをランダム切片・係数モデル(RIC)と呼ぶことにする。

ric.model <- lmer(normexam ~ standLRT + (standLRT + 1 | school), data=data)
summary(ric.model)
## Linear mixed model fit by REML ['lmerMod']
## Formula: normexam ~ standLRT + (standLRT + 1 | school)
##    Data: data
## 
## REML criterion at convergence: 9328
## 
## Scaled residuals: 
##    Min     1Q Median     3Q    Max 
## -3.832 -0.632  0.034  0.683  3.456 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev. Corr
##  school   (Intercept) 0.0921   0.304        
##           standLRT    0.0150   0.122    0.49
##  Residual             0.5536   0.744        
## Number of obs: 4059, groups: school, 65
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept)  -0.0116     0.0401   -0.29
## standLRT      0.5565     0.0201   27.67
## 
## Correlation of Fixed Effects:
##          (Intr)
## standLRT 0.365

最小モデル(null)、ランダム切片モデル(RI)、ランダム切片・係数モデル(RIC)を比較してみる。

anova(null.model, ri.model, ric.model)
## refitting model(s) with ML (instead of REML)
## Data: data
## Models:
## null.model: normexam ~ 1 + (1 | school)
## ri.model: normexam ~ standLRT + (1 | school)
## ric.model: normexam ~ standLRT + (standLRT + 1 | school)
##            Df   AIC   BIC logLik deviance  Chisq Chi Df Pr(>Chisq)    
## null.model  3 11017 11036  -5505    11011                             
## ri.model    4  9365  9390  -4679     9357 1653.4      1    < 2e-16 ***
## ric.model   6  9329  9367  -4658     9317   40.4      2    1.7e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

AIC, BIC ともに RIC モデルがいいみたい。

変量効果を考えることでよりよいモデル構築ができた。

(続きはあとで書く)