【統計検定１級対策】一般化線形モデル（GLM）の具体例

前回記事に続いてGLMについて備忘録まとめていきます。

前回の記事はこちら。

【統計検定１級対策】一般化線形モデル（GLM）についてのまとめ

指数型分布族についてGLMの一般式を前回書きましたが、本当にその通りになっているのか、みていこうと思います。

GLMの具体例１：正規分布を用いた線形回帰モデル

モデルの定義

観測値 \(y_i\) は平均 \(\mu_i\)、分散 \(\sigma^2\) の正規分布に従う：

\[ y_i \sim \mathcal{N}(\mu_i, \sigma^2) \]

正規分布の確率密度関数は以下のように書ける：

\[ f(y_i; \mu_i, \sigma^2) = \frac{1}{\sqrt{2\pi\sigma^2}} \exp\{ -\frac{(y_i – \mu_i)^2}{2\sigma^2} \} \]

リンク関数は恒等関数（\(g(\mu_i) = \mu_i\)）であり、線形予測子は：

\[ \mu_i = \eta_i = x_i^\top \beta \]

ここで

\(x_i = (x_{1i}, x_{2i})^\top\)
\(\beta = (\beta_1, \beta_2)^\top\)

としておきます。

対数尤度関数の導出

上記を代入して、対数尤度関数 \(\log L\) を求めると：

\[ \log L = -\frac{n}{2} \log(2\pi\sigma^2) – \sum_{i=1}^n \frac{(y_i – x_i^\top \beta)^2}{2\sigma^2} \]

尤度の微分と最尤推定

パラメータ \(\beta\) に関して微分：

連鎖律を使って分解：

\[ \frac{\partial \log L}{\partial \beta} = \sum_{i=1}^n \frac{(y_i – x_i^\top \beta)}{\sigma^2} \cdot x_i \]

ここでいったん前回みたGLMのスコア関数の一般式をみてみますと

\[ \sum_{i=1}^n (y_i – \mu_i) \cdot \text{Var}^{-1}(y_i) \cdot \frac{\partial \mu_i}{\partial \eta_i} \cdot \frac{\partial \eta_i}{\partial \beta} = 0 \]

ここで\(\frac{\partial \eta_i}{\partial \beta}=x_i\)であり、\(\frac{\partial \mu_i}{\partial \eta_i}=1\)であることからキチンと合致していることが分かります。

リンク関数と線形予測子が単純な構造だからこそ、このような式になるのがわかります。

あとはスコア関数＝０としておけば、βの推定ができます。

\[ \sum_{i=1}^n (y_i – x_i^\top \beta) x_i = 0 \]

正規方程式の導出

上記を行列で表してみると正規方程式が導出できます。まず以下の様に行列を定義しましょう。

\(X\) は説明変数を並べた計画行列（デザイン行列）を

\[
X = \begin{pmatrix}
x_{1,1} & x_{2,1} \\
x_{1,2} & x_{2,2} \\
\vdots & \vdots
\end{pmatrix}
\]

\(y\) は観測データベクトル：

\[
y = \begin{pmatrix}
y_1 \\
y_2 \\
\vdots
\end{pmatrix}
\]

そうすると

\[ X^\top (y – X\beta) = 0 \]

整理して

\[ X^\top y = X^\top X \beta \]

両辺に \((X^\top X)^{-1}\) をかけることで、\(\beta\) の推定値を得ることができます。

\[ \hat{\beta} = (X^\top X)^{-1} X^\top y \]

この手順は通常の線形回帰で、GLMの中では特例として解釈できます
リンク関数が恒等関数（\(g(\mu) = \mu\)）で、分布が正規分布である場合に可能です

というのも前回のGLMのスコア関数の一般式を再度みてみますと

\[ \sum_{i=1}^n (y_i – \mu_i) \cdot \text{Var}^{-1}(y_i) \cdot \frac{\partial \mu_i}{\partial \eta_i} \cdot \frac{\partial \eta_i}{\partial \beta} = 0 \]

これが恒等リンクなので、\(\mu = \eta\) であり、微分した部分は１となります。そのため解析的に解けるようになっています。

例２：Poisson分布

確率密度関数

\[ f(y \mid \mu) = \frac{e^{-\mu} \mu^y}{y!} \]

対数尤度関数（複数観測値の場合）

\[ \log L(\mu_i) = \log \prod_{i=1}^n \frac{e^{-\mu_i} \mu_i^{y_i}}{y_i!} = \sum_{i=1}^n \left[ -\mu_i + y_i \log \mu_i – \log \Gamma(y_i + 1) \right] \]

ここで \(\mu_i > 0\) なので、リンク関数として

\[ g(\mu_i) = \log \mu_i = \eta_i \]

とするのが自然と言えます。

線形予測子との関係は

\[ \eta_i = \mathbf{x}_i^\top \beta \Rightarrow \mu_i = \exp(\mathbf{x}_i^\top \beta) \]

先ほどのGLMの一般式を思い出すと、これは微分するのが辛そうなのがわかります。

これを用いた対数尤度は：

\[ \log L(\mu_i) = \sum_{i=1}^n \left[ -\exp(\mathbf{x}_i^\top \beta) y_i \mathbf{x}_i^\top \beta – \log \Gamma(y_i + 1) \right] \]

（最後の項は \(\beta\) に依存しないので微分には影響しない）

スコア関数の導出

対数尤度を \(\beta_j\) で微分すると：

\[ \frac{\partial \log L(\mu_i)}{\partial \beta_j} = \frac{\partial \log L}{\partial \mu_i} \cdot \frac{\partial \mu_i}{\partial \eta_i} \cdot \frac{\partial \eta_i}{\partial \beta_j} \]

展開すると：

\[ = \sum_{i=1}^n \left( -1 + \frac{y_i}{\mu_i} \right) x_{ij} \cdot \frac{\partial \mu_i}{\partial \eta_i} = \sum_{i=1}^n \frac{(y_i – \mu_i)}{\mu_i} x_{ij} \cdot \frac{\partial \mu_i}{\partial \eta_i} \]

ポアソン分布では平均と分散は一致しており、いずれも \(\mu\) ですので、GLM のスコア関数の一般形：

\[ \sum_i (y_i – \mu_i) \cdot \text{Var}^{-1}(y_i) \cdot \frac{\partial \mu_i}{\partial \eta_i} \cdot \frac{\partial \eta_i}{\partial \beta} \]

と一致していることが分かります。