1 Example: 例題

1.1 Question

MASS パッケージの nlschools データを使い、子供をレベル1、クラスをレベル2とみなし、子供の言語 (lang) の成績を従属変数として、以下の要領で複数のモデルを推定し、その結果を比較しなさい。

  1. 独立変数として IQ と SES を投入し、Model 1: ランダム切片モデル、Model 2: IQ のランダム傾きを仮定したモデル、Model 3: IQ を全体平均センタリングしてランダム傾きを仮定したモデル、の3つをそれぞれ推定し、Model 4: グループ平均センタリングしてランダム傾きを仮定したモデル、 のフィッティングとランダム効果の分散を比較せよ。
  2. さらに Model 5: Model 3 に SES のクラス平均を独立変数として追加したモデル、Model 6: Model 4 に SES のクラス平均を独立変数として追加したモデルを推定し、グループ間効果とグループ内効果の大きさを比較検討しなさい。

1.2 Answer

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 で少し異なるが、この程度の誤差は分布の歪みやグループサイズのばらつきによって生じうるということであろう。

2 Practice: 練習問題

  1. faraway パッケージの jsp データを使い、パーソン・イヤーがレベル1、学校をレベル2とみなし、子供の数学の成績 (math) を従属変数として、以下の要領で複数のモデルを推定し、それらの適合度を比較しなさい。

  2. gender, social, raven を独立変数として raven についてランダム傾きを仮定し、例題と同じ要領でモデルを作って比較しなさい。 2 練習問題と同じ要領で、raven を全体平均センタリングしたモデルとグループ平均センタリングしたモデルを推定し、グループ間効果とグループ間効果を比較検討しなさい。

inserted by FC2 system