REML の導出
線型混合モデルにおける REML の導出方法は二種類あり、一つは誤差対比を用いて尤度関数を分解する方法で、もう一つはベイズ法的な周辺尤度関数に基づく方法である。 本稿では Harville (1974) による周辺尤度関数を使った導出についてまとめる。
記法と最尤推定量 (ML)
期待値ベクトルが $\mathbf{X}\mathrm{\beta}$、分散共分散行列が $\mathbf{V}$ の多変量正規分布に従う観測値確率変数ベクトルを $\mathbf{Y}$、パラメータベクトルを $\mathrm{\theta} = (\mathrm{\beta}', \mathrm{vec}(\mathbf{V})')'$、尤度関数を $L(\mathrm{\theta}) = f(\mathbf{Y}; \mathrm{\theta})$、対数尤度関数を $l(\mathrm{\theta})=\log L(\mathrm{\theta})$ とする。 ただし、$\mathbf{X}$ は $n \times p$ 次元、$\mathrm{\beta}$ は $p \times 1$ 次元、$\mathbf{V}$ は $n \times n$ 次元、$\mathbf{Y}$ は $n \times 1$ 次元とし、$\mathrm{\beta}$ のパラメータ空間を $\mathrm{\Theta}_{\mathrm{\beta}}$ とおく。
$\mathrm{\beta}$ の最尤推定量については、尤度方程式から
$\hat{\mathrm{\beta}}=(\mathbf{X}'\mathbf{V}^{-1}\mathbf{X})^{-1}\mathbf{X}'\mathbf{V}^{-1}\mathbf{Y}$
が得られる。
一方で $\mathbf{V}$ の最尤推定量については、通常はプロファイル尤度法が用いられる。
$\mathbf{V}$ を既知として得られる $\mathrm{\beta}$ の最尤推定量を代入した対数プロファイル尤度関数は
\[
pl(\mathbf{V} \mid \hat{\mathrm{\beta}}(\mathbf{V})) =
-{\dfrac{n}{2}} \log{2\pi} - {\dfrac{1}{2}} \log{\mathbf{V}} - {\dfrac{1}{2}} \,
\mathbf{r}'\mathbf{V}^{-1}\mathbf{r},
\]
であり、これを最大化するような $\hat{\mathbf{V}}$ を数値最適化によって求める。
ただし、$\mathbf{r} = \mathbf{Y} -
\mathbf{X}(\mathbf{X}'\mathbf{V}^{-1}\mathbf{X})^{-1}\mathbf{X}'\mathbf{V}^{-1}\mathbf{Y}$ である。
Harville (1974)
対数プロファイル尤度を用いて推定する部分では $\mathbf{V}$ の推定のみに興味がある。
そのため、$\mathrm{\beta}$ を周辺化した周辺尤度関数 $\int_{\mathrm{\beta} \in \mathrm{\Theta}_{\mathrm{\beta}}}
L(\mathrm{\theta}) \mathrm{d}\mathrm{\beta}$ を用いて REML を導出する。
導出では
\[
(\mathbf{Y}-\mathbf{X}\mathrm{\beta})'\mathbf{V}^{-1}(\mathbf{Y}-\mathbf{X}\mathrm{\beta})=
(\mathbf{Y}-\mathbf{X}\hat{\mathrm{\beta}})'\mathbf{V}^{-1}(\mathbf{Y}-\mathbf{X}\hat{\mathrm{\beta}})+
(\mathrm{\beta}-\hat{\mathrm{\beta}})'(\mathbf{X}'\mathbf{V}^{-1}\mathbf{X})(\mathrm{\beta}-\hat{\mathrm{\beta}})
\]
を用いる。
Proof. \[ \begin{aligned} (\mathbf{Y}-\mathbf{X}\mathrm{\beta})'\mathbf{V}^{-1}(\mathbf{Y}-\mathbf{X}\mathrm{\beta})&= (\mathbf{Y}-\mathbf{X}\hat{\mathrm{\beta}}+\mathbf{X}\hat{\mathrm{\beta}}-\mathbf{X}\mathrm{\beta})'\mathbf{V}^{-1}(\mathbf{Y}-\mathbf{X}\hat{\mathrm{\beta}}+\mathbf{X}\hat{\mathrm{\beta}}-\mathbf{X}\mathrm{\beta}) \\&=(\mathbf{Y}-\mathbf{X}\hat{\mathrm{\beta}})'\mathbf{V}^{-1}(\mathbf{Y}-\mathbf{X}\hat{\mathrm{\beta}})+2(\mathbf{Y}-\mathbf{X}\hat{\mathrm{\beta}})'\mathbf{V}^{-1}(\mathbf{X}\hat{\mathrm{\beta}}-\mathbf{X}\mathrm{\beta})+(\mathbf{X}\mathrm{\beta}-\mathbf{X}\hat{\mathrm{\beta}})'\mathbf{V}^{-1}(\mathbf{X}\mathrm{\beta}-\mathbf{X}\hat{\mathrm{\beta}}) \\&=(\mathbf{Y}-\mathbf{X}\hat{\mathrm{\beta}})'\mathbf{V}^{-1}(\mathbf{Y}-\mathbf{X}\hat{\mathrm{\beta}})+2(\mathbf{Y}'-\mathbf{Y}'\mathbf{V}^{-1}\mathbf{X}(\mathbf{X}'\mathbf{V}^{-1}\mathbf{X})^{-1}\mathbf{X}')\mathbf{V}^{-1}\mathbf{X}(\hat{\mathrm{\beta}}-\mathrm{\beta})+(\mathrm{\beta}-\hat{\mathrm{\beta}})'(\mathbf{X}'\mathbf{V}^{-1}\mathbf{X})(\mathrm{\beta}-\hat{\mathrm{\beta}}) \\&=(\mathbf{Y}-\mathbf{X}\hat{\mathrm{\beta}})'\mathbf{V}^{-1}(\mathbf{Y}-\mathbf{X}\hat{\mathrm{\beta}})+2(\mathbf{Y}'\mathbf{V}^{-1}\mathbf{X}-\mathbf{Y}'\mathbf{V}^{-1}\mathbf{X}(\mathbf{X}'\mathbf{V}^{-1}\mathbf{X})^{-1}\mathbf{X}'\mathbf{V}^{-1}\mathbf{X})(\hat{\mathrm{\beta}}-\mathrm{\beta})+(\mathrm{\beta}-\hat{\mathrm{\beta}})'(\mathbf{X}'\mathbf{V}^{-1}\mathbf{X})(\mathrm{\beta}-\hat{\mathrm{\beta}}) \\&= (\mathbf{Y}-\mathbf{X}\hat{\mathrm{\beta}})'\mathbf{V}^{-1}(\mathbf{Y}-\mathbf{X}\hat{\mathrm{\beta}})+0+(\mathrm{\beta}-\hat{\mathrm{\beta}})'(\mathbf{X}'\mathbf{V}^{-1}\mathbf{X})(\mathrm{\beta}-\hat{\mathrm{\beta}}) \\&= (\mathbf{Y}-\mathbf{X}\hat{\mathrm{\beta}})'\mathbf{V}^{-1}(\mathbf{Y}-\mathbf{X}\hat{\mathrm{\beta}})+(\mathrm{\beta}-\hat{\mathrm{\beta}})'(\mathbf{X}'\mathbf{V}^{-1}\mathbf{X})(\mathrm{\beta}-\hat{\mathrm{\beta}}). \end{aligned} \] $\blacksquare$
上述の関係式を用いれば、周辺尤度関数は
\[
\begin{aligned}
\int_{\mathrm{\beta} \in \mathrm{\Theta}_{\mathrm{\beta}}} L(\mathrm{\theta}) \mathrm{d}\mathrm{\beta}&=
{\dfrac{1}{(2\pi)^{n/2} |\mathbf{V}|^{1/2}}} \int_{\mathrm{\beta} \in \mathrm{\Theta}_{\mathrm{\beta}}}
\exp\left\{-{\dfrac{1}{2}} \,
(\mathbf{Y}-\mathbf{X}\mathrm{\beta})'\mathbf{V}^{-1}(\mathbf{Y}-\mathbf{X}\mathrm{\beta})\right\}
\mathrm{d}\mathrm{\beta}
\\&=
{\dfrac{1}{(2\pi)^{n/2} |\mathbf{V}|^{1/2}}} \exp\left\{-{\dfrac{1}{2}} \,
(\mathbf{Y}-\mathbf{X}\hat{\mathrm{\beta}})'\mathbf{V}^{-1}(\mathbf{Y}-\mathbf{X}\hat{\mathrm{\beta}})\right\}
\int_{\mathrm{\beta} \in \mathrm{\Theta}_{\mathrm{\beta}}} \exp\left\{-{\dfrac{1}{2}} \,
(\mathrm{\beta}-\hat{\mathrm{\beta}})'(\mathbf{X}'\mathbf{V}^{-1}\mathbf{X})(\mathrm{\beta}-\hat{\mathrm{\beta}})\right\}
\mathrm{d}\mathrm{\beta},
\end{aligned}
\]
とかける。
ここで $\mathbf{W}^{-1}=\mathbf{X}'\mathbf{V}^{-1}\mathbf{X}$ ($p \times p$ 次元) とおけば、
\[
\begin{aligned}
\int_{\mathrm{\beta} \in \mathrm{\Theta}_{\mathrm{\beta}}} \exp\left\{-{\dfrac{1}{2}} \,
(\mathrm{\beta}-\hat{\mathrm{\beta}})'\mathbf{W}^{-1}(\mathrm{\beta}-\hat{\mathrm{\beta}})\right\}
\mathrm{d}\mathrm{\beta}
&=
(2\pi)^{p/2} |\mathbf{W}|^{1/2} \int_{\mathrm{\beta} \in \mathrm{\Theta}_{\mathrm{\beta}}}
{\dfrac{1}{(2\pi)^{p/2} |\mathbf{W}|^{1/2}}} \exp\left\{-{\dfrac{1}{2}} \,
(\mathrm{\beta}-\hat{\mathrm{\beta}})'\mathbf{W}^{-1}(\mathrm{\beta}-\hat{\mathrm{\beta}})\right\}
\mathrm{d}\mathrm{\beta}
\\&=
(2\pi)^{p/2} |\mathbf{W}|^{1/2},
\end{aligned}
\]
である (一行目右辺の積分は多変量正規分布の密度関数の全領域と一致するため $1$ になる)。
したがって、
\[
\begin{aligned}
\int_{\mathrm{\beta} \in \mathrm{\Theta}_{\mathrm{\beta}}} L(\mathrm{\theta}) \mathrm{d}\mathrm{\beta}
&=
{\dfrac{(2\pi)^{p/2} |\mathbf{W}|^{1/2}}{(2\pi)^{n/2} |\mathbf{V}|^{1/2}}} \exp\left\{-{\dfrac{1}{2}} \,
(\mathbf{Y}-\mathbf{X}\hat{\mathrm{\beta}})'\mathbf{V}^{-1}(\mathbf{Y}-\mathbf{X}\hat{\mathrm{\beta}})\right\}
\\&=
{\dfrac{|(\mathbf{X}'\mathbf{V}^{-1}\mathbf{X})^{-1}|^{1/2}}{(2\pi)^{(n-p)/2} |\mathbf{V}|^{1/2}}}
\exp\left\{-{\dfrac{1}{2}} \,
(\mathbf{Y}-\mathbf{X}\hat{\mathrm{\beta}})'\mathbf{V}^{-1}(\mathbf{Y}-\mathbf{X}\hat{\mathrm{\beta}})\right\},
\end{aligned}
\]
である。
この周辺尤度関数について対数をとれば、
\[
\log \int_{\mathrm{\beta} \in \mathrm{\Theta}_{\mathrm{\beta}}} L(\mathrm{\theta})
\mathrm{d}\mathrm{\beta}=
-{\dfrac{n-p}{2}} \log{2\pi} - {\dfrac{1}{2}} \log|\mathbf{V}| - {\dfrac{1}{2}}
\log|\mathbf{X}'\mathbf{V}^{-1}\mathbf{X}| - {\dfrac{1}{2}} \mathbf{r}'\mathbf{V}^{-1}\mathbf{r},
\]
が得られる。
これは、誤差対比に基づいて導出した REML の対数プロファイル尤度に等しい。
個人的には周辺尤度に基づいても REML が導かれるというのはなかなか面白く、納得感もあります。
参考文献
- Harville D. Bayesian inference for variance components using only error contrasts. Biometrika 1974; 61(2): 383–385.
履歴
- 2016/05/27 公開