MASS パッケージの nlschools データを使い、子供をレベル1、クラスをレベル2とみなし、子供の言語 (lang) の成績を従属変数として、以下の要領で複数のモデルを推定し、その結果を比較しなさい。
library(MASS) # nlme パッケージの呼び出し
head(nlschools) # データの中身を少しだけ確認
## lang IQ class GS SES COMB
## 1 46 15.0 180 29 23 0
## 2 45 14.5 180 29 10 0
## 3 33 9.5 180 29 15 0
## 4 46 11.0 180 29 23 0
## 5 20 8.0 180 29 10 0
## 6 30 9.5 180 29 10 0
library(lme4)
## Loading required package: Matrix
m1 <- lmer(lang ~ IQ + SES + (1 | class), data = nlschools)
m2 <- lmer(lang ~ IQ + SES + (IQ | class), data = nlschools)
iq.c <- nlschools$IQ - mean(nlschools$IQ) # 全体平均センタリング
m3 <- lmer(lang ~ iq.c + SES + (iq.c | class), data = nlschools)
iq.m <- ave(nlschools$IQ, nlschools$class) # グループ平均
iq.gc <- nlschools$IQ - iq.m
m4 <- lmer(lang ~ iq.gc + SES + (iq.gc | class), data = nlschools)
library(texreg)
## Warning: package 'texreg' was built under R version 3.3.1
## Version: 1.36.7
## Date: 2016-06-21
## Author: Philip Leifeld (University of Glasgow)
##
## Please cite the JSS article in your publications -- see citation("texreg").
screenreg(list(m1, m2, m3, m4))
##
## ====================================================================================
## Model 1 Model 2 Model 3 Model 4
## ------------------------------------------------------------------------------------
## (Intercept) 9.39 *** 9.01 *** 36.25 *** 35.47 ***
## (0.87) (1.10) (0.50) (0.57)
## IQ 2.25 *** 2.30 ***
## (0.07) (0.08)
## SES 0.17 *** 0.16 *** 0.16 *** 0.18 ***
## (0.01) (0.01) (0.01) (0.01)
## iq.c 2.30 ***
## (0.08)
## iq.gc 2.23 ***
## (0.09)
## ------------------------------------------------------------------------------------
## AIC 15150.67 15136.40 15136.40 15218.28
## BIC 15179.35 15176.54 15176.54 15258.42
## Log Likelihood -7570.34 -7561.20 -7561.20 -7602.14
## Num. obs. 2287 2287 2287 2287
## Num. groups: class 133 133 133 133
## Var: class (Intercept) 9.16 63.95 9.16 18.16
## Var: Residual 40.03 39.24 39.24 39.30
## Var: class IQ 0.22
## Cov: class (Intercept) IQ -3.61
## Var: class iq.c 0.22
## Cov: class (Intercept) iq.c -1.02
## Var: class iq.gc 0.22
## Cov: class (Intercept) iq.gc -1.10
## ====================================================================================
## *** p < 0.001, ** p < 0.01, * p < 0.05
Model 1 (ランダム切片モデル)よりも、Model 2 (ランダム傾きモデル)のほうが AIC も BIC も小さく、若干フィッティングがよい。切片の分散 “Var: class (Intercept)” をみると、9.16 から 63.95 に増加しており、IQ をセンタリングしていないことに起因していると考えられる。Model 2 と Model 3 はフィッティングは同じであるが、切片の分散が異なり、Model 3 のほうは Model 1 と同じであり、こちらのほうがクラスによる平均的な成績の違いを適切に反映していると思われる。Model 4 では IQ をグループ平均センタリングしたためにグループ間の違いが IQ では説明されていないために、その分、切片の分散が大きくなっている。クラスによる IQ の傾きの分散は Model 2 や Model 3 と同じである。
m5 <- update(m3, ~ . + iq.m)
m6 <- update(m4, ~ . + iq.m)
screenreg(list(m5, m6), custom.model.names = c("Model 5", "Model 6"))
##
## ========================================================
## Model 5 Model 6
## --------------------------------------------------------
## (Intercept) 23.71 *** -3.58
## (3.81) (3.50)
## iq.c 2.26 ***
## (0.08)
## SES 0.16 *** 0.16 ***
## (0.01) (0.01)
## iq.m 1.07 *** 3.38 ***
## (0.32) (0.30)
## iq.gc 2.26 ***
## (0.09)
## --------------------------------------------------------
## AIC 15128.65 15126.64
## BIC 15174.53 15172.52
## Log Likelihood -7556.32 -7555.32
## Num. obs. 2287 2287
## Num. groups: class 133 133
## Var: class (Intercept) 8.36 8.27
## Var: class iq.c 0.21
## Cov: class (Intercept) iq.c -0.80
## Var: Residual 39.17 39.05
## Var: class iq.gc 0.27
## Cov: class (Intercept) iq.gc -0.82
## ========================================================
## *** p < 0.001, ** p < 0.01, * p < 0.05
Model 6 の推定結果を見ると、IQ のグループ間効果が 3.38 で、グループ内効果は 2.26 で、グループ間効果のほうが大きい。Model 5 の iq.m がグループ間効果とグループ内効果の差を示しており、有意である。ただし、Model 5 から計算したグループ間効果が \(2.26 + 1.07 = 3.33\) であるのに対して Model 6 では 3.38 であり、若干のずれがあるのがわかる。AIC 等の値も Model 5 と Model 6 で少し異なるが、この程度の誤差は分布の歪みやグループサイズのばらつきによって生じうるということであろう。
faraway パッケージの jsp データを使い、パーソン・イヤーがレベル1、学校をレベル2とみなし、子供の数学の成績 (math) を従属変数として、以下の要領で複数のモデルを推定し、それらの適合度を比較しなさい。
gender, social, raven を独立変数として raven についてランダム傾きを仮定し、例題と同じ要領でモデルを作って比較しなさい。 2 練習問題と同じ要領で、raven を全体平均センタリングしたモデルとグループ平均センタリングしたモデルを推定し、グループ間効果とグループ間効果を比較検討しなさい。