ある保険会社が、契約者の年齢・運転歴・車種から、年間の事故件数を予測したいとします。事故件数は0, 1, 2, … という非負の整数値であり、通常の重回帰分析を適用すると、予測値が負になったり、小数点のある値になったりして、カウントデータとしての性質が反映されません。
あるいは、臨床試験で薬の投与量から副作用の有無（0か1か）を予測する問題を考えてください。目的変数が二値データの場合、正規分布の仮定は明らかに不適切です。目的変数は0と1しか取らず、対称でもなく、分散も平均に依存します。
こうした「正規分布に従わない目的変数」に対して、回帰分析を自然に拡張する統一的な枠組みが一般化線形モデル（GLM: Generalized Linear Model）です。ロジスティック回帰、ポアソン回帰、ガンマ回帰などは、すべてGLMの特殊ケースとして理解できます。GLMの理論は1972年にNelder & Wedderburnが体系化し、統計学に革命的な影響を与えました。
GLMを理解すると、以下のような応用が開けます。
- 保険数理: 事故件数のモデリング（ポアソン回帰）、保険金額のモデリング（ガンマ回帰）
- 生態学: 種の個体数（ポアソン回帰）や出現確率（ロジスティック回帰）のモデリング
- 医学: 二値の応答変数（発症の有無）、生存時間データ、用量-反応関係
- 品質管理: 不良品の発生率、欠陥数のモデリング
- 社会科学: 投票行動（二値）、犯罪件数（カウント）の分析
本記事では、GLMの3つの構成要素（確率分布、リンク関数、線形予測子）を定式化し、指数型分布族の性質を利用してパラメータの最尤推定を導出します。IRLSアルゴリズムの詳細な導出を行い、Pythonでポアソン回帰を実装して理論を確認します。
本記事の内容
- GLMの3つの構成要素
- 指数型分布族の定義と性質（期待値・分散の導出）
- リンク関数と正準リンク
- 最尤推定とスコア方程式の導出
- IRLS（反復重み付き最小二乗法）の導出
- 代表的なGLM（正規、ポアソン、ロジスティック、ガンマ）
- 偏差とモデルの適合度
- Pythonでの実装と可視化
前提知識
この記事を読む前に、以下の記事を読んでおくと理解が深まります。
GLMの3つの構成要素
なぜ統一的な枠組みが必要か
重回帰分析、ロジスティック回帰、ポアソン回帰、ガンマ回帰。これらは一見すると別々の手法に見えますが、実はすべて同じ構造を共有しています。GLMが登場する以前は、各手法が個別に開発・教育されていました。GLMの枠組みは、これらを1つの統一理論としてまとめ、新しい応用への拡張を容易にしました。
GLMは、以下の3つの要素で特徴づけられます。
1. 確率分布（Random Component）: 目的変数 $y_i$ が従う確率分布。正規分布、ポアソン分布、二項分布、ガンマ分布などの指数型分布族から選びます。データの性質（連続・離散・非負・0-1など）に応じて適切な分布を選択します。
2. 線形予測子（Systematic Component）: 説明変数の線形結合 $\eta_i = \bm{x}_i^\top\bm{\beta}$。通常の重回帰と同じ形です。GLMの「線形」は、パラメータ $\bm{\beta}$ に関して線形であることを指し、$\bm{x}$ の非線形変換（多項式や交互作用項）を含めることもできます。
3. リンク関数（Link Function）: 目的変数の期待値 $\mu_i = E[y_i]$ と線形予測子 $\eta_i$ を結びつける単調な微分可能関数 $g$。
$$ \begin{equation} g(\mu_i) = \eta_i = \bm{x}_i^\top\bm{\beta} \end{equation} $$
通常の線形回帰は $g(\mu) = \mu$（恒等リンク）、正規分布を持つGLMの特殊ケースです。ロジスティック回帰は $g(\mu) = \ln(\mu/(1-\mu))$（ロジットリンク）、二項分布を持つGLMです。
リンク関数の直感的な理解
なぜリンク関数が必要なのでしょうか。この疑問に対して、具体例から直感的に理解しましょう。
線形予測子 $\eta = \bm{x}^\top\bm{\beta}$ は $-\infty$ から $+\infty$ まで任意の値を取ります。しかし、期待値 $\mu$ には制約がある場合が多いです。
- ポアソン分布: 平均 $\mu > 0$（事故件数の期待値は正でなければならない）
- 二項分布: 確率 $\mu \in (0, 1)$（発症確率は0から1の間）
- ガンマ分布: 平均 $\mu > 0$（保険金額の期待値は正）
リンク関数は、$\mu$ の制約のある空間から $\eta$ の制約のない空間への「変換」の役割を果たします。逆リンク関数 $g^{-1}$ が $\eta$ を $\mu$ に変換します。たとえばポアソン回帰では $g(\mu) = \ln\mu$ なので $\mu = e^{\eta} > 0$ となり、平均が正であることが自動的に保証されます。
この3つの構成要素を変えるだけで、さまざまなタイプのデータに対応できるのがGLMの強みです。次に、GLMで使う分布族の数学的性質を見ていきましょう。
指数型分布族
定義
GLMで使う確率分布は指数型分布族（exponential family）に属する分布です。指数型分布族とは、密度関数（または質量関数）が次の標準形で書ける分布の総称です。
$$ \begin{equation} f(y | \theta, \phi) = \exp\left(\frac{y\theta – b(\theta)}{a(\phi)} + c(y, \phi)\right) \end{equation} $$
ここで各項の意味は次の通りです。
- $\theta$: 自然パラメータ（natural parameter） — 分布の位置を決めるパラメータ
- $\phi$: 分散パラメータ（dispersion parameter） — 分布の広がりを制御するパラメータ
- $a(\phi)$: 分散パラメータの関数。多くの場合 $a(\phi) = \phi/w_i$（$w_i$ は既知の重み）
- $b(\theta)$: 累積母関数（cumulant function）と呼ばれ、分布の正規化を保証する関数
- $c(y, \phi)$: $\theta$ に依存しない正規化定数
一見すると複雑な形に見えますが、この形に統一することで、期待値・分散・最尤推定の理論を分布ごとに個別に導出する必要がなくなります。
期待値と分散の導出
指数型分布族の最も便利な性質は、$b(\theta)$ の微分から期待値と分散が直接得られることです。これを導出しましょう。
密度関数の積分は1なので、次が成り立ちます。
$$ \int f(y|\theta, \phi)\,dy = 1 $$
両辺を $\theta$ で微分します（積分と微分の交換を仮定）。
$$ \int \frac{\partial}{\partial\theta}f(y|\theta,\phi)\,dy = 0 $$
$f$ の中身を微分すると次のようになります。
$$ \frac{\partial}{\partial\theta}\ln f = \frac{y – b'(\theta)}{a(\phi)} $$
したがって $\frac{\partial f}{\partial\theta} = f \cdot \frac{y – b'(\theta)}{a(\phi)}$ です。これを積分に代入すると次のようになります。
$$ \int \frac{y – b'(\theta)}{a(\phi)} f(y)\,dy = 0 $$
$a(\phi) \neq 0$ なので次が得られます。
$$ E[y] – b'(\theta) = 0 $$
以上から、期待値は次のように表されます。
$$ \begin{equation} E[y] = \mu = b'(\theta) \end{equation} $$
同様に $\theta$ で2回微分すると、分散が得られます。
$$ \begin{equation} \text{Var}(y) = a(\phi) \cdot b”(\theta) = a(\phi) \cdot V(\mu) \end{equation} $$
ここで $V(\mu) = b”(\theta)$ は分散関数と呼ばれ、平均 $\mu$ の関数として分散を表します。分散が平均に依存する（等分散でない）ことが、正規分布の場合と大きく異なる点です。
代表的な分布の指数型表現
主要な分布を指数型分布族の形に書き下し、それぞれの構成要素を確認しましょう。
正規分布 $N(\mu, \sigma^2)$:
$$ f(y) = \frac{1}{\sqrt{2\pi\sigma^2}}\exp\left(-\frac{(y-\mu)^2}{2\sigma^2}\right) = \exp\left(\frac{y\mu – \mu^2/2}{\sigma^2} – \frac{y^2}{2\sigma^2} – \frac{1}{2}\ln(2\pi\sigma^2)\right) $$
したがって $\theta = \mu$, $b(\theta) = \theta^2/2$, $a(\phi) = \sigma^2$ です。$V(\mu) = b”(\theta) = 1$（分散が平均に依存しない）。
ポアソン分布 $\text{Poi}(\mu)$:
$$ f(y) = \frac{\mu^y e^{-\mu}}{y!} = \exp(y\ln\mu – \mu – \ln y!) $$
$\theta = \ln\mu$, $b(\theta) = e^\theta = \mu$, $a(\phi) = 1$ です。$V(\mu) = b”(\theta) = e^\theta = \mu$（分散が平均に等しい）。
二項分布 $\text{Bin}(n, p)$ を $y/n$ で書くと:
$\theta = \ln(p/(1-p))$（ロジット）, $b(\theta) = \ln(1+e^\theta)$, $V(\mu) = \mu(1-\mu)$ です。
ガンマ分布 $\text{Gamma}(\alpha, \beta)$:
$\theta = -1/\mu$, $b(\theta) = -\ln(-\theta)$, $V(\mu) = \mu^2$ です。
これらをまとめると次の表のようになります。
| 分布 | $\theta$ | $b(\theta)$ | $V(\mu)$ | 正準リンク |
|---|---|---|---|---|
| 正規 | $\mu$ | $\theta^2/2$ | $1$ | 恒等 $g(\mu)=\mu$ |
| ポアソン | $\ln\mu$ | $e^\theta$ | $\mu$ | 対数 $g(\mu)=\ln\mu$ |
| 二項 | $\ln(\mu/(1-\mu))$ | $\ln(1+e^\theta)$ | $\mu(1-\mu)$ | ロジット $g(\mu)=\ln(\mu/(1-\mu))$ |
| ガンマ | $-1/\mu$ | $-\ln(-\theta)$ | $\mu^2$ | 逆数 $g(\mu)=1/\mu$ |
正準リンク関数
正準リンク関数は $g(\mu) = \theta$ と定義される、自然パラメータと線形予測子を直接等しくするリンク関数です。$\theta = g(\mu) = \eta = \bm{x}^\top\bm{\beta}$ となるため、自然パラメータが説明変数の線形関数になります。
正準リンクを使うと、以下の利点があります。
- スコア方程式が簡潔な形: $\bm{X}^\top(\bm{y} – \bm{\mu}) = \bm{0}$ という美しい形になります
- 十分統計量の存在: $\bm{X}^\top\bm{y}$ が $\bm{\beta}$ の十分統計量になります
- 凹性の保証: 対数尤度が $\bm{\beta}$ に関して凹であることが保証されます
正準リンクを使う必要はありませんが、理論的な簡潔さと計算の安定性から、実用上よく使われます。
最尤推定の理論に進みましょう。
最尤推定
対数尤度関数
$n$ 個の独立な観測 $(y_i, \bm{x}_i)$ に対する対数尤度関数は次のようになります。
$$ \begin{equation} \ell(\bm{\beta}) = \sum_{i=1}^{n}\left(\frac{y_i\theta_i – b(\theta_i)}{a(\phi)} + c(y_i, \phi)\right) \end{equation} $$
$\theta_i$ は $\bm{x}_i^\top\bm{\beta}$ とリンク関数を通じて次のように決まります。
$$ \eta_i = \bm{x}_i^\top\bm{\beta} \to \mu_i = g^{-1}(\eta_i) \to \theta_i = (b’)^{-1}(\mu_i) $$
スコア方程式の導出
対数尤度を $\beta_j$ で微分します。連鎖律を使って次のように書きます。
$$ \frac{\partial\ell}{\partial\beta_j} = \sum_{i=1}^{n}\frac{\partial\ell_i}{\partial\theta_i}\cdot\frac{\partial\theta_i}{\partial\mu_i}\cdot\frac{\partial\mu_i}{\partial\eta_i}\cdot\frac{\partial\eta_i}{\partial\beta_j} $$
各項を計算します。
$\frac{\partial\ell_i}{\partial\theta_i} = \frac{y_i – b'(\theta_i)}{a(\phi)} = \frac{y_i – \mu_i}{a(\phi)}$（$b'(\theta) = \mu$ を使用）
$\frac{\partial\theta_i}{\partial\mu_i} = \frac{1}{b”(\theta_i)} = \frac{1}{V(\mu_i)}$（$b'(\theta) = \mu$ の逆関数の微分）
$\frac{\partial\mu_i}{\partial\eta_i} = \frac{1}{g'(\mu_i)}$（リンク関数の逆関数の微分）
$\frac{\partial\eta_i}{\partial\beta_j} = x_{ij}$
これらを合わせると、スコア方程式が得られます。
$$ \begin{equation} \frac{\partial \ell}{\partial \beta_j} = \sum_{i=1}^{n}\frac{(y_i – \mu_i)x_{ij}}{a(\phi)V(\mu_i)g'(\mu_i)} = 0 \end{equation} $$
ベクトル形式で書くと次のようになります。
$$ \frac{\partial\ell}{\partial\bm{\beta}} = \frac{1}{a(\phi)}\bm{X}^\top\bm{W}_s\bm{\Delta}(\bm{y} – \bm{\mu}) = \bm{0} $$
ここで $\bm{W}_s = \text{diag}(1/V(\mu_i))$、$\bm{\Delta} = \text{diag}(1/g'(\mu_i))$ です。
正準リンクの場合: $g'(\mu) = 1/V(\mu)$ となるため、$V(\mu_i) \cdot g'(\mu_i) = 1$ が成り立ち、スコア方程式は次のように簡略化されます。
$$ \sum_{i=1}^{n}(y_i – \mu_i)\bm{x}_i = \bm{0} \quad \Longleftrightarrow \quad \bm{X}^\top(\bm{y} – \bm{\mu}) = \bm{0} $$
これはロジスティック回帰で見たスコア方程式と同じ形です。正規回帰の正規方程式 $\bm{X}^\top(\bm{y} – \bm{X}\bm{\beta}) = \bm{0}$ も同じ形であり、GLMの統一的な枠組みの美しさがここに現れています。
フィッシャー情報行列
ヘッセ行列の期待値（フィッシャー情報行列）は次のようになります。
$$ \begin{equation} \bm{I}(\bm{\beta}) = \frac{1}{a(\phi)}\bm{X}^\top\bm{W}\bm{X} \end{equation} $$
ここで $\bm{W} = \text{diag}(w_i)$ は対角重み行列で、$w_i = 1/(V(\mu_i)(g'(\mu_i))^2)$ です。
この情報行列は推定量の漸近的な精度を決定し、IRLSアルゴリズムの各反復で使われます。
スコア方程式は $\bm{\mu}$ を通じて $\bm{\beta}$ に非線形に依存するため、反復法で解く必要があります。次にIRLSアルゴリズムを導出しましょう。
IRLS（反復重み付き最小二乗法）
アルゴリズムの導出
スコア方程式をニュートン法（フィッシャー・スコアリング法）で解きます。更新式は次の通りです。
$$ \bm{\beta}^{(t+1)} = \bm{\beta}^{(t)} + [\bm{I}(\bm{\beta}^{(t)})]^{-1}\bm{U}(\bm{\beta}^{(t)}) $$
ここで $\bm{U} = \partial\ell/\partial\bm{\beta}$ はスコアベクトル、$\bm{I}$ はフィッシャー情報行列です。
右辺を整理すると、次の形になります。
$$ \bm{\beta}^{(t+1)} = (\bm{X}^\top\bm{W}^{(t)}\bm{X})^{-1}\bm{X}^\top\bm{W}^{(t)}\bm{z}^{(t)} $$
ここで $\bm{z}^{(t)}$ は作業応答変量（working response）と呼ばれ、次のように定義されます。
$$ \begin{equation} z_i^{(t)} = \eta_i^{(t)} + (y_i – \mu_i^{(t)})g'(\mu_i^{(t)}) \end{equation} $$
$z_i$ は、現在の線形予測子 $\eta_i$ に残差 $y_i – \mu_i$ のリンク関数スケールでの補正を加えたものです。
この更新式は重み付き最小二乗法の正規方程式と全く同じ形をしています。つまり、各反復で「重み $\bm{W}^{(t)}$ で重み付けした作業応答変量 $\bm{z}^{(t)}$ に対する重回帰分析」を繰り返しているのです。これが反復重み付き最小二乗法（IRLS: Iteratively Reweighted Least Squares）という名前の由来です。
代表的なGLMのIRLS
ポアソン回帰（対数リンク）の場合、IRLSの各量は次のようになります。
- 逆リンク: $\mu_i = e^{\eta_i}$
- 重み: $w_i = \mu_i$（分散関数 $V(\mu) = \mu$、$g'(\mu) = 1/\mu$ より $w = 1/(V \cdot (g’)^2) = \mu$）
- 作業応答変量: $z_i = \eta_i + (y_i – \mu_i)/\mu_i$
ロジスティック回帰（ロジットリンク）の場合:
- 逆リンク: $\mu_i = 1/(1 + e^{-\eta_i})$
- 重み: $w_i = \mu_i(1 – \mu_i)$
- 作業応答変量: $z_i = \eta_i + (y_i – \mu_i)/(\mu_i(1-\mu_i))$
収束と初期値
IRLSの収束は通常高速（5-10回の反復で収束）ですが、初期値の選択が重要です。一般的な初期値の選択法は次の通りです。
- $\hat{\mu}_i^{(0)} = y_i + \delta$（観測値に小さな定数を加える）
- $\hat{\eta}_i^{(0)} = g(\hat{\mu}_i^{(0)})$
- $\hat{\bm{\beta}}^{(0)} = (\bm{X}^\top\bm{X})^{-1}\bm{X}^\top\hat{\bm{\eta}}^{(0)}$
正準リンクを使う場合、対数尤度の凹性が保証されるため、任意の初期値からの収束が保証されます。
偏差とモデルの適合度
偏差（Deviance）
GLMにおける残差二乗和に相当する量が偏差（Deviance）です。偏差は、推定モデルの対数尤度と飽和モデル（saturated model、パラメータ数 = データ数で完全にあてはめたモデル）の対数尤度の差として定義されます。
$$ \begin{equation} D = 2[\ell(\hat{\bm{\theta}}_{\text{sat}}) – \ell(\hat{\bm{\theta}}_{\text{model}})] \end{equation} $$
正規線形モデルでは、偏差は残差二乗和に一致します。ポアソン回帰では次のようになります。
$$ D = 2\sum_{i=1}^{n}\left[y_i\ln\frac{y_i}{\hat{\mu}_i} – (y_i – \hat{\mu}_i)\right] $$
偏差残差
偏差から各観測の寄与を分離した偏差残差は次のように定義されます。
$$ \begin{equation} d_i = \text{sign}(y_i – \hat{\mu}_i)\sqrt{d_i^2} \end{equation} $$
ここで $d_i^2$ は第 $i$ 観測の偏差への寄与です。偏差残差は正規線形モデルの通常の残差のGLM版であり、Q-Qプロットなどのモデル診断に使用されます。
偏差残差がおおむね正規分布に従い、系統的なパターンがなければ、モデルの適合は良好と判断できます。
Pythonでポアソン回帰を実装しましょう。
Pythonによる実装と可視化
ポアソン回帰のスクラッチ実装
IRLSアルゴリズムを用いてポアソン回帰をスクラッチで実装します。
import numpy as np
import matplotlib.pyplot as plt
from scipy import stats
np.random.seed(42)
# --- ポアソン回帰のデータ生成 ---
n = 200
x = np.random.uniform(0, 3, n)
X = np.column_stack([np.ones(n), x])
beta_true = np.array([0.5, 0.8]) # 真のパラメータ
# ポアソン分布からサンプリング
eta = X @ beta_true
mu = np.exp(eta) # 逆リンク関数: exp
y = np.random.poisson(mu)
# --- IRLS によるポアソン回帰 ---
beta = np.zeros(2)
log_lik_history = []
for iteration in range(20):
eta_hat = X @ beta
mu_hat = np.exp(eta_hat)
# 重み行列（対角成分）
w = mu_hat
# 作業応答変量
z = eta_hat + (y - mu_hat) / mu_hat
W = np.diag(w)
beta_new = np.linalg.solve(X.T @ W @ X, X.T @ W @ z)
# 対数尤度
ll = np.sum(y * (X @ beta_new) - np.exp(X @ beta_new))
log_lik_history.append(ll)
if np.max(np.abs(beta_new - beta)) < 1e-8:
print(f"収束: {iteration + 1} 回で収束")
beta = beta_new
break
beta = beta_new
print("=== ポアソン回帰の結果 ===")
print(f"真のパラメータ: beta = {beta_true}")
print(f"推定値: beta = {beta}")
# --- 比較: 線形回帰（不適切） ---
beta_ols = np.linalg.lstsq(X, y, rcond=None)[0]
print(f"OLS推定値: beta = {beta_ols}")
IRLSはわずか数回の反復で収束しています。これはニュートン法に基づく二次収束（各反復で精度が2乗で改善される）の恩恵です。推定値が真のパラメータに近いことも確認できます。
GLMの可視化
ポアソン回帰のフィット、対数スケールでの線形性、偏差残差のQ-Qプロットを可視化します。
import numpy as np
import matplotlib.pyplot as plt
from scipy import stats
np.random.seed(42)
n = 200
x = np.random.uniform(0, 3, n)
X = np.column_stack([np.ones(n), x])
beta_true = np.array([0.5, 0.8])
eta = X @ beta_true
mu = np.exp(eta)
y = np.random.poisson(mu)
# IRLS
beta = np.zeros(2)
for _ in range(20):
eta_hat = X @ beta
mu_hat = np.exp(eta_hat)
w = mu_hat
z = eta_hat + (y - mu_hat) / mu_hat
beta_new = np.linalg.solve(X.T @ np.diag(w) @ X, X.T @ (w * z))
if np.max(np.abs(beta_new - beta)) < 1e-8:
beta = beta_new
break
beta = beta_new
beta_ols = np.linalg.lstsq(X, y, rcond=None)[0]
fig, axes = plt.subplots(1, 3, figsize=(16, 5))
# (a) データとフィット
ax = axes[0]
x_plot = np.linspace(0, 3, 200)
X_plot = np.column_stack([np.ones(200), x_plot])
mu_plot = np.exp(X_plot @ beta)
mu_ols = X_plot @ beta_ols
ax.scatter(x, y, alpha=0.3, s=15, color="steelblue", label="Data")
ax.plot(x_plot, mu_plot, "r-", linewidth=2, label="Poisson GLM")
ax.plot(x_plot, mu_ols, "g--", linewidth=1.5, label="OLS (incorrect)")
ax.set_xlabel("$x$", fontsize=11)
ax.set_ylabel("$y$ (count)", fontsize=11)
ax.set_title("(a) Poisson regression vs OLS", fontsize=12)
ax.legend(fontsize=9)
ax.grid(True, alpha=0.3)
# (b) 対数スケールでの線形性
ax = axes[1]
bins = np.linspace(0, 3, 15)
bin_centers = 0.5 * (bins[:-1] + bins[1:])
bin_means = []
for i in range(len(bins) - 1):
mask = (x >= bins[i]) & (x < bins[i+1])
if mask.sum() > 0:
bin_means.append(np.mean(y[mask]))
else:
bin_means.append(np.nan)
ax.scatter(bin_centers, np.log(np.array(bin_means) + 0.1), color="steelblue",
s=50, label="Binned means (log)")
ax.plot(x_plot, X_plot @ beta, "r-", linewidth=2, label="Linear predictor $\\eta$")
ax.set_xlabel("$x$", fontsize=11)
ax.set_ylabel("$\\ln(\\bar{y})$", fontsize=11)
ax.set_title("(b) Log-linear relationship", fontsize=12)
ax.legend(fontsize=9)
ax.grid(True, alpha=0.3)
# (c) 残差のQ-Qプロット（偏差残差）
mu_fitted = np.exp(X @ beta)
sign = np.sign(y - mu_fitted)
deviance_resid = sign * np.sqrt(2 * np.where(y > 0, y * np.log(y / mu_fitted), 0)
- 2 * (y - mu_fitted))
res_sorted = np.sort(deviance_resid)
theoretical = stats.norm.ppf(np.arange(1, n+1) / (n+1))
ax = axes[2]
ax.scatter(theoretical, res_sorted, alpha=0.5, s=15, c="steelblue")
ax.plot([-3, 3], [-3, 3], "r--", linewidth=1.5)
ax.set_xlabel("Theoretical quantiles", fontsize=11)
ax.set_ylabel("Deviance residuals", fontsize=11)
ax.set_title("(c) Q-Q plot of deviance residuals", fontsize=12)
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig("glm_poisson.png", dpi=150, bbox_inches="tight")
plt.show()
この結果から、GLMの適切さが確認できます。
-
データとフィット（図a）: ポアソンGLM（赤実線）はデータの非線形なパターン（指数的増加）をよく捉えています。一方、通常のOLS（緑破線）は直線を当てはめるため、$x$ が小さい領域では負の予測値が出てしまい、カウントデータには不適切です。GLMはリンク関数を通じて予測値が常に非負であることを保証しています
-
対数スケールでの線形性（図b）: $\ln(\bar{y})$ と $x$ の関係がほぼ線形であり、対数リンク関数の選択が適切であることを確認しています。GLMはこの対数スケールでの線形関係をモデル化しており、元のスケールでは指数的な関係になります
-
偏差残差のQ-Qプロット（図c）: 偏差残差がおおむね直線上に並んでおり、モデルの適合が良好であることを示しています。裾の部分で若干のずれが見られますが、ポアソン分布の離散性を考慮すると許容範囲です
IRLSの収束過程の可視化
IRLSの各反復でのパラメータ値と対数尤度の変化を追跡し、収束の速さを可視化します。
import numpy as np
import matplotlib.pyplot as plt
np.random.seed(42)
n = 200
x = np.random.uniform(0, 3, n)
X = np.column_stack([np.ones(n), x])
beta_true = np.array([0.5, 0.8])
y = np.random.poisson(np.exp(X @ beta_true))
beta = np.zeros(2)
beta_history = [beta.copy()]
ll_history = []
for iteration in range(15):
eta_hat = X @ beta
mu_hat = np.exp(eta_hat)
w = mu_hat
z = eta_hat + (y - mu_hat) / mu_hat
beta = np.linalg.solve(X.T @ np.diag(w) @ X, X.T @ (w * z))
beta_history.append(beta.copy())
ll = np.sum(y * (X @ beta) - np.exp(X @ beta))
ll_history.append(ll)
beta_history = np.array(beta_history)
fig, axes = plt.subplots(1, 2, figsize=(13, 5))
ax = axes[0]
ax.plot(range(len(beta_history)), beta_history[:, 0], "o-", label="$\\beta_0$", color="steelblue")
ax.plot(range(len(beta_history)), beta_history[:, 1], "s-", label="$\\beta_1$", color="coral")
ax.axhline(beta_true[0], color="steelblue", linestyle="--", alpha=0.5)
ax.axhline(beta_true[1], color="coral", linestyle="--", alpha=0.5)
ax.set_xlabel("Iteration", fontsize=11)
ax.set_ylabel("Parameter value", fontsize=11)
ax.set_title("(a) IRLS parameter convergence", fontsize=12)
ax.legend(fontsize=10)
ax.grid(True, alpha=0.3)
ax = axes[1]
ax.plot(range(1, len(ll_history)+1), ll_history, "o-", color="steelblue")
ax.set_xlabel("Iteration", fontsize=11)
ax.set_ylabel("Log-likelihood", fontsize=11)
ax.set_title("(b) Log-likelihood convergence", fontsize=12)
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig("glm_irls_convergence.png", dpi=150, bbox_inches="tight")
plt.show()
IRLSの収束特性がこの図から明確に読み取れます。
-
パラメータの収束（図a）: $\beta_0$ と $\beta_1$ の両方が、わずか4-5回の反復で真の値の近くに収束しています。破線は真のパラメータ値を示しており、推定値がそれに接近していく様子が見えます。初期値（ゼロベクトル）から出発しても、迅速に収束する頑健さがIRLSの強みです
-
対数尤度の単調増加（図b）: 対数尤度は各反復で単調に増加し、急速に最大値に到達しています。フィッシャー・スコアリング法は、対数尤度を単調に増加させる性質を持つため、発散の心配がありません
まとめ
本記事では、一般化線形モデル（GLM）の理論を解説しました。
- 3つの構成要素: 確率分布（指数型分布族）、線形予測子（$\eta = \bm{x}^\top\bm{\beta}$）、リンク関数（$g(\mu) = \eta$）。この3要素を組み合わせることで、さまざまなデータ型に対応する回帰モデルが構成できます
- 指数型分布族: $f(y|\theta,\phi) = \exp((y\theta – b(\theta))/a(\phi) + c(y,\phi))$ の形の分布族。$b'(\theta) = \mu$, $b”(\theta) = V(\mu)$ から期待値と分散が直接求まります
- 正準リンク: 自然パラメータ $\theta$ と線形予測子を直接等しくするリンク関数。スコア方程式が $\bm{X}^\top(\bm{y} – \bm{\mu}) = \bm{0}$ という簡潔な形になります
- IRLS: ニュートン法に基づく反復アルゴリズムで最尤推定量を計算。各反復が重み付き最小二乗問題として表現されます
- 偏差: モデルの適合度を測る量で、正規線形モデルの残差二乗和を一般化したもの
- 統一的な枠組み: 正規回帰、ロジスティック回帰、ポアソン回帰、ガンマ回帰はすべてGLMの特殊ケースであり、同じ理論・同じアルゴリズムで扱えます
次のステップとして、以下の記事も参考にしてください。
- ロジスティック回帰 — GLMの二項分布版
- 混合効果モデル — GLMにランダム効果を追加した拡張
- リッジ回帰とLASSO — 正則化回帰との関連