nlme パッケージの MathAchieve データを使い、以下の要領でランダム傾きモデルを推定し、結果を解釈しなさい。
データの概要はランダム切片モデルの練習で見ているので、それらは省略して学校別のデータの分布をまず概観する。学校の数が多く、すべての学校について分布を見るのは現実的でないので、16校を無作為に選んで、数学の成績とその他の変数の関連を見てみよう。
library(nlme) # nlme パッケージの呼び出し
head(MathAchieve)
## 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
# 学校を 16 無作為抽出
id.school <- levels(MathAchieve$School) # 学校のIDの一覧を作る。levels()で因子の値リストを返す
id.school[1:40] # 160校すべてのIDを表示するのは冗長なので最初の40校だけ見る
## [1] "8367" "8854" "4458" "5762" "6990" "5815" "7172" "4868" "7341" "1358" "4383" "2305" "8800" "3088" "8775" "7890"
## [17] "6144" "6443" "5192" "6808" "2818" "9340" "4523" "6816" "2277" "8009" "5783" "3013" "7101" "4530" "9021" "4511"
## [33] "2639" "3377" "6578" "9347" "3705" "3533" "1296" "4350"
sample.school <- sample(id.school, 16) # 16校を無作為に非復元抽出。sample(ベクトル, 抽出数) でベクトルの要素を非復元抽出
sample.school
## [1] "4223" "6089" "3332" "2995" "1317" "9586" "8874" "5640" "1477" "5783" "4383" "8775" "2818" "7734" "4931" "9225"
d0 <- subset(MathAchieve, School %in% sample.school) #subset(データフレーム, 条件式) で条件に合った行だけ抽出
summary(d0$School) # "ベクトル1 %in% ベクトル2" でベクトル1の要素がベクトル2に含まるれかどうかを TRUE or FALSE で返す
## 4383 8775 2818 5783 1317 4223 4931 2995 7734 8874 9225 5640 6089 1477 3332 9586
## 25 48 42 29 48 45 58 46 22 36 36 57 33 62 38 59
library(ggplot2)
p <- ggplot(d0, aes(x = SES, y = MathAch, colour = Sex))
p <- p + geom_point() + facet_wrap(~ School)
p
図1 を見ると学校ごとに SES の傾きが違うことがわかる。また運が良ければ男子校や女子高も含まれていることがわかるだろう。校内の男女差ははっきりしない。同様にして、人種による成績の差を箱ひげ図で示したのが図2 である。
p <- ggplot(d0, aes(x = Minority, y = MathAch))
p <- p + geom_boxplot()
p <- p + facet_wrap(~ School)
p
人種による成績の違いははっきりしているが、人種間の差はある程度学校によって異なることもわかるだろう。また、やはり運が良ければ、マイノリティだけ、マジョリティだけの学校もあることがわかる。次にランダム切片モデルといくつかのランダム傾きモデルを推定した結果が以下である。
library(lme4)
## Loading required package: Matrix
##
## Attaching package: 'lme4'
## The following object is masked from 'package:nlme':
##
## lmList
m1 <- lmer(MathAch ~ Sex + Minority + SES + (1 | School), data = MathAchieve) # ランダム切片
m2 <- lmer(MathAch ~ Sex + Minority + SES + (SES | School), data = MathAchieve) # SES にランダム傾きを仮定
m3 <- lmer(MathAch ~ Sex + Minority + SES + (Sex | School), data = MathAchieve) # Sex にランダム傾きを仮定
m4 <- lmer(MathAch ~ Sex + Minority + SES + (Minority | School), data = MathAchieve)
m5 <- lmer(MathAch ~ Sex + Minority + SES + (SES + Sex | School), data = MathAchieve) # SES と Sex にランダム傾きを仮定
m6 <- lmer(MathAch ~ Sex + Minority + SES + (SES + Minority | School), data = MathAchieve)
m7 <- lmer(MathAch ~ Sex + Minority + SES + (Sex + Minority | School), data = MathAchieve)
AIC(m1, m2, m3, m4, m5, m6, m7) # AIC の計算
## df AIC
## m1 6 46406.38
## m2 8 46405.12
## m3 8 46403.81
## m4 8 46396.15
## m5 11 46403.95
## m6 11 46396.66
## m7 11 46394.56
BIC(m1, m2, m3, m4, m5, m6, m7) # BIC の計算
## df BIC
## m1 6 46447.66
## m2 8 46460.16
## m3 8 46458.85
## m4 8 46451.19
## m5 11 46479.62
## m6 11 46472.34
## m7 11 46470.24
# ランダム切片モデルとそこそこAIC が小さかったモデルをまとめて表示
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").
"表1: ランダム切片モデルとランダム傾きモデルの推定結果"
## [1] "表1: ランダム切片モデルとランダム傾きモデルの推定結果"
screenreg(list(m1, m3, m4, m7))
##
## ===============================================================================================
## Model 1 Model 2 Model 3 Model 4
## -----------------------------------------------------------------------------------------------
## (Intercept) 14.11 *** 14.11 *** 14.14 *** 14.13 ***
## (0.20) (0.21) (0.19) (0.20)
## SexFemale -1.23 *** -1.28 *** -1.25 *** -1.29 ***
## (0.16) (0.18) (0.16) (0.18)
## MinorityYes -2.96 *** -2.96 *** -3.04 *** -3.04 ***
## (0.21) (0.20) (0.25) (0.25)
## SES 2.09 *** 2.10 *** 2.06 *** 2.07 ***
## (0.11) (0.11) (0.11) (0.11)
## -----------------------------------------------------------------------------------------------
## AIC 46406.38 46403.81 46396.15 46394.56
## BIC 46447.66 46458.85 46451.19 46470.24
## Log Likelihood -23197.19 -23193.90 -23190.07 -23186.28
## Num. obs. 7185 7185 7185 7185
## Num. groups: School 160 160 160 160
## Var: School (Intercept) 3.67 4.48 3.10 3.96
## Var: Residual 35.91 35.76 35.70 35.53
## Var: School SexFemale 0.98 1.03
## Cov: School (Intercept) SexFemale -1.19 -1.24
## Var: School MinorityYes 1.81 1.86
## Cov: School (Intercept) MinorityYes 0.75 0.57
## Cov: School SexFemale MinorityYes 0.20
## ===============================================================================================
## *** p < 0.001, ** p < 0.01, * p < 0.05
# confint(m7) # 最後のモデルの係数の区間推定。時間がかかるので実行しない
AIC が最も小さかったのは、表1の Model 4 で、性別と人種にランダム傾きを仮定したモデルであった。SES にランダム傾きを仮定したモデルは一貫して当てはまりが悪かった。切片と性別のランダム効果の共分散を見ると負の値を示しており、切片が大きい学校ほど性別の係数が小さくなる(つまり、この場合は男女差が大きくなる)傾向があることがわかる。いっぽう切片と人種のランダム効果の共分散は正の値を示しており、切片が大きいほど人種の係数は大きくなる(この場合は人種による差が小さくなる)傾向がある。このように性別と人種でかなり傾向が異なる点は興味深い。性別と人種のランダム効果の共分散は正なので、性別の格差が大きい学校では人種による格差も大きい傾向があるということである。なお、BIC はランダム切片モデルが一番小さいので、モデルの選択は慎重におこなったほうがよかろう。
faraway パッケージの jsp データを使い、以下の要領でランダム傾きモデルを推定し、結果を解釈しなさい。
jsp はロング形式のパネルデータのようであり、同じ個人を 2 ~ 3 年間追跡した結果と思われる。つまり、おそらくこのデータの 1行は個人ではなく、ある年のある個人(パーソン・イヤー)に対応する。それゆえ、本来であれば、上の例題とは少し異なる分析をするのが適切であるが、手間がかかるので、ここではパーソン・イヤーを個人とみなして、例題と同じように分析してよい。