同じ人から複数回測定したデータを、すべて独立したデータとして解析すると間違った結論を導きます。混合効果モデルはこの問題に対処するための手法ですが、数式だけ見ていてもなかなか理解しにくいです。この記事では、シミュレーションデータを自分で作成し、単純回帰と混合効果モデルの結果を比較することで、混合効果モデルが何をしているのかを体感します。 ## 混合効果モデルとは 混合効果モデル（Mixed-effects model）は、データに階層構造（クラスター構造）がある場合に使用する統計モデルです。「固定効果」と「変量効果（ランダム効果）」の両方を含むため、この名前がついています。 - **固定効果（Fixed effects）**: すべての対象に共通する効果。推定したい母集団のパラメータ - **変量効果（Random effects）**: 個人やグループごとに異なる効果。ある分布からランダムに抽出されたと仮定 典型的な例として以下のようなデータがあります。 - 同一患者の繰り返し測定データ（縦断研究） - 学校内の生徒データ（生徒は学校にネストされている） - 施設ごとの患者データ（多施設共同研究） ## なぜ混合効果モデルが必要か ### 問題：独立性の仮定違反 通常の回帰分析は、観測値が互いに独立であることを仮定しています。しかし、同じ人から複数回測定したデータでは、同一人物のデータは互いに相関しています。 例えば、血圧を3回測定した場合を考えます。 - 高血圧の人は3回とも高い値を示す傾向がある - 低血圧の人は3回とも低い値を示す傾向がある このような相関を無視して通常の回帰分析を行うと、標準誤差が過小評価され、誤った結論を導く可能性があります。 ### 解決策：分散の分離 混合効果モデルでは、総変動を以下の2つに分離します。 $ \text{総変動} = \text{個人間変動} + \text{個人内変動} $ ![[mixed_model_concept.png]] 上図は混合効果モデルの構造を示しています。まず、大きな青い分布（個人間分布）から各個人の「真の血圧」が決まります。次に、その真の値を中心とした小さな分布（個人内分布）から実際の測定値が生成されます。点線は各個人の真の値を、下の点は実際の観測データを表しています。 数式で表すと、観測値 $y_{ij}$（個人 $i$ の $j$ 回目の測定）は次のようにモデル化されます。 $ y_{ij} = \mu + u_i + e_{ij} $ - $\mu$：全体平均（固定効果） - $u_i \sim N(0, \sigma^2_{\text{between}})$：個人 $i$ の変量効果（個人間変動） - $e_{ij} \sim N(0, \sigma^2_{\text{within}})$：残差（個人内変動） ## シミュレーションデータで理解する 実際にデータを生成して、混合効果モデルの動作を確認します。 ### データ生成の設定 以下の条件でシミュレーションデータを作成します。 - 30人の被験者、各3回測定（計90観測） - 全体平均血圧：120 mmHg - 個人間SD：15 mmHg（人によるばらつきが大きい） - 個人内SD：3 mmHg（同じ人の測定誤差は小さい） ```stata clear all set seed 12345 * パラメータ設定 local n_subjects = 30 local n_obs = 3 local grand_mean = 120 local sigma_between = 15 local sigma_within = 3 * データ生成 set obs `=`n_subjects' * `n_obs'' gen id = ceil(_n / `n_obs') bysort id: gen time = _n * 個人間変動（idごとに固定の値） gen u_between = . forvalues i = 1/`n_subjects' { local random_effect = rnormal(0, `sigma_between') replace u_between = `random_effect' if id == `i' } * 真の血圧値と観測値 gen bp_true = `grand_mean' + u_between gen e_within = rnormal(0, `sigma_within') gen bp = bp_true + e_within ``` ### 生成されたデータの構造 ``` +--------------------------------------------------------------------+ | id time grand_~n u_between bp_true e_within bp | |--------------------------------------------------------------------| 1. | 1 1 120 14.388 134.388 -2.732075 131.6559 | 2. | 1 2 120 14.388 134.388 3.314332 137.7023 | 3. | 1 3 120 14.388 134.388 1.445811 135.8338 | |--------------------------------------------------------------------| 4. | 2 1 120 -.0270755 119.9729 1.863381 121.8363 | 5. | 2 2 120 -.0270755 119.9729 .2746946 120.2476 | 6. | 2 3 120 -.0270755 119.9729 .7975975 120.7705 | |--------------------------------------------------------------------| 7. | 3 1 120 8.161156 128.1612 -1.891799 126.2694 | 8. | 3 2 120 8.161156 128.1612 1.739509 129.9007 | 9. | 3 3 120 8.161156 128.1612 -.0756916 128.0855 | +--------------------------------------------------------------------+ ``` このデータの特徴を確認します。 - `u_between`（個人間変動）は同じIDでは同じ値 - `e_within`（個人内変動）は測定ごとに異なる - 個人1は真の血圧が134.4 mmHgで、3回の測定値は131.7、137.7、135.8とほぼその周辺に分布 - 個人2は真の血圧が120.0 mmHgで、測定値は121.8、120.2、120.8とやはり真の値の周辺に分布 ## 分析1：単純回帰（各人1回目のデータのみ） まず、各人から1回目の測定値のみを使って単純回帰を行います。 ```stata preserve keep if time == 1 regress bp restore ``` **出力結果：** ``` Source | SS df MS Number of obs = 30 -------------+---------------------------------- F(0, 29) = 0.00 Model | 0 0 . Prob > F = . Residual | 5669.20596 29 195.489861 R-squared = 0.0000 -------------+---------------------------------- Adj R-squared = 0.0000 Total | 5669.20596 29 195.489861 Root MSE = 13.982 ------------------------------------------------------------------------------ bp | Coefficient Std. err. t P>|t| [95% conf. interval] -------------+---------------------------------------------------------------- _cons | 125.6904 2.55271 49.24 0.000 120.4696 130.9113 ------------------------------------------------------------------------------ ``` - 推定された平均血圧：125.7 mmHg - 標準誤差：2.55 - 理論上の標準誤差：$15 / \sqrt{30} = 2.74$ 30人の独立したサンプルからの推定なので、標準誤差は個人間SDを$\sqrt{n}$で割った値に近くなります。 ## 分析2：混合効果モデル（全データ） 次に、全90観測を使って混合効果モデルを推定します。 ```stata mixed bp || id:, stddeviations ``` **出力結果：** ``` Mixed-effects ML regression Number of obs = 90 Group variable: id Number of groups = 30 Obs per group: min = 3 avg = 3.0 max = 3 Wald chi2(0) = . Log likelihood = -296.01862 Prob > chi2 = . ------------------------------------------------------------------------------ bp | Coefficient Std. err. z P>|z| [95% conf. interval] -------------+---------------------------------------------------------------- _cons | 126.2381 2.497834 50.54 0.000 121.3424 131.1338 ------------------------------------------------------------------------------ ------------------------------------------------------------------------------ Random-effects parameters | Estimate Std. err. [95% conf. interval] -----------------------------+------------------------------------------------ id: Identity | sd(_cons) | 13.53998 1.784845 10.45713 17.53169 -----------------------------+------------------------------------------------ sd(Residual) | 3.395909 .3100026 2.839569 4.06125 ------------------------------------------------------------------------------ LR test vs. linear model: chibar2(01) = 137.88 Prob >= chibar2 = 0.0000 ``` ### 結果の解釈 **固定効果部分：** - 推定された平均血圧：126.2 mmHg - 標準誤差：2.50（単純回帰とほぼ同じ） **変量効果部分：** - 個人間SD（`sd(_cons)`）：13.5（設定値15に近い） - 個人内SD（`sd(Residual)`）：3.4（設定値3に近い） ### 分散成分の比較 | パラメータ | 設定値 | 推定値 | |:-----------|:------:|:------:| | 個人間SD | 15.0 | 13.5 | | 個人内SD | 3.0 | 3.4 | 混合効果モデルは、設定したパラメータをほぼ正確に推定できています。 ## ICC（級内相関係数） ICC（Intraclass Correlation Coefficient）は、総変動のうち個人間変動が占める割合を示します。 $ \text{ICC} = \frac{\sigma^2_{\text{between}}}{\sigma^2_{\text{between}} + \sigma^2_{\text{within}}} $ ```stata estat icc ``` **出力結果：** ``` Intraclass correlation ------------------------------------------------------------------------------ Level | ICC Std. err. [95% conf. interval] -----------------------------+------------------------------------------------ id | .940819 .0179763 .8941036 .9676714 ------------------------------------------------------------------------------ ``` ### ICCの比較 | | 理論値 | 推定値 | |:--|:------:|:------:| | ICC | 0.962 | 0.941 | ICCが0.94ということは、観測された変動の94%が個人間の差異に由来し、同一人物内の測定誤差は6%に過ぎないことを意味します。 ### ICCの解釈の目安 | ICC | 解釈 | |:----|:-----| | < 0.40 | 低い相関 | | 0.40 - 0.75 | 中程度の相関 | | > 0.75 | 高い相関 | 今回のデータはICC = 0.94と非常に高く、同じ人を繰り返し測定しても新しい情報はほとんど得られないことを示しています。 ### ICCが小さい例 比較のため、個人間変動と個人内変動が同程度（どちらもSD=10）の場合を見てみます。 ```stata * パラメータを変更 local sigma_between = 10 local sigma_within = 10 ``` この設定でシミュレーションを行うと、以下のようなデータになります。 ``` +---------------------------------------------------------+ | id time u_between bp_true e_within bp | |---------------------------------------------------------| 1. | 1 1 9.592 129.592 -9.106916 120.4851 | 2. | 1 2 9.592 129.592 11.04777 140.6398 | 3. | 1 3 9.592 129.592 4.819369 134.4114 | |---------------------------------------------------------| 4. | 2 1 -.0180503 119.9819 6.21127 126.1932 | 5. | 2 2 -.0180503 119.9819 .9156487 120.8976 | 6. | 2 3 -.0180503 119.9819 2.658658 122.6406 | +---------------------------------------------------------+ ``` 同じ人の中でも測定値のばらつきが大きくなっています（個人1：120.5, 140.6, 134.4）。 混合効果モデルの結果： ``` ------------------------------------------------------------------------------ Random-effects parameters | Estimate Std. err. [95% conf. interval] -----------------------------+------------------------------------------------ id: Identity | sd(_cons) | 6.495376 1.791062 3.783471 11.15111 -----------------------------+------------------------------------------------ sd(Residual) | 11.3197 1.033343 9.465231 13.53751 ------------------------------------------------------------------------------ ``` ``` Intraclass correlation ------------------------------------------------------------------------------ Level | ICC Std. err. [95% conf. interval] -----------------------------+------------------------------------------------ id | .2477016 .1185843 .0864188 .5340348 ------------------------------------------------------------------------------ ``` ICCは0.25と低く、観測値の変動の約25%しか個人間の差異で説明できません。このような場合、同じ人を繰り返し測定することで、より多くの情報が得られます。 ## 重要なポイント ### 1. 標準誤差の比較 | 分析 | サンプルサイズ | 標準誤差 | |:-----|:-------------:|:--------:| | 単純回帰（1回目のみ） | 30 | 2.55 | | 混合効果モデル（全データ） | 90 | 2.50 | サンプルサイズが30から90に3倍になっても、標準誤差はほとんど変わりません。 ### 2. なぜ精度が向上しないのか もし90の観測が独立であれば、標準誤差は次のようになるはずです。 $ SE = \frac{SD}{\sqrt{n}} = \frac{14}{\sqrt{90}} \approx 1.48 $ しかし、実際には同一人物内のデータは相関しているため、「有効サンプルサイズ」は30に近い値になります。 ### 3. 混合効果モデルの利点 平均値の推定精度向上には限界がありますが、混合効果モデルには以下の利点があります。 - **分散の分離**: 個人間変動と個人内変動を分離して推定できる - **ICCの算出**: クラスター内の相関の程度を定量化できる - **正しい推論**: クラスター構造を考慮した正しい標準誤差を得られる - **欠損値への対応**: 測定回数が個人ごとに異なっても対応可能 - **自由度の節約**: 個人を固定効果として扱うと人数分のパラメータが必要だが、変量効果なら分散パラメータ1つで済む ## Stataでの混合効果モデルの基本構文 ```stata mixed 従属変数 独立変数 || グループ変数: ``` ### オプション | オプション | 説明 | |:-----------|:-----| | `stddeviations` | 分散ではなく標準偏差で結果を表示 | | `variance` | 分散で結果を表示（デフォルト） | | `mle` | 最尤法（デフォルト） | | `reml` | 制限付き最尤法 | ### 事後コマンド | コマンド | 説明 | |:---------|:-----| | `estat icc` | ICCを計算 | | `predict, fitted` | 予測値を計算 | | `predict, reffects` | 変量効果（BLUPs）を計算 | ## 参考文献 - Rabe-Hesketh S, Skrondal A. Multilevel and Longitudinal Modeling Using Stata. 3rd ed. Stata Press; 2012. - Stata Manual: `help mixed` ## 完全なコード 以下のコードをコピーしてdo-fileとして保存すれば、この記事の内容を再現できます。 ```stata *------------------------------------------------------------------------------- * 混合効果モデル学習用シミュレーション *------------------------------------------------------------------------------- clear all set more off * パラメータ設定 local n_subjects = 30 local n_obs = 3 local grand_mean = 120 local sigma_between = 15 local sigma_within = 3 * 再現性のためseed設定 set seed 12345 *------------------------------------------------------------------------------- * 1. データ生成 *------------------------------------------------------------------------------- set obs `=`n_subjects' * `n_obs'' * ID（個人識別子）と時点を生成 gen id = ceil(_n / `n_obs') bysort id: gen time = _n * 個人間変動（idごとに固定） gen u_between = . forvalues i = 1/`n_subjects' { local random_effect = rnormal(0, `sigma_between') replace u_between = `random_effect' if id == `i' } * 真の血圧値と観測値 gen bp_true = `grand_mean' + u_between gen e_within = rnormal(0, `sigma_within') gen bp = bp_true + e_within * ラベル付け label variable id "個人ID" label variable time "測定時点" label variable bp_true "真の血圧 (mmHg)" label variable bp "観測血圧 (mmHg)" * データ確認 list id time u_between bp_true e_within bp in 1/9, sepby(id) *------------------------------------------------------------------------------- * 2. 単純回帰（各人1回目のデータのみ） *------------------------------------------------------------------------------- display _n "=== 単純回帰（n=30） ===" preserve keep if time == 1 regress bp restore *------------------------------------------------------------------------------- * 3. 混合効果モデル（全データ） *------------------------------------------------------------------------------- display _n "=== 混合効果モデル（n=90） ===" mixed bp || id:, stddeviations *------------------------------------------------------------------------------- * 4. ICC *------------------------------------------------------------------------------- display _n "=== ICC ===" estat icc display _n "理論上のICC: " %5.3f `sigma_between'^2 / (`sigma_between'^2 + `sigma_within'^2) ```