ランダム切片モデルの練習
太郎丸 博
2016年10月28日
1 Example: 例題
1.1 Question
nlme パッケージの MathAchieve データを使い、以下の要領でランダム切片モデルを推定し、結果を解釈しなさい。
- 子供がレベル1、学校をレベル2とみなし、
- 子供の数学の成績を従属変数、その他の変数を独立変数とみなす。
1.2 Answer
1.2.1 Summary of Data: データの概要
用いるデータの要約統計量を示したのが以下の表1 であり、従属変数である数学の成績 (MathAch) のヒストグラム、数学の成績と Minority (マイノリティかどうか) 、Sex (性別)、SES (社会経済的地位) の関係を図示したのが、以下の図1 である。
library(nlme) # nlme パッケージの呼び出し
head(MathAchieve) # データの最初の6行だけながめる## Grouped Data: MathAch ~ SES | School
## School Minority Sex SES MathAch MEANSES
## 1 1224 No Female -1.528 5.876 -0.428
## 2 1224 No Female -0.588 19.708 -0.428
## 3 1224 No Male -0.528 20.349 -0.428
## 4 1224 No Male -0.668 8.781 -0.428
## 5 1224 No Male -0.158 17.898 -0.428
## 6 1224 No Male 0.022 4.583 -0.428# ?MathAchieve # ? でデータの解説が見られる。 "#" はつけずに "?" から実行せよ
library(knitr) # 以下の kable() を使わないなら knitr はロードしなくていい
kable(summary(MathAchieve, digits =2), caption = "表1: 用いるデータの記述統計") # kable()は使わなくてもよい| School | Minority | Sex | SES | MathAch | MEANSES | |
|---|---|---|---|---|---|---|
| 2305 : 67 | No :5211 | Male :3390 | Min. :-3.75800 | Min. :-2.8 | Min. :-1.1880 | |
| 5619 : 66 | Yes:1974 | Female:3795 | 1st Qu.:-0.53800 | 1st Qu.: 7.3 | 1st Qu.:-0.3170 | |
| 4292 : 65 | NA | NA | Median : 0.00200 | Median :13.1 | Median : 0.0380 | |
| 8857 : 64 | NA | NA | Mean : 0.00014 | Mean :12.7 | Mean : 0.0061 | |
| 4042 : 64 | NA | NA | 3rd Qu.: 0.60200 | 3rd Qu.:18.3 | 3rd Qu.: 0.3330 | |
| 3610 : 64 | NA | NA | Max. : 2.69200 | Max. :25.0 | Max. : 0.8310 | |
| (Other):6795 | NA | NA | NA | NA | NA |
par(mfrow =c(2, 2)) # 描画領域を2行3列に分割
hist(MathAchieve$MathAch) # MathAch のヒストグラムの描画
plot(MathAch ~ Minority, data = MathAchieve) # ~ の前が数値、後ろが因子だと箱ひげ図を描画
plot(MathAch ~ I(Minority : Sex), data = MathAchieve) # I(Minority : Sex) は二つの変数を組み合わせてグループ分け
plot(MathAch ~ SES, data = MathAchieve) # ~ の前後がどちらも数値だと散布図
図1: 数学の成績のヒストグラムと独立変数との関連を示す箱ひげ図と散布図
数学の成績はマイノリティよりマジョリティのほうが高く、女性より男性のほうが高い傾向が見られる。また、平均的に見れば、SES が高いほど数学の成績もよい。
1.2.2 Random Intercept Models
次に、ランダム切片モデルを推定したのが、以下の表である。
library(lme4)## Loading required package: Matrix##
## Attaching package: 'lme4'## The following object is masked from 'package:nlme':
##
## lmListm0 <- lmer(MathAch ~ 1 + (1 | School), data = MathAchieve)
m1 <- lmer(MathAch ~ Sex + Minority + SES + (1 | School), data = MathAchieve)
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")."表2: ランダム切片モデルの推定結果"## [1] "表2: ランダム切片モデルの推定結果"screenreg(list(m0, m1))##
## =====================================================
## Model 1 Model 2
## -----------------------------------------------------
## (Intercept) 12.64 *** 14.11 ***
## (0.24) (0.20)
## SexFemale -1.23 ***
## (0.16)
## MinorityYes -2.96 ***
## (0.21)
## SES 2.09 ***
## (0.11)
## -----------------------------------------------------
## AIC 47122.79 46406.38
## BIC 47143.43 46447.66
## Log Likelihood -23558.40 -23197.19
## Num. obs. 7185 7185
## Num. groups: School 160 160
## Var: School (Intercept) 8.61 3.67
## Var: Residual 39.15 35.91
## =====================================================
## *** p < 0.001, ** p < 0.01, * p < 0.058.61 / (8.61 + 39.15) # Model 1 の ICC (Intra-Class Correlation)## [1] 0.1802764切片のみを推定した Model 1 の級内相関係数は 0.18 なので、ある程度、学校によって数学の成績の平均に違いがあることがわかる。個人レベルの独立変数を投入した Model 2 では、すべての独立変数が統計的に有意であり、すでに図 1 で確認したのと同じ傾向を示している。グループによる切片のばらつき (Var: School (Intercept)) を Model 1 と Model 2 で比較すると、8.61 から 3.67 に減少しており、学校によって SESの 平均値や Minority の比率が異なることが、学校間の数学の成績の違いを生み出していると考えられる。
2 Practice: 練習問題
faraway パッケージの jsp データを使い、以下の要領でランダム切片モデルを推定し、結果を解釈しなさい。
- パーソン・イヤー（以下の解説を参照）がレベル1、学校をレベル2とみなし、
- 子供の数学の成績 (math) を従属変数、性別 (gender)、出身階級 (class)、の変数を独立変数とする。
2.0.1 Note: 解説
jsp はロング形式のパネルデータのようであり、同じ個人を 2 ～ 3 年間追跡した結果と思われる。つまり、おそらくこのデータの 1行は個人ではなく、ある年のある個人（パーソン・イヤー）に対応する。それゆえ、本来であれば、上の例題とは少し異なる分析をするのが適切であるが、手間がかかるので、ここではパーソン・イヤーを個人とみなして、例題と同じように分析してよい。