[Stata][Multileve]Stataでのマルチレベル分析について

マルチレベルMultilevel分析

色々な資料はインターネット上に転がっているが、ここではPrinceton大学の講義資料に沿って学習する。

https://www.princeton.edu/~otorres/Multilevel101.pdf

Stataのxtmixedについてのリンクはこちら。

まとめ

マルチレベル解析はデータに二つ以上の階層構造があるときに用いる。例えば多数の病院での臨床試験や、国ごとの国際アンケート調査などのデータがこれに当たる。マルチレベル分析を用いることにより、

• グループや単位ごとに異なる効果を検証できる。
• グループレベルでの平均を推定できる

ことが挙げられる。また、利点としては

• 通常の回帰分析では無視されてしまうグループや単位間での平均の違いを検出できる。
• グループの数が多くなると回帰分析ではサンプル数が足りなくなったり、一般化ができなくなるという問題を回避できる。

があるだろう。

実際のデータを使って手を動かす。

use https://dss.princeton.edu/training/schools.dta, replace

でデータをimportする。データは下のようになる。

変数としては、”school”, “student”, “y”, “x1”, “x2”, “x3″の\(6\)つがある。学生の数は\(4,059\)人。”school”は学校のIDで全部で\(65\)校、”student”はそのまま個々の生徒のIDであるが学校ごとにそれぞれ振られている（つまり、学校\(1\)のID\(1\)と、学校\(2\)のID\(1\)が存在する）。”y”は点数（マイナスもあるよう）、”x1″はReadingの試験の点数（これもマイナスがある）、”x2″は性別で女性が\(1\)、最後の”x3″は学校のタイプで”Mixed”, “Boys only”, “Girls only”の\(3\)種類。以下のコマンドでまず学校ごとの平均点の変数を生成する。

bysort school: egen y_mean = mean(y)

以下のコマンドで個々の学生のスコア\((y)\)と学校ごとのスコアの平均\((y_{mean})\)をグラフにできる。

twoway scatter y school, msize(tiny) || connected y_mean school, connect(L) clwidth(thick) clcolor(black) mcolor(black) msymbol(none) || , ytitle(y)

通常の回帰分析

まずは学校ごとの回帰分析、つまり各学校に属する個人のReadingテストの成績がy(Score)とどのような線形関係にあるのかをグラフにしてみる。

statsby inter = _b[_cons] slope = _b[x1], by(school) saving(ols, replace): regress y x1
sort school
merge school using ols
drop _merge
gen yhat_ols = inter + slope*x1
sort school x1
seperate y, by(school)
seperate yhat_ols, by(school)
twoway connected yhat_ols1-yhat_ols65 x1 || lfit y x1, clwidth(thick) clcolor(black) legend(off) ytitle(y) 

出力されるグラフがこちら。

通常の回帰分析では階層構造が消失してしまい、単なる学校ごとの分析になってしまっている。そこで以下でマルチレベル分析を行う。

変動切片モデル Varying-intercept model (null)

最も単純なモデルで、\(y_i = {\alpha}_{j[i]} + \epsilon_{i}\)というモデルで表される。つまり、傾きは持たず、群の間で切片だけが異なるという推定モデルになる。

xtmixed y || school:, mle nolog

“mle”はmaximum likelihood estimationの略になる。Stataのmultilevel分析はdefaultでmleになっている。制限付き最尤推定法を指定するときは”reml”にする。結果は以下。”nolog”オプションでlogの出力を抑制できる。

上の結果で\(y\)の下にある_consが全体の切片である。schoolにあるsd(_cons)が学校レベル（レベル\(2\)）の標準偏差、sd(Residual)が個人レベル（レベル\(1\)）の標準偏差である。”Prob >= chibar2 = 0.0000″がRandom-effects = 0という帰無仮説に対する検定結果になる。つまり、\(\alpha = 0.05\)とすれば、この帰無仮説は棄却されるので、ランダム効果ありという解釈すれば良い。

級内相関Intraclass correlation (ICC)は以下の式で表示される。$$\begin{eqnarray}ICC & = &\frac{({\text{sigma_u}})^2}{(\text{sigma_u})^2 + (\text{sigma_e})^2} \\ & = & \frac{{\text{sd(_cons)}}^2}{{\text{sd(_cons)}}^2+{\text{sd(residual)}}^2} \\ & = & \frac{{4.11}^2}{{4.11}^2+{9.21}^2}\\ & = & 0.17\end{eqnarray}$$もしもICCが\(0\)に近づくと、集団ごとのグルーピングは無意味になっていく。つまり、単純な回帰を行ったほうが良いということになる。逆に、ICCが\(1\)に近づくと、群内の個人個人が全員同一で、個人レベルでの差異を説明するものはないということになる。

StataでICCを計算するときは、”xtmixed”コマンドの後に、

estat icc

とすると良い。

なお、もとのpdfではStata 13からxtmixedのコマンドはmixedに変わったとあるが、mixedで解析を行うと標準偏差ではなく分散が表示される。

変動切片モデル Varying-intercept model (one level-1 predictor)

これは群間で切片だけが異なり、傾きは同一とみなす推定モデルになる。式で表すと\(y_i = {\alpha}_{j[i]} + \beta x_i + {\epsilon}_i\)となる。

xtmixed y x1 || school:, mle nolog

\(x1\)の部分の係数(Coef.)が全体の共通の傾き\(\beta\)である。その他の結果の解釈については上と同一になる。また、$$\begin{eqnarray}ICC & = &\frac{({\text{sigma_u}})^2}{(\text{sigma_u})^2 + ({\text{sigma_e}})^2} \\ & = & \frac{{\text{sd(_cons)}}^2}{{\text{sd(_cons)}}^2+{\text{sd(residual)}}^2} \\ & = & \frac{{3.03}^2}{{3.03}^2+{7.52}^2}\\ & = & 0.14\end{eqnarray}$$である。

変動切片・傾きモデル Varying-intercept, varying-coefficient model

いわゆる混合モデル（Miexed model）である。群間で切片も傾きも異なるとみなし推定する。式は\(y_i = {\alpha}_{j[i]} +{\beta}_{j[i]}x_i + {\epsilon}_i\)となる。

xtmixed y x1 || school: x1, mle nolog covariance(unstructure)

$$\begin{eqnarray}ICC & = &\frac{({\text{sigma_u}})^2}{(\text{sigma_u})^2 + ({\text{sigma_e}})^2} \\ & = & \frac{{\text{sd(_cons)}}^2+sd(x1)}{{\text{sd(_cons)}}^2+sd(x1)+{\text{sd(residual)}}^2} \\ & = & \frac{{0.12}^2+{3.01}^2}{{0.12}^2+{3.01}^2+{7.44}^2}\\ & = & 0.14\end{eqnarray}$$である。

covarianceオプションは変数として”ind(ependent)”, “ex(changeable)”, “id(dntity)”, “un(structured)”の\(4\)つがある。defaultでは”independent”が採用されている。このオプションはランダム効果の共分散行列の構造を指定し、ランダム効果の式ごとに指定することができる。各ランダム効果式に指定することができる。”independent”オプションの共分散構造では、ランダム効果式内の各ランダム効果に個別の分散を許容し、すべての共分散が\(0\)であることを仮定する。”exchangeable”オプションの構造では、すべてのランダム効果に対して 1 つの共通分散と 1 つの共通ペア共分散が指定される。”identity”は”multiple of the indetity”の省略であり、つまりすべての分散が等しく，すべての共分散がゼロであることである。”unstructured”では、すべての分散と共分散を区別することができる。方程式が\(p\)個のランダム効果項で構成されている場合，”unstructured”の共分散行列は\(\displaystyle \frac{p(p+1)}{2}\)個の一意なパラメータを持つ。

変動傾きモデル Varying-slope model

上のVarying-intercept, varying-coefficient modelとは違い、切片は群間で異なるとはしない。式で書くと\(y_i = \alpha + {\beta}_{j[i]}x_i + {\epsilon}_i\)である。

xtmixed y x1 || _all: R.x1, mle nolog

Stataでは”_all”は予約語で、すべての変数を意味する。”R.”は変数x1に対してランダム効果を課している。R.を付けないと結果は同じであるが、分散が計算されない（ランダム効果を想定していないので当然）。

Likelihood ratioを用いたモデルの比較

尤度比を用いた検定(lrtest)でmleでfitしたモデルの比較を行う。この検定は二つのモデルが優位に異なるかどうか、モデルの対数尤度を比較する。

*Fitting random intercepts and storing results
quietly xtmixed y x1 || school:, mle nolog
est(imates) store  ri

*Fitting random coefficients and storing results
quietly xtmixed y x1 || school: x1, mle nolog covariance(unstructure)
est(imates) store rc

*Running the likelihood-ratio test to compare
lrtest ri rc, stats

この結果の解釈であるが、帰無仮説は二つのモデルに差はない、である。”Prob > chi2″が\(0.05\)よりも小さいときに帰無仮説を棄却する。この場合、rcモデルの方が対数尤度が小さく、すなわちAICもBICも小さくなっており、よりfitしたモデルと言える。

これまで出てきた\(4\)つのモデル全てについて、lrtestを行うと、rcモデルが最もfitしたモデルであることが言える。

Varying-intercept, varying-coefficient modelを使ったPost estimation

xtmixed y x1 || school: x1, mle nolog covariance(unstructure) variance

とすると、標準偏差ではなく分散が表示される。この分散を使って分散共分散行列を表示することができる。

estat recovariance

相関を見るにはcorrelationオプションを付ける。

estat recovariance, corr(elation)

切片と\(x1\)の相関を見ると、\(y\)と\(x1\)の平均の”close relationship”を表している。

ランダム効果の推定

ランダム効果\(u\)を推定するには、predictコマンドとreffectsオプションを用いる。こうすることでランダム効果のBest Linear unbiased predictions (BLUPs)を与えることができる。これは、切片と傾きの変化の変動量を表す。数式では、$$y_i = {\alpha}_{j[i]} +{\beta}_{j[i]}x_i + {\epsilon}_i$$を$$y_i = {\alpha}_{j[i]} + {\beta}_{j[i]}x_i + u_{{\alpha}_i} + u_{{\beta}_{j[i]}} + {\epsilon}_i$$と、Fixed-effectsである\({\alpha}_{j[i]} + {\beta}_{j[i]}x_i\)およびRandom-effectsである\(u_{{\alpha}_i} + u_{{\beta}_{j[i]}} + {\epsilon}_i\)に分解することに相当する。

xtmixed y x1 || school: x1, mle nolog covariance(unstructure)
predict u*, reffects

とすると、新しい変数\(u1, u2\)が生成される。\(u1\)は\(u_{\beta}\)、\(u2\)は\(u_{\alpha}\)である。具体的な数値を入れると、$$y_i = -0.12 + 0.56×1 + u_{\alpha} + u_{\beta}$$となる。

bysort school: gen(erate) groups = (_n == 1)
list school u2 u1 if school <= 10 & groups

とすることで学校ごとの\(u2, u1\)をリスト表示できる。

例えば学校\(1\)の場合、$$\begin{eqnarray}y_1 & = & -0.12 + 0.56×1 + 3.75 + 0.12×1 \\ & = & (-0.12 + 3.75) + (0.56 + 0.12)x1 \\ & = & 3.63 + 0.68 x1 \end{eqnarray}$$と推定式を作ることができる。

これを学校ごとにすべて行うには、

gen intercept = _b[_cons] + u2
gen slope = _b[x1] + u1
list school intercept slope if school <= 10 & groups

となる。学校ごとの推定値は

gen yhat = intercept + (slope*x1)

とするか、xtmixedの後に

predict yhat_fit, fitted

とすると良い。

list school yhat yhat_fit if school <= 10 & groups

plotも簡単。

twoway connected yhat_fit x1 if school <= 10, connect(L)

以下のコマンドで残差プロットも描ける（図は省略）。xtmixedの後に行うこと。

predict resid, residuals
predict resid_std, rstandart
qnorm resid_std