モデルが推定のための様々な仮定をどの程度満たしていたのか、チェックする作業をモデルの診断 (diagnostics) という。もしも満たしていなければ、推定結果はゆがんでいる可能性がある。多少なら仮定を満たしていなくても問題ない場合もあるが、モデルの仮定と現実の乖離が大きくなるほど、推定値は信用できなくなっていく。
独立変数間の関連が強すぎることを多重共線性 (multicolinearity) という。多重共線性が強すぎると、標準誤差が非常に大きくなり、推定値が不安定になる(ちょっとしたモデルの変更で推定値が大きく変化する)。多重共線性のチェックにはさまざまな統計指標が用いられるが、VIF (Variance Inflation Factor) が最もよく知られている。VIF は大きいほど多重共線性が強いが、どの程度大きければ、モデルの信頼性に問題が生じるかは、サンプル・サイズや説明変数の分布の偏りなどにも影響を受けるので、一概には言えない。VIF だけでなく、標準誤差や推定値の不安定性を見て判断するのがよいだろう。よく VIF が 2 以上だとダメとか 5 以上だとダメといった基準を聞くが、サンプルサイズが大きく、説明変数や残差がおかしな分布をしていなければ、VIFが 10 以上あっても安定した通常の推定が可能である。
R でマルチレベル分析をする場合、merMod オブジェクトに対応する VIF の関数がないので、lm() で同じ固定効果のモデルを推定し、それを引数として car パッケージの vif() 関数で計算すればよい。
vif(lmオブジェクト)
VIF は、独立変数(固定効果)の相関行列から計算されるので、固定効果部分のモデルが同じならば、回帰分析でもロジスティック回帰分析でも、マルチレベル分析でも、すべて同じになる。以下は計算例。
head(data1)
## id.firm id.worker firm.size educ wage educ.c size.c
## 1 a 1 10 11 242 -2.104 -24.0204
## 2 a 2 10 7 267 -6.104 -24.0204
## 3 a 3 10 11 449 -2.104 -24.0204
## 4 a 4 10 7 238 -6.104 -24.0204
## 5 a 5 10 13 345 -0.104 -24.0204
## 6 b 6 209 7 259 -6.104 -22.0304
library(lme4)
rand.slope1 <- lmer(wage ~ educ.c * size.c + (educ.c | id.firm), data = data1)
lm1 <- lm(wage ~ educ.c*size.c, data = data1) # rand.slope2 と同じ固定効果を OLS で推定
library(car)
vif(lm1)
## educ.c size.c educ.c:size.c
## 1.319346 1.312923 1.006803
マルチレベル分析では、ランダム効果は
という仮定が置かれている。これらが満たされていない場合、推定結果が歪んでいる可能性があるので、ランダム効果がどう分布しているかはチェックしたほうが良い。
R の場合、レベル1のランダム効果 (\(R_{ij}\))の分布は、
plot(merModオブジェクト)
で、簡単に図示できる。以下は計算例。
plot(rand.slope1)
横軸 (fitted(.)) は従属変数のモデルからの予測値で、縦軸がレベル1 の残差(ランダム効果)である。\(U_j\) と \(R_{ij}\) が無相関で、平均ゼロで分散が一定なら、残差も平均ゼロで分散は一定の正規分布になるはずである。上の例の場合モデルからの予測値が 300 あたりで、やや分散が大きくなる傾向が見られるので分散の不均一性を疑ったほうがいいかもしれない。しかし、横軸 (fitted(.)) の値にかかわらず、おおむね残差の平均はゼロなので、独立変数に関しては非線形の関係はあまり考えなくてもよさそうである。
残差が正規分布しているかどうかは Quantile-Quantile プロットでみることができる。R の場合 lattice パッケージの qqmath() 関数を使うのが簡単である。
qqmath(merModオブジェクト)
以下は計算例。
library(lattice)
qqmath(rand.slope1)
これが一直線にならんでいれば完全な正規分布である。数が 125 しかないので、この程度の乖離ならばおおむね正規分布とみなしてよいと思う。
R の場合、グループレベルのランダム効果の推定値は
ranef(merModオブジェクト)
で取り出すことができる。ただし返り値はリスト形式なので注意。以下は計算例。
ranef2 <- ranef(rand.slope1)
ranef2
## $id.firm
## (Intercept) educ.c
## a 4.48821582 1.36350738
## b -3.04302314 -0.92446189
## c 5.24529903 1.59350713
## d 6.01374435 1.82695867
## e -5.81032381 -1.76516007
## f -8.25850307 -2.50891007
## g 4.47507934 1.35951655
## h -2.04504067 -0.62127762
## i 15.86470259 4.81965215
## j 14.46980731 4.39588687
## k -0.04806699 -0.01460262
## l 8.64200340 2.62541639
## m -12.79028445 -3.88565254
## n -11.22093654 -3.40888905
## o -4.44533902 -1.35048152
## p -13.35619723 -4.05757526
## q -3.05044305 -0.92671604
## r -17.58625965 -5.34265635
## s -2.03208628 -0.61734211
## t -12.49735331 -3.79666088
## u -9.26246510 -2.81391092
## v 0.94110592 0.28590534
## w 23.97756541 7.28431712
## x 16.05921135 4.87874336
## y 5.26958779 1.60088599
あとは、グループレベルの変数とこれらのランダム効果の推定値の散布図を作ればよい。
firm.size <- tapply(data1$firm.size, data1$id.firm, mean) # グループ単位のデータを作る
par(mfrow=c(1,2))
plot(ranef2$id.firm[,1] ~ firm.size) # 切片のグループレベルのランダム効果
plot(ranef2$id.firm[,2] ~ firm.size) # 教育の傾きのランダム効果
企業規模が大きい方がランダム効果の分散が大きくなっているように見える。
ランダム効果の Quantile-Quantile プロットも lattice パッケージの qqmath() 関数で作れる。以下は計算例。
qqmath(ranef2)
## $id.firm
これぐらいならおおむね正規分布と言えるだろう。感覚的でわからないという人は同じサイズの正規分布する乱数を何度か作って Quantile-Quantile プロットしてみればよい。この程度直線的なら正規分布する乱数とみなせる。ただし、数が多くなるほど一直線に近くなるので、有効サンプル・サイズと同じ数だけ乱数を発生させるのがポイントである。
set.seed(1986)
qqmath(rnorm(125))