GLMM の推定アルゴリズムに Penalized quasi-likelihood (PQL) という方法があるのですが、整理のためメモを書くことにしました。
PQL
混合効果モデルのパラメータ推定において、周辺尤度の積分を必要とする方法があることはよく知られていると思います。
GLMM ではリンク関数が非線形なものを扱うが故に、この積分が非常に大変でした (現代の計算機では苦労せず解けるものも多い)。
一昔前はこの積分が大変だったため、非線形な関数を線形な関数に近似してあげて、LMM の範疇で扱えるようにした方法が PQL です。
GLMM として、確率変数 $\mathbf{Y}_i$ の期待値が、変量効果 $\mathbf{b}_i$ を与えた下で、
\[
\mathrm{E}[\mathbf{Y}_i|\mathbf{b}_i]=g^{-1}(\mathbf{X}_i\boldsymbol{\beta}+\mathbf{Z}_i\textbf{b}_i)=g^{-1}(\boldsymbol{\xi}_i)=\boldsymbol{\mu}_i
\] とモデル化された一般的なものを考えます ($\mathbf{X}_i$, $\mathbf{Z}_i$ はデザイン行列、$\boldsymbol{\beta}$ は固定効果のパラメータベクトル、$\boldsymbol{\xi}_i = \mathbf{X}_i\boldsymbol{\beta}+\mathbf{Z}_i\textbf{b}_i$、$g^{-1}$ はリンク関数の逆関数)。
PQL では $\boldsymbol{\xi}_i$ の関数である $\boldsymbol{\mu}_i$ に対して $\boldsymbol{\beta}=\hat{\boldsymbol{\beta}}$, $\mathbf{b}_i=\hat{\mathbf{b}_i}$ の周りの一次の Taylor 展開、
\[
g^{-1}(\mathbf{X}_i\boldsymbol{\beta}+\mathbf{Z}_i\textbf{b}_i) \approx g^{-1}(\mathbf{X}_i\hat{\boldsymbol{\beta}}+\mathbf{Z}_i\hat{\textbf{b}_i}) + \hat{\Delta}_i\mathbf{X}_i(\boldsymbol{\beta}-\hat{\boldsymbol{\beta}}) + \hat{\Delta}_i\mathbf{Z}_i(\textbf{b}_i-\hat{\textbf{b}_i})
\] を考えます、ただし、
\[
\hat{\Delta}_i= \mathrm{diag} \left( \frac{\partial g^{-1}(\xi_{ij})}{\partial \xi_{ij}} \Big|_{\boldsymbol{\xi}_i=\hat{\boldsymbol{\xi}}_i}\right)_{1 \leq j \leq n_i}
\] とします ($n_i \times n_i$ 次元行列)。
これを整理すると、
\[
\hat{\Delta}_i^{-1}(\boldsymbol{\mu}-g^{-1}(\mathbf{X}_i\hat{\boldsymbol{\beta}}+\mathbf{Z}_i\hat{\textbf{b}_i}))+ \mathbf{X}_i\hat{\boldsymbol{\beta}}+ \mathbf{Z}_i\hat{\textbf{b}_i} \approx \mathbf{X}_i\boldsymbol{\beta}+ \mathbf{Z}_i\textbf{b}_i
\] が得られます。
左辺をよく見ると、これは、$\hat{\Delta}_i^{-1}(\mathbf{Y}_i-g^{-1}(\mathbf{X}_i\hat{\boldsymbol{\beta}}+\mathbf{Z}_i\hat{\textbf{b}_i}))+ \mathbf{X}_i\hat{\boldsymbol{\beta}}+ \mathbf{Z}_i\hat{\textbf{b}_i} \equiv \mathbf{P}_i$ の期待値になっていることが分かります。
確率変数 $\mathbf{Y}_i$ 以外は全て定数とみなして新しい作業確率変数として $\mathbf{P}_i$ を定義します。
また、同じく $\mathbf{Y}_i$ 以外が全て定数とすれば、$\mathrm{Var}(\mathbf{P}_i\mid\mathbf{b}_i)= \hat{\Delta}_i^{-1} \mathrm{Var}(\mathbf{Y}_i\mid\mathbf{b}_i) \hat{\Delta}_i^{-1}$ が簡単に分かります。
そして、右辺からは非線形なリンク関数がうまく消去できていて、線形化することができました。
以上の近似を使って、LMM、
\[
\mathbf{P}_i=\mathbf{X}_i\boldsymbol{\beta}+\mathbf{Z}_i\mathbf{b}_i+\boldsymbol{\epsilon}^*_i,~
\mathrm{Var}(\boldsymbol{\epsilon}^*_i)= \mathrm{Var}(\mathbf{P}_i \mid \textbf{b}_i)
\] を当てはめたものが、PQL と呼ばれる方法です。
また、類似の Marginal quasi-likelihood (MQL) という方法では、変量効果の方を母集団平均で置き換えるため、$\mathbf{b}_i=\mathbf{0}$ を用います。
ただし、MQL は変量効果が非常に小さい場合でないとうまくいかないことが知られています。
また、アウトカムが二項分布に従う場合は、PQL と MQL による近似は特にうまくいかないケースがあることが知られています。
感想
とまあ、一次の Taylor 展開を使ったシンプルな方法だという事が分かります。
計算は早くなりますが、デメリットもあります。
作業確率変数 $\mathbf{P}_i$ に置き換えるため、$\mathbf{y}_i$ に基づく尤度が定義できない事はデメリットの一つでしょう。
どこかに quasi-likelihood を使う方法だから良くないと書いてあった気がしますが、そうではないんです。
本質的な問題点は作業確率変数を使った近似の方ですよね。
よって、利用可能であれば、適応型ガウス–エルミート求積法などの周辺尤度を直接計算するような方法が望ましいでしょう。
ちなみに SAS 9.3 の GLIMMIX Procedure は METHOD = RSPL (PQL で REML 推定) がデフォルトです。
