線形混合効果モデルの同時・周辺尤度関数

目的

本稿は，線形混合モデル \[ \g{Y}_i=\g{X}_i\gx{\beta}+\g{Z}_i\g{b}_i+\gx{\epsilon}_i, \] \[ \begin{pmatrix} \g{b}_i \\ \gx{\epsilon}_i \end{pmatrix} \sim N\left( \begin{pmatrix} \g{0} \\ \g{0} \end{pmatrix}, \begin{pmatrix} \g{G} & \g{O}\\ \g{O} & \g{R}_i \end{pmatrix} \right), \] における同時尤度関数と周辺尤度関数の関係について簡単にメモしたものである． $\g{Y}_i$は$n$次元ベクトル，$\g{X}_i$は$n\times p$行列，$\gx{\beta}$は$p$次元ベクトル，$\g{Z}_i$は$n\times q$行列，$\g{b}_i$は$q$次元ベクトル，$\gx{\epsilon}_i$は$n$次元ベクトルとする．

条件付き分布

条件付き分布は \[ f(\g{Y}_i; \gx{\beta}, \g{R}_i \mid \g{b}_i)= \frac{f(\g{Y}_i,\g{b}_i; \gx{\beta}, \g{R}_i, \g{G})}{f(\g{b}_i; \g{G})}, \] で定義されている．ここでは，線形混合モデルのモデル式を密度関数に代入したときの，$f(\g{Y}_i; \gx{\beta}, \g{R}_i \mid \g{b}_i)$の計算結果を示しておく．

$\g{Y}_i$と$\g{b}_i$の同時分布と$\g{b}_i$の分布

線形混合モデルのもとでは， \[ \begin{pmatrix} \g{b}_i \\ \g{Y}_i \end{pmatrix} \sim N\left( \begin{pmatrix} \g{0} \\ \g{X}_i\gx{\beta} \end{pmatrix}, \g{\Sigma}_i \right), \] \[ \g{\Sigma}_i= \begin{pmatrix} \g{G} & \g{G}\g{Z}_i'\\ \g{Z}_i\g{G} & \g{V}_i \end{pmatrix}, \g{V}_i=\g{R}_i+\g{Z}_i\g{G}\g{Z}_i', \] が成り立つ．

以下では，この導出について示す． $\g{b}_i$は定義から$N(\g{0}, \g{G})$に従う． $\g{Y}_i=\g{X}_i\gx{\beta}+\gx{\epsilon}_i^*,~ \gx{\epsilon}_i^*=\g{Z}_i\g{b}_i+\gx{\epsilon}_i$より，$\g{Z}_i\g{b}_i$と$\gx{\epsilon}_i$は独立に多変量正規分布にしたがう確率変数の線形結合であるため，$\g{Y}_i \sim N(\g{X}_i\gx{\beta}, \g{V}_i)$である．あとは$\Cov [\g{b}_i, \g{Y}_i]$を求めればよく， \[ \E\left[ \begin{pmatrix} \g{b}_i \\ \g{Y}_i - \g{X}_i\gx{\beta} \end{pmatrix} \begin{pmatrix} \g{b}_i' & (\g{Y}_i - \g{X}_i\gx{\beta})' \end{pmatrix} \right] = \E\left[ \begin{pmatrix} \g{b}_i\g{b}_i' & \g{b}_i(\g{Y}_i - \g{X}_i\gx{\beta})' \\ (\g{Y}_i - \g{X}_i\gx{\beta})\g{b}_i' & (\g{Y}_i - \g{X}_i\gx{\beta})(\g{Y}_i - \g{X}_i\gx{\beta})' \end{pmatrix} \right] \] となる． $\E[\g{b}_i\g{b}_i']=\g{G}$および$\E[(\g{Y}_i - \g{X}_i\gx{\beta})(\g{Y}_i - \g{X}_i\gx{\beta})']=\g{V}_i$は定義より明らかである．あとは$\E[\g{b}_i\g{Y}_i' - \g{b}_i\gx{\beta}'\g{X}_i']$を求めればよく，$\E[\g{b}_i]=\g{0}$および$\g{b}_i$と$\gx{\epsilon}_i$の独立性から， \[ \E[\g{b}_i\g{Y}_i'-\g{b}_i\gx{\beta}'\g{X}_i']= \E[\g{b}_i\gx{\beta}'\g{X}_i'+\g{b}_i\g{b}_i'\g{Z}_i'+\g{b}_i\gx{\epsilon}_i'-\g{b}_i\gx{\beta}'\g{X}_i']= \E[\g{b}_i\g{b}_i'\g{Z}_i']=\g{G}\g{Z}_i', \] となる．よって，$\Cov [\g{b}_i, \g{Y}_i]=\g{\Sigma}_i$となることを示すことができる．

各密度関数は \[ f(\g{Y}_i,\g{b}_i; \gx{\beta}, \g{R}_i, \g{G})=\frac{1}{(2\pi)^{(n+q)/2}|\g{\Sigma}_i|^{1/2}} \exp\left\{ -\frac{1}{2} \begin{pmatrix} \g{b}_i \\ \g{Y}_i-\g{X}_i\gx{\beta} \end{pmatrix}' \g{\Sigma}_i^{-1} \begin{pmatrix} \g{b}_i \\ \g{Y}_i-\g{X}_i\gx{\beta} \end{pmatrix} \right\}, \] \[ f(\g{b}_i; \g{G})=\frac{1}{(2\pi)^{q/2}|\g{G}|^{1/2}} \exp\left( -\frac{1}{2} \g{b}_i'\g{G}^{-1}\g{b}_i \right), \] で定義される．

条件付き分布$\g{Y}_i|\g{b}_i$

以下では，条件付き分布 \[ \begin{aligned} f(\g{Y}_i; \gx{\beta}, \g{R}_i \mid \g{b}_i)&= \frac{f(\g{Y}_i,\g{b}_i; \gx{\beta}, \g{R}_i, \g{G})}{f(\g{b}_i; \g{G})} \\&= \frac{1}{(2\pi)^{n/2}|\g{R}_i|^{1/2}} \exp\left\{ -\frac{1}{2} \begin{pmatrix} \g{Y}_i-\g{X}_i\gx{\beta}-\g{Z}_i\g{b}_i \end{pmatrix}' \g{R}_i^{-1} \begin{pmatrix} \g{Y}_i-\g{X}_i\gx{\beta}-\g{Z}_i\g{b}_i \end{pmatrix} \right\}, \end{aligned} \] を計算する．

ブロック行列の行列式の性質から， \[ |\g{\Sigma}_i|= |\g{G}||\g{V}_i-\g{Z}_i\g{G}\g{G}^{-1}\g{G}\g{Z}_i'|= |\g{G}||\g{R}_i+\g{Z}_i\g{G}\g{Z}_i'-\g{Z}_i\g{G}\g{Z}_i'|= |\g{G}| |\g{R}_i|, \] である．よって， \[ \frac{(2\pi)^{q/2}|\g{G}|^{1/2}}{(2\pi)^{(n+q)/2}|\g{\Sigma}_i|^{1/2}}= \frac{1}{(2\pi)^{n/2}|\g{R}_i|^{1/2}}, \] が成り立つ．

次に， \[ \begin{pmatrix} \g{b}_i \\ \g{Y}_i-\g{X}_i\gx{\beta} \end{pmatrix}' \g{\Sigma}_i^{-1} \begin{pmatrix} \g{b}_i \\ \g{Y}_i-\g{X}_i\gx{\beta} \end{pmatrix}- \g{b}_i'\g{G}^{-1}\g{b}_i = (\g{Y}_i-\g{X}_i\gx{\beta}-\g{Z}_i\g{b}_i)' \g{R}_i^{-1} (\g{Y}_i-\g{X}_i\gx{\beta}-\g{Z}_i\g{b}_i), \tag{1} \] を示す．ブロック行列の逆行列の性質から， \[ \g{\Sigma}_i^{-1}= \begin{pmatrix} \g{G}^{-1}+\g{G}^{-1}\g{G}\g{Z}_i'\g{R}_i^{-1}\g{Z}_i\g{G}\g{G}^{-1} & -\g{G}^{-1}\g{G}\g{Z}_i'\g{R}_i^{-1} \\ -\g{R}_i^{-1}\g{Z}_i\g{G}\g{G}^{-1} & \g{R}_i^{-1} \end{pmatrix} = \begin{pmatrix} \g{G}^{-1}+\g{Z}_i'\g{R}_i^{-1}\g{Z}_i & -\g{Z}_i'\g{R}_i^{-1} \\ -\g{R}_i^{-1}\g{Z}_i & \g{R}_i^{-1} \end{pmatrix}, \] である．したがって，(1)の左辺第1項にこれを代入すると， \[ \begin{aligned} &\begin{pmatrix} \g{b}_i' & \g{Y}_i'-\gx{\beta}'\g{X}_i' \end{pmatrix} \begin{pmatrix} \g{G}^{-1}+\g{Z}_i'\g{R}_i^{-1}\g{Z}_i & -\g{Z}_i'\g{R}_i^{-1} \\ -\g{R}_i^{-1}\g{Z}_i & \g{R}_i \end{pmatrix} \begin{pmatrix} \g{b}_i \\ \g{Y}_i-\g{X}_i\gx{\beta} \end{pmatrix} \\&\quad= \g{b}_i'(\g{G}^{-1}+\g{Z}_i'\g{R}_i^{-1}\g{Z}_i)\g{b}_i -(\g{Y}_i'-\gx{\beta}'\g{X}_i')\g{R}_i^{-1}\g{Z}_i\g{b}_i -\g{b}_i'\g{Z}_i'\g{R}_i^{-1}(\g{Y}_i-\g{X}_i\gx{\beta}) +(\g{Y}_i'-\gx{\beta}'\g{X}_i')\g{R}_i^{-1}(\g{Y}_i-\g{X}_i\gx{\beta}) \\&\quad= \g{b}_i'\g{G}^{-1}\g{b}_i+ \g{b}_i'\g{Z}_i'\g{R}_i^{-1}\g{Z}_i\g{b}_i -(\g{Y}_i'-\gx{\beta}'\g{X}_i')\g{R}_i^{-1}\g{Z}_i\g{b}_i -\g{b}_i'\g{Z}_i'\g{R}_i^{-1}(\g{Y}_i-\g{X}_i\gx{\beta}) +(\g{Y}_i'-\gx{\beta}'\g{X}_i')\g{R}_i^{-1}(\g{Y}_i-\g{X}_i\gx{\beta}), \end{aligned} \tag{2} \] となる．したがって， \[ \begin{aligned} &\begin{pmatrix} \g{b}_i \\ \g{Y}_i-\g{X}_i\gx{\beta} \end{pmatrix}' \g{\Sigma}_i^{-1} \begin{pmatrix} \g{b}_i \\ \g{Y}_i-\g{X}_i\gx{\beta} \end{pmatrix}- \g{b}_i'\g{G}^{-1}\g{b}_i \\&\quad= \g{b}_i'\g{Z}_i'\g{R}_i^{-1}\g{Z}_i\g{b}_i -(\g{Y}_i'-\gx{\beta}'\g{X}_i')\g{R}_i^{-1}\g{Z}_i\g{b}_i -\g{b}_i'\g{Z}_i'\g{R}_i^{-1}(\g{Y}_i-\g{X}_i\gx{\beta}) +(\g{Y}_i'-\gx{\beta}'\g{X}_i')\g{R}_i^{-1}(\g{Y}_i-\g{X}_i\gx{\beta}) \\&\quad= (\g{Y}_i-\g{X}_i\gx{\beta}-\g{Z}_i\g{b}_i)'\g{R}_i^{-1}(\g{Y}_i-\g{X}_i\gx{\beta}-\g{Z}_i\g{b}_i), \end{aligned} \] となる．以上より， \[ f(\g{Y}_i; \gx{\beta}, \g{R}_i \mid \g{b}_i)= \frac{1}{(2\pi)^{n/2}|\g{R}_i|^{1/2}} \exp\left\{ -\frac{1}{2} \begin{pmatrix} \g{Y}_i-\g{X}_i\gx{\beta}-\g{Z}_i\g{b}_i \end{pmatrix}' \g{R}_i^{-1} \begin{pmatrix} \g{Y}_i-\g{X}_i\gx{\beta}-\g{Z}_i\g{b}_i \end{pmatrix} \right\}, \] である．

同時尤度関数と周辺尤度関数

線形混合モデルにおける同時尤度関数とは， \[ L_J(\gx{\beta},\g{R},\g{G}; \g{Y},\g{b})= \prod_{i=1}^{n} f(\g{Y}_i; \gx{\beta}, \g{R}_i \mid \g{b}_i) f(\g{b}_i; \g{G}), \] のことである．ただし，$\g{R}$は$\g{R}_i$が含むパラメータ全体を行列で表したものとし，$\g{Y}=(\g{Y}_1', \ldots, \g{Y}_n')'$, $\g{b}=(\g{b}_1', \ldots, \g{b}_n')'$とする．

一般に$\g{b}$は観測不能なため，これを同時尤度関数から積分消去した周辺尤度関数 \[ L_M(\gx{\beta},\g{R},\g{G}; \g{Y})= \prod_{i=1}^{n} \int_{\g{b}_i \in \mathbb{R}^q} f(\g{Y}_i; \gx{\beta}, \g{R}_i \mid \g{b}_i) f(\g{b}_i; \g{G}) \idif \g{b}_i, \] を使った推測手法が発展してきた．この積分の結果は，定義から明らかに$\g{Y}_i \sim N(\g{X}_i\gx{\beta},\g{V}_i)$から作った尤度関数に一致するが，同時尤度関数を積分して求めることができるため，これを示しておく．

行列和$\g{V}_i=\g{R}_i+\g{Z}_i\g{G}\g{Z}_i'$の性質を用いれば， $\g{Y}_i$についての周辺尤度関数は \[ \begin{aligned} f(&\g{Y}_i; \gx{\beta},\g{R},\g{G})= \frac{1}{(2\pi)^{n/2}|\g{V}_i|^{1/2}} \exp\left\{ -\frac{1}{2} (\g{Y}_i-\g{X}_i\gx{\beta})'\g{V}_i^{-1}(\g{Y}_i-\g{X}_i\gx{\beta}) \right\} \\&= \frac{1}{(2\pi)^{n/2}(|\g{R}_i||\g{G}||\g{Q}_i|)^{1/2}} \exp\left\{ -\frac{1}{2} (\g{Y}_i'-\gx{\beta}'\g{X}_i')(\g{R}_i^{-1}-\g{R}_i^{-1}\g{Z}_i\g{Q}_i^{-1}\g{Z}_i'\g{R}_i^{-1})(\g{Y}_i-\g{X}_i\gx{\beta}) \right\}, \end{aligned} \tag{3} \] と変形できる．ただし，$\g{Q}_i=\g{G}^{-1}+\g{Z}_i'\g{R}_i^{-1}\g{Z}_i$の$q \times q$行列である．また，同時密度関数の一部である(2)が \[ \begin{aligned} \g{b}_i'\g{Q}_i\g{b}_i &-(\g{Y}_i'-\gx{\beta}'\g{X}_i')\g{R}_i^{-1}\g{Z}_i\g{b}_i -\g{b}_i'\g{Z}_i'\g{R}_i^{-1}(\g{Y}_i-\g{X}_i\gx{\beta}) +(\g{Y}_i'-\gx{\beta}'\g{X}_i')\g{R}_i^{-1}(\g{Y}_i-\g{X}_i\gx{\beta}) \\&= \g{b}_i'\g{Q}_i\g{b}_i -(\g{Y}_i'-\gx{\beta}'\g{X}_i')\g{R}_i^{-1}\g{Z}_i\g{Q}_i^{-1}\g{Q}_i\g{b}_i -\g{b}_i'\g{Q}_i\g{Q}_i^{-1}\g{Z}_i'\g{R}_i^{-1}(\g{Y}_i-\g{X}_i\gx{\beta}) \\&\quad +(\g{Y}_i'-\gx{\beta}'\g{X}_i')(\g{R}_i^{-1}\g{Z}_i\g{Q}_i^{-1}\g{Q}_i\g{Q}_i^{-1}\g{Z}_i'\g{R}_i^{-1})(\g{Y}_i-\g{X}_i\gx{\beta}) \\&\quad -(\g{Y}_i'-\gx{\beta}'\g{X}_i')(\g{R}_i^{-1}\g{Z}_i\g{Q}_i^{-1}\g{Q}_i\g{Q}_i^{-1}\g{Z}_i'\g{R}_i^{-1})(\g{Y}_i-\g{X}_i\gx{\beta}) +(\g{Y}_i'-\gx{\beta}'\g{X}_i')\g{R}_i^{-1}(\g{Y}_i-\g{X}_i\gx{\beta}) \\&= (\g{b}_i'-\g{Y}_i'\g{R}_i^{-1}\g{Z}_i\g{Q}_i^{-1})\g{Q}_i(\g{b}_i-\g{Q}_i^{-1}\g{Z}_i'\g{R}_i^{-1}\g{Y}_i)+(\g{Y}_i'-\gx{\beta}'\g{X}_i')(\g{R}_i^{-1}-\g{R}_i^{-1}\g{Z}_i\g{Q}_i^{-1}\g{Z}_i'\g{R}_i^{-1})(\g{Y}_i-\g{X}_i\gx{\beta}), \end{aligned} \] と変形できることを用いれば，$\g{Y}_i$についての同時尤度関数の積分が \[ \begin{aligned} \int_{\g{b}_i \in \mathbb{R}^q} &f(\g{Y}_i; \gx{\beta}, \g{R}_i \mid \g{b}_i) f(\g{b}_i; \g{G}) \idif \g{b}_i = \int_{\g{b}_i \in \mathbb{R}^q} f(\g{Y}_i,\g{b}_i; \gx{\beta}, \g{R}_i, \g{G}) \idif \g{b}_i \\&= \frac{1}{(2\pi)^{n/2}|\g{R}_i|^{1/2}|\g{G}|^{1/2}} \exp\left\{ -\frac{1}{2} (\g{Y}_i'-\gx{\beta}'\g{X}_i')(\g{R}_i^{-1}-\g{R}_i^{-1}\g{Z}_i\g{Q}_i^{-1}\g{Z}_i'\g{R}_i^{-1})(\g{Y}_i-\g{X}_i\gx{\beta}) \right\} \\&\qquad \times \frac{1}{|\g{Q}|^{1/2}}\int_{\g{b}_i \in \mathbb{R}^q} \frac{1}{(2\pi)^{q/2}|\g{Q}^{-1}|^{1/2}} \exp\left\{ -\frac{1}{2} (\g{b}_i'-\g{Y}_i'\g{R}_i^{-1}\g{Z}_i\g{Q}_i^{-1})(\g{Q}_i^{-1})^{-1}(\g{b}_i-\g{Q}_i^{-1}\g{Z}_i'\g{R}_i^{-1}\g{Y}_i) \right\} \idif \g{b}_i \\&= \frac{1}{(2\pi)^{n/2}(|\g{R}_i||\g{G}||\g{Q}_i|)^{1/2}} \exp\left\{ -\frac{1}{2} (\g{Y}_i'-\gx{\beta}'\g{X}_i')(\g{R}_i^{-1}-\g{R}_i^{-1}\g{Z}_i\g{Q}_i^{-1}\g{Z}_i'\g{R}_i^{-1})(\g{Y}_i-\g{X}_i\gx{\beta}) \right\} \times 1, \end{aligned} \] となり，(3)に一致することを確認できる．ただし，上式の3行目の積分値は，平均ベクトル$\g{Q}_i^{-1}\g{Z}_i'\g{R}_i^{-1}\g{Y}_i$，分散共分散行列$\g{Q}_i^{-1}$の正規分布の密度関数の全領域を積分したものに等しいため$1$となる(これは単に積分値を求めているだけなので，実際の確率変数とその平均ベクトル・分散共分散行列とは無関係に成立する)．

同時尤度を少し変形したLee & Nelder (1996)のh-likelihood (hierarchical log-likelihood)に基づく方法もある．これは変量効果を積分消去しないで反復計算や近似計算などで対処するアプローチである(線形モデルのときはあまり意味はないが)．

補遺

行列$\g{A}$を$n\times n$，$\g{B}$を$n\times m$，$\g{C}$を$m\times n$，$\g{D}$を$m\times m$で，$\g{A}$が正則とする．

ブロック行列 \[ \g{F}=\begin{pmatrix} \g{A} & \g{B} \\ \g{C} & \g{D} \\ \end{pmatrix} \] の行列式は， \[ |\g{F}|=|\g{A}| |\g{D} - \g{C}\g{A}^{-1}\g{B}|, \] である．また，$\g{S}=\g{D}-\g{C}\g{A}^{-1}\g{B}$とおくと，逆行列は \[ \g{F}^{-1}=\begin{pmatrix} \g{A}^{-1} + \g{A}^{-1}\g{B}\g{S}^{-1}\g{C}\g{A}^{-1} & -\g{A}^{-1}\g{B}\g{S}^{-1} \\ -\g{S}^{-1}\g{C}\g{A}^{-1} & \g{S}^{-1} \\ \end{pmatrix}, \] である．

行列$\g{D}$も正則のとき，行列和 \[ \g{A}+\g{B}\g{D}\g{C} \] について， \[ |\g{A}+\g{B}\g{D}\g{C}|= |\g{A}| |\g{D}| |\g{D}^{-1}+\g{C}\g{A}^{-1}\g{B}|, \] および， \[ (\g{A}+\g{B}\g{D}\g{C})^{-1}= \g{A}^{-1}- \g{A}^{-1}\g{B}(\g{D}^{-1}+\g{C}\g{A}^{-1}\g{B})^{-1}\g{C}\g{A}^{-1}, \] が成り立つ．

参考文献

Lee Y, Nelder JA. Hierarchial generalized linear models (with Discussion). Journal of the Royal Statistical Society: Series B 1996; 58(4): 619–656. DOI: 10.1111/j.2517-6161.1996.tb02105.x

履歴

2019/09/19 公開