社会科学では、個人と集団といったレベルの異なる対象を同時に扱うことがしばしばある。例えば、労働者の賃金は、個人の年齢や学歴だけでなく、彼女の勤め先の規模や経営状況の影響も受けるかもしれない。このような複数のレベルの独立変数が、個人のレベルの従属変数に及ぼす影響を分析するための分析法がマルチレベル分析 (Multilevel Analysis) である。マルチレベル分析は、
など、いくつかの呼び名で呼ばれているが、数学的には同じものである(使われる文脈とニュアンスは多少違う)。
マルチレベル分析は回帰分析の発展型の一つと位置づけられるが、複数のレベルの残差(マルチレベル分析ではランダム効果と呼ばれる)を同時に扱える点が、通常の最小二乗法 (Ordinal Least Square OLS) での回帰分析とは異なる点である。それによって標準誤差を適切に推定することができる。
なお、以下では二つのレベルのうち、
と呼ぶ。レベル1の単位は個人である必要はなく、時点でも職業でも文献でも、何でもよい。ただし、以下ではグループ数が最低でも10以上はあり、各グループ内にミクロ・レベルのケースが2つ以上含まれている(少数のグループでケースが1 でも可)をマルチレベル・データと呼ぶ。
架空のデータでマルチレベル分析がどのようなものか、例示しよう。下の図は、25 の企業を無作為抽出し、その従業員を 5 人ずつサンプリングした結果である(あくまで架空)。
# data1 という架空のデータの data frame をあらかじめ作ってある
# データは左から順に、企業ID, 個人ID, 企業規模(従業員数), 教育年数、賃金
head(data1, 12)
## id.firm id.worker firm.size educ wage
## 1 a 1 10 11 242
## 2 a 2 10 7 267
## 3 a 3 10 11 449
## 4 a 4 10 7 238
## 5 a 5 10 13 345
## 6 b 6 209 7 259
## 7 b 7 209 9 304
## 8 b 8 209 6 205
## 9 b 9 209 14 294
## 10 b 10 209 9 321
## 11 c 11 412 8 244
## 12 c 12 412 15 374
summary(data1)
## id.firm id.worker firm.size educ wage
## a : 5 Min. : 1 Min. : 10 Min. : 6.0 Min. :188.0
## b : 5 1st Qu.: 32 1st Qu.:1214 1st Qu.:11.0 1st Qu.:310.0
## c : 5 Median : 63 Median :2421 Median :13.0 Median :389.0
## d : 5 Mean : 63 Mean :2412 Mean :13.1 Mean :404.7
## e : 5 3rd Qu.: 94 3rd Qu.:3612 3rd Qu.:16.0 3rd Qu.:480.0
## f : 5 Max. :125 Max. :4810 Max. :20.0 Max. :817.0
## (Other):95
# 企業(規模)別に教育年数と賃金の散布図を作成
library(ggplot2)
g <- ggplot(data1, aes(x = educ, y = wage)) # グラフのデータを指定
g <- g + geom_point() # グラフの種類を指定
g <- g + facet_wrap(~firm.size) # 企業規模別にグラフを作るよう指定
g
企業によって傾向は様々だが、おおむね教育年数が上がるほど賃金も上がる傾向が見られる。また、従業員数の多い企業ほど賃金が高い傾向があることがわかるだろう。
このような傾向を踏まえると、回帰分析の知識があれば、以下のような回帰式が考えられるだろう。 \[ Y_{ij} = \beta_0 + \beta_1 X_{ij} + \beta_2 Z_{j} + R_{ij}. \qquad (1.1) \] ただし、\(Y_{ij}\) は \(j\) という企業で働く \(i\) という個人の賃金、\(X_{ij}\) は \(j\) という企業で働く \(i\) という個人の教育年数、\(Z_j\) は \(j\) という企業の規模、\(R_{ij}\) は残差である。これは、OLS で推定可能な通常の回帰モデルである。
しかし、このモデルがこのデータに適当かどうかは疑ってかかるべきである。なぜなら、このデータは、個人が企業にネスト (nest) されているからである。OLS では残差は独立に分布しているという仮定があるが、同じ企業に属している個人の残差は相関している可能性がある。例えば、賃金の高い産業に属する企業の労働者は全員、上のモデルの予測値よりも高めの賃金を得ているかもしれない。そのような場合、残差は相関している、という。このように残差が相関している場合に無理に OLS でパラメータを推定すると、標準誤差を過小に(場合によっては過大に)推定し、第一種の(あるいは第二種の)過誤を犯しやすくなってしまう。また、OLS では、係数の検定を t 分布を使って行うのが通例であるが、その際の t 分布の自由度は、 \[ N - \mathrm{推定したパラメータの数} \qquad (1.2) \] である。ただし、\(N\) はサンプル・サイズ。しかし、マルチレベル・データの場合、サンプル・サイズは、個人レベルとグループレベルで異なるので、自由度をどう見積もればいいのかも判然としない。もしも OLS の際の自由度で検定を行えば、やはり \(p\) 値を過小に推定し、第一種の過誤を犯しやすくなってしまう。
このような場合に用いられるのが、マルチレベル分析である。マルチレベル分析と通常の回帰分析の一番大きな違いは、マルチレベル分析では、個人レベルだけでなく、グループレベルの残差もモデルに組み込むという点にある。あるいは、個人レベルとグループレベルの回帰式を同時に推定する、といってもいい。すなわち、 \[ \begin{array}{lll} Y_{ij} & = \beta_{0j} + \beta_1 X_{ij} + R_{ij} & (1.3) \\ \beta_{0j} & =\gamma_{00} + \gamma_{01} Z_{j} + U_j & (1.4) \end{array} \] が、最も単純なマルチレベル・モデルで、ランダム切片モデルと呼ばれるものである。(1.3) 式が個人レベル、(1.4) 式がグループレベルの回帰式である。個人レベルは、教育年数に賃金を回帰させるモデルになっているが、通常の回帰分析との違いは、切片の添字に \(j\) が付いている点である。すなわち、切片の位置がグループによって異なると仮定されている。
次にそれぞれのグループの切片を企業規模で予測しているのが、下の方のグループレベルの回帰式である。この式には、\(U_j\) というグループレベルの残差が付いている。すなわち、各企業の切片の位置は、企業規模である程度予測できるかもしれないが、企業規模とは異なる要因によっても影響を受けているかもしれない、ということを仮定したモデルである。この2つの式を同時推定するのがマルチレベル・モデルである。
(1.3) 式に (1.4) 式を代入すると、 \[ \begin{array}{lll} Y_{ij} & = \ \gamma_{00} + \gamma_{01} Z_{j} + U_j + \beta_1 X_{ij} + R_{ij} & \\ & = \ \gamma_{00} + \gamma_{01} Z_{j}+ \beta_1 X_{ij} + U_j + R_{ij} & (1.5) \end{array} \] である。(1.3), (1.4) 式と、(1.5) 式は、表現の仕方が違うだけで、どちらも同じランダム切片モデルを表している。\(U_j\) と \(R_{ij}\) は、測定されないランダム変数なので、ランダム効果と呼ばれる。このモデルがランダム切片モデルと呼ばれるのも、(1.4) 式で示されているように、切片の位置がランダム効果の影響を受けて変化すると仮定されているからである。ランダム効果に対して、独立変数が持つ効果(つまり \(\gamma_{01}\) や \(\beta_1\))を固定効果 (fixed effect) と呼ぶ。
マルチレベル・モデルである (1.5) 式と通常の回帰モデルである (1.1) 式の違いは、係数の表記法の違いを除けば、実質的にはグループレベルの残差である \(U_j\) をモデルに含んでいるかどうかという点である。この \(U_j\) の分散を推定することによって、標準誤差を適切に推定することができるようになる。
マルチレベル分析では、以下のような仮定が置かれる。 \[ \begin{array}{llll} \mathrm{E}(U_j) & = \ \mathrm{E}(R_{ij}) & = \ 0 & (1.6)\\ \mathrm{cov}(X_{ij}, \; R_{ij}) &= \ \mathrm{cov}(Z_{j}, \; R_{ij}) &= \ 0. & (1.7) \\ \mathrm{cov}(X_{ij}, \; U_{j}) &= \ \mathrm{cov}(Z_{j}, \; U_{j}) & = \ 0. & (1.8) \\ \mathrm{cov}(U_j, \; R_{ij}) & = \ 0. & & (1.9) \end{array} \] ただし
である。(1.6) 式は残差の期待値はゼロであり、(1.7) 式と (1.8) 式は独立変数と残差の共分散がゼロであることを示しているが、これらの仮定は OLS の際にも置かれるものなので、特に目新しいものではない。マルチレベル分析で置かれる実際に新しい仮定は (1.9) 式で、レベルの異なるランダム効果どうしの共分散はゼロであると仮定されている。これらの仮定が正しくなければ、分析結果は歪んだものになる可能性がある。
それでは、上の架空のデータを (1.1) 式で表される通常の回帰モデルと、 (1.5) 式のランダム切片モデルの両方で推定してみよう。r のマルチレベル・モデルのための関数は、lme4 パッケージの lmer() 関数を使う。ランダム切片モデルの場合は、
lmer(モデル式 + (1 | グループID変数), data = データフレーム名) # ランダム切片モデルの場合
とすればよい。すなわち、ほとんど lm() と同じように引数を指定すればよいが、グループを識別するための変数を指定する必要がある。独立変数は個人レベルでもグループレベルでも、同じようにモデル式の中で指定すればよい。つまり、(1.5) 式に近い書き方になる。以下が分析例である。なお、切片のみのモデルと、教育年数だけを説明変数として投入したモデルもあわせて推定する。
ols1 <- lm(wage ~ educ + I(firm.size/100), data = data1) # 従業員数を100で割って単位を変換(係数が見やすく)
library(lme4)
mlm0 <- lmer(wage ~ 1 + (1 | id.firm), data = data1)
mlm1 <- lmer(wage ~ educ + (1 | id.firm), data = data1)
mlm2 <- lmer(wage ~ educ + I(firm.size/100) + (1 | id.firm), data = data1)
library(texreg)
screenreg(list(mlm0, mlm1, mlm2, ols1)) # screenreg() は複数の回帰分析の結果をまとめて表示
##
## ===========================================================================
## Model 1 Model 2 Model 3 Model 4
## ---------------------------------------------------------------------------
## (Intercept) 404.69 *** 51.59 11.68 13.65
## (18.36) (27.08) (24.77) (23.83)
## educ 26.95 *** 24.79 *** 24.60 ***
## (1.91) (1.93) (2.02)
## I(firm.size/100) 2.83 *** 2.84 ***
## (0.54) (0.44)
## ---------------------------------------------------------------------------
## AIC 1517.18 1396.74 1378.75
## BIC 1525.66 1408.05 1392.89
## Log Likelihood -755.59 -694.37 -684.38
## Num. obs. 125 125 125 125
## Num. groups: id.firm 25 25 25
## Var: id.firm (Intercept) 6827.53 2114.70 687.61
## Var: Residual 8012.45 3209.57 3181.18
## R^2 0.74
## Adj. R^2 0.74
## RMSE 61.86
## ===========================================================================
## *** p < 0.001, ** p < 0.01, * p < 0.05
上の Model 1 ~ 3 がランダム切片モデルで、Model 4 が OLS での推定結果である。マルチレベル・モデルでは、グループレベルと個人レベルのランダム効果の分散が重要な意味を持つので、しばしば分析結果をまとめる際にも表記される。上の表では “Var: id.firm (Intercept)” がグループレベルのランダム効果 (\(U_j\)) の分散で、“Var: Residual” が個人レベルのランダム効果 (\(R_{ij}\)) の分散である。
Model 3 と 4 を比べると、切片や独立変数の係数は、どちらのモデルでもほぼ同じであることがわかるだろう。しかし、標準誤差は違っており、教育年数 educ の標準誤差はランダム切片モデルのほうが小さく、企業規模 firm.size の標準誤差は、ランダム切片モデルのほうが大きくなっていることがわかるだろう。このようにして、マルチレベル・モデルのほうが正確な推定が可能になることが知られている。
上の Model 1, Model 2 と Model 3 の結果を比較することで、分析結果への理解が深まる。特にランダム効果の分散に注意が必要である。以下では個人レベルと集団レベルの分散の関係について解説しよう。
とすると、一般的に \[ \mathrm{var}(Y_{ij}) = \tau_{00} + \sigma^2 \qquad (1.10) \] という関係が成り立つ。(1.10) 式の証明は割愛するが、比較的簡単な計算で証明できるし、これがわかると構造方程式モデルなどの理解も深まるので、中級者への道を目指す人は挑戦してほしい。
(1.10) 式は、従属変数の分散をグループ内平均の分散とグループ内偏差の分散に分割できるということを意味する。それゆえ、\(\tau_{00} / \mathrm{var}(Y_{ij})\) が大きいほどグループによって従属変数の値が異なるということである。逆に言えば、これが小さいほど、グループによる違いは小さくなる。
実は、上で推定した 切片のみのランダム切片モデル(Model 1)は、まさにこの \(\tau_{00}\) と \(\sigma^2\) を推定するモデルであり、一元配置の分散分析とよく似たモデルである。切片のみのランダム切片モデルを式で表すと、 \[ \begin{array}{lll} Y_{ij} & = \ \beta_{0j} + R_{ij} & (1.11)\\ \beta_{0j} & = \ \gamma_{00} + U_{j} & (1.12) \end{array} \] であり、(1.11) 式に (1.12) 式を代入すると、 \[ Y_{ij} = \gamma_{00} + U_{j} + R_{ij} \qquad (1.13) \] である。\(\gamma_{00}\) は \(Y_{ij}\) のサンプル全体の平均 (\(\bar Y_{\cdot\cdot}\) と表記する) の推定値であり、\(U_{j}\)は 各グループの平均の全体平均からの偏差 (\(\bar Y_{\cdot j} - \bar Y_{\cdot\cdot}\)) の推定値であり、\(R_{ij}\) は、グループ内偏差(同じ \(R_{ij}\) で表記)の推定値である。それゆえ、\(U_{j}\) の分散は \(\tau_{00}\) と一致し、\(R_{ij}\) の分散も \(\sigma^2\) と一致する。これらも証明を省略するが、簡単な計算で証明できる。
このような知識をもってさきほどの推定結果の “Var: id.firm (Intercept)” と “Var: Residual” を見ると、Model 1 はグループレベルの分散が 9095.61 で、個人レベルの分散が 7337.20 であるから、従属変数の分散のうち \[ 9095.61 / (9095.61 + 7337.20) =0.55 \] つまり、55% はグループレベルの分散であることがわかる。
Model 2 では、個人レベルの変数である教育年数を独立変数として投入しているが、Model 1 と比べると、“Var: Residual” だけでなく、“Var: id.firm (Intercept)” もかなり減少しており、企業間の平均的な賃金の違いは、かなり企業による教育年数の違いで説明できることがわかる。Model 3 ではグループレベルの変数である企業規模を投入しているが、グループレベルの残差の分散 “Var: id.firm (Intercept)” はそれによってさらに小さくなったものの、個人レベルの残差の分散は若干の減少があるもののほとんど変化がないことがわかる。
最初の図1.1 を見ると、教育年数の傾きは、規模の小さな企業ではほとんどゼロだが、規模の大きな企業ほどはっきりした正の値をとっているように見えないだろうか。つまり、企業規模と教育年数の交互作用効果があるようにも見える。このような関係は、OLS による回帰分析の場合、 \[ Y_{ij} = \beta_0 + \beta_1 X_{ij} + \beta_2 Z_{j} + \beta_3 X_{ij} Z_{j} + R_{ij}. \qquad (1.14) \] と表せる。なお、以下では \(X_{ij}\) も、\(Z_{j}\) も、それぞれのサンプル全体の平均でセンタリングしてあると仮定する。交互作用効果を推定する場合、センタリングしていないと、係数の解釈が著しく難しくなることが多いからである。さて、このモデルのような個人レベルの変数とグループレベルの変数の交互作用効果を、交差レベル交互作用効果 (cross-level interaction effect) と呼ぶことがある。交差レベル交互作用を仮定する場合、上記のような OLS での推定は不適切である。なぜなら、規模以外の企業の特性も教育年数の傾きに影響を与えている可能性があるからである。このような企業レベルの効果を考慮しないと、やはり標準誤差や p 値の推定が歪んでしまう。そこで、このような状況をマルチレベル・モデルで表現すると、 \[ \begin{array}{lll} Y_{ij} & = \ \beta_{0j} + \beta_{1j} X_{ij} + R_{ij} & (1.15) \\ \beta_{0j} & = \ \gamma_{00} + \gamma_{01} Z_{j} + U_{0j} & (1.16) \\ \beta_{1j} & = \ \gamma_{10} + \gamma_{11} Z_{j} + U_{1j} & (1.17) \end{array} \] である。このようなモデルは、ランダム傾きモデルと呼ばれる。(1.15) 式は (1.3) 式とよく似ているが、教育年数の係数の添字に \(j\) が加わっている。つまり、教育年数の傾きがグループによって異なることを示している。(1.16) 式は (1.4) 式と同じで、 (1.15) 式の切片が、企業規模とグループレベルのランダム効果によって影響を受けるということを示している。(1.17) 式は、教育年数の傾きが企業規模とグループレベルのランダム効果によって影響を受けることを示している。
(1.15) 式に (1.16) 式と (1.17) 式を代入すると、 \[ \begin{array}{lll} Y_{ij} & = \ ( \gamma_{00} + \gamma_{01} Z_{j} + U_{0j}) + (\gamma_{10} + \gamma_{11} Z_{j} + U_{1j}) X_{ij} + R_{ij} & \\ & = \ \gamma_{00} + \gamma_{01} Z_{j} + \gamma_{10}X_{ij} + \gamma_{11} Z_{j} X_{ij} + U_{0j} + U_{1j} X_{ij} + R_{ij} & (1.18) \end{array} \]
ランダム傾きモデルでも、ランダム切片モデルのと同じように、ランダム効果に関して、以下のような仮定が置かれている。
ただし、\(U_{0j}\) と \(U_{1j}\) は相関が仮定され、データから共分散が推定される。このような仮定のもとにランダム傾きモデルを推定してみよう。
R のスクリプトは、
lmer(モデル式 + (ランダム傾きを仮定する変数 | グループID変数), data = データフレーム名)
で、以下が、分析例。
data1$educ.c <- data1$educ - mean(data1$educ) # 教育年数のセンタリング
data1$size.c <- (data1$firm.size - mean(data1$firm.size))/100 # 企業規模をセンタリングして100人単位に
r.itcpt <- lmer(wage ~ educ.c + size.c + (1 | id.firm), data =data1)
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)
screenreg(list(r.itcpt, r.slope1, r.slope2), custom.model.names = c("Model 3", "Model 5", "Model 6"))
##
## ======================================================================
## Model 3 Model 5 Model 6
## ----------------------------------------------------------------------
## (Intercept) 404.69 *** 395.50 *** 391.76 ***
## (7.28) (6.60) (6.11)
## educ.c 24.79 *** 24.70 *** 25.05 ***
## (1.93) (2.63) (1.98)
## size.c 2.83 *** 1.98 *** 2.80 ***
## (0.54) (0.43) (0.44)
## educ.c:size.c 0.56 ***
## (0.12)
## ----------------------------------------------------------------------
## AIC 1378.75 1370.04 1359.43
## BIC 1392.89 1389.84 1382.06
## Log Likelihood -684.38 -678.02 -671.72
## Num. obs. 125 125 125
## Num. groups: id.firm 25 25 25
## Var: id.firm (Intercept) 687.61 435.85 262.00
## Var: Residual 3181.18 2600.97 2598.79
## Var: id.firm educ.c 97.14 24.18
## Cov: id.firm (Intercept) educ.c 205.76 79.59
## ======================================================================
## *** p < 0.001, ** p < 0.01, * p < 0.05
Model 5 は (1.17) 式から \(\gamma_{11} Z_{j}\) を取り除いたモデルで、教育年数の傾きは、平均的な傾きである \(\gamma_{10}\) と、そこからの企業による偏差である \(U_{1j}\) だけで予測されている。Model 5 の推定結果で、まず注目すべきは \(U_{1j}\) の分散を示す “Var: id.firm educ.c” である。これを見ると、97.14 あり、標準偏差は\(\sqrt{97.14}= 9.9\) である。それゆえ、\(U_{1j}\) が正規分布するならば、教育年数の傾きは企業によって異なるが、95% は \(24.70 - 9.9 * 2 = 4.90\) から \(24.70 + 9.9 * 2 = = 44.5\) のあいだで分布していることがわかる。また、AIC や BIC をランダム切片モデルである Model 3 と比較すると、Model 5 のほうが小さいことがわかるだろう。これは Model 5 のほうが当てはまりが良いことを示している(詳しくは後述)。
Model 6 は Model 5 に交差レベル交互作用を加えたもので、有意な交互作用が見られる。これによってグループレベルの分散もかなり減少している。教育年数の傾きに対するランダム効果(\(U_{1j}\))の分散を示す “Var: id.firm educ.c” は、モデル5では、97.14 あったが、モデル6では 24.18 まで減少しており、教育年数の傾きのバラつきのうち、\(1 - 24.18/97.14 =0.75\) は、企業規模で説明できていることになる。
ランダム切片モデルともランダム傾きモデルも、グループによる従属変数の違いを完全に統制しているわけではない。それゆえ、マルチレベル分析をしても、従属変数とモデルに含まれる独立変数の両方と相関する変数がモデルに投入されていなければ、その投入されていない変数がグループ・レベルであれ、個人レヴェルであれ、係数の推定値は歪む。
上の例で言えば、モデル2 と モデル3 の教育年数の係数が違うことに注目すべきである。モデル2 はランダム切片モデルで、グループ・レベルの独立変数は投入されておらず、モデル3 は企業規模がグループ・レベルの独立変数として投入されている。教育年数の係数は、26.95 から 24.79 まで減少しており、企業規模を統制するかどうかで、教育年数の係数が変わってくることがわかる。これは教育年数と企業規模が相関していることが原因である。
以上のように、マルチレベル分析は回帰分析の発展型の一つであり、通常の回帰分析では1つしかないランダム効果を、複数に分割し、複数の係数のグループによるバラつきを推定することを可能にしたものである。これによってより正確な標準誤差と p値の推定が可能になった。また、ランダム効果の分散がモデルによってどう異なるのかを検討することで、独立変数がグループ間のバラつきに及ぼす影響と、グループ内のバラつきに及ぼす影響を別々に検討できるというメリットもあるのである。