📝

一般化線形モデル入門 第7章 解答例

に公開
統計学
一般化線形モデル
ロジスティック回帰モデル
irls
idea

https://zenn.dev/honwakasan/articles/b9ae08fd8a694e の第7章の解答例です.

演習問題 7.1

説明変数を放射線量, 目的変数を白血病での死亡率 ( = 白血病死亡者数 / 合計死亡者数 ) として, 単変量ロジスティック回帰モデル

\begin{aligned} \pi_i = \frac{\exp(\beta_0 + \beta_1 x_i)}{1 + \exp(\beta_0 + \beta_1 x_i)} \end{aligned}

つまり,

\begin{aligned} \log \left(\frac{\pi_i}{1 - \pi_i} \right) = \beta_0 + \beta_1 x_i \end{aligned}

を構築する. 元のデータセットは以下で作成した.

df = pd.DataFrame(
    [
        [13, 5, 5, 3, 4, 18],
        [378, 200, 151, 47, 31, 33]
    ]
)
df.columns = [f"放射線量_{v}" for v in ["0", "1-9", "10-49", "50-99", "100-199", "200+"]]
df.index = ["白血病", "その他のがん"]
df.index.name = "死因"
死因 放射線量_0 放射線量_1-9 放射線量_10-49 放射線量_50-99 放射線量_100-199 放射線量_200+
白血病 13 5 5 3 4 18
その他のがん 378 200 151 47 31 33
合計 391 205 156 50 35 51

説明変数と目的変数の設定

なお、ここでは各放射線量ごとの死亡率をそのまま使用し, 「1-9と10-49をまとめて1-49の1つのbinにする」といった処理はせずにモデル構築することを考える.[1] 実際には, 被曝放射線量と白血病死亡率の間には, 「多量の放射線を浴びたかどうか」が関連していると思われるが, そういったことは今は考慮せずにモデルを作ってみる. 以下のコードで, 説明変数と目的変数に分割する.

df_tmp = df.T
y = (df_tmp["白血病"] / df_tmp["合計"]).values
X = np.stack([np.ones(len(y)), list(range(1, 7))])

次に, 対数尤度関数のスコア及び情報行列を求める.

スコアの算出

各被曝放射線量における白血病死亡者数を Y_i, 死亡者数を n_i, 白血病死亡率を \Pi_i = Y_i / n_i とする. また, Y_i が二項分布 Bin(n_i, \pi_i) に従うと仮定する. なお, サンプル数 N は N = 6 であることに注意する. すると対数尤度関数は

\begin{aligned} l(\pi_1, \ldots, \pi_N; y_1, \ldots, y_N) &= \sum_{i= 1}^N \left(y_i \log \left(\frac{\pi_i}{1 - \pi_i} \right) + n_i \log(1-\pi_i) + \log {}_{{n_i}}\mathrm{C}_{{y_i}} \right) \\ &= \sum_{i= 1}^N \left( y_i(\beta_0 + \beta_1 x_i) - n_i\log(1 + \exp(\beta_0 + \beta_1 x_i)) + \log {}_{{n_i}}\mathrm{C}_{{y_i}} \right) \\ \end{aligned}

となる. 二行目の展開は, ロジスティック回帰モデルであることより従う. パラメータ \beta_0, \beta_1 で対数尤度関数を偏微分することで, スコアを得る.

\begin{aligned} U_0 &= \frac{\partial l}{\partial \beta_0} \\ &= \sum_{i = 1}^N \left(y_i - n_i \left(\frac{\exp(\beta_0 + \beta_1 x_i)}{1 + \exp(\beta_0 + \beta_1 x_i)} \right) \right) \\ &= \sum_{i = 1}^N \left(y_i - n_i \pi_i \right) \\[2ex] \end{aligned}
\begin{aligned} U_1 &= \frac{\partial l}{\partial \beta_1} \\ &= \sum_{i = 1}^N \left(y_i x_i - n_i x_i \left(\frac{\exp(\beta_0 + \beta_1 x_i)}{1 + \exp(\beta_0 + \beta_1 x_i)} \right) \right) \\ &= \sum_{i = 1}^N \left(y_i - n_i \pi_i \right) x_i \\ \end{aligned}

となる.

情報行列の算出

情報行列 \mathfrak{J} とは, その第 (j,k) 成分が

\begin{aligned} \mathfrak{J}_{jk} &= E[U_j U_k] \\ &= \sum_{i = 1}^N \frac{x_{ij} x_{ik}}{Var(Y_i)} \left( \frac{\partial \mu_i}{\partial \eta_i} \right)^2 \end{aligned}

であるような行列であった.[1:1] ここでは, 少し注意が必要である. 本文では, 上記のように記載されているのだが, これはあくまで一般化線形モデルの条件を満たす確率変数 Y_i を対象とした場合である.

今回の問題のように白血病死亡者数 Y_i ではなく, 白血病死亡率 \Pi_i を対象とした場合, \Pi_i = Y_i / n_i を一般化線形モデルの仮定を満たす独立な確率変数とみなしている.

つまり,

\begin{aligned} E[\Pi_i] = E \left[\frac{Y_i}{n_i} \right] = \frac{\mu_i}{n_i}= \pi_i \end{aligned}

であるため, リンク関数と期待値の関係も

\begin{aligned} \eta_i = g(\pi_i) = \log \left(\frac{\pi_i}{1 - \pi_i} \right) \end{aligned}

となる. 故に, 上記の情報行列も

\begin{aligned} \mathfrak{J}_{jk} &= E[U_j U_k] \\ &= \sum_{i = 1}^N \frac{x_{ij} x_{ik}}{Var(\Pi_i)} \left( \frac{\partial \pi_i}{\partial \eta_i} \right)^2 \end{aligned}

のように修正される. 偏微分の項については

\begin{aligned} \frac{\partial \eta_i}{\partial \pi_i} &= \frac{1}{\pi_i (1 - \pi_i)} \end{aligned}

となり, 分散の項は, Y_i \sim Bin(n_i, \pi_i) であることから,

\begin{aligned} Var(\Pi_i) &= \frac{Var(Y_i)}{n_i^2} \\ &= \frac{n_i \pi_i (1 - \pi_i)}{n_i^2} = \frac{\pi_i (1 - \pi_i)}{n_i} \end{aligned}

を得る. 以上を修正された情報行列に代入することで,

\begin{aligned} \mathfrak{J}_{jk} &= \sum_{i = 1}^N n_i x_{ij} x_{ik} \pi_i (1 - \pi_i) \\ \end{aligned}

を得る. 線形成分は定数項を含むので, x_{i0} = 1, x_{i1} = x_i のように読み替えることで, 次の情報行列を得る.

\mathfrak{J} = \left[ \begin{matrix} \displaystyle \sum_{i = 1}^N n_i \pi_i (1 - \pi_i) & \displaystyle \sum_{i = 1}^N n_i x_i \pi_i (1 - \pi_i) \\ \displaystyle \sum_{i = 1}^N n_i x_i \pi_i (1 - \pi_i) & \displaystyle \sum_{i = 1}^N n_i x_i^2 \pi_i (1 - \pi_i) \end{matrix} \right]

更新式の定義

パラメータの更新式として, 次の(4.21)を採用する.

\begin{aligned} \mathbb{b}^{(m)} &= \mathbb{b}^{(m - 1)} + (\mathfrak{J}^{(m - 1)})^{-1}U^{(m - 1)} \\ \end{aligned}

Pythonでのパラメータ最適化

以上で準備が整ったので, 以下のコードで係数を推定する.

def pred_pi(beta, X):
    z = beta @ X
    return np.exp(z) / (1 + np.exp(z))

# 初期値
beta = np.array([0.0, 0.0])

for _ in range(10):
    # piの推定
    pi = pred_pi(beta, X)  # shape (N,)
    
    # スコアの推定 (2,)
    U = np.sum((y - n*pi) * X, axis=1)  # Xの各行で重み付き和
    
    # 情報行列の推定 (2x2)
    W = n * pi * (1 - pi)  # shape (N,)
    J = X @ (W[:, None] * X.T)  # 2xN @ N×2 = 2x2
    beta = beta + np.linalg.inv(J)@U

betaの出力結果は以下になった.

beta_0 beta_1
-10.81132773 0.97888537

係数値が正の値を取ることから, 放射線被曝量と白血病死亡率との間には関連がありそう.[2]

脚注

  1. https://routledge.com/downloads/C9500/C9500_Solutions.pdf には3つの被曝放射線量帯を定義してロジスティック回帰を実行しているため, ここで載せている結果と異なっている. ↩︎ ↩︎

  2. 本当は標準誤差を求めて係数が0かどうかの検定もやった方が良いのだが, 面倒なので省略 ↩︎

Discussion