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

2025/11/09に公開
 演習問題4.1https://zenn.dev/honwakasan/articles/b9ae08fd8a694e の第4章の解答例です.
この章では, 数値データを用いた解析があるが, pythonを使って実装する.
表4.5を作成しておく.
df_45 = pd.DataFrame(
    {
        "year":[1984, 1985, 1986, 1987, 1988],
        "四半期_1":[1, 27, 43, 88, 110],
        "四半期_2":[6, 39, 51, 97, 113],
        "四半期_3":[16, 31, 63, 91, 149],
        "四半期_4":[23, 30, 70, 104, 159]
    }
)

 (a)以下のコードで, 期間と患者数を散布図で描く.
array_45 = df_45.drop(["year"], axis = 1).values.reshape(-1)
plt.scatter(
    list(range(1, len(array_45) + 1)),
    array_45
)

plt.xlabel("期間")     # x軸ラベル
plt.ylabel("患者数")   # y軸ラベル

plt.xticks(range(1, len(array_45) + 1))
plt.show()
出力結果は以下であった.


 (b)期間および患者数を対数変換し散布図を描く.
plt.scatter(
    np.log(list(range(1, len(array_45) + 1))),
    np.log(array_45)
)

plt.xlabel("期間（対数表示）")     # x軸ラベル
plt.ylabel("患者数（対数表示）")   # y軸ラベル

plt.show()
出力結果は以下となった.


 (c)IRLSをする際に必要な W, b, z を求め, 最終的な \beta を推定する. 本問における一般化線形モデルは次のようなものであった.
Y_i を各期間 i における患者数とし, i = 1, \ldots, 20 とする. パラメータ \lambda_i のポアソン分布に従うと仮定する. サンプル i の説明変数を x_i とすると

\begin{aligned}
E[Y_i] &= \lambda_i \\[6pt]
g(\lambda_i) &= \log(\lambda_i) = \beta_1 + \beta_2 x_i  \\[6pt]
\eta_i &= g(\lambda_i) \\[6pt]
\end{aligned}
である.

 \beta の初期値を設定初期値 \beta_1^{(0)} = 0, \beta_2^{(0)} = 0 と設定する.

 Wの値を求める一般化線形モデルの式から
 \eta_i^{(m-1)} = \beta_1^{(m-1)} + \beta_2^{(m - 1)} x_i を得る. また,
\eta_i^{(m-1)} = \log(\lambda_i^{(m-1)}) より \lambda_i^{(m-1)} で微分すると

\begin{aligned}
\left( \frac{\partial \eta_i^{(m-1)}}{\partial \lambda_i^{(m-1)}} \right) = \frac{1}{\lambda_i^{(m-1)}} \\[6pt]
\end{aligned}
となる. さらに, Y_i はポアソン分布に従うので, 期待値と分散が等しい. 以上から, (4.23)より,

\begin{aligned}
w_{ii}^{(m-1)} &= \frac{1}{Var(Y_i)} \left(\frac{\partial \lambda_i^{(m-1)}}{\partial \eta_i^{(m-1)}}  \right)^2 \\
  &= \frac{1}{\lambda_i^{(m-1)}} (\lambda_i^{(m-1)})^2 \\
  &= \lambda_i^{(m-1)}
\end{aligned}
よって,
 W^{(m-1)} = \mathrm{diag}(w_{11}^{(m-1)}, ..., w_{nn}^{(m-1)}) で表せる.

 z を求める(4.24)より,

\begin{aligned}
z_i^{(m - 1)} 
    &= \sum_{k = 1}^p x_{ik} b_k^{(m - 1)} + (y_i - \lambda_i^{(m-1)}) \left( \frac{\partial \eta_i^{(m-1)}}{\partial \lambda_i^{(m-1)}}\right) \\
    &= \beta_1^{(m-1)} + \beta_2^{(m-1)} x_i + (y_i - \lambda_i^{(m-1)}) \frac{1}{\lambda_i^{(m-1)}} \\
    &= \eta_i^{(m-1)} +  \frac{y_i - \lambda_i^{(m-1)}}{\lambda_i^{(m-1)}} \\
\end{aligned}
である[1]. よってベクトル表記すると,

  
\mathrm{z}^{(m - 1)} =
\begin{bmatrix}
z_1^{(m -1)} \\ \vdots \\ z_n^{(m -1)}
\end{bmatrix}

である.

 bの推定値を更新(4.25)より

\begin{aligned}
\mathbf{b}^{(m)} = \left( X^T W^{(m-1)} X  \right)^{-1} X^T W^{(m - 1)}  \mathbf{z} ^{(m -1)}
\end{aligned}
に対し, これまで計算してきた W, z を代入して b を更新すれば良い.

 Pythonによる実装以上を踏まえ, 以下のコードでIRLSを実行した.

# 1列目は定数部分なので1を入れる.
# x_iはlog(i)とみなすので対数変換も行う.
X = np.stack(
    (np.ones(len(array_45)),
    np.log(list(range(1, len(array_45) + 1)))), 
    axis = 1
)

y = array_45.copy()

# パラメータの初期値
beta = np.ones(X.shape[1])

# 繰り返し回数
n_iters = 10

list_beta = []
for t in range(n_iters):

    # etaの値を推定 
    eta_t = X @ beta

    # 連結関数が対数関数であったことに注意
    lambda_t = np.exp(eta_t)
    W_t = np.diag(lambda_t)

    # (4.24)より 
    z_t = eta_t + ( y - lambda_t ) / lambda_t        

    # (4.25)より以下でbetaを推定
    beta = np.linalg.inv((X.T)@W_t@X)@(X.T)@W_t@z_t

    list_beta.append(beta)

df_beta = pd.DataFrame(list_beta, columns = ["beta_0", "beta_1"])
なお, 出力結果は以下である. 数回程度で収束しているのがわかる.


回数
beta_0
beta_1


0
0.652354
1.65199

1
0.841857
1.42955

2
0.98454
1.33373

3
0.995952
1.32664

4
0.995998
1.32661

5
0.995998
1.32661

6
0.995998
1.32661

7
0.995998
1.32661

8
0.995998
1.32661

9
0.995998
1.32661


 (d)statsmodelsライブラリを用いて, 同様のポアソン回帰を実行し, (c)の結果を比較する.
import statsmodels.api as sm
import statsmodels.formula.api as smf

data = pd.DataFrame(
    np.stack((X[: , 1], y),axis = 1),
    columns = ["x", "y"]
)
model = smf.glm(formula='y ~ x', data=data, family=sm.families.Poisson())
result = model.fit()

print(pd.DataFrame(result.params.values.reshape(1, 2), columns = ["beta_0", "beta_1"]))
出力結果は以下となった. (c)の結果と一致していることが確認できた.



beta_0
beta_1


0
0.995998
1.32661


 演習問題4.2問題文の表4.6のデータを以下のコードで作成する.
import pandas as pd
import numpy as np
import japanize_matplotlib
import matplotlib.pyplot as plt
import statsmodels.api as sm
import statsmodels.formula.api as smf

# データ
x = [3.36, 2.88, 3.63, 3.41, 3.78, 4.02, 4.00, 4.23, 3.73,
     3.85, 3.97, 4.51, 4.54, 5.00, 5.00, 4.72, 5.00]

y = [65, 156, 100, 134, 16, 108, 121, 4, 39,
     143, 56, 26, 22, 1, 1, 5, 65]

# DataFrame 作成
df = pd.DataFrame({
    "x": x,
    "y": y
})

 (a)散布図を以下のコードで作成する.
plt.title("xとyの散布図")
plt.scatter(df["x"], df["y"])
plt.xlabel("x（初期の白血球数）")
plt.ylabel("y（生存時間;週）")
plt.show()
出力結果は以下であった. 初期白血球数が多くなると, 生存時間が減少していくように見える.


 (b)
E[Y_i] = \exp(\beta_1 + \beta_2 x_i)
として表されるとき, 連結関数として g(x) = \log x を用いれば

\log E[Y_i] = \beta_1 + \beta_2 x_i
と表せるので, 対数関数が適切である.

 (c)
\begin{aligned}
E[Y] 
 &= \int_0^\infty y e^{-y \theta} dy \\
 &= \frac{1}{\theta} \\

E[Y^2] 
 &= \int_0^\infty y^2 e^{-y \theta} dy \\
 &= \frac{2}{\theta^2}
\end{aligned}
となるので

\begin{aligned}
Var[Y] &= E[Y^2] - E[Y]^2 \\
    &= \frac{1}{\theta^2}
\end{aligned}
また, statsmodelsでモデルを当てはめる.
import statsmodels.api as sm
import statsmodels.formula.api as smf

model_exp = smf.glm(
    formula="y ~ x",
    data=df,
    family=sm.families.Gamma(link=sm.families.links.log())
)
results_exp = model_exp.fit()
結果は以下であった. xの係数が負の値になっており, 初期白血球数の増加が生存時間に負の影響を与えていることがわかる. これは, (a)の散布図とも整合的な結果である.


item
coef
std err
z
p-value
下側2.5%点
上側2.5%点


Intercept
8.4775
1.603
5.287
0.000
5.335
11.620

x
-1.1093
0.387
-2.865
0.004
-1.868
-0.350


 (c) スクラッチの場合Y_i \sim \operatorname{Ex}(\theta_i), E[Y_i] = \exp(\beta_1 + \beta_2 x_i) = \mu_i とする. 連結関数を対数関数とすれば

U_1 = \sum_{i = 1}^N \left(\frac{Y_i - \mu_i}{\mu_i} \right),\quad  U_2 = \sum_{i = 1}^N \left(\frac{Y_i - \mu_i}{\mu_i} \right)x_i 
である. また, 情報行列は

\begin{aligned}
\mathfrak{J} = 
\begin{pmatrix} 
N & \displaystyle \sum_{i = 1}^N x_i \\
\displaystyle
 \sum_{i = 1}^N x_i & \displaystyle
 \sum_{i = 1}^N x_i^2\\
\end{pmatrix}
\end{aligned}
であるから, 更新式

\boldsymbol{b}^{(t+1)} = \boldsymbol{b}^{(t)} + (\mathfrak{J}^{(t)})^{-1} \boldsymbol{U}
により更新すれば良い. 更新は以下のPythonコードで実行する.
import numpy as np
import pandas as pd
from tqdm import tqdm
from scipy.stats import norm

# 情報行列
J = np.array(
    [
        [len(df), df["x"].sum()],
        [df["x"].sum(), (df["x"]**2).sum()]
    ]
)

beta = np.zeros(2)
tol = 1e-8
max_iter = 1000
for j in tqdm(range(max_iter)):
    
    mu_i = np.exp(beta[0] + beta[1] * df["x"])
    
    U = np.array(
        [
            ((df["y"] - mu_i) / mu_i).sum(),
            (((df["y"] - mu_i) / mu_i) * df["x"]).sum()
        ]
    )

    delta = np.linalg.inv(J) @ U
    beta_new = beta + delta

    # 収束判定
    if np.linalg.norm(beta_new - beta, ord=2) < tol:
        beta = beta_new
        print(f"Converged at iteration {j+1}")
        break

    beta = beta_new.copy()
得られた結果は以下である. なお, 標準誤差は \mathfrak{J} の対角成分の平方根を採用した.


item
coef
std err


\beta_1
8.47749765
1.65480771

\beta_2
-1.10929792
0.39965447


 (d)以下のコードで標準化残差を算出する.
# 当てはめ値
df["pred"] = np.exp(results_exp.params["x"] * df["x"] + results_exp.params["Intercept"])

# 標準化残差の算出
df["std_resid"] = (df["y"] - df["pred"]) / df["pred"]
標準化残差として以下を得る. 最もズレが大きいのは, 16行目のx = 65, y = 5.00のとき, r = 2.467310 となるレコードである.


x
y
pred
r


3.36
65
115.613954
-0.437784

2.88
156
196.904661
-0.207738

3.63
100
85.690865
0.166986

3.41
134
109.376030
0.225131

3.78
16
72.555482
-0.779479

4.02
108
55.596488
0.942569

4.00
121
56.843734
1.128643

4.23
4
44.043046
-0.909180

3.73
39
76.693459
-0.491482

3.85
143
67.134675
1.130047

3.97
56
58.767262
-0.047088

4.51
26
32.283752
-0.194641

4.54
22
31.227063
-0.295483

5.00
1
18.746519
-0.946657

5.00
1
18.746519
-0.946657

4.72
5
25.574902
-0.804496

5.00
65
18.746519
2.467310

残差をプロットすると以下のようになっており, 大多数は小さい値であることがわかる.


 演習問題4.3
 最尤推定量の導出Y_1, \ldots , Y_N \sim N(\log \beta, \sigma^2) なので対数尤度関数は

\begin{aligned}
l(\beta; y_1, \ldots, y_N) &= - \frac{n}{2} \log 2\pi \sigma^2 - \frac{1}{2\sigma^2} \sum_{i = 1}^N (y_i - \log \beta)^2 \\
\end{aligned}
であり \beta で偏微分すると

\begin{aligned}
U 
 &= \frac{\partial l(\beta; y_1, \ldots, y_N)}{\partial \beta} \\
 &= \sum_{i = 1}^N \frac{y_i - \log \beta}{\sigma^2} \frac{1}{\beta} \\
\end{aligned}
となる. U = 0 とおくことで, 最尤推定量 \beta_{MLE} = \exp(\bar{y}) を得る.

 (4.18)式のチェックこの問題設定において, 本文(4.18)

 
\begin{aligned}
U_j = \sum_{i = 1}^N \frac{y_i - \mu_i}{Var(Y_i)}x_{ij} \left(\frac{\partial \mu_i}{\partial \eta_i} \right)
\end{aligned}

が成り立つことを確認する.
 
\begin{aligned}
E[Y_i] &= \mu_i = \log \beta \\
Var(Y_i) &= \sigma^2 \\
\frac{\partial \mu_i}{\partial \eta_i}x_{ij} &= \frac{\partial \mu_i}{\partial \beta_j} = \frac{\partial \log\beta}{\partial \beta} = \frac{1}{\beta}
\end{aligned}
より, これらを(4.18)式の右辺に代入することで, U と同じ値を得る.

 (4.25)式のチェック最尤推定量 \beta_{MLE} が(4.25)式

 
\begin{aligned}
X^T W X b^{(m)} = X^T W z
\end{aligned}

を満たしているか確認する.[2] b^{(m)} に \beta_{MLE} を代入し, X^T W X に 情報量 \mathfrak{J} を代入すると, (4.25)式の左辺は
 
\begin{aligned}
X^T W X b^{(m)} &= \mathfrak{J} \beta_{MLE} \\
    &= \sum_{i = 1}^N \frac{x_{ik}^2}{Var(Y_i)} \left(\frac{\partial \mu_i}{\partial \eta_i} \right)^2 \\
    &= \frac{N}{\beta_{MLE}^2 \sigma^2}
\end{aligned}
となる.[3] 一方, (4.25)式の右辺は
 
\begin{aligned}
X^T W z 
 &= \sum_{i = 1}^N \frac{x_{ij}^2}{Var(Y_i)} \left(\frac{\partial \mu_i}{\partial \eta_i} \right)^2 
  + \sum_{i = 1}^N \frac{y_i - \mu_i}{Var(Y_i)}x_{ij} \left(\frac{\partial \mu_i}{\partial \eta_i} \right) \\
 &= \frac{N}{\beta_{MLE}^2 \sigma^2} + U \\
 &= \frac{N}{\beta_{MLE}^2 \sigma^2}

\end{aligned}
となる. ただし, 右辺第2式において, U = 0 であることを用いた.

脚注
説明変数を表す行列 X の1列目は, 線形モデルの定数項に相当するので, x_{ik} のように記載されていても, 値は1である. なお, 定数項を除外したモデルの場合はその限りではない. ↩︎
\beta の最尤推定量 \beta_{MLE} が(4.25)式を満たすかどうかをチェックしろ、というのが著者の意図だと理解したのでこのような証明を行っている. 本文の演習問題には, 「・・・また, この場合について式（4.18）と(4.25)が成り立つことを確かめよ.」と記載されており, "この場合" が指している事象がやや曖昧なのだが上記ように解釈した. ↩︎
本文では, パラメータ数がp個ある場合を想定しているため, x_{ij} や x_{ik} のように記載されているが, 今回はパラメータは \beta の1個だけなので, このように記載している. 本当は, x_i のように記載した方が見通しやすいかもしれない. ↩︎
一般化線形モデル入門第4章解答例

演習問題4.1

(a)

(b)

(c)

$\beta$ の初期値を設定

Wの値を求める

$z$ を求める

bの推定値を更新

Pythonによる実装

(d)

演習問題4.2

(a)

(b)

(c)

(c) スクラッチの場合

(d)

演習問題4.3

最尤推定量の導出

(4.18)式のチェック

(4.25)式のチェック

Discussion

回数	beta_0	beta_1
0	0.652354	1.65199
1	0.841857	1.42955
2	0.98454	1.33373
3	0.995952	1.32664
4	0.995998	1.32661
5	0.995998	1.32661
6	0.995998	1.32661
7	0.995998	1.32661
8	0.995998	1.32661
9	0.995998	1.32661
item	coef	std err	z	p-value	下側2.5%点	上側2.5%点
Intercept	8.4775	1.603	5.287	0.000	5.335	11.620
x	-1.1093	0.387	-2.865	0.004	-1.868	-0.350
x	y	pred	r
3.36	65	115.613954	-0.437784
2.88	156	196.904661	-0.207738
3.63	100	85.690865	0.166986
3.41	134	109.376030	0.225131
3.78	16	72.555482	-0.779479
4.02	108	55.596488	0.942569
4.00	121	56.843734	1.128643
4.23	4	44.043046	-0.909180
3.73	39	76.693459	-0.491482
3.85	143	67.134675	1.130047
3.97	56	58.767262	-0.047088
4.51	26	32.283752	-0.194641
4.54	22	31.227063	-0.295483
5.00	1	18.746519	-0.946657
5.00	1	18.746519	-0.946657
4.72	5	25.574902	-0.804496
5.00	65	18.746519	2.467310