マルチレベルデータの解析方法(2)：混合効果モデル（lme4パッケージ）

マルチレベルデータとは
- マルチレベルデータの解析
混合効果モデルの概要
- 混合効果モデル = 固定効果 + 変量効果
- クラスター効果をモデルに含める
lme4パッケージで混合効果モデルを使う
- サンプルデータの準備
- glmer( )を使ってモデルを当てはめる
おわりに
参考資料

マルチレベルデータとは

「通常の」解析*1では、「ある患者のアウトカムは、他の患者のアウトカムから影響を受けない」と考えて解析している。もう少し正確に言えば、アウトカムの変化のうち、説明変数で説明できない誤差の部分に関しては「互いに独立である（mutually independent）」という前提で解析している。

では、「アウトカムの誤差が互いに独立である」とは考えにくい状況とはどういったものだろうか。

1つ目は、アウトカムが類似する傾向がある亜集団（=クラスター）が存在する場合である。例えば、研究対象集団がいくつかので世帯で構成されているとき、食事や生活習慣などは世帯内で類似する傾向があるだろう。このようなとき、世帯が異なる研究参加者同士は独立と考えてよいが、同一世帯内の参加者を独立と考えるのは無理がある。世帯の他にも、学校や地域などがクラスターになりうる。

2つ目は、同一対象者において反復してアウトカムが測定される場合である。この場合のデータの単位は「人」ではなく「評価時点」になる。同じ人から測定されたアウトカムは類似した傾向を持つと考える方が普通だろう。眼や四肢・指趾といった「臓器単位」あるいは個々の腫瘍といった「病変単位」でアウトカムをみるような研究の場合も、疾患のメカニズムによっては反復測定と同様に考える必要がある。各測定に対して、人がクラスター単位となっていると考えることもできる。

アウトカムの誤差の出どころが、「測定単位（観測単位）」とその上層の「クラスター単位」の複数箇所にわたる、という意味でマルチレベルデータと呼ばれることがある。解析手法の観点からは、以下の用語はいずれも同じものを指していると考えて差し支えない。

クラスターデータ（clustered data）
反復測定データ（repeated-measurement data）= パネルデータ（panel data）
マルチレベルデータ（multi-level data）

マルチレベルデータの解析

マルチレベルデータを解析する際には、誤差の相関構造や階層性を考慮する必要である。それは以下の理由のため。

データの情報量を過大評価しないため：
アウトカムの誤差に相関があるということは、観測単位間で情報の一部が重複しているということ。これを考慮しないと情報量を水増ししていることになってしまう。

クラスターを考慮した適切なモデルを使うことで検出力アップ：
観測値の誤差の中でクラスター効果として説明できることがあれば、クラスター単位での異質性を考慮することで検出力が上がるかもしれない。

解析する方法としては大きく以下の2種類がある。

一般化推定方程式（Generalized Estimating Equation, GEE）
混合効果モデル（Mixed-effects model, MEM）

前者は「誤差間の相関構造を定義して解析する方法」で、個々の観測単位（の誤差）間にある相関構造行列を与えることで対応するアプローチ。
後者は「誤差の発生源を階層化する方法」で、上層（例：世帯）が与えられた条件下では下層（例：世帯構成員）の誤差は独立であると考えるアプローチ。
この記事では後者のみ解説する。

混合効果モデルの概要

混合効果モデルといってもモデルの基本は一般化線形モデル（Generalized Linear Model, GLM）なので、分布族（family）やリンク関数（link function）についての説明は割愛。

混合効果モデル = 固定効果 + 変量効果

固定効果（fixed effect, FE）と変量効果（random effect, RE）の両方を含んだモデルを混合効果モデルと呼ぶ。現実的には固定効果を含まないモデルはないので、変量効果を含むモデルを混合効果モデルと呼んでいる。

では固定効果と変量効果の違いは何か。実は統一された定義はないみたい。なので固定効果・変量効果の話をするときはお互いにどういう意味で使っているのか注意が必要（特に計量経済学分野と疫学分野は違う気がする）。

Gelman先生によると、次のような定義・解釈があるとのこと。

FEは対象者間で共通の効果、REは対象者間で変動する
FEは興味ある効果、REはそれ自体に興味がないもの
集団の大部分に関連するならFE、無視できるくらいごく一部にしか関連しないならRE
効果が変数の実現値だと仮定できるのであればRE
FEは最小二乗法・最尤法で推定されるが、REはshrinkageで推定される

この後の混合効果モデルの説明を読むときには、定義1をイメージするのが分かりやすくなっていいと思う。あと定義4に関連して、「FEのパラメータは確かにモデルの係数（coefficient）だが、REは厳密に言えばパラメータではなく観察されていない変数（unobserved random variable）である」という点は重要。

クラスター効果をモデルに含める

傾向が類似する亜集団をここではクラスターと呼んでいるので、その亜集団の違いがもたらすアウトカムへの影響の違いを「クラスター効果」と呼ぶことにする。*2

同一対象者に複数回アウトカムを測定したデータを想定する（対象者=クラスター）。対象者 $i$ の時点 $t$ における説明変数を $X_{it}$ 、アウトカム変数を $Y_{it}$ として、次の式のような回帰モデルを当てはめる。

$Y_{it} = \beta_0 + \beta_1 X_{it} + \varepsilon_{it}$

ここで、 $\beta_0$ は「切片」、 $\beta_1$ は「傾き」、 $\varepsilon_{it}$ はモデルから期待されるアウトカム $E[Y_{i}$ ]と測定値 $Y_i$ の間の「誤差」を表している。このモデルでは、集団全体におけるXとYの関係性を１つの直線で表している。

この場合は個人間の差は無視して、集団全体で画一的な切片と傾きを設定して推定した。でも、個人毎にスタート地点（=切片）は異なっているし、増加速度（=傾き）にも緩急がある。これを考慮してやる方が、適切な推定値がえられるかもしれないし、推定誤差が小さくて済むかもしれない。

このクラスターによる影響を変量効果としてモデルに組み込んだものが、切片だけにクラスター効果を考慮した「ランダム切片モデル（random intercept model）」や、切片と傾きの両方にクラスター効果を考慮した「ランダム切片+傾きモデル（random intercept and slope model）」である。

クラスター効果は変量効果に含めるものとして説明してきたが、必ずしも決まっているわけではない。固定効果としてもモデルに組み入れることができる。例えば、切片のみにクラスター効果を含めたモデルを考える。

固定効果として含めたときのモデル式は次のように書ける。

$\begin{align} Y_{it} &= (\beta_0 + \beta_{0,i}^{'}) + \beta_1 X_{it} + \varepsilon_{it} \\\ \varepsilon_{it} &\sim N(0, \sigma^2) \end{align}$

このモデルでは、

n人の対象者のクラスター効果はn-1 個の $\beta_{0,i}^{'}$ として導入されていて*3、 $\beta_0$ と合わせてn個分の切片を別々に設定している
誤差は1種類だけ

ある対象者のクラスター効果 $\beta_{0,i}^{'}$ は、その対象者の測定値のみから推定される。他の対象者のデータに制約を受けないが、新しく加わった対象者については全く予想のしようがない。

これに対して、変量効果として含めたときのモデル式は次のように書ける。

$\begin{align} Y_{it} &= (\beta_0 + b_{i}) + \beta_1 X_{it} + \varepsilon_{it} \\\ b_i &\sim N(0, \sigma_b^2) \\\ \varepsilon_{it}|b_i &\sim N(0, \sigma^2) \end{align}$

このモデルでは、

n人の対象者のクラスター効果がn個の $b_i$ として導入されている*4
誤差を「対象者間誤差 $b_i$ 」と「対象者内誤差 $\varepsilon_{it}$ 」に分解している。誤差の合計が最小になるような値を推定する。
クラスター効果が与えられた（=同一対象者内である）条件下では、対象者内誤差は互いに独立

前述の固定効果モデルと違い、ある対象者のクラスター効果 $b_i$ は他の対象者のデータからも緩やかに制約を受ける。他の対象者の情報を適度に借りているとも考えられる。また、新しく加わった対象者についても、 $\sigma_b^{2}$ のお陰で多少予想が立つ。

lme4パッケージで混合効果モデルを使う

サンプルデータの準備

ChickWeightというデータセットを使う。ヒヨコの餌と体重変化の関係性をみたデータらしい（ヒヨコ50匹、測定数578回）。

> head(ChickWeight)
Grouped Data: weight ~ Time | Chick
  weight Time Chick Diet
1     42    0     1    1
2     51    2     1    1
3     59    4     1    1
4     64    6     1    1
5     76    8     1    1
6     93   10     1    1

ヒヨコ毎に体重は類似していると考える方が自然だろう。なのでChickについてクラスター効果を考える。

glmer( )を使ってモデルを当てはめる

lme4パッケージのglmer( )を使う。「GLMじゃなくてLMでいいならlmer( )を使ってよ！」とエラーメッセージが出るが、無視してここではglmer( )で統一しておく。

変量効果については、指定したいクラスターを|の後ろに書いて指定する。どこにクラスター効果を設定するかについて3パターンだけ紹介する（複数のクラスター効果の階層化についてはまた別の機会に）。固定効果の部分の書き方はglm( )と同じだが、切片があることを意識して1を明示的に書いた。

ランダム切片モデル（random intercept model）

切片部分にだけクラスター効果を入れたければ、固定効果の式の後に"1|クラスター"を追加する。固定効果と区別するために変量効果は"()"で囲む方がよい。

library(lme4)
fit_1 <- glmer(weight ~ 1+Time + (1|Chick),
             data = ChickWeight,
             family = gaussian(link="identity"))

結果は次のとおり。

> summary(fit_1)
Linear mixed model fit by REML ['lmerMod']
Formula: weight ~ 1 + Time + (1 | Chick)
   Data: ChickWeight

REML criterion at convergence: 5619.4

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-3.0086 -0.5528 -0.0888  0.4975  3.4588 

Random effects:
 Groups   Name        Variance Std.Dev.
 Chick    (Intercept) 717.9    26.79   
 Residual             799.4    28.27   
Number of obs: 578, groups:  Chick, 50

Fixed effects:
            Estimate Std. Error t value
(Intercept)  27.8451     4.3877   6.346
Time          8.7261     0.1755  49.716

Correlation of Fixed Effects:
     (Intr)
Time -0.422

まず、当てはめられたモデルや推定に関する概要が表示される。次に、Random effects:で変量効果についての結果（説明は次のモデルで）、Fixed effects:で固定効果についての結果が表示されている。

各クラスターに割り当てられた変量効果を確認したいときはranef( )を使う。

> head(ranef(fit_1)$Chick)
   (Intercept)
18   0.2754554
16 -26.3026906
15 -25.2830406
13 -50.5775531
9  -38.3765089
20 -40.8929742

ランダム切片＋傾きモデル, REに相関あり（random intercept and slope model, RE correlated）

Timeの係数についてもクラスター効果を入れたければ、変量効果のカッコの中に1とTimeを入れる。

fit_2 <- glmer(weight ~ 1+Time + (1+Time|Chick),
               data = ChickWeight,
               family = gaussian(link="identity"))

> summary(fit_2)
Linear mixed model fit by REML ['lmerMod']
Formula: weight ~ 1 + Time + (1 + Time | Chick)
   Data: ChickWeight

REML criterion at convergence: 4827.5

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-2.6731 -0.5563 -0.0277  0.5010  3.4939 

Random effects:
 Groups   Name        Variance Std.Dev. Corr 
 Chick    (Intercept) 140.54   11.855        
          Time         14.14    3.761   -0.95
 Residual             163.51   12.787        
Number of obs: 578, groups:  Chick, 50

Fixed effects:
            Estimate Std. Error t value
(Intercept)  29.1780     1.9573   14.91
Time          8.4531     0.5408   15.63

Correlation of Fixed Effects:
     (Intr)
Time -0.871

変量効果の部分にCorr（相関係数）が表示されている。これは「切片に対するRE」と「Timeの傾きに対するRE」の相関係数で、強い負の相関がある（つまり、「ベースラインが高い方が時間による増加が緩やか」）ことがわかる。

ResidualのStd.Dev.は、対象者が特定された条件下で、各測定値がどれくらいバラついているかを示している。

ランダム切片＋傾きモデル, REに相関なし（random intercept and slope model, RE uncorrelated）

1つ前のモデルでは、切片と傾きのREの間に相関があることを想定したが、クラスター効果の間に関連がないことを想定したい場合は次のいずれかの方法で書く。

誤：(1 | Chick) + (Time | Chick)
正：(1 | Chick) + (Time-1 | Chick)
正：(1 | Chick) + (0+Time | Chick)

切片は含まれることがデフォルトなので、切片を含まないことを明示しなければいけない。

fit_3 <- glmer(weight ~ 1+Time + (1|Chick)+(0+Time|Chick),
             data = ChickWeight,
             family = gaussian(link="identity"))

> summary(fit_3)
Linear mixed model fit by REML ['lmerMod']
Formula: weight ~ 1 + Time + (1 | Chick) + (0 + Time | Chick)
   Data: ChickWeight

REML criterion at convergence: 4890.5

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-2.5601 -0.5614 -0.0059  0.4823  3.4968 

Random effects:
 Groups   Name        Variance Std.Dev.
 Chick    (Intercept) 114.95   10.722  
 Chick.1  Time         12.29    3.506  
 Residual             166.06   12.886  
Number of obs: 578, groups:  Chick, 50

Fixed effects:
            Estimate Std. Error t value
(Intercept)   29.048      1.825   15.92
Time           8.466      0.507   16.70

Correlation of Fixed Effects:
     (Intr)
Time -0.079