1.と 2. をまとめて行った結果が以下の表1である。
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.428fixed1 <- lm(MathAch ~ -1 + School + Sex + Minority + SES, data = MathAchieve) # -1 はエラーが出ないおまじない。グループ変数を最初に指定すること
# summary(fixed1) 学校ごとの切片がすべて表示されるので省略
library(lme4)## Loading required package: Matrix##
## Attaching package: 'lme4'## The following object is masked from 'package:nlme':
##
## lmListrand.REML <- lmer(MathAch ~ Sex + Minority + SES + (1 | School), data = MathAchieve)
rand.ML <- lmer(MathAch ~ Sex + Minority + SES + (1 | School), data = MathAchieve, REML = FALSE)
library(texreg)## 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 固定効果モデル、ランダム切片モデル（REML と ML）の推定結果"## [1] "表1 固定効果モデル、ランダム切片モデル（REML と ML）の推定結果"screenreg(list(fixed1, rand.REML, rand.ML),
custom.model.names = c("Fixed", "REML", "ML"), # 表頭に示される各モデルの名前を指定
omit.coef = "School") # 固定効果モデルの学校ごとの切片の表示を省略##
## ==================================================================
## Fixed REML ML
## ------------------------------------------------------------------
## SexFemale -1.16 *** -1.23 *** -1.23 ***
## (0.17) (0.16) (0.16)
## MinorityYes -2.92 *** -2.96 *** -2.96 ***
## (0.22) (0.21) (0.21)
## SES 1.91 *** 2.09 *** 2.09 ***
## (0.11) (0.11) (0.11)
## (Intercept) 14.11 *** 14.11 ***
## (0.20) (0.20)
## ------------------------------------------------------------------
## R^2 0.83
## Adj. R^2 0.83
## Num. obs. 7185 7185 7185
## RMSE 5.99
## AIC 46406.38 46398.84
## BIC 46447.66 46440.11
## Log Likelihood -23197.19 -23193.42
## Num. groups: School 160 160
## Var: School (Intercept) 3.67 3.64
## Var: Residual 35.91 35.90
## ==================================================================
## *** p < 0.001, ** p < 0.01, * p < 0.05固定効果モデルとランダム切片モデルの係数を比較すると、性別は 0.07、マイノリティのダミーは 0,04 の差で、それぞれの標準誤差よりもずっと小さく、ほとんど同じ値であると考えられる。SES の係数は 0.18 の違いで、標準誤差よりもずっと大きいので、学校によって SES が異なる（つまり両者に相関がある）と考えたほうがよいかもしれない。
REML と ML の結果を比較すると、係数もその標準誤差も同じであるが、REML のほうが ML よりもランダム効果の分散を少しだけ大きく推定していることがわかる。
3.次に Raudenbush and Bryk (2002) の概算法で自由度を計算するために、個人の数 (N) とグループの数 (M)、そして推定したパラメータの数 (k) をまず求める。
N <- nobs(rand.REML) # レベル1のサンプル・サイズ
M <- ngrps(rand.REML) # レベル2のサンプル・サイズ
k <- length(fixef(rand.REML)) # coef() で固定効果の値のベクトルを返す。length(ベクトル) でベクトルの要素の数を返す今、検定しようとしているモデルに投入している独立変数は、すべて個人レベルの変数の主効果だから、自由度は \(N - k\) （ちなみに切片はグループレベルの変数なので自由度は\(M - 1\)）である。それゆえ、独立変数の t検定の p値は以下のようになる。
coef.REML <- summary(rand.REML)$coefficients
coef.REML## Estimate Std. Error t value
## (Intercept) 14.114511 0.1970283 71.636967
## SexFemale -1.229794 0.1627085 -7.558267
## MinorityYes -2.961472 0.2057554 -14.393164
## SES 2.089424 0.1057058 19.766408t.REML <- coef.REML[-1, 3] # t値だけとってくる
p.REML <- pt(abs(t.REML), df = N - k) # abs() で絶対値を返す。
2 * (1 - p.REML) # pt() は 帰無仮説のもとで t値以下の値をとる確率を返すので、両側検定のために 1 から引いて2倍## SexFemale MinorityYes SES
## 4.596323e-14 0.000000e+00 0.000000e+00すべての変数が 0.1% 水準で有意である。
次に一般化 Satterthwaite 法を使って p値を計算すると以下のようになる。
library(lmerTest)##
## Attaching package: 'lmerTest'## The following object is masked from 'package:lme4':
##
## lmer## The following object is masked from 'package:stats':
##
## steprand.REML2 <- lmer(MathAch ~ Sex + Minority + SES + (1 | School), data = MathAchieve)
summary(rand.REML2)## Linear mixed model fit by REML
## t-tests use Satterthwaite approximations to degrees of freedom ['lmerMod']
## Formula: MathAch ~ Sex + Minority + SES + (1 | School)
## Data: MathAchieve
##
## REML criterion at convergence: 46394.4
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.2427 -0.7216 0.0342 0.7620 2.8631
##
## Random effects:
## Groups Name Variance Std.Dev.
## School (Intercept) 3.674 1.917
## Residual 35.909 5.992
## Number of obs: 7185, groups: School, 160
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 14.1145 0.1970 263.0000 71.637 < 2e-16 ***
## SexFemale -1.2298 0.1627 6641.0000 -7.558 4.62e-14 ***
## MinorityYes -2.9615 0.2058 4830.0000 -14.393 < 2e-16 ***
## SES 2.0894 0.1057 6832.0000 19.766 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) SexFml MnrtyY
## SexFemale -0.434
## MinorityYes -0.292 0.014
## SES -0.078 0.058 0.195さらに Wald 検定すると以下のようになる。
library(aod)
wald.test(b = fixef(rand.REML), Sigma = vcov(rand.REML), Terms = 2) # 女性の係数の Wald 検定## Wald test:
## ----------
##
## Chi-squared test:
## X2 = 57.1, df = 1, P(> X2) = 4.1e-14wald.test(b = fixef(rand.REML), Sigma = vcov(rand.REML), Terms = 3) # マイノリティの係数の Wald 検定## Wald test:
## ----------
##
## Chi-squared test:
## X2 = 207.2, df = 1, P(> X2) = 0.0wald.test(b = fixef(rand.REML), Sigma = vcov(rand.REML), Terms = 4) # SESの係数の Wald 検定## Wald test:
## ----------
##
## Chi-squared test:
## X2 = 390.7, df = 1, P(> X2) = 0.0最後に切片の分散の区間推定を行った結果は以下の通りである。
# プロファイル区間推定を使った片側検定
confint(rand.REML, parm = 1, level = .9, oldNames = FALSE) # 片側 5% 水準で検定したい場合## Computing profile confidence intervals ...## 5 % 95 %
## sd_(Intercept)|School 1.697002 2.146152confint(rand.REML, parm = 1, level = .98, oldNames = FALSE) # 片側 1% 水準で検定したい場合## Computing profile confidence intervals ...## 1 % 99 %
## sd_(Intercept)|School 1.61723 2.255471confint(rand.REML, parm = 1, level = .998, oldNames = FALSE) # 片側 0.1% 水準で検定したい場合## Computing profile confidence intervals ...## 0.1 % 99.9 %
## sd_(Intercept)|School 1.532205 2.38615# パラメトリック・ブートストラップ区間推定を用いた片側検定
confint(rand.REML, parm = 1, level = .9, oldNames = FALSE, method ="boot") # 片側 5% 水準で検定したい場合## Computing bootstrap confidence intervals ...## 5 % 95 %
## sd_(Intercept)|School 1.695102 2.137305## 以下同様だが、計算に時間がかかるので実行しない
# confint(rand.REML, parm = 1, level = .98, oldNames = FALSE, method ="boot") # 片側 1% 水準で検定したい場合
# confint(rand.REML, parm = 1, level = .998, oldNames = FALSE, method ="boot") # 片側 0.1% 水準で検定したい場合
# セミパラメトリック・ブートストラップ
confint(rand.REML, parm = 1, level = .9, oldNames = FALSE,
method ="boot", type = "semiparametric", use.u = TRUE) # 片側 5% 水準で検定したい場合## Computing bootstrap confidence intervals ...## Warning in optwrap(object@optinfo$optimizer, ff, x0, lower = lower, control = control$optCtrl, : convergence code 3 from
## bobyqa: bobyqa -- a trust region step failed to reduce q## 5 % 95 %
## sd_(Intercept)|School 1.544439 1.80447## 以下同様だが、計算に時間がかかるので実行しない
# confint(rand.REML, parm = 1, level = .98, oldNames = FALSE,
# method ="boot", type = "semiparametric", use.u = TRUE) # 片側 1% 水準で検定したい場合
# confint(rand.REML, parm = 1, level = .998, oldNames = FALSE,
# method ="boot", type = "semiparametric", use.u = TRUE) # 片側 0.1% 水準で検定したい場合いずれも信頼区間の下限値が 0 よりも大きいので、切片の分散は統計的に有意である。