Rを通じて統計学を学ぶ備忘録ブログ

SPSSからRに移行したい私のような人向けのR解説ブログ兼学習用備忘録。

第59回 マルチレベル分析〜R実践〜

第59回はRでマルチレベル分析を行っていきます。用いるデータはマルチレベル分析の生みの親でもあるRaudenbushが実際に使った高校のデータを使っていきます。データはここからダウンロードできます。 

□hsb12の内容

hsb12の内容は、高校をサンプリングした後に生徒をサンプリングする2段抽出法であり、まさにマルチモデル分析の格好の例です。データはこんな感じです。

f:id:teruaki-sugiura:20151003014420p:plain

級内相関係数を計算してみましょう。級内相関係数は、集団間分散が大きくなると係数は高くなり、集団内分散が大きくなると、係数は小さくなります。−1から1までの範囲を取ります。級内相関が高い場合、1つの集団から発生したサンプルではなく、異なる集団から発生したサンプルと考える必要があります。Rで計算してみると(パッケージpsychのICC関数)だいたい、0.27程度です。つまり、sesは学校間で異なるということがわかります。

 □マルチレベル分析

前回の記事でも説明したように、今回の最終目標は以下の式に辿りつくことです。

f:id:teruaki-sugiura:20151003014525p:plain

この式が意味するところは、切片と回帰係数に関して、固定・変量効果を同時に推定することです。この時の回帰係数のことをランダム係数と呼ぶこともあります。以下の図のようになります。あくまでイメージです。

f:id:teruaki-sugiura:20151003014732p:plain 

傾き、切片の固定効果(平均)があり、その平均から集団間の影響、つまり確率的なばらつき(変量効果)があることを意味しています。

 □マルチレベル分析と中心化

マルチレベル分析を実際に使っていくとなると、説明変数を中心化する必要があります。なぜなら、級内相関がある変数に関しては集団レベルと個人レベルの効果がごちゃまぜになっているからです。つまり、上記の分析を行う際に「sesが高い生徒のテストの成績がよいのか」「sesが高い学校のテストの成績がよいのか」のどちらが正しい結果(効果)なのかを明確にできません。それゆえ「集団平均中心化(centering within cluster: cwc)」と「全体平均中心化(Grand mean centering:gmc)」を行うことで効果を分離する必要があります。「集団平均中心化」では、集団の平均値から得点を引くことで、集団レベルの効果を除去し、純粋な個人レベルでの効果に変換します。次に、「全体平均中心化」ですが、これは全体の平均値から得点を引くことになりますが、これが意味するところはなんでしょうか。さきほどの集団平均中心化をすることで、「集団間の効果」がなくなってしまっているので、これを表すために全体平均中心化は行われます。つまり、Level2の変数に関しては全体平均中心化を行い、Level1の変数については集団平均中心化をする必要があります(おそらく一般的な方法)。また、このように変数を中心化することでマルチコが起こりにくくなります。 

□Rでの実践

パッケージLmerTestのインストールは忘れずに。lme4パッケージでもマルチレベル分析はできますが、有意性検定が行われないのでLmerTestのほうがいいかも。

hsb12 <- read.csv("hsb12.csv",header=TRUE)
hsb12$ ses.group.mean <- tapply(hsb12$ses,hsb12$school,mean,na.rm=TRUE[as.character(hsb12$school)]
hsb12$ses.cwc <- hsb12$ses-hsb12$ses.group.mean
hsb12$ses.gmc<- hsb12$ses.group.mean-mean(hsb12$ses.group.mean, na.rm=TRUE)
result.null <- lmer(mathach~1+(1|school),data=hsb12,REML=FALSE)
result.clossinter <- lmer(mathach~ses.cwc+ses.gmc+sector+ses.gmc:ses.cwc+ses.cwc:sector+(ses.cwc|school),data=hsb12,REML=FALSE)

f:id:teruaki-sugiura:20151003020802p:plain

f:id:teruaki-sugiura:20151003020809p:plain

f:id:teruaki-sugiura:20151003021203p:plain

数学の成績は集団間(学校間)で異なり、数学の成績に対するsesの影響にも集団間(学校間)で異なることが示されたことになります。

 

 □単純主効果の検定

重回帰分析の交互作用同様に、PreacherのTwo-Way Interaction Effects in HLMのオンラインツールを使います。共分散行列はvcov関数で得られます。共分散行列と回帰分析のサマリーを使って、表を埋めていきます。

vcov.merMod(result.clossinter)

f:id:teruaki-sugiura:20151003015846p:plainf:id:teruaki-sugiura:20151003025341p:plain

xx <- c(-0.6606,0.6606) # <-- change to alter plot dims
yy <- c(9.3762,14.2156) # <-- change to alter plot dims
leg <- c(-0.6606,9.9811) # <-- change to alter legend location
x <- c(-0.6606,0.6606) # <-- x-coords for lines
y1 <- c(10.1827,14.0747)
y3 <- c(12.4956,14.2156)
plot(xx,yy,type='n',font=2,font.lab=2,xlab='SES',ylab='Mathach',main='H LM 2-Way Interaction Plot')
lines(x,y1,lwd=3,lty=1,col=1)
lines(x,y3,lwd=3,lty=6,col=2)
points(x,y1,col=1,pch=16)
points(x,y3,col=2,pch=16)
legend("bottomleft",legend=c('sector=0','sector=1'),lwd=c(3,3),lty=c(1,6),col=c(1,2))

f:id:teruaki-sugiura:20151003024832p:plain

SESが高くなるとsectorにあまり関係なく、数学の成績が高くなっていることがわかります。 以上でおしまい。