線形混合効果モデルメモ
線形混合効果モデルについて調べてみたメモ。
線形混合効果モデルは分野によって呼び方が異なるらしい。
- 教育学
- 階層線形モデル(Hierarchical linear model)
- 計量経済学
- 線形混合効果モデル(Linear mixed-effects model)
- 生物統計学
- ランダム係数回帰モデル(Random-coefficient regression model)
- 社会学
- マルチレベルリニアモデル(Multilevel linear model)
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 155normexam が試験の得点を標準化した数値で、standLRT が入学時に実施された読解力テストの得点を標準化したもの。
この2つの関係を見てみると、かなりきれいな相関が見える。
data <- Exam
plot(normexam ~ standLRT, data)
これに通常の線形モデル
\[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)
よくわからないので、級内相関係数 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 ' ' 1AIC, BIC ともに RIC モデルがいいみたい。
変量効果を考えることでよりよいモデル構築ができた。
（続きはあとで書く）