#はじめに
線形混合効果モデルにおける固定効果βの最尤推定量とその分散の導出について、基本中の基本だが、参考書にはその導出が載っていないし、日本語のウェブ上でも見当たらなかったのでまとめました（英語ならたくさんあります）。

基本的には自分自身の勉強のためであり、定義があいまい部分があります。

線形混合効果モデル

##相関を考慮できる
線形混合効果モデルは階層構造をもつ結果変数Yの相関を考慮して解析するためのモデルの一つである。ここでの階層構造とは、例えば複数施設からそれぞれ患者をサンプリングする場合や、経時データにおける患者の繰り返し測定を行うような状況を指す。前者は施設間内で相関が発生し、後者では個人内で相関が発生するため、両者とも結果変数Yが独立であると仮定している通常の回帰分析や分散分析では不適切であるといえる。
##モデルの設定（notation）
線形混合効果モデルを次のように設定
$\mathbf{Y_i = X_i β + Z_ib_i+ ε　}$
${\mathbf{Y_i} \sim MVN( \mathbf{X_i β, ZGZ'+R_i })　　\mathbf{b_i} \sim MVN(0,\mathbf{G})　　\mathbf{ε} \sim MVN(0,\mathbf{R_i})}$
分散成分：$\mathbf{V_i(\alpha)} = \mathbf{ZGZ'+R_i} とおく（ \alphaはパラメータ）$

$\mathbf{Y_i}：結果変数ベクトル(T_i×1）$
$\mathbf{X_i}：固定効果のデザイン行列(T_i×p)$
$\mathbf{β}：固定効果(p×1)$
$\mathbf{Z_i}：変量効果のデザイン行列(T_i×q)$
$\mathbf{b_i}：変量効果(q×1)$
$\mathbf{ε}：誤差(T_i×1)$
$\mathbf{G}：変量効果の分散共分散行列(q×q)$
$\mathbf{R}：誤差の分散共分散行列(T_i×T_i)$

$\mathbf{X_i}の転置を\mathbf{X_i'}で表現$

εなど、計算では関係ありませんので気にしないでください。

最尤推定量の導出

尤度関数
$$L( \mathbf{β,θ;Y} )=\prod_{i=1}^n \frac{1}{(2\pi)^{T_i}}
\frac{1}{|\mathbf{V_i(\alpha)}|^\frac{1}{2}}
\exp{ \left[-\frac{1}{2} (\mathbf{Y_i - X_i β})' \mathbf{V_i ^{-1} (\alpha) }(\mathbf{Y_i - X_i β}) \right]}$$
対数尤度関数
$$\log{ L(\mathbf{β,θ;Y}) } =
\sum_{i=1}^n \left[ \log{\frac{1}{(2\pi)^{T_i}} } +
\log{\frac{1}{|\mathbf{V_i(\alpha)}|^\frac{1}{2}}}
-\frac{1}{2} (\mathbf{Y_i - X_i β})' \mathbf{V_i ^{-1} (\alpha)} (\mathbf{Y_i - X_i β})\right]$$

$ここで、\mathbf{V_i(\alpha)}のパラメータ\alpha=\alpha_0で固定して
　\mathbf{W_i} = \mathbf{V_i ^{-1} (\alpha_0) }とおく。$

$βと関係ない部分を消去して$

$$\frac{\partial}{\mathbf{\partial β}}\log{ L(\mathbf{β,θ;Y}) } =
\frac{\partial}{\partial β}\sum_{i=1}^n \left[ -\frac{1}{2} (\mathbf{Y_i - X_i β})' \mathbf{W_i} (\mathbf{Y_i - X_i β})\right]$$

行列部分を展開すると

$$\frac{\partial}{\mathbf{\partial β}}\log{ L(\mathbf{β,θ;Y}) }=\frac{\partial}{\mathbf{\partial β}} \sum_{i=1}^n \left[ -\frac{1}{2} (\mathbf{Y_i' - β' X_i'}) \mathbf{W_i} (\mathbf{Y_i - X_i β})\right]$$

$$=\frac{\partial}{\mathbf{\partial β}}\sum_{i=1}^n \left[ -\frac{1}{2} (\mathbf{Y_i' W_i Y_i - β ' X_i W_i Y_i - Y_i' W_i X_i β + β' X_i' W_i X_i β})\right]$$

$\mathbf{β}のベクトル微分をすると$

$$\frac{\partial}{\partial β}\log{ L(\mathbf{β,θ;Y}) }=\sum_{i=1}^n \left[ -\frac{1}{2} (\mathbf{-X_i' W_i Y_i - X_i' W_i Y_i + 2 X_i'W_i X_i β})\right]$$

$$=\sum_{i=1}^n \left[\mathbf{ X_i' W_i Y_i - X_i'W_i X_i β}\right]$$

$\sum_{i=1}^n [ \mathbf{X_i' W_i Y_i - X_i'W_i X_i β}] = \mathbf{0} とおくと$

$$\sum_{i=1}^n \left( \mathbf{X_i'W_i X_i} \right) \mathbf{β} = \sum_{i=1}^n \mathbf{X_i' W_i Y_i} $$

$両辺に(\sum_{i=1}^n \mathbf{X_i'W_i X_i} )^{-1}を左から掛けると $

$$ \mathbf{\widehat β} = \left(\sum_{i=1}^n \mathbf{X_i'W_i X_i}\right)^{-1} \sum_{i=1}^n \mathbf{X_i' W_i Y_i}$$

最尤推定量の分散の導出

$ 最尤推定量 \mathbf{\widehat β} の分散Var\left[\mathbf{\widehat β}\right]$がフィッシャー情報量の逆数に一致することを利用する。
フィッシャー情報量は対数尤度の一回微分のマイナスをとった期待値であるから

$$E\left[-\frac{\partial^2}{\mathbf{\partial β^2}}\log{ L(\mathbf{β,θ;Y}) } \right]
= \sum_{i=1}^n \mathbf{X_i' W_i Y_i}$$

$$Var\left[ \mathbf{ \widehatβ } \right] = \left( \sum_{i=1}^n \mathbf{X_i' W_i Y_i } \right) ^{-1}$$

#まとめ
線形混合効果モデルにおいて
最尤推定量
$$\mathbf{\widehat β} = \left(\sum_{i=1}^n \mathbf{X_i'W_i X_i}\right)^{-1} \sum_{i=1}^n \mathbf{X_i' W_i Y_i}$$

その分散
$$Var \left[ \mathbf{ \widehatβ } \right] = \left( \sum_{i=1}^n \mathbf{X_i' W_i Y_i } \right) ^{-1}$$

$Var\left[\mathbf{\widehat β}\right]$の導出においてはそれがフィッシャー情報量の逆数に一致することを利用した。