マルチレベル分析
具体的なsoftwareのcodeではなく、極めて簡単な部分からマルチレベル解析について学習する。ただし、簡単なStataのcodeは掲載する。
Multilevel分析とは階層性のあるデータの解析に用いられる。紛らわしいのは、Hierarchical linear model, nested model, mixed model（混合モデル）、そしてrandom modelなど様々な名称が混在していることで、専門分野によって呼び方が異なるのであるが、これらはすべてマルチレベル分析と同一のものを指していると考えられる。
なお、マルチレベル分析は階層ベイズモデルに内包されると考えて良い。
どのようなデータセットに対してマルチレベル分析が用いられるか
研究参加者のデータが階層になっている場合がマルチレベル分析の適応である。例えば、全国規模での臨床試験を考えてみる。複数の病院である薬剤の効果を調べたとき、病院Aの患者は似たような属性を有している（レベル\(2\)）と考えられる。もちろん患者個人の薬剤に対する反応（レベル\(1\)）が重要であるのだが、このように入れ子構造になっているデータに対して、マルチレベル分析は威力を発揮する。入れ子構造としてはいくつあっても構わない。個人の在住する県や東北、関東地方などの地域、さらに国レベルなどと入れ子を拡大していくことができる。
また、パネルデータなど、同じ個人に対して複数のアンケートを行った状況を考えてみる。例えば、3年連続で1年おきに、年収や世帯状況などを調査するときの年収を考えてみる。個人の年収は（平均して）年齢とともに上昇していくことが期待される。この場合、\(3\)ポイントの時系列データがあるので、個人のデータをバラバラに扱うのでなく、反復データとしてまとめて扱いたいときに、同じくマルチレベル分析の出番になる。
通常の回帰モデルとの違いはなにか
通常の回帰分析で用いられるのはOrdinary Least Square (OLS)、つまり最小二乗法推定による計算である。\(X, Y\)などはベクトルとして、次の単純な式で表される。$$Y = {\beta}_0 + {\beta}_1 X+\epsilon$$この式は各々のデータが互いに独立であるならば正しいのであるが、例えば上の図で同じ病院Aの患者たちは同じような治療を受け、相関correlationがあると考えられる。この群ごとの相関のことを級内相関intra-class correlation (ICC)という。逆に言えば、ICCが存在しないときは、通常の回帰分析で十分であるということでもある。
群の数が\(2, 3\)個であれば層別化して回帰分析を行えばいいが、例えば世界\(160\)か国の調査などでは、回帰式を\(160\)個立てることになり、多重検定の問題が生じる（参考）。国ごとのダミー変数を投入しても、そもそも被説明変数の違いをもたらすものはなにかという疑問に対して答えられない。
マルチレベルモデルの導入
説明変数を持たないモデル
もっとも単純なモデルは、説明変数がなく、以下のように表される。$$\begin{eqnarray}Y & = & {\beta}_0 + r\\ {\beta}_0 & = & {\gamma}_{00} + u_0\end{eqnarray}$$ここで、\({\beta}_0\)は\(Y\)の平均値、\(r\)は個人の差異に由来する誤差である。これをレベル\(1\)の誤差という。\({\gamma}_{00}\)は真の平均、\(u_0\)はグループ・群ごとの差異に由来する誤差である。これがレベル\(2\)の誤差である。\(r\)の分散を\({\sigma}^2\)、\(u_0\)の分散を\({\tau}^2\)とすると、級内相関ICCは以下の式で表される。$$ICC = \frac{{\tau}^2}{{\tau}^2+{\sigma}^2}$$
同じエクスポージャーでもアウトカムが異なるとき、その差が個人だけの差に由来する部分と、集団の差に由来する部分に分解される、というのが重要で、上の式はすべての誤差に対する各集団ごとの分散の割合である。定義からICCは\(0\leq ICC\leq 1\)となる。ICCが\(0\)に近い（\(\tau\)が小さい）ということは、級内相関がないということで、つまり集団ごとに分ける意味はない。この場合は単純な回帰分析を行えば良いだろう。逆にICCが\(1\)に近い場合、\(\tau\)は大きく、\(\sigma\)は\(\tau\)に比べて無視できる。個人レベルでの差異を説明するものはなく、各群に属する個人は同一だと考えて良いということだろう。
Stataのcodeはこちらを参照。
説明変数を投入したモデル
同様に\(X, Y\)をベクトルとして、$$Y = {\beta}_0 + {\beta}_1 X + r \tag{a}\label{a}$$とする。ここで、$$\begin{eqnarray}{\beta}_0 & = & {\gamma}_{00} + u_0\\ {\beta}_1 & = & {\gamma}_{10} + u_1\end{eqnarray} \tag{b}\label{b}$$とすると、回帰式の切片にも傾きにもレベル\(2\)の分散が導入されることになる。\(u_0, u_1\)の分散が\(0\)に近い場合、群間の差というのはほとんどないので、わざわざマルチレベル分析を行う価値はないだろう。ちなみに、この場合は変動切片・傾きモデルとも呼ばれる。切片にはレベル\(2\)の分散を導入しない場合（変動傾きモデル）や、傾きにはレベル\(2\)の分散を導入しない場合（変動切片モデル）なども考えられることはすぐに分かる。
\eqref{b}を\eqref{a}に代入すると、$$Y = {\gamma}_{00}+u_0 + ({\gamma}_{10} + u_1)X + r$$となるが、これを書き直すと$$Y = {\gamma}_{00}+ {\gamma}_{10}X + (u_0 + u_1X + r) \tag{c}\label{c}$$となる。\eqref{c}で括弧外の部分が固定効果、括弧内の部分がランダム効果（個人の差異によらない部分）となる。マルチレベル回帰式をこのような形で表すことを混合表記という。
固定効果、というのは有限の要素によってその効果が決定されることをいう。ランダム効果というのは、母集団のとり方（群の分け方）によって効果が異なる、サンプリングの方法によって誤差があるということの言明である。こうやってみると、通常の回帰分析との違いがよくわかる。通常の回帰分析では切片と傾きはすべて固定効果によって決定されるが、マルチレベル分析では切片や傾きに対してグループレベルの効果（ランダム効果）を導入しているというわけである。
Stataでのコード
Stataでのxtmixedコマンドは、デフォルトでmle（最尤法）オプションになっている。
まず、説明変数を代入しないモデルが以下。式は\(y_i = a_{j[i]} + {\epsilon}_i\)となる。
xtmixed y || group(level 2):, mle nolog
estat iccまた、ランダム切片モデルはこちら。式は\(y_i = a_{j[i]} + \beta x_i + {\epsilon}_i\)となる。
xtmixed y x1 || group:, mle nolog
estat iccランダム切片・傾きモデルはこちら。式は\(y_i = a_{j[i]} + \beta_{j[i]} x_i + {\epsilon}_i\)となる。傾きにランダム効果を導入したい個体レベルの変数をgroup:の後に記載していく。
xtmixed y x1 || group: x1, mle nolog covariance(unstructure)
estat iccランダム傾きモデルはこちら。式は\(y_i = a + \beta_{j[i]} x_i + {\epsilon}_i\)となる。これは結構式が難しい。
xtmixed y x1 || _all: R.x1, mle nolog
estat iccモデル比較を行う場合は、解析を行った後適当に名前をつけて保存する。
est(imates) store (name)
est(imates) stats *アステリスクを入力すると、メモリに保存されたモデルをすべて比較できる。
xtreg (fe, be, re)との違い
Stataでパネルデータを解析するとき、”xtreg”というコマンドもある。このコマンドとxtmixedの違いに混乱している人も多いだろう（自分も）。ここでは以下のUCLAのwebsiteを参考に考えてみる。
データのimportはこちら。
use https://stats.idre.ucla.edu/stat/stata/faq/spf24, replaceこれは\(8\)人の被験者に各々\(4\)回の調査を行ったものである。\(y\)はスコア、\(a\)は\(1\)か\(2\)でグループを表す。8人の被験者が\(4\)人\(4\)人のグループに分けられている。\(b\)は各々の調査で、\(1, 2, 3, 4\)まである。\(s\)は個人IDになる。デザインは、被験者内因子(within)と被験者間因子(between)の両方を持つ混合モデルである。被験者内因子\((b)\)は\(4\)水準で、被験者間因子\((a)\)は\(2\)水準である。分析を簡単にするために、ここでは\(a\cdot b\)の交互作用項を考慮しない。
within-subject
まずはwithin-subjectをみるため、”xtreg, fe”コマンドを用いる。”fe”オプションはfixed-effectsを表し、これは被験者内因子と同一である。コードはこちら。
xtreg y i.a i.b, i(s) fe
固定効果モデルでは時間によって変化しない変数（例えば性別や出生時の体重など）は推定されない。ここでも、グループ変数\(a\)は時間によって変化していないので、Coef. は0 (omitted)になっている。また、グループ変数\(a\)は時間によって変化せず、誤差との相関がないという仮定なので、R-sqのbetweenの部分は空白になっている。
between-subject
つぎはbeオプション。
xtreg y i.a i.b, i(s) be
both within and between
これはxtreg, reコマンドになる。グループ\(a\)と個人の各計測の両方を推定している。この結果が”xtmixed”の結果と一致する。
xtreg y i.a i.b, i(s) re
xtmixed y i.a i.b || s:
なお、この例ではxtreg, feもreもbeも同じ結果（同じ係数）になるが、これはwithin-effectsとbetween-effectsが直行するためであり、一般的には結果は異なる。
関連記事
参考記事
参考にした「理論と方法」のpdf。以下のblogも役に立つ。
コメント