はじめに
千葉大学・株式会社Nospareの川久保です．
以前の記事で，線形モデルと線形混合モデルにおけるAIC（赤池情報量規準，Akaike information criterion）型の情報量規準について紹介しました．そのときは，線形混合モデルの説明は簡単にしかできませんでしたので，この記事で応用例も交えながら解説したいと思います．
線形混合モデルとは
モデルのかたち
線形混合モデルは一般に以下のように書けます:
$$
y = X\beta + Zb + \varepsilon, \quad b \sim (0,G), \quad \varepsilon \sim (0,R), \tag{1}
$$
$X\beta$の項を固定効果（fixed effect），$Zb$の項を変量効果（random effect）と呼び，2つの効果が"混ざった"モデルということで，線形混合モデルと呼ばれています．$G$と$R$はそれぞれ，$b$と$\varepsilon$の分散共分散行列です．$b$や$\varepsilon$は，しばしば正規分布の仮定がなされます．$Z$のデザインによって，個人効果や空間効果などを自由にモデリングできますが，個別のモデル例は後の「応用例」のsectionで紹介します．
BLUP
線形モデルにおいて「OLSがBLUEである」ことは，米倉先生の記事でも解説されていますが，線形混合モデルにおいてはBLUEではなくBLUP（ブラップ）という概念が出てきます．適当な既知の列ベクトル$c$と$d$に対して，混合効果 $\mu = c^\top\beta + d^\top b$ を予測したいとき，線形で不偏なクラスの予測量で最良な（分散の小さい）予測量のことをBLUP (Best Linear Unbiased Predictor, BLUP) といい，以下のかたちになることが知られています:
$$
\hat{\mu} = c^\top \hat{\beta} + d^\top \hat{b},
$$
ただし，
\begin{align}
&\hat{\beta} = (X^\top \Sigma^{-1}X)^{-1} X^\top \Sigma^{-1} y, \quad \hat{b} = GZ^\top \Sigma^{-1}(y - X\hat{\beta}), \\
&\Sigma = Var(y) = ZGZ^\top + R,
\end{align}
です．なぜBLUE（推定量）ではなくBLUP（予測量）なのかというと，混合効果$\mu$には変量効果（すなわち確率変数）である$b$が含まれているため，$\mu$自体が確率変数だからです．確率変数は通常，「推定する」とは言わず「予測する」と言います．実用上は，$G$や$R$のなかに含まれている分散パラメータが未知であるため，その推定値を$\hat{\mu}$にプラグインして用います．これをEBLUP（empirical BLUP，イーブラップ）といいます．EBLUPは$b$や$\varepsilon$に分布の仮定をおかなくても導かれますが，もし$b$や$\varepsilon$に正規分布を仮定すると，二乗損失のもとでの$\mu$のempirical best predictor（$E[\mu \mid y]$に$\beta$と分散成分の推定値をプラグインしたもの）として同じ形が導かれます．
分散成分の推定については，変量効果を積分消去した$y$の周辺分布$y \sim (X\beta, \Sigma)$にもとづいた最尤（Maximum Likelihood, ML）法，および不偏に近い推定量が得られるように修正した尤度にもとづいた制限最尤（Restricted ML, REML）法がよく用いられます．
応用例
経時データ・パネルデータ
複数の個人に対して，ある変数を時間経過とともに繰り返し観測したデータのことを，経時データやパネルデータと言います．主に生物統計・医療統計系では経時データ，計量経済学ではパネルデータという用語が使われている印象です．経時データ・パネルデータは，しばしば線形混合モデルでモデリングされます．
$y_{it} \ (i=1,\dots,n; \ t=1,\dots,T)$ を個人$i$の$t$時点での観測だとします．共変量として$x_{it}$が利用可能であるとき，
$$
y_{it} = \beta_0 + x_{it} \beta_1 + b_i + \varepsilon_{it}, \quad b_i \overset{\text{iid}}{\sim} (0,\tau^2), \quad \varepsilon_{it} \overset{\text{iid}}{\sim} (0,\sigma^2), \tag{2}
$$
というモデルを考えます．$b_i$は$x_{it}$で捉えられない個人$i$特有の効果（個別効果）と解釈でき，例えば個人$i$の観測できない体質だったり能力だったりなどが想定されます．このモデルは線形混合モデルの1つで，行列表記で(1)式のように表現することが可能です．このモデルでは，$V(y_{it}) = \tau^2 + \sigma^2, \ Cov(y_{it},y_{it'}) = \tau^2 \ (t \not= t')$であり，同一個体の異時点のデータに相関を持たせることができます．なお，計量経済学では，個別効果を変量効果ではなく固定効果としてモデリングすることが多いです．
(2)式のモデルは切片のみがrandom（変量効果が入っている）ということで，random intercept modelとも呼ばれますが，その発展として$x_{it}$にかかる係数（傾き）もrandomにした
\begin{align}
&y_{it} = (\beta_0 + b_{0i}) + x_{it}(\beta_1 + b_{1i}) + \varepsilon_{it}, \tag{3}\\
&b_{0i} \sim (0,\tau_0^2), \quad b_{1i} \sim (0,\tau_1^2), \quad \varepsilon_{it} \sim (0,\sigma^2),
\end{align}
といったrandom coefficient modelも考えられます．また個人$i$の効果だけでなく，時点$t$の効果も入れたい場合，
$$
y_{it} = \beta_0 + x_{it}\beta_1 + b_i + v_t + \varepsilon_{it}, \quad b_i \sim (0,\tau_b^2), \quad v_t \sim (0,\tau_v^2), \quad \varepsilon \sim (0,\sigma^2)
$$
といったモデルも考えられます．すべて線形混合モデル(1)式で行列表記可能です．
クラスターデータ
経時データ・パネルデータにおける観測$y_{it}$の添字を$t$から$j$に変えた$y_{ij} \ (i=1\dots,n; \ j=1,\dots,n_i)$を考え，$i$番目の地域（クラスター）における$j$番目の個人（家計）の観測だとします．例えば，総務省統計局が行う家計調査では，$y_{ij}$として地域$i$の$j$番目の家計の所得が観測されます．このような観測に対して，(2)式と同様のモデル
$$
y_{ij} = \beta_0 + x_{ij} \beta_1 + b_i + \varepsilon_{ij}, \quad b_i \sim (0,\tau^2), \quad \varepsilon_{ij} \sim (0,\sigma^2), \tag{4}
$$
を考えると，$b_i$は地域効果（クラスター効果）と解釈できます．標本調査において，全国では十分なサンプルが得られていても，地域ごとのサンプルサイズ$n_i$が小さい場合には，地域ごとの標本平均$n_i^{-1}\sum_{j=1}^{n_i}y_{ij}$では地域平均（地域ごとの平均所得など）の推定が不安定となります．そこで(4)式のようなモデルを考え，EBLUPにより地域平均を予測します．線形混合モデルなどの統計モデルの力を借りて，地域（クラスター）ごとの平均などの母集団特性を，小標本から安定的に推定することを目指す問題を小地域推定と言います．この話題については，また別の記事で詳しく解説したいと思います．
Rによる解析例
Rで線形混合モデルを用いた解析を行いたい場合，lme4というパッケージが便利です．cAICの解説でも利用したlme4内のsleepstudyというデータセットの解析例をいくつか見ていきたいと思います．
library(lme4)
fm1 <- lmer(Reaction ~ Days + (1 | Subject), sleepstudy)
睡眠不足の日数Daysが反応時間Reactionに与える影響を調べています．(2)式のモデルにおいて，$y_{it}$を$i$番目の個人Subjectの$t$期目の反応時間，$x_{it} = t-1$（1期目を睡眠不足0日目と数えているため）としたモデルです．Rのformulaにおける(1 | Subject)でrandom intercept $b_i$の項を指定しています．fm1に解析結果が格納されていますが，その要約は以下です．
> summary(fm1)
Linear mixed model fit by REML ['lmerMod']
Formula: Reaction ~ Days + (1 | Subject)
Data: sleepstudy
REML criterion at convergence: 1786.5
Scaled residuals:
Min 1Q Median 3Q Max
-3.2257 -0.5529 0.0109 0.5188 4.2506
Random effects:
Groups Name Variance Std.Dev.
Subject (Intercept) 1378.2 37.12
Residual 960.5 30.99
Number of obs: 180, groups: Subject, 18
Fixed effects:
Estimate Std. Error t value
(Intercept) 251.4051 9.7467 25.79
Days 10.4673 0.8042 13.02
Correlation of Fixed Effects:
(Intr)
Days -0.371
結果のRandom effectsの項目から，(2)式における$\tau^2$の推定値が1378.2，$\sigma^2$の推定値が960.5であることが読み取れます．これらの分散パラメータは，lmer()関数のデフォルトではREML法により推定されます．またFixed effectsの項目から，(2)式における$\beta_1$の推定値が10.4673であることがわかり，睡眠不足の日数が1日伸びると，平均的に反応時間が10.4673（ミリ秒）遅くなると解釈されます．変量効果の予測値$\hat{b}_i$を見たい場合，ranef()関数を用います．
> ranef(fm1)
$Subject
(Intercept)
308 40.783710
309 -77.849554
310 -63.108567
330 4.406442
331 10.216189
332 8.221238
333 16.500494
334 -2.996981
335 -45.282127
337 72.182686
349 -21.196249
350 14.111363
351 -7.862221
352 36.378425
369 7.036381
370 -6.362703
371 -3.294273
372 18.115747
308や309といった数字はSubject変数の値で，経時データにおける個人のIDみたいなものです．308のSubjectの$b_i$の予測値は40.783710，309のSubjectの$b_i$の予測値は-77.849554であり，308のSubjectは反応時間が長く，309のSubjectは反応時間が短いという個体差があることが分かります．
切片だけでなく傾きにも変量効果を入れたrandom coefficient model (3)式を当てはめたい場合は，
fm2 <- lmer(Reaction ~ Days + (Days | Subject), sleepstudy)
とします．ただし，この書き方だと(3)式における$b_{0i}$と$b_{1i}$の相関を仮定しているモデルになるので，独立だと仮定したモデルを当てはめたい場合は，
fm3 <- lmer(Reaction ~ Days + (1 | Subject) + (0 + Days | Subject), sleepstudy)
または
fm4 <- lmer(Reaction ~ Days + (Days || Subject), sleepstudy)
と書いてあげます．fm3のformulaにおける(0 + Days | Subject)という書き方は，「切片には変量効果を入れず，Daysの係数にのみ変量効果を入れる」という書き方です．この書き方を使えば，傾きにのみ変量効果を入れた
fm5 <- lmer(Reaction ~ Days + (0 + Days | Subject), sleepstudy)
という推定もできます．
おわりに
株式会社Nospareには，統計学の様々な分野を専門とする研究者が所属しております．統計アドバイザリーやビジネスデータの分析につきましては株式会社Nospare までお問い合わせください．
