GLMをもう少し理解したい③ - 統計コンサルの議事メモ

前回の記事において、GLMでは以下の方程式を用いてパラメータベクトルを推定するという話をしました：

$\mathbf{b}^ {(m)} = \mathbf{b}^ {(m-1)} + [\mathfrak{J}^ {(m-1)}]^ {-1} \mathbf{U}^ {(m-1)}$

ushi-goroshi.hatenablog.com

今回はその続きです。

※ 1/25 記事を修正しました

最尤推定

上の式には情報行列 $\mathfrak{J}$ の逆行列が入っているので、前から情報行列を乗じます：

$\mathfrak{J}^ {(m-1)}\mathbf{b}^ {(m)} = \mathfrak{J}^ {(m-1)}\mathbf{b}^ {(m-1)} + \mathbf{U}^ {(m-1)} \tag{1}$

ところで、（ここもまた天下りですが）スコア $U$ は以下の式で表されますが*1、

$\mathbf{U_{j}} = \Sigma_{i=1}^ {N} [ \frac{(y_{i} - \mu_{i})}{var(Y_{i})} x_{ij} ( \frac{\partial{\mu_{i}}}{\partial{\eta_{i}}} ) ] \tag{2}$

$\mathfrak{J}$ はスコア $U$ の分散共分散行列であり、上記の式から

$\begin{align} \mathfrak{J}_{jk} &= E\left\{\Sigma_{i=1}^ {N}[ \frac{(Y_{i} - \mu_{i})}{var(Y_{i})} x_{ij} (\frac{\partial{\mu_{i}}}{\partial{\eta_{i}}}) ] \Sigma_{l=1}^ {N}[ \frac{(Y_{l} - \mu_{l})}{var(Y_{l})} x_{lk} (\frac{\partial{\mu_{l}}}{\partial{\eta_{l}}}) ]\right\} \\ &= \Sigma_{i=1}^ {N} \frac{E[ (Y_{i} - \mu_{i})^ {2} ] x_{ij}x_{ik}} {[var(Y_{i})]^ 2} (\frac{\partial{\mu_{i}}}{\partial{\eta_{i}}})^ {2} \end{align}$

となり、さらに $E[ (Y_{i} - \mu_{i})^ {2}] = var(Y_{i})$ から以下のようになります：

$\mathfrak{J}_{jk} = \Sigma_{i=1}^ {N} \frac{x_{ij}x_{ik}} {var(Y_{i})} (\frac{\partial{\mu_{i}}}{\partial{\eta_{i}}})^ {2}$

これを行列で表記すると：

$\mathfrak{J} = \mathbf{X^ {T} WX}$

となります。ただし $\mathbf{W}$ は

$w_{ii} = \frac{1}{var(Y_{i})} (\frac{\partial{\mu_{i}}}{\partial{\eta_{i}}})^ {2}$

となる対角行列です。よって冒頭の式(1)の左辺は $\mathbf{X^ {T} WXb}^ {(m)}$ と書けます。

次に右辺ですが、式(2)から以下のように書けます：

$\Sigma_{k=1}^ {p} \Sigma_{i=1}^ {N} \frac{x_{ij}x_{ik}}{var(Y_{i})} (\frac{\partial{\mu_{i}}}{\partial{\eta_{i}}} )^ {2}b_{k}^ {(m-1)} + \Sigma_{i=1}^ {N} \frac{(y_{i}-\mu_{i})x_{ij}}{var(Y_{i})} (\frac{\partial{\mu_{i}}}{\partial{\eta_{i}}} )$

このとき1項目と2項目を $\mathbf{X^ {T}W}$ でくくり、残りを $\mathbf{z}$ とすると

$\mathbf{X^ {T}WXb}^ {(m)} = \mathbf{X^ {T}Wz}$

が得られます。これは線形モデルに対する重み付き最小二乗法を適用して得られる正規方程式と同じ形となりますが、 $\mathbf{z}$ や $\mathbf{W}$ は $\mathbf{b}$ に依存するため、反復的に解く必要があります。
なお $\mathbf{X^ {T}W}$ でくくる際、 $(\frac{\partial{\mu_{i}}}{\partial{\eta_{i}}})^ {2}$ を使っているため二項目にはこの逆数である $(\frac{\partial{\eta_{i}}}{\partial{\mu_{i}}})$ が残り、 $\mathbf{z}$ にはそれが渡ります。すなわち：

$\mathbf{z_{i}} = \Sigma_{k=1}^ {p} x_{ik}b_{k}^ {(m-1)} + (y_{i} - \mu_{i})( \frac{\partial{\eta_{i}}}{\partial{\mu_{i}}} )$

です。

※ 以下のセクションで $\eta$ と $\mu$ の関係を取り違えて説明していたので修正しました

さて、この $\mathbf{z}$ に含まれる $(\frac{\partial{\eta_{i}}}{\partial{\mu_{i}}})$ はどのような形になるでしょうか？GLMにおいて $\eta$ は~~線形予測子 $\mathbf{X^ {T}b}$ をリンク関数で変換したもの~~線形予測子で説明されるもので、 $\mu$ はそれを逆リンク関数で変換したものでした。つまり、

$\eta = g(\mu) \\ \mu = g^ {-1}(\eta)$

という関係であることを思い出すと、 $(\frac{\partial{\eta_{i}}}{\partial{\mu_{i}}})$ は分析者が事前に設定したリンク関数に依存することになります。例えばリンク関数としてlogを指定していれば、

$(\frac{\partial{\eta_{i}}}{\partial{\mu_{i}}}) = (\frac{\partial{log(\mu_{i})}}{\partial{\mu_{i}}}) = \frac{1}{\mu_{i}}$

となります。identityなら1になるでしょう。

同様に $w_{ii}$ に含まれる $(\frac{\partial{\mu_{i}}}{\partial{\eta_{i}}})$ は、 $\mu$ が $g^ {-1}(\eta)$ であることから逆リンク関数に依存します。例えば逆リンク関数がexpなら

$(\frac{\partial{\mu_{i}}}{\partial{\eta_{i}}}) = (\frac{\partial{exp(\eta_{i})}}{\partial{\eta_{i}}}) = exp(\eta_{i}) = \mu_{i}$

となるでしょう。

以上から、 $\mathbf{X}$ 、 $\mathbf{W}$ 、 $\mathbf{z}$ が分かれば $\mathbf{b}$ を更新することができ、またそれぞれについても具体的に計算できそうです。GLMではこの反復重み付き最小二乗法（IRLS, Iteratively Reweighted Least Squares）によって最尤推定量を求めます。

それでは次回は、具体的な数値を使って計算してみましょう。

*1:Dobsonの式(4.18)より。この式の導出がちょっと煩雑だったので割愛。