そもそも私たちはどのような状況下でマルチレベル・モデルを使い、どのような状況下で使うべきではないのか。一概には言えないが、幾つかの原則をのべておこう。
グループ数が少ない場合、グループ・レベルのランダム効果の分散の推定が不正確になる。それゆえ、特にグループ・レベルのランダム効果そのものに注目するならば、グループ数を十分に確保する必要がある。どのぐらいあれば十分かは、以下のような諸要因によっても変わってくるので、一概には言えない。
ただし、各グループに最低でも数十人以上の個人が含まれ、バランスが取れたデータで著しい独立変数の偏りや多重共線性はないという前提で考えると、
ただし、これらは絶対的基準ではないので、ケースバイケースで柔軟に対応すべきであろう。
もしも、単にグループの効果を統制して、個人レベルの変数の固定効果を推定／検定したいというだけならば、マルチレベルモデルではなく、固定効果モデルのほうが適当な場合も多い。固定効果モデルとは、パネル・データ分析の用語法であるが、パネル・データにかぎらず、マルチレベル・データの多くに応用が利く。固定効果モデルとは、 \[ Y_{ij} = \beta_{0j} + \beta_1 X_{ij} + R_{ij} \] というモデルをたてるが、この点では、ランダム切片モデルと同じである。ランダム切片モデルと違うのは、\(\beta_{0j}\)はグループ・レベルのランダム効果によって変動するのではなく、\(\beta_1\)とおなじ固定効果であり、その他の説明変数と相関していると仮定されている点である（ランダム効果モデルでは独立と仮定する）。そのため、グループの効果を完全に統制できる。
固定効果モデルを推定するには、単純にグループのダミー変数をモデルに投入して OLS で推定すればよい。25 の企業における賃金と教育年数の例で言えば、企業ダミーを 24 作って モデルに投入すればよい。
fixed1 <- lm(wage ~ id.firm + educ, data = data1)
summary(fixed1)##
## Call:
## lm(formula = wage ~ id.firm + educ, data = data1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -139.399 -36.802 -1.402 30.805 135.601
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 63.1874 32.0933 1.969 0.051764 .
## id.firmb -11.5990 35.7230 -0.325 0.746100
## id.firmc -36.4008 35.7070 -1.019 0.310484
## id.firmd -22.0026 35.9151 -0.613 0.541527
## id.firme -38.4041 36.2690 -1.059 0.292236
## id.firmf -62.2036 36.1333 -1.722 0.088283 .
## id.firmg -10.8023 35.8717 -0.301 0.763941
## id.firmh -26.8018 35.7986 -0.749 0.455823
## id.firmi 2.5931 37.3217 0.069 0.944749
## id.firmj 19.9951 36.5053 0.548 0.585109
## id.firmk -8.2031 36.0152 -0.228 0.820297
## id.firml 0.7972 35.9629 0.022 0.982360
## id.firmm 10.3987 35.7437 0.291 0.771717
## id.firmn -35.4049 36.5053 -0.970 0.334482
## id.firmo 6.5954 36.4222 0.181 0.856674
## id.firmp -5.8062 36.9844 -0.157 0.875573
## id.firmq -2.0041 36.2690 -0.055 0.956045
## id.firmr -5.2075 37.5666 -0.139 0.890032
## id.firms 39.1938 36.9844 1.060 0.291842
## id.firmt 7.3954 36.4222 0.203 0.839516
## id.firmu 40.9931 37.3217 1.098 0.274706
## id.firmv 68.7936 37.0928 1.855 0.066625 .
## id.firmw 154.5925 37.5666 4.115 8.01e-05 ***
## id.firmx 122.5943 36.7802 3.333 0.001209 **
## id.firmy 137.9915 38.1034 3.622 0.000464 ***
## educ 25.0013 2.0235 12.356 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 56.43 on 99 degrees of freedom
## Multiple R-squared: 0.8261, Adjusted R-squared: 0.7822
## F-statistic: 18.82 on 25 and 99 DF, p-value: < 2.2e-16library(car) # 多重共線性だけを病的に気にする社会学者が多いので念のため
vif(fixed1) # VIF をチェック## GVIF Df GVIF^(1/(2*Df))
## id.firm 1.589053 24 1.009695
## educ 1.589053 1 1.260577ただし、グループ数が多いと推定に著しく時間がかかる場合もあるので、その場合はグループ内センタリングした変数を使った推定が必要になる。これについてはセンタリングの節でもう少し詳しく説明する。
国際比較分析が典型的であるが、グループが無作為抽出されていないデータをマルチレベル分析したいという場合はある。グループも個人も無作為抽出されていないとマルチレベル分析してはいけないという説を耳にすることもあるが、この説は必ずしも正しくはない。
もしも母集団に属するすべてのグループのすべての個人に分析結果を一般化したいならば、グループも個人も無作為抽出しなければならない。しかし、サンプルとして得られたグループ以外への一般化をしないならば、グループが無作為に抽出されている必要はない。例えば、OECD に加盟する 35カ国のデータがあり、その外に一般化しないならば、このデータを使ってマルチレベル分析するのは何の問題もない。
従属変数がグループから一切影響を受けない場合、例えデータがマルチレベル構造をしていたとしても、通常の回帰分析で問題なく推定できる。多段抽出の場合、グループの影響を受けなかったとしても、単純無作為抽出に比べて標準誤差が大きくなるが、グループ数が非常に多い場合、慣習的には単純無作為抽出とみなして分析されることが多い。
しかし、厳密には分析してみなければ従属変数がグループから影響を受けるかどうかわからない場合も多い。それゆえ、理論上（あるいは先行研究や常識などから）、グループが従属変数に何らかの影響を持っていると考えられる場合は、マルチレベル分析を行うことになる。
一つの目安として、従属変数の級内相関係数 (Intra Class Correlation: ICC) を計算することもある。級内相関係数 (ICC) とは、従属変数の分散を \(\mathrm{var}(Y_{ij})\)、従属変数の各グループの平均値 (\(\bar Y_{\cdot j}\)) の分散を \(\tau_{00}\)、とすると、 \[ \mathrm{ICC} = \frac{\tau_{00}}{\mathrm{var}(Y_{ij})} \] である。これがほとんどゼロの場合、グループによる従属変数の違いはあまり無い場合が多いので、マルチレベル分析をする必要性を判断する基準の一つとして用いられている。しかし、級内相関係数がゼロでも、他の独立変数を投入すると、グループによる違いが出ることもある。つまり一種の擬似無相関が生じていることもあるので、級内相関係数はあくまで参考程度にとどめ、理論的な予測を優先すべきであろう。
R でこれを計算する場合、切片のみのランダム切片モデルから、\(\tau_{00}\)と\(\sigma^2\)をもとめ、式 (1.10) のように \(\mathrm{var}(Y_{ij}) = \tau_{00} + \sigma^2\) なので、これらから級内相関係数を計算するのが、意外と簡単かもしれない。例えば、以下のように。
library(lme4)
rand.itcpt0 <- lmer(wage ~ 1 + (1 | id.firm), data = data1)
rand.itcpt0## Linear mixed model fit by REML ['lmerMod']
## Formula: wage ~ 1 + (1 | id.firm)
## Data: data1
## REML criterion at convergence: 1511.176
## Random effects:
## Groups Name Std.Dev.
## id.firm (Intercept) 82.63
## Residual 89.51
## Number of obs: 125, groups: id.firm, 25
## Fixed Effects:
## (Intercept)
## 404.782.63^2 / (82.63^2 + 89.51^2) # 級内相関係数## [1] 0.4600963おなじことだが、以下のようにすれば丸め誤差が減るので、さらに正確であろう。
tau00 <- as.numeric(summary(rand.itcpt0)$varcor[[1]]) # レベル2のランダム効果の分散
sigma2 <- summary(rand.itcpt0)$sigma^2 # レベル1のランダム効果の分散
tau00 / (tau00 + sigma2) # 級内相関係数## [1] 0.4600767また、級内相関係数専用の関数もいくつかある。multilevel パッケージの ICC1() 関数は一元配置の分散分析の結果を引数として、級内相関係数を返す。まず一元配置の分散分析は
aov(従属変数 ~ グループを示す因子, data = データフレーム名)
で、この結果 (aov オブジェクトという) を引数として、
ICC1(aovオブジェクト)
とすればやはり級内相関係数を返す。以下は計算例。
aov0 <- aov(wage ~ id.firm, data = data1)
summary(aov0)## Df Sum Sq Mean Sq F value Pr(>F)
## id.firm 24 1011602 42150 5.261 1.53e-09 ***
## Residuals 100 801245 8012
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1library(multilevel)
ICC1(aov0)## [1] 0.4600767このセクションの議論をまとめると、
場合にマルチレベル分析を行うべきであろう。
マルチレベル・モデルの推定には、主に最尤法 (Maximum Likelihood: ML) と制限最尤法 (REstricted Maximum Likelihood: REML) の二種類が主に使われている（ベイズ推定もあるが、計算に時間がかかるらしい）。ML はすべてのパラメータを最尤法で推定するが、REML はまず一般化最小二乗法で固定効果を推定し、その推定値を尤度関数に代入して、ランダム効果の分散と共分散を最尤推定する。第一段階で固定効果の推定値は決定され、第二段階の最尤推定の際には、ランダム効果の分散と共分散だけが推定される。
どちらの推定法を使うべきか？ 結論から言えば、殆どの場合で REML を使うべきである。なぜなら REML のほうが固定効果の標準誤差が小さく、ランダム効果の分散・共分散の推定もより正確であると言われている。特にグループ数が少ない場合、 ML はグループ・レベルの分散・共分散を過小に推定すると言われており、REML の利用が推奨される。
それでは例外的に ML が推奨されるのはどのような場合か？ ML の場合、モデルを比較する際に尤度比検定や AIC を使うことができるが、REML の場合はそれらを用いることはしばしば無意味である。それゆえ、ML のほうがモデルの比較がやりやすいというメリットがある。ところで、 グループ数が非常に多い場合、ML の推定結果と REML の推定結果はほとんど同じになる。それゆえ、グループ数が多く ML と REML の推定結果がほぼおなじでなおかつ尤度比検定のような尤度を使った計算がしたいならば、ML を行うことが推奨される。
とはいえ、ほとんどの場合、REML で十分と思われる。社会学に限って言えば、確かに AIC のようなモデルの適合度が使えたほうがいいのだが、比較するモデルが階層的な関係にあれば Wald 検定やブートストラップ検定をすればいいし、おおまかに適合度がわかりさえすればいいなら決定係数を使うことができる。さきほど REML の場合は尤度を使ったモデルの比較は無意味と述べたが、このあたりの議論は論争が続いており、あまりはっきりした結論が出ているわけではない。それゆえ、あえて ML を使う必要はほとんどないと思われる。
R の lmer() 関数の場合、REML がデフォルトであるが、ML で推定する場合、
REML = FALSE
という引数を追加で指定すればよい。賃金と教育年数の例で計算してみよう。
data1$educ.c <- data1$educ - mean(data1$educ) # 教育年数のセンタリング
data1$size.c <- (data1$firm.size - mean(data1$firm.size))/100 # 企業規模をセンタリングして100人単位に
r.slope1 <- lmer(wage ~ educ.c * size.c + (educ.c | id.firm), data =data1)
r.slope2 <- lmer(wage ~ educ.c * size.c + (educ.c | id.firm), data =data1, REML = FALSE)
library(texreg)
screenreg(list(r.slope1, r.slope2), custom.model.names = c("Model 2.1 REML", "Model 2.1 ML"))##
## =============================================================
## Model 2.1 REML Model 2.1 ML
## -------------------------------------------------------------
## (Intercept) 391.76 *** 391.79 ***
## (6.11) (5.96)
## educ.c 25.05 *** 25.08 ***
## (1.98) (1.92)
## size.c 2.80 *** 2.81 ***
## (0.44) (0.43)
## educ.c:size.c 0.56 *** 0.56 ***
## (0.12) (0.12)
## -------------------------------------------------------------
## AIC 1359.43 1365.33
## BIC 1382.06 1387.95
## Log Likelihood -671.72 -674.66
## Num. obs. 125 125
## Num. groups: id.firm 25 25
## Var: id.firm (Intercept) 262.00 231.64
## Var: id.firm educ.c 24.18 19.71
## Cov: id.firm (Intercept) educ.c 79.59 67.56
## Var: Residual 2598.79 2545.25
## =============================================================
## *** p < 0.001, ** p < 0.01, * p < 0.05固定効果の推定値はほとんど同じだが、ランダム効果の分散が全般に REML のほうが大きくなっている。 そのため係数の標準誤差も REML のほうがやや大きい。REML のほうが保守的な（つまり、有意な結果が出にくい）分析法と言える。
結論だけ言えば、自分が使っている統計ソフトのデフォルトの検定結果を使っておけばほぼ問題ない。しかし、マルチレベル分析における検定に関しては、専門家の間で論争があり、単純な t 検定では批判を受ける可能性もないとは言えない。
すでに述べたように、マルチレベル分析では、個人を単位として考えるか、グループを単位として考えるかで、サンプルサイズが違ってくる。そのため、モデル（の残差）の自由度がはっきりしない。通常の回帰分析では、自由度は、 \[ N - \mathrm{推定したパラメータ数} \] であった。それゆえ、個人とグループのどちらを単位とするかで \(N\) の値が変わってくるというわけである。
しかし、もしもグループ数が 100 を超えていれば、どちらを単位にしても 自由度は 100 に近いので、t 値は正規分布におおむね近似すると考えられる（ただし、推定するパラメータが多いほど正規分布からは離れる）。それゆえ、グループ数が 100 を超えるような場合は、
といった基準で考えればよい。つまり、グループ数が多ければ、自由度がわからなくても係数が有意かどうかは判断できる。
しかし、実際にはグループ数が少ないことも珍しくない。そういう場合は上のような推論は使えないが、例えば t 値が 5 を超えていれば、自由度が 2 であっても両側 5% 水準で有意なので、自由度をどう考えるにせよ、ほとんどの場合で有意な結果だと考えて差し支え無いだろう（ただし、グループ数が少ないのにたくさんのグループレベルの独立変数を投入しているような場合は除く）。それゆえ、自由度について深く考えなくてもおおむね判断できる場合も多い。
実際、t値のみを表記してアスタリスクはつけず、何パーセント水準で有意かは明示しない、というマルチレベル分析のテキストは多い。近年、有意性検定に固執することは危険であるという批判が声高に叫ばれているので、t 値のみを表記するというやり方でも、私は構わないと思う。
しかし、学界の慣習上、何パーセント水準で有意かはっきりさせたい場合もある。社会学でも係数が両側 5% 水準で有意かどうかは重要な意味を持ち続けているので、p 値を計算しないのは勇気がいるだろう。そこで、自由度を概算する方法がいくつか提唱されており、それらのうちのどれかを用いて自由度と p 値を計算すればよい。一番簡単なのは、Reudenbush and Bryk (2002) が提唱している自由度の計算法で、個人レベルの変数に関する t 値の自由度は、 \[ N - \mathrm{切片を含む推定したすべての固定効果の数} \] で、これは OLS での自由度の数え方と同じである。いっぽうグループレベルの変数に関する t 値の自由度は、グループ数を \(M\) とすると、 \[ M - \mathrm{切片を含む推定したグループレベルの固定効果の数} \] である。どの変数が個人レベルでどの変数がグループレベルかについては、
と分類する。例えば、Model 2.1 の係数の自由度と p 値は、以下のように計算できる。
# さきほど推定したランダム係数モデルの推定結果、r.slpe1 を使う
M <- ngrps(r.slope1) # ngrps(MerModオブジェクト) で推定に用いたグループ数を返す
N <- nobs(r.slope1) # ngrps(MerModオブジェクト) で推定に用いた観察数を返す
fixef(r.slope1) # fixef() は固定効果の推定値だけを返す## (Intercept) educ.c size.c educ.c:size.c
## 391.7609482 25.0549777 2.8020757 0.5630951この例では、“educ.c” の係数だけが個人レベルで、その他はすべてグループレベルなので、以下のように計算できる。
df1 <- c(M - 3, N - 4, M - 3, M - 3) # 自由度
df1## id.firm id.firm id.firm
## 22 121 22 22(coef.tab1 <- summary(r.slope1)$coefficients) # 係数の表だけをとってくる## Estimate Std. Error t value
## (Intercept) 391.7609482 6.1146332 64.069411
## educ.c 25.0549777 1.9798499 12.654988
## size.c 2.8020757 0.4407357 6.357723
## educ.c:size.c 0.5630951 0.1242731 4.531109(t1 <- coef.tab1[, "t value"]) # t値だけとってくる## (Intercept) educ.c size.c educ.c:size.c
## 64.069411 12.654988 6.357723 4.5311092 * (1 - pt(t1, df = df1)) # 両側検定の p 値。 pt(t値, 自由度) で t値に対応する累積確率を返す## (Intercept) educ.c size.c educ.c:size.c
## 0.000000e+00 0.000000e+00 2.132121e-06 1.648607e-04他にも自由度の概算法はいくつかある。これまで texreg パッケージの screenreg() 関数で作った表でも検定結果が表示されてきたが、これらは自由度が無限大とみなしている、つまり正規分布に近似するという仮定のもとに p 値が計算されている。これはいささか有意な結果が出やすすぎるかもしれない。
また、lmerTest パッケージの lmer() 関数を使うと、一般化 Satterthwaite 法で自由度を計算する。引数の指定は、lme4 パッケージの lmer() とほとんど同じである。
library(lmerTest)
r.slope3 <- lmer(wage ~ educ.c * size.c + (educ.c | id.firm), data =data1)
summary(r.slope3)## Linear mixed model fit by REML t-tests use Satterthwaite approximations
## to degrees of freedom [lmerMod]
## Formula: wage ~ educ.c * size.c + (educ.c | id.firm)
## Data: data1
##
## REML criterion at convergence: 1343.4
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.61452 -0.74551 -0.05109 0.52242 2.88719
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## id.firm (Intercept) 262.00 16.186
## educ.c 24.18 4.917 1.00
## Residual 2598.79 50.978
## Number of obs: 125, groups: id.firm, 25
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 391.7609 6.1146 36.0200 64.069 < 2e-16 ***
## educ.c 25.0550 1.9798 33.2600 12.655 2.89e-14 ***
## size.c 2.8021 0.4407 41.1500 6.358 1.32e-07 ***
## educ.c:size.c 0.5631 0.1243 24.3100 4.531 0.000133 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) educ.c size.c
## educ.c 0.266
## size.c -0.039 -0.356
## educ.c:sz.c -0.317 0.045 0.321Type III の分散分析で検定するという方法もある。その場合、Kenward-Roger 法を選べる。
anova(r.slope3, ddf="Kenward-Roger")## Analysis of Variance Table of type III with Kenward-Roger
## approximation for degrees of freedom
## Sum Sq Mean Sq NumDF DenDF F.value Pr(>F)
## educ.c 400191 400191 1 31.061 153.991 1.420e-13 ***
## size.c 101035 101035 1 32.185 38.878 5.371e-07 ***
## educ.c:size.c 48701 48701 1 19.604 18.740 0.0003395 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1ランダム効果の分散（\(\tau_{00}\)）が有意にゼロから離れているかどうかが問題になる場合もある。その場合、\(\tau_{00} = 0\) を帰無仮説として検定したい。もしもランダム効果の分散がゼロから十分に離れていれば、\(\tau_{00}\) は正規分布することが知られている。それゆえ、標準誤差さえ求められれば、いちおう検定はできる。とはいえ、ゼロから十分に離れているかどうかわからないから検定するのであるから、この方法は簡便で計算も早いものの、使っていいのかどうか確信の持てない場合も多いだろう。
REMLで推定し、グループ数が40以上あるならば、ランダム効果をゼロに固定したモデルを作り、尤度比検定するという方法がまず考えられる。ランダム切片モデルなら、同じ固定効果を仮定したモデルを OLS で推定すれば、切片のランダム効果をゼロに固定したモデルが得られる。傾きのランダム効果を尤度比検定したい場合は、同じ固定効果を仮定したランダム傾きモデルとランダム切片モデルを推定し、両者の尤度比を検定すればよい。
尤度比検定の基本的な考え方を解説しよう。今、モデル1 に含まれるパラメータ \(\beta_1, \ \beta_2 \ldots, \beta_k\) （上の例なら\(\tau_{00}\)や\(\tau_{01}\)）がゼロかどうか検定したいとしよう。ベータは一つだけでもいいし、複数あってもよい。つまり帰無仮説は、 \[ \beta_1 = \ \beta_2 = \ldots, = \beta_k = 0 \] である。言い換えればこれらはすべてゼロだと仮定されている。それゆえ対立仮説は \(\beta_1, \ \beta_2 \ldots, \beta_k\) のうちどれか一つ以上はゼロではない、ということになる。 これらのパラメータをゼロに固定したモデルをモデル0 と呼ぶことにする。モデル0、モデル1の尤度をそれぞれ\(L_0\)、\(L_1\)とすると、 \[ -2 (\log L_0 - \log L_1) = -2 log\frac{L_0}{L_1} \] は帰無仮説のもとでカイ二乗分布に近似する。自由度はゼロに固定したパラメータの数である。上の式のように二つのモデルの尤度の比をもとにして統計量を計算するので、尤度比検定と呼ばれる。尤度そのものもよりも対数尤度を計算に使うほうが多い。
以下は、ランダム切片モデルの切片のランダム効果の尤度比検定の例。
r.itcpt1 <- lmer(wage ~ educ + (1 | id.firm), data = data1) # モデル1（ランダム切片モデル）
ols1 <- lm(wage ~ educ, data = data1) # モデル0（切片のランダム効果の分散をゼロと仮定)
L1 <- logLik(r.itcpt1) # logLik() でモデルの対数尤度を返す
L1 # deviance() は -2*対数尤度を返すはずだが、## 'log Lik.' -694.368 (df=4)L0 <- logLik(ols1) # lm lmer で計算の仕方が違うのでここでは使っていない
L0## 'log Lik.' -709.9103 (df=3)# -2 対数尤度
dif01 <- -2 * as.numeric(L0 - L1) # 変な属性がついててうっとうしいので
dif01 # as.numeric() でただの数値に変換## [1] 31.08477# -2対数尤度に対応する p値の計算
1 - pchisq(dif01, df = 1) # pchisq() で指定したカイ二乗値と自由度に対応する下側の累積確率を返す## [1] 2.470022e-08p値が 0.05 未満なので切片のランダム効果の分散はゼロより大きいと判断できる。
以下は、ランダム切片モデルとランダム傾きモデルを尤度比検定をして、傾きのランダム効果の分散、および傾きと切片のランダム効果の共分散がゼロかどうか検定した例。
# 以下のモデル2と上の例のモデル1（ランダム切片モデル）を尤度比検定
r.slope1 <- lmer(wage ~ educ + (educ | id.firm), data = data1) # モデル2（ランダム傾きモデル）
L2 <- logLik(r.slope1)
L2## 'log Lik.' -684.3534 (df=6)dif12 <- -2 * as.numeric(L1 - L2)
dif12## [1] 20.029191 - pchisq(dif12, 2) # 傾きの分散と傾きと切片の共分散をゼロに固定しているので自由度は 2## [1] 4.47422e-05# 後で述べるようにMerMod オブジェクトどうしならば以下のようにしても同じことができる
anova(r.itcpt1, r.slope1, refit = FALSE)## Data: data1
## Models:
## object: wage ~ educ + (1 | id.firm)
## ..1: wage ~ educ + (educ | id.firm)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## object 4 1396.7 1408.0 -694.37 1388.7
## ..1 6 1380.7 1397.7 -684.35 1368.7 20.029 2 4.474e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1なお、尤度比検定は、固定効果の検定にも応用できるが、グループレベルの固定効果の場合は ML でグループ数が十分にたくさん（40以上）あることが前提である。
\(-2\log\frac{L_0}{L_1}\)がカイ二乗分布に近似することを利用するのが、通常の尤度比検定であるが、グループ数が少ない（目安は40以下）場合は、このような近似が得られない。そこで正確検定 (Exact Test) を行うことが考えられる。REMLで推定したランダム切片モデルならば、RLRsim パッケージの exactRLRT() 関数で尤度比の正確検定ができる。
exactRLRT(REMLで推定したランダム切片モデルのMerModオブジェクト)
以下は利用例。
library(RLRsim)
exactRLRT(r.itcpt1)##
## simulated finite sample distribution of RLRT.
##
## (p-value based on 10000 simulated values)
##
## data:
## RLRT = 22.299, p-value < 2.2e-16スクリプトを書く手間だけ考えれば、尤度比検定するよりも区間推定するのが簡単である（が、計算には時間がかかる）。R の場合、プロファイル区間推定とブートストラップ区間推定（パラメトリックとセミパラメトリック）が可能で、ブートストラップの場合は上記の自由度／グループ数の問題に拘泥する必要が無いので、ランダム効果の分散だけでなく、固定効果の検定にも応用可能である。
区間推定して信頼区間の中にゼロが含まれないならば（つまり信頼区間の下限値がゼロより大きければ）、ランダム効果の分散は有意にゼロよりも大きいと考えることもできよう。以下で紹介するプロファイル区間推定やブートストラップ区間推定は厳密には有意性検定とは異なるロジックに従っているので、通常の有意性検定のロジックにこだわるならば、用いるべきではない。しかし、あまりこだわる必要はないと私は思う。
プロファイル区間推定とは、尤度比検定を使った区間推定であり、95% 信頼区間の上限値ならば以下のように計算を行う。
同様の作業を行えば、下限値も計算できる。
分散はゼロ以上の値をとるので、分散\(=0\)を帰無仮説とする検定は必然的に片側検定になる。区間推定を検定の代用として用いる場合、\(1 - 限界値\times 2\) を水準として設定すべきであろう。例えば、5% 水準で分散がゼロかどうか検定したいならば、\(1 - 0.05 \times 2 = 0.9\)、つまり 90% 水準の区間推定をし、その下限値がゼロより大きければ帰無仮説を棄却できる。
R でランダム効果の検定のために用いるならば、lme4 パッケージの confint() 関数で以下のように引数を指定すればよい。
confint(MerModオブジェクト, parm = “theta_”, level = 0.9, oldNames =FALSE)
confint(r.slope1, parm = "theta_", level= .9, oldNames = FALSE)## 5 % 95 %
## sd_(Intercept)|id.firm 34.95670 135.97380
## cor_educ.(Intercept)|id.firm -1.00000 1.00000
## sd_educ|id.firm 5.59940 13.62416
## sigma 45.79964 57.72044パラメトリック・ブートストラップとは、推定したモデルが正しいという前提で、シミュレーションによって新しいサンプル（慣例に従ってこれをリサンプル (resample) と呼ぶ）を作り、そのリサンプルで再び同じモデルを推定し、パラメータの推定値をえる。このようなリサンプリングを数百～数千回行い、得られたパラメータの分布を、サンプリングの際の分布であるとみなして、区間推定する方法である。
パラメトリック・ブートストラップで区間推定する場合、
confint(MerModオブジェクト, level =.9, oldNames = FALSE, method = “boot”)
とすればよく、セミパラメトリック・ブートストラップで区間推定する場合、
confint(MerModオブジェクト, level =.9, oldNames = FALSE, method = “boot”, type = “semiparametric”, use.u = TRUE)
とすればよい。以下は計算例である。
confint(r.slope1, oldNames = FALSE, method = "boot", level =.9)## 5 % 95 %
## sd_(Intercept)|id.firm 37.223277 135.6108171
## cor_educ.(Intercept)|id.firm -1.000000 -0.9544231
## sd_educ|id.firm 5.348729 13.3718929
## sigma 43.979245 56.4292335
## (Intercept) -1.859716 86.4595510
## educ 22.902488 30.9256346confint(r.slope1, oldNames = FALSE, method = "boot", level =.9,
type = "semiparametric", use.u = TRUE)## 5 % 95 %
## sd_(Intercept)|id.firm 49.065958 127.2781579
## cor_educ.(Intercept)|id.firm -1.000000 -0.9626437
## sd_educ|id.firm 6.563282 12.0753708
## sigma 39.902004 51.8120321
## (Intercept) 2.069745 68.7699599
## educ 25.047087 30.1571950セミパラメトリック区間推定はランダム効果が正規分布しているかどうか疑わしい場合に特に利用が推奨される。
モデルの適合度を比較したい場合、決定係数を計算するというアプローチと、尤度を使う（尤度比検定や AIC, BIC）というアプローチがあるが、どちらも十分に満足の行く方法とはいえない。もしも最尤法 (ML) で推定していれば、尤度を使うことにあまり大きな問題はない。しかし、制限最尤法（REML）で推定している場合、尤度比検定はランダム効果のみが異なるモデル間の比較にしか使えないし、AIC を使うことに否定的な議論も根強い。これは REML が最尤推定ではないからである。一般にモデルの適合度は、モデルの推定において最大化（あるいは最小化）したものを基準に計算すべきである、という考え方があるようで、ロジスティック回帰分析で通常の残差をもとにした決定係数が計算されず、尤度をもとにした擬似決定係数が用いられるのも、この考え方に従ったものと思われる（ロジスティック回帰分析には最尤推定が用いられ、 OLS は用いられないから）。
いっぽう決定係数を計算するという方法にも当然批判がある。ML でも REML でも残差平方和を最小化しているわけではないからである。REML は一般化最小二乗法と最尤法の二段階推定なので、単純に1つの統計量を最適化しているわけではない。そのため、何を使っても批判はありうるというわけである。それゆえ、グループ数が十分に多く ML でも REML でも推定結果がほぼ同じで、ぜひモデルの適合度を比較したいという時には ML で推定すればよい。
REML を使ったほうが推定が正確になると考えられるような状況では、どの指標を使ってもそれなりに批判がありうるということをふまえておけば、どの適合度指標を使っても構わないと私は思う。社会学の論文でランダム効果モデルが用いられる場合、適合度が表示されないことも珍しくない。
AIC は R の場合 \[ -2 \times \mathrm{対数尤度} + 2 \times \mathrm{推定したパラメータ数} \] で定義される（多少方言もある）。\(-2 \times 対数尤度\) は逸脱度 (deviance) とも呼ばれる。「推定したパラメータ数」とは、切片を含む固定効果すべての数とランダム効果の分散と共分散すべてである。例えば上の Model 2.4 は、固定効果を３つ、ランダム効果の分散を２つ推定しているので、「推定したパラメータ数」は５つである。
BICは \[ -2 \times \mathrm{対数尤度} + \log N \times \mathrm{推定したパラメータ数} \] で定義されるが、\(N\) をどうとるのか、という問題を惹起してしまうので、AIC よりもさらに面倒かもしれない。
R の場合、
AIC(MerModオブジェクト)
BIC(MerModオブジェクト)
で、AIC/BIC のみを返す。以下は使用例。
data1$educ.gm <- ave(data1$educ, data1$id.firm) # グループ平均
rand.itcpt4 <- lmer(wage ~ educ.c + educ.gm + size.c + (1 | id.firm), data = data1)
AIC(rand.itcpt4)## [1] 1374.986BIC(rand.itcpt4)## [1] 1391.955AIC も BIC も値が小さいほどモデルの適合度がよい。ただし、被説明変数が違っていたり、異なるサンプルに対する推定結果を比較するのに尤度や逸脱度や AIC や BIC を使っても無意味である。
尤度比検定とは、階層的な関係にある二つのモデルの比較に用いる。階層的な関係とは、いっぽうのモデル（モデル 1）は他方のモデル（モデル 2）のパラメータのうちの幾つかをゼロに固定した（独立変数をモデルから除外すればそのパラメータをゼロに固定したのと同じことになる）ものである場合をいう。モデル 2 のパラメータの幾つかをゼロに固定したモデルをモデル 1 と呼び、モデル1、モデル2の逸脱度をそれぞれ \(D_1\), \(D_2\), それぞれの推定したパラメータ数を \(k_1\), \(k_2\) とすると、帰無仮説（固定したパラメータはすべて本当にゼロである）のもとで、\(D_1 - D_2\) は、自由度が \(k_2 - k_1\) のカイ二乗分布をする。これが有意になれば、帰無仮説は棄却され、モデル1 で固定したパラメータのうち1つ以上はゼロではない（つまり、モデル 1 よりは、モデル 2 のほうがデータに適合したモデルである）、と考えられる。
R で固定効果が異なる複数のモデルに関して尤度比検定をする場合、
anova(merMod オブジェクト, merMod オブジェクト, …)
とすればよい。以下は使用例。
rand.itcpt5 <- lmer(wage ~ educ.c + (1 | id.firm), data = data1, REML = FALSE)
rand.itcpt6 <- lmer(wage ~ educ.c + size.c + (1 | id.firm), data = data1, REML = FALSE)
anova(rand.itcpt5, rand.itcpt6)## Data: data1
## Models:
## object: wage ~ educ.c + (1 | id.firm)
## ..1: wage ~ educ.c + size.c + (1 | id.firm)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## object 4 1406.4 1417.7 -699.19 1398.4
## ..1 5 1388.1 1402.2 -689.04 1378.1 20.291 1 6.652e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1二つのモデルの逸脱度の差が \(1398.4 - 1378.1 = 20.291\) 、自由度の差が \(5 - 4 =1\)、p値が \(6.652\times 10^{-6} = 0.000006652\) なので、0.1% 水準で後のモデルのほうが当てはまりがよい。
ランダム効果の部分だけが異なるモデル間の比較なら REML でも尤度比検定できる。その場合、
anova(merMod オブジェクト, merMod オブジェクト, …, refit = FALSE)
とする。anova() はデフォルトでは REML で推定していても自動的に ML で推定し直してから、尤度比検定するので、REML のままで尤度比検定したい場合は、" refit = FALSE" と指定する。以下は使用例
rand.itcpt7 <- lmer(wage ~ educ.c + size.c + (1 | id.firm), data = data1)
rand.slope1 <- lmer(wage ~ educ.c + size.c + (educ.c | id.firm), data = data1)
anova(rand.itcpt7, rand.slope1, refit=FALSE)## Data: data1
## Models:
## object: wage ~ educ.c + size.c + (1 | id.firm)
## ..1: wage ~ educ.c + size.c + (educ.c | id.firm)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## object 5 1378.8 1392.9 -684.38 1368.8
## ..1 7 1370.0 1389.8 -678.02 1356.0 12.71 2 0.001738 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1ふたつ目のモデルは、一つ目のモデルと固定効果部分は同じだが、ランダム傾きを仮定している点で異なる。ランダム傾きを仮定したことで、\(U_{1j}\) の分散、\(U_{0j}\) と \(U_{1j}\) の共分散をさらに推定することになり、パラメータの数が二つ増えている。それゆえ、逸脱度の差の自由度は 2 で逸脱度の差が 20.291 だから 0.1% 水準で有意な結果であり、ふたつ目のモデル（つまりランダム傾きを仮定したモデル）のほうが採択される。
ランダム切片モデルの場合、決定係数を計算したいモデル（モデル 1 ）の個人レベルと集団レベルのランダム効果の分散 を\(\sigma^2_1\), \(\tau_{001}\) とする。切片のみのランダム切片モデル（モデル 0 ）の個人、集団レベルのランダム効果の分散をそれぞれ \(\sigma^2_0\), \(\tau_{000}\) とすると、このモデル 1 の決定係数 (\(R^2_1\)）は \[ R^2_1 = 1 - \frac{\sigma^2_1 + \tau_{001}}{\sigma^2_0 + \tau_{000}} \] で計算できる。これは固定効果で従属変数の分散がどれだけ説明できているかを示す。ただし、ランダム傾きモデルにこれを応用するのは難しい。以下は R での計算例。
VarCorr(rand.itcpt7)## Groups Name Std.Dev.
## id.firm (Intercept) 26.222
## Residual 56.402rvar1 <-26.222^2 + 56.402^2
rand.itcpt.null <- lmer(wage ~ 1 + (1 | id.firm), data = data1)
VarCorr(rand.itcpt.null)## Groups Name Std.Dev.
## id.firm (Intercept) 82.629
## Residual 89.512rvar0 <- 82.629^2 + 89.512^2
1 - rvar1 / rvar0 # 決定係数## [1] 0.7392997ランダム傾きモデルの場合、ベースライン・モデルと比べて、それぞれの分散がどれだけ減少したか、それぞれ計算するというやり方も考えられる。例えば、ミクロレベルの独立変数とランダム効果だけをモデルに投入し、グループレベルの変数は無いモデルを推定し、その後、グループレベルの変数を投入してどれだけグループレベルのランダム効果の分散が減少したか、見ることもできよう。以下は計算例。
rand.slope0 <- lmer(wage ~ educ.c + (educ.c | id.firm), data = data1)
rand.slope2 <- lmer(wage ~ educ.c*size.c + (educ.c | id.firm), data = data1)
screenreg(list(rand.slope0, rand.slope2), custom.model.names = c("Model 2.6", "Model 2.7"))##
## =========================================================
## Model 2.6 Model 2.7
## ---------------------------------------------------------
## (Intercept) 394.61 *** 391.76 ***
## (9.41) (6.11)
## educ.c 26.56 *** 25.05 ***
## (2.47) (1.98)
## size.c 2.80 ***
## (0.44)
## educ.c:size.c 0.56 ***
## (0.12)
## ---------------------------------------------------------
## AIC 1380.71 1359.43
## BIC 1397.68 1382.06
## Log Likelihood -684.35 -671.72
## Num. obs. 125 125
## Num. groups: id.firm 25 25
## Var: id.firm (Intercept) 1591.37 262.00
## Var: id.firm educ.c 85.36 24.18
## Cov: id.firm (Intercept) educ.c 368.57 79.59
## Var: Residual 2646.89 2598.79
## =========================================================
## *** p < 0.001, ** p < 0.01, * p < 0.051 - 262.00 / 1591.37 # 各企業間の賃金の差がどの程度説明されているか## [1] 0.8353621 - 24.18 / 85.36 # 教育年数の係数バラつきがどの程度説明されているか## [1] 0.71672911 - 2598.79/2646.89 # 企業内の賃金の差が、企業規模による教育の効果の違いを考慮することでどの程度説明されたか## [1] 0.01817227このような決定係数の計算は、ベンチマークのモデルを何に設定するかで変わってくるので、注意が必要であろう。