最近はAIを使った予測が流行っていますが、予測よりも現象を解釈したい場合もあるかと思います。
行った施策に対しての効果を分析したり、事象が起こった後に要因を解釈したり...
そんな時に使うのがマルチレベル分析です。現実の仮定を反映しやすく、様々な仮説について解析をすることができます。
今回はそんなマルチレベル分析を導入する必要性や、解析、結果の解釈方法などをまとめていきます!
Rでの実行方法は補足として記事中に記載してあるので、参考にしてみてください!
マルチレベル分析とは?
マルチレベル分析とは、複数の階層性をもつデータを解析する時に用いられる分析手法です。例えばテストの点数と勉強時間との関連を調べる時に対象となる学校を選び、そこから数人の生徒を選んだ場合や、残業時間と従業員の満足度を調べる際に、企業を選び、従業員を選んだ場合、季節ごとに商品の売り上げを調査した場合など、二段階以上のレベルからサンプルを抽出したデータを階層性のあるデータと呼びます。
なぜマルチレベル分析が必要か?
ここでは、テストの点数と勉強時間の関連性を調べることにしましょう。
以下の二つの疑問を考えてみます。
- 学校間で平均点数に違いがあるか
- 個人レベルで勉強時間と点数に関連はあるのか
とりあえず、二つ目の疑問から考えてみます。
全部で7校の学校から生徒を数人ずつ選んだと仮定して、テストの点数と勉強時間から回帰直線を求め、散布図に重ねてみます。
X軸が勉強時間、Y軸が点数を表しています。ここからは勉強時間が1時間長いと点数が19点高くなることが分かります。
(上記のデータセットは以下のように作成しました)
library(dplyr) #データセットの作成 set.seed(123) mu <- list(20,20,30,30,40,50,50) N_st <- list(10,5,40,25,5,15,10) ID <- map(N_st,function(n)seq(1:n)) %>% flatten_dbl() X <- map2(N_st,mu,rnorm,sd=10) %>% flatten_dbl() sctype <- sample(c(0,1),replace = T,7) school <- rep(1:7, times=unlist(N_st)) a0 <- 1000 b0 <- 5 a <- a0 + rnorm(length(N_st), mean=0, sd=300) %>% sort(decreasing = T) b <- b0 + 10*sctype+rnorm(length(N_st), mean=0, sd=10) sctype <- sctype[order(b)] b <- b[order(b)] data_frame(X=X,school = as.factor(school), a=a[school],b=b[school],sctype=sctype[school],ID=ID) %>% mutate(Y=rnorm(length(X),a+b*X,100)) %>% select(X,ID,school,sctype,Y)-> st_df
Xは勉強時間、IDは学校内での生徒番号、schoolは学校識別ID、sctypeは公立か私立か(後から使います)、Yはテストの点数を表しています。
X | ID | school | sctype | Y |
---|---|---|---|---|
14.39524 | 1 | 1 | 0 | 1453.686 |
17.69823 | 2 | 1 | 0 | 1341.678 |
35.58708 | 3 | 1 | 0 | 1430.010 |
20.70508 | 4 | 1 | 0 | 1476.366 |
21.29288 | 5 | 1 | 0 | 1520.955 |
37.15065 | 6 | 1 | 0 | 1333.399 |
今度は、学校ごとに色付けをしてしてみましょう。
どうでしょうか?どうやら学校ごとに勉強時間や点数に傾向があり、単純に個人レベルで勉強時間が長いほど点数が伸びるとは言えなさそうです。
さらに、仮に以下のような関係であれば、個人の勉強時間と点数に関連がなくても、学校全体で見ると関連があるように見えてしまいます。
これがマルチレベル分析を導入するべき理由です。
マルチレベル分析を利用すると、個人レベルと学校レベルの関連を分離させてデータの解釈をすることができます。
マルチレベル分析の適応
さて、マルチレベル分析を適応すべきかどうかはどのように判断すれば良いでしょうか?もちろんデータが収集された背景を考えて決めても良いですが、何か指標が欲しいところですよね。
ところで、目的変数である点数Yのバラツキには、学校間でのバラツキと個人間でのバラツキが含まれています。
マルチレベル分析は、バラツキを学校レベルと個人レベルに分離させて解析するというのが大まかな目的となります。
そこで級内相関というものを使います。級内相関とは同じ群内にある値がどれぐらい似ているかという指標になるのですが、ここでは学校間の点数のバラツキが学校内のバラツキに比べてどれぐらい大きいかを表します。
つまり、個人よりも学校による点数の違いが大きい場合は、マルチレベル分析を適応するべきだということになります。
この級内相関を求める際には、今回は傾きを0とした線形混合モデルの中のランダム化切片モデルというものを使います。難しそうな名前ですが、各学校の点数の平均値(切片)が全体平均を中心とした正規分布に従い、学校内の生徒の点数は学校の平均値を中心とした正規分布に従うことを前提とするモデルのことを指します。
傾きを0としているので、一元配置分散分析と似たようなイメージです。
下の図はX軸を各学校、Y軸を点数としたグラフです。
赤色の正規分布の分散が、学校間の分散、黒色が学校内の分散を表しています。
(ちなみに正規分布を仮定しないようなデータは線形混合モデルでは対応できないので、ベイズ推定による階層ベイズモデルを使います。)
難しい方はざっくりとイメージだけでもつかめればOKです。
式で表すと
となります。
ここから級内相関は、
として求めることができます。
これは、群内のバラツキと比べて、群間のバラツキがどれぐらい大きいか、を表す指標になり、まさにマルチレベル分析を適応すべきかどうかを判断できる指標ですね。
これらのパラメータを求めてみると、
となってICCは98742/(98743+28116)=約78%となります。
---補足---
学校内分散や学校間分散の求め方についてご質問頂きましたので、コメントにて回答させて頂きました。
----------
明確な基準はないですが、回帰分析は全ての観測値が独立であることを前提としているので、級内相関が0.1(10%)でもマルチレベル分析を適応すべきではないかとされています。
この級内相関が意味するところは、点数変動の78%は学校間変動によって説明されるという意味になります。
ところで、学校間の違いによる影響が大きいということは、学校間で平均点に大きな違いがあるということです。
つまり、学校間で平均点数に違いがあり、それは切片平均(1617点)を中心に314点(98742の平方根)ぐらいのバラツキがあるということで、第一の疑問に対して答えることができました!
Rで実行するには{lme4}パッケージを使います。{lme4}の使い方についてはこちらも参考にしてみてください。
library(lme4) st_df %>% lme4::lmer(data=.,Y~(1|school),REML = F) %>% summary() >Linear mixed model fit by maximum likelihood ['lmerMod'] Formula: Y ~ (1 | school) Data: . AIC BIC logLik deviance df.resid 1471.5 1479.6 -732.7 1465.5 107 Scaled residuals: Min 1Q Median 3Q Max -2.77716 -0.57499 -0.02915 0.56146 3.01817 Random effects: Groups Name Variance Std.Dev. school (Intercept) 98742 314.2 Residual 28116 167.7 Number of obs: 110, groups: school, 7 Fixed effects: Estimate Std. Error t value (Intercept) 1617.6 120.5 13.42
REML=Fはパラメーターの推定方法に制限つき最尤法(REsidual Maximum. Likelihood estimation)を使うか、通常の最尤法を使うかを表します。今回は通常の最尤法を使うためFALSEとしました。
結果の見方ですが、Random effectsのschool(intercept)が学校間分散、Residualが全分散から学校間分散を引いた学校内分散、Fixed effectsの(intecept)が切片平均を表しています。
> #ICC > 98742/(98743+28116) [1] 0.7783602
試しにランダムなデータで分散分析モデルを試してみると
set.seed(123) d <- data_frame(x=c(rnorm(70,50,10)), school = as.factor(rep(c(1:7),rep(10,7)))) d %>% lmer(data=.,x~(1|school),REML = F) %>% summary() singular fit Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's method [lmerModLmerTest] Formula: x ~ (1 | school) Data: . AIC BIC logLik deviance df.resid 512.4 519.2 -253.2 506.4 67 Scaled residuals: Min 1Q Median 3Q Max -2.2648 -0.6313 -0.0134 0.6535 2.3255 Random effects: Groups Name Variance Std.Dev. school (Intercept) 0.00 0.000 Residual 81.17 9.009 Number of obs: 70, groups: school, 7 Fixed effects: Estimate Std. Error df t value Pr(>|t|) (Intercept) 50.738 1.077 70.000 47.12 <2e-16 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 convergence code: 0 singular fit
singular fitと出ます。切片の分散(学校間分散)がとても小さいため推定がうまくできていません。つまり、学校間の違いはほとんどないことを表しています。
ランダム化切片モデル
分析を進めていきます。先ほど学校間でテストの平均点に差があることが分かりましたので、仮説を少し変えて以下の疑問に答えていきましょう。
- 勉強時間の影響を除いた上で、学校間で平均点数に違いがあるか
- 個人レベルでも勉強時間と点数に関連はあるのか
まずは、勉強時間の影響を除いた上で学校間で平均点に差があるかどうかを明らかにしていきます。
どいうことかというと、テストの点数だけではなく、勉強時間にも学校間で差があることが考えられるため、仮に勉強時間が同じ生徒を比べた時にも、学校間で差があるかどうかを検証します。
イメージは以下のようになります。
この切片分散が0でなければ勉強時間の影響を除いた後でも、学校ごとに点数に違いがあることになります。ここでは傾きが学校間で一緒であるという前提のもとに解析するため、共分散分析のイメージに近いですね。
式で書いてみると
これを求めてみると、
となりました。
ここから、学校間分散が0ではないため勉強時間の影響を除いても学校間には183点(33510の平方根)ぐらいの差があると言えそうです。
さらに、級内相関は33510/(33510+16318)=0.67(67%)となり、少し小さくなりました。学校内分散の減少率よりも学校間分散の減少率が大きい、つまり勉強時間は学校レベルの違いに大きく影響したことが分かります。
st_df %>% lme4::lmer(data=.,Y~X+(1|school),REML = F) %>% summary() Linear mixed model fit by maximum likelihood [lmerMod] Formula: Y ~ X + (1 | school) Data: . AIC BIC logLik deviance df.resid 1410.0 1420.8 -701.0 1402.0 106 Scaled residuals: Min 1Q Median 3Q Max -2.16891 -0.65933 -0.02279 0.63100 2.67048 Random effects: Groups Name Variance Std.Dev. school (Intercept) 33510 183.1 Residual 16318 127.7 Number of obs: 110, groups: school, 7 Fixed effects: Estimate Std. Error t value (Intercept) 1171.88 85.43 13.717 X 12.74 1.36 9.365 Correlation of Fixed Effects: (Intr) X -0.558 #説明変数を入れないモデル Random effects: Groups Name Variance Std.Dev. school (Intercept) 98742 314.2 Residual 28116 167.7 Number of obs: 110, groups: school, 7
ここでは、学校間の違いに興味があるためRandom effectsのみ記載しました。
説明変数の中心化
さて、次の疑問である個人レベルでも勉強時間と点数に関連はあるのかに答えるために、中心化の概念を説明しておきます。中心化とは説明変数のバラツキを個人レベルと集団レベルに分離して解析するための手法です。
先ほどのモデルの傾きb0を確認してみると16となっていて、これは勉強時間が1時間長いと点数が16点高くなることを示しています。
しかし、勉強時間には学校レベルと個人レベルの影響が混在しているため、個人の関係を表しているのか、そもそも勉強時間が平均的に長い学校の点数が高いのか分かりません。
そのため、各個人の勉強時間から学校ごとの平均勉強時間を引くことで、学校による勉強時間の違いを除くことができます。 これを集団平均中心化(Centering Within Cluster:CWC)と呼びます。 以下はCWCを行う前とCWCを行った後の勉強と点数の関係を表したグラフです。
下のグラフでは、勉強時間の平均値が0になっており、学校による勉強時間の違いが排除されていることがわかります。
このCWCを行った後の勉強時間を説明変数とすることで、純粋な個人レベルの勉強時間と点数の関係を推定することができます。
また、全体平均中心化(centerring at the grand mean:CGM)という手法もあります。これは、各個人の勉強時間から全体の勉強時間平均を引く手法です。
さて、説明変数にCWC(集団平均中心化)を行った勉強時間をモデルに加えても傾きが有意であれば、個人レベルで勉強時間と点数に関係があると言えそうです。
式で表すと
となります。
これを求めてみると、
となりました。ここから傾き12、標準誤差1.39なので、個人レベルでも勉強時間と点数に関係があることが分かりました!
説明変数を中心化したため、切片平均はそれぞれの学校で勉強時間が学校平均の人の点数の平均を表し、学校間分散はそのバラツキ具合となります。(ややこしい....)
#集団平均中心化のXを作る st_df %>% group_by(school) %>% mutate(X_sc=mean(X),X_cwc=X-X_sc) -> st_df2 st_df2 %>% lme4::lmer(data=.,Y~X_cwc+(1|school),REML = F) Linear mixed model fit by maximum likelihood ['lmerMod'] Formula: Y ~ X_cwc + (1 | school) Data: . AIC BIC logLik deviance df.resid 1417.2 1428.0 -704.6 1409.2 106 Scaled residuals: Min 1Q Median 3Q Max -2.26259 -0.65575 -0.01902 0.65463 2.67923 Random effects: Groups Name Variance Std.Dev. school (Intercept) 99957 316.2 Residual 16278 127.6 Number of obs: 110, groups: school, 7 Fixed effects: Estimate Std. Error t value (Intercept) 1617.087 120.504 13.419 X_cwc 12.063 1.394 8.655 Correlation of Fixed Effects: (Intr) X_cwc 0.000
Fixed effectsの(intercept)とX_CWCが切片と傾きになります。
---2月14日追記---
学校レベルの影響を追加
さて、中心化を行うと学校レベルの影響が排除されることが分かりました。するとこのモデルでは、学校間での比較をすることができなくなります。
そこで、マルチレベル分析において集団平均中心化を用いた際は、その集団平均を集団レベルの変数に加えることが推奨されています。
つまり、平均的に勉強時間が長い学校は点数が高いのかどうか、という学校レベルでの議論できるようになります。
また、集団レベルの変数を投入する際には全体平均中心化(CGM)をすることで、今後の解析が行いやすくなります。CGMを行う詳しい理由については、下記の書籍を参考にしてみてください!
Rで学ぶ マルチレベルモデル[入門編]: 基本モデルの考え方と分析
- 作者: 尾崎幸謙,川端一光,山田剛史
- 出版社/メーカー: 朝倉書店
- 発売日: 2018/09/12
- メディア: 単行本
- この商品を含むブログを見る
以上を式で表すと
ここでは個人レベルの切片に関わるパラメーターにはaをつけておきます。 少しややこしですが、これを求めてみると
となります。個人レベルでの傾きは12点と変化していません。しかし、学校レベルでの平均勉強時間の傾きが推定されていて25となっています。つまり、勉強時間が全体平均よりも長い学校にいる生徒は点数が高くなるということが分かります。
さらに学校種(公立、私立)も変数に加えてみましょう。ここでは、学校種によって点数に差があるかどうかを検討できます。
これを求めると
となりました。学校種の勉強時間の影響は、公立(0)に比べて私立(1)になると83点高くなります。
と言いたいところですが標準誤差が非常に大きく傾きが有意ではないため、学校種によって点数に違いがあるとは言えなさそうです。
これは学校種をモデルに加えても、学校間分散や学校内分散に大きな変化がみられなかったことからも推測できます。
このように、同じモデルの中で学校レベルと個人レベルの二つのレベルでデータを解釈することができます!
#勉強時間の学校平均に全体平均中心化を行う st_df2 %>% ungroup(school) %>% mutate(X_sc_cgm = X_sc - mean(X))->st_df3 #実行する st_df3 %>% lme4::lmer(data=.,Y~X_cwc+X_sc_cgm+(1|school),REML = F) %>% summary() Linear mixed model fit by maximum likelihood [ lmerMod] Formula: Y ~ X_cwc + X_sc_cgm + (1 | school) Data: . AIC BIC logLik deviance df.resid 1404.5 1418.0 -697.3 1394.5 105 Scaled residuals: Min 1Q Median 3Q Max -2.33690 -0.65150 -0.03606 0.65453 2.82370 Random effects: Groups Name Variance Std.Dev. school (Intercept) 10723 103.5 Residual 16297 127.7 Number of obs: 110, groups: school, 7 Fixed effects: Estimate Std. Error t value (Intercept) 1595.514 42.135 37.867 X_cwc 12.063 1.395 8.650 X_sc_cgm 25.636 3.611 7.099 Correlation of Fixed Effects: (Intr) X_cwc X_cwc 0.000 X_sc_cgm -0.085 0.000 #学校種追加 st_df3 %>% lme4::lmer(data=.,Y~X_cwc+X_sc_cgm+sctype+(1|school),REML = F) %>% summary() Linear mixed model fit by maximum likelihood [ lmerMod] Formula: Y ~ X_cwc + X_sc_cgm + sctype + (1 | school) Data: . AIC BIC logLik deviance df.resid 1406.4 1422.6 -697.2 1394.4 104 Scaled residuals: Min 1Q Median 3Q Max -2.3264 -0.6497 -0.0284 0.6601 2.8163 Random effects: Groups Name Variance Std.Dev. school (Intercept) 10247 101.2 Residual 16315 127.7 Number of obs: 110, groups: school, 7 Fixed effects: Estimate Std. Error t value (Intercept) 1564.128 86.921 17.995 X_cwc 12.063 1.395 8.645 X_sc_cgm 22.429 8.560 2.620 sctype 83.016 201.773 0.411 Correlation of Fixed Effects: (Intr) X_cwc X_sc_c X_cwc 0.000 X_sc_cgm 0.784 0.000 sctype -0.880 0.000 -0.910
まとめ
今回はマルチレベル分析の必要性とその基本であるランダム化切片モデルをまとめました。
マルチレベル分析は、説明変数の中心化や結果の解釈などが少し複雑ですが、適切な解析をするには必須の分析法になると思います。
今回は扱わなかった、傾きにもランダム性を仮定したランダム切片傾きモデルや、複数のモデルから適切なモデルを選択する方法などもあるので、次回まとめていきます!
※本記事は筆者が個人的に学んだことをまとめた記事なります。数学の記法や詳細な理論、筆者の勘違い等で誤りがあった際はご指摘頂けると幸いです。
参考
こちらの内容を参考にしました。マルチレベル分析の基礎から具体的な使い方がまとめてあり、とても分かりやすい内容になってます。実践的に使いたい方にはおすすめです!
Rで学ぶ マルチレベルモデル[入門編]: 基本モデルの考え方と分析
- 作者: 尾崎幸謙,川端一光,山田剛史
- 出版社/メーカー: 朝倉書店
- 発売日: 2018/09/12
- メディア: 単行本
- この商品を含むブログを見る
こちらも同じくマルチレベル分析の方法について書かれていますが、今回扱った線形混合モデルに加え構造方程式モデルなども解説されていて、心理学、社会学をテーマにした内容になっています。
- 作者: 清水裕士
- 出版社/メーカー: ナカニシヤ出版
- 発売日: 2014/10/01
- メディア: 単行本
- この商品を含むブログを見る