nlme パッケージの MathAchieve データを使い、子供がレベル1、学校をレベル2とみなし、子供の数学の成績を従属変数として、以下の要領で複数のモデルを推定し、それらの適合度を比較しなさい。
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.428m0 <- lm(MathAch ~ Sex + Minority + SES, data = MathAchieve) # モデル 0
library(lme4)## Loading required package: Matrix##
## Attaching package: 'lme4'## The following object is masked from 'package:nlme':
##
## lmListm1 <- lmer(MathAch ~ Sex + Minority + SES + (1 | School), data = MathAchieve)
# 愚直な尤度比検定（グループ数が160で違うのはランダム効果の部分だけなので尤度比検定OK）
ll0 <- logLik(m0) # それぞれのモデルの対数尤度
ll1 <- logLik(m1)
ll0; ll1 # ; で区切ると二つのコマンドを同じ行に書ける## 'log Lik.' -23374.79 (df=5)## 'log Lik.' -23197.19 (df=6)dif01 <- -2 * (ll0 - ll1)
dif01 <- as.numeric(dif01)
dif01## [1] 355.19361 - pchisq(dif01, df = 6 - 5) # 尤度比検定の p値## [1] 0# グループ数が160なので必要はないが、この場合は RLRsim パッケージの exactRLRT() も使える
library(RLRsim)## Warning: package 'RLRsim' was built under R version 3.3.2exactRLRT(m1)##
## simulated finite sample distribution of RLRT.
##
## (p-value based on 10000 simulated values)
##
## data:
## RLRT = 365.1, p-value < 2.2e-16#AIC と BIC
AIC(m0, m1) ## df AIC
## m0 5 46759.58
## m1 6 46406.38-2 * c(ll0, ll1) + 2 * (5:6) # これでも上と同じ結果が得られる。## [1] 46759.58 46406.38BIC(m0, m1)## df BIC
## m0 5 46793.98
## m1 6 46447.66-2 * c(ll0, ll1) + log(nobs(ll0)) * (5:6) # これでも上と同じ結果が得られる。## [1] 46793.98 46447.66# 決定係数の比較（この場合同じになるので必要ないが、お遊びでやってみた）
## モデル 0 の決定係数
(r2.m0 <- summary(m0)$r.squared) # カッコで全体をくくると代入した値を表示## [1] 0.171313## モデル 1 の決定係数
### ベースラインモデル（ふつうは切片のみのモデル）のランダム効果の分散の総和
m1.baseline <- lmer(MathAch ~ 1 + (1 | School), data = MathAchieve)
VarCorr(m1.baseline)## Groups Name Std.Dev.
## School (Intercept) 2.9350
## Residual 6.2569(rvar0 <- 2.9350 ^ 2 + 6.2569 ^ 2)## [1] 47.76302var(MathAchieve$MathAch) # 切片のみのモデルのランダム効果の分散の総和は従属変数の分散に一致## [1] 47.31026### モデル 1 のランダム効果の分散の総和
VarCorr(m1)## Groups Name Std.Dev.
## School (Intercept) 1.9167
## Residual 5.9924(rvar1 <- 1.9167 ^ 2 + 5.9924 ^ 2)## [1] 39.5826r2.m1 <- 1 - rvar1 / rvar0 # モデル 1 の決定係数
r2.m0; r2.m1 # 推定法の違いがあるのでぴったりとは一致しない## [1] 0.171313## [1] 0.1712711尤度比検定でも AIC でも BIC でもモデル 0 よりもモデル 1 のほうが当てはまりがよい。固定効果の部分はまったく同じなので、決定係数で比較しても無意味である。
m2 <- update(m1, ~. - SES) # update() はモデルを少しだけ修正するのに使う。
m2 # この場合 m1 の独立変数から SES を取り除いたモデルを推定## Linear mixed model fit by REML ['lmerMod']
## Formula: MathAch ~ Sex + Minority + (1 | School)
## Data: MathAchieve
## REML criterion at convergence: 46759.32
## Random effects:
## Groups Name Std.Dev.
## School (Intercept) 2.464
## Residual 6.121
## Number of obs: 7185, groups: School, 160
## Fixed Effects:
## (Intercept) SexFemale MinorityYes
## 14.382 -1.389 -3.701# REMLで固定効果の有無を検討する際には意味がないがいちおう尤度比検定
anova(m2, m1, refit=FALSE)## Data: MathAchieve
## Models:
## m2: MathAch ~ Sex + Minority + (1 | School)
## m1: MathAch ~ Sex + Minority + SES + (1 | School)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## m2 5 46769 46804 -23380 46759
## m1 6 46406 46448 -23197 46394 364.94 1 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Kenward Roger の自由度で F 検定（REML でグループ数が少ないときにも使える！）
library(pbkrtest)
KRmodcomp(m1, m2)## F-test with Kenward-Roger approximation; computing time: 0.99 sec.
## large : MathAch ~ Sex + Minority + SES + (1 | School)
## small : MathAch ~ Sex + Minority + (1 | School)
## stat ndf ddf F.scaling p.value
## Ftest 390.19 1.00 6852.19 1 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# AIC と BIC（REML なのでこれを使うことに批判もあろうが）
AIC(m2, m1); BIC(m2, m1)## df AIC
## m2 5 46769.32
## m1 6 46406.38## df BIC
## m2 5 46803.72
## m1 6 46447.66# 決定係数
VarCorr(m2)## Groups Name Std.Dev.
## School (Intercept) 2.4637
## Residual 6.1211rvar2 <- 2.4637 ^ 2 + 6.1211 ^ 2
r2.m2 <- 1 - rvar2 / rvar0
r2.m2; r2.m1## [1] 0.08846466## [1] 0.1712711どの基準で比較してもモデル 2 よりもモデル 1 のほうが当てはまりがよい。
m3 <- lmer(MathAch ~ Sex + Minority + SES + (Minority | School), data = MathAchieve)
# 尤度比検定
anova(m1, m3, refit = FALSE)## Data: MathAchieve
## Models:
## m1: MathAch ~ Sex + Minority + SES + (1 | School)
## m3: MathAch ~ Sex + Minority + SES + (Minority | School)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## m1 6 46406 46448 -23197 46394
## m3 8 46396 46451 -23190 46380 14.235 2 0.0008106 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# AIC, BIC
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, m3))##
## =================================================================
## Model 1 Model 2
## -----------------------------------------------------------------
## (Intercept) 14.11 *** 14.14 ***
## (0.20) (0.19)
## SexFemale -1.23 *** -1.25 ***
## (0.16) (0.16)
## MinorityYes -2.96 *** -3.04 ***
## (0.21) (0.25)
## SES 2.09 *** 2.06 ***
## (0.11) (0.11)
## -----------------------------------------------------------------
## AIC 46406.38 46396.15
## BIC 46447.66 46451.19
## Log Likelihood -23197.19 -23190.07
## Num. obs. 7185 7185
## Num. groups: School 160 160
## Var: School (Intercept) 3.67 3.10
## Var: Residual 35.91 35.70
## Var: School MinorityYes 1.81
## Cov: School (Intercept) MinorityYes 0.75
## =================================================================
## *** p < 0.001, ** p < 0.01, * p < 0.05# ふつうの決定係数を比較しても無意味なので、個人レベルのランダム効果の分散を比較
1 - 35.91 / 6.2569 ^ 2 # モデル 1 の個人レベルの決定係数## [1] 0.082730451 - 35.70 / 6.2569 ^ 2 # モデル 1 の個人レベルの決定係数## [1] 0.0880946BIC ではモデル 1 のほうがあてはまりがよいが、その他の指標ではモデル 3 のほうがよい。
faraway パッケージの jsp データを使い、パーソン・イヤーがレベル1、学校をレベル2とみなし、子供の数学の成績 (math) を従属変数として、以下の要領で複数のモデルを推定し、それらの適合度を比較しなさい。
モデル 1 と モデル 1 に raven のランダム傾きを仮定したモデル（モデル 3）の比較。