🔥

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

に公開
2件
統計学
ポアソン回帰
一般化線形モデル
idea

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

演習問題9.2

(a)

各カテゴリーにおける請求率を計算し、AGE, CAR, DISTとの関係をプロットする. サンプルデータは以下で作成した.

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
from scipy.stats import chi2

plt.style.use("ggplot")

df = pd.DataFrame(
    {
        "CAR":[1]*4 + [2]*4+ [3]*4+[4]*4,
        "AGE":[1, 2, 3, 4]*4,
        "y_d0":[65, 65, 52, 310, 98, 159, 175, 877, 41, 117, 137, 477, 11, 35, 39, 167],
        "n_d0":[317, 476, 486, 3259, 486, 1004, 1355, 7660, 223, 539, 697, 3442, 40, 148, 214, 1019],
        "y_d1":[2, 5, 4, 36, 7, 10, 22, 102, 5, 7, 16, 63, 0, 6, 8, 33],
        "n_d1":[20, 33, 40, 316, 31, 81, 122, 724, 18, 39, 68, 344, 3, 16, 25, 114]
        
    }
)

プロットの仕方はいくつか考えられるが, AGE, CAR, DIST と請求率の折れ線グラフを描いてみる.

# 合計
df["n"] = df["n_d0"] + df["n_d1"]
df["y"] = df["y_d0"] + df["y_d1"]

fig, axes = plt.subplots(1, 2, figsize=(14, 6), sharey=True)

# AGEごとに請求率をプロット
for age, group in df.groupby("AGE"):
    group["ratio"] = group["y"] / group["n"]
    axes[0].scatter(group["AGE"], group["ratio"], marker='o', label=f"AGE = {age}")

axes[0].set_xlabel("AGE")
axes[0].set_ylabel("claim_ratio")
axes[0].set_title("Claim Ratio by AGE")
axes[0].legend(title="AGE")

# CARごとに
for car, group in df.groupby("CAR"):
    group["ratio"] = group["y"] / group["n"]
    axes[1].scatter(group["CAR"], group["ratio"], marker='o', label=f"CAR = {car}")

axes[1].set_xlabel("CAR")
axes[1].set_title("Claim Ratio by CAR")
axes[1].legend(title="CAR")

plt.tight_layout()
plt.show()

AGEとCARの可視化結果は以下となった. AGEが高くなるほど（ = 契約者年齢が高くなるほど）請求率は下がり, CARが高くなるほど（ = 乗用車の保険区分が高くなるほど）請求率は上がった.

また, DISTと請求率の関係は以下のコードで可視化した.

plt.figure(figsize = (8, 5))
plt.plot(df.index,df["y_d0"] / df["n_d0"], marker='o', label = "DIST = 0")
plt.plot(df.index,df["y_d1"] / df["n_d1"], marker='o', label = "DIST = 1")
plt.ylabel("claim_ratio")
plt.legend()
plt.show()

可視化結果は以下となった. DIST = 1（ = 住居区分がロンドンなどの主要都市）の方が請求率が高い傾向にある.

(b)

ポアソン回帰を適用して, AGE, CAR, DISTの3変数の主効果と交互作用を求める. 問題文中にある通り, AGEやCARを指示変数に変換する必要があるため, pandas.get_dummies などで変換する.

car_age_cols = ["CAR", "AGE"]

# DIST=0 と DIST=1 のデータを作成
df_arrange = pd.concat([
    df[car_age_cols + ["y_d0", "n_d0"]]
        .rename(columns={"y_d0": "y", "n_d0": "n"})
        .assign(DIST=0),
    
    df[car_age_cols + ["y_d1", "n_d1"]]
        .rename(columns={"y_d1": "y", "n_d1": "n"})
        .assign(DIST=1)
]).reset_index(drop=True)

# CARとAGEをダミー変数に変換
df_dummies = pd.get_dummies(df_arrange, columns = ["CAR", "AGE"], dtype = int)

# 交互作用項の追加
for i in range(2, 5):
    df_dummies[f"CAR_{i}*DIST"] = df_dummies[f"CAR_{i}"] * df_dummies["DIST"]
    df_dummies[f"AGE_{i}*DIST"] = df_dummies[f"AGE_{i}"] * df_dummies["DIST"]
    for j in range(2, 5):
        df_dummies[f"AGE_{i}*CAR_{j}"] = df_dummies[f"AGE_{i}"]*df_dummies[f"CAR_{j}"]
        df_dummies[f"AGE_{i}*CAR_{j}*DIST"] = df_dummies[f"AGE_{i}"]*df_dummies[f"CAR_{j}"]*df_dummies["DIST"]

# カラムの並び替え
sort_cols = [
    'CAR_2', 'CAR_3', 'CAR_4', 
    'AGE_2', 'AGE_3', 'AGE_4', 'DIST',

    'AGE_2*DIST', 'AGE_3*DIST', 'AGE_4*DIST', 
    'CAR_2*DIST', 'CAR_3*DIST', 'CAR_4*DIST', 
    
    'AGE_2*CAR_2', 'AGE_2*CAR_3', 'AGE_2*CAR_4',
    'AGE_3*CAR_2', 'AGE_3*CAR_3', 'AGE_3*CAR_4',
    'AGE_4*CAR_2', 'AGE_4*CAR_3', 'AGE_4*CAR_4',
    
    'AGE_2*CAR_2*DIST', 'AGE_2*CAR_3*DIST', 'AGE_2*CAR_4*DIST',
    'AGE_3*CAR_2*DIST', 'AGE_3*CAR_3*DIST', 'AGE_3*CAR_4*DIST',
    'AGE_4*CAR_2*DIST', 'AGE_4*CAR_3*DIST', 'AGE_4*CAR_4*DIST',
    'y', 'n'
]
df_dummies = df_dummies[sort_cols]

によりデータセットを整形する. df_arrange の先頭5行を表示すると以下である.

CAR_2 CAR_3 CAR_4 AGE_2 AGE_3 AGE_4 DIST CAR_2 * DIST CAR_3*DIST CAR_4*DIST ... AGE_2 * CAR_3 * DIST AGE_2 * CAR_4 * DIST AGE_3 * CAR_2 * DIST AGE_3 * CAR_3 * DIST AGE_3 * CAR_4 * DIST AGE_4 * CAR_2 * DIST AGE_4 * CAR_3 * DIST AGE_4 * CAR_4 * DIST y n
0 0 0 0 0 0 0 0 0 0 ... 0 0 0 0 0 0 0 0 65 317
0 0 0 1 0 0 0 0 0 0 ... 0 0 0 0 0 0 0 0 65 476
0 0 0 0 1 0 0 0 0 0 ... 0 0 0 0 0 0 0 0 52 486
0 0 0 0 0 1 0 0 0 0 ... 0 0 0 0 0 0 0 0 310 3259
1 0 0 0 0 0 0 0 0 0 ... 0 0 0 0 0 0 0 0 98 486

次に, ポアソン回帰の回帰式を以下の式で定式化する. ここで, i は各共変量パターンを表し, 本文では, i = 1, \cdots, 32 である. また, AGEやCARは4種類のカテゴリーが存在するが, AGE_1 及び CAR_1 を基準カテゴリーとみなしてモデリングすることにした.

\begin{aligned} \log(claim\_rate_i) &= \log (n_i) + \beta_0 \\ &+ \sum_{k = 2}^4 \beta^{age}_k (AGE\_k)_i + \sum_{k = 2}^4 \beta^{car}_k (CAR\_k)_i + \beta^{dist} DIST_i \\ &+ \sum_{k = 2}^4 \beta^{age*dist}_k (AGE\_k * DIST)_i + \sum_{k = 2}^4 \beta^{car*dist}_k (CAR\_k * DIST)_i \\ &+ \sum_{k ,\ l = 2}^4 \beta^{age*car}_{k, l} ( AGE\_k * CAR\_l )_i \\ &+ \sum_{k ,\ l = 2}^4 \beta^{age*car*dist}_{k, l} (AGE\_k * CAR\_l * DIST )_i \end{aligned}

次のコードで, Newton-Raphson法による回帰係数の推定を行う.

# Newton-Raphson法
def NewtonRaphson(n_iters, beta, df, expvar):

    df = df.copy()
    
    # 使用変数
    X = np.hstack([np.ones(len(df)).reshape(len(df), 1), df[expvar].values])
    
    y = df["y"].values
    n = df["n"].values

    # ポアソン回帰モデルにおける期待値
    def _mu(X, beta, n):
        return n*np.exp(np.sum(X*beta, axis = 1))
        
    # Newton-Raphsonによる反復更新
    for _ in range(n_iters):

        # スコアUの計算
        mu_ = _mu(X, beta, n)
        U = (X*(y - mu_).reshape(len(df), 1)).sum(axis = 0)
    
        # 情報行列Jの計算
        list_J = []
        for j in range(X.shape[1]):
            X_j = X[:, j]
            list_J_jk = []
            for k in range(X.shape[1]):
                X_k = X[:, k]
                list_J_jk.append(np.sum(mu_*X_j*X_k))
            list_J.append(list_J_jk)
        J = np.array(list_J)
    
        # 係数の更新
        beta = beta + np.linalg.inv(J)@U

    # 標準誤差の推定
    std_error = np.sqrt(np.diag(np.linalg.inv(J)))
    
    return beta, J, std_error

今回は, 交互作用項を含めた形で推定するため, 以下で係数を求める.

# 使用する説明変数
expvar = [v for v in df_dummies.columns if v not in ["y", "n"]]

# 係数の初期値
beta = np.zeros(len(expvar)+1)

# 反復回数
n_iters = 30

# Newton-Raphson法による推定
beta_b, J_b, std_error_b = NewtonRaphson(n_iters, beta, df_dummies, expvar)

df_coef_b = pd.DataFrame({"ポアソン回帰の係数推定値": beta_b, "標準誤差":std_error}, index = ["定数"] + expvar)

df_coef_bの出力結果は以下となる.[1]

係数 推定値 標準誤差
定数 -1.584515 1.240347e-01
CAR_2 -0.016727 1.599647e-01
CAR_3 -0.109085 1.994364e-01
CAR_4 0.293530 3.260272e-01
AGE_2 -0.406516 1.754116e-01
AGE_3 -0.650450 1.860521e-01
AGE_4 -0.768089 1.364200e-01
DIST -0.718071 7.179029e-01
AGE_2*DIST 0.822032 8.548504e-01
AGE_3*DIST 0.650450 8.857852e-01
AGE_4*DIST 0.898451 7.391808e-01
CAR_2*DIST 0.831235 8.175854e-01
CAR_3*DIST 1.130736 8.601017e-01
CAR_4*DIST -27.977383 1.136949e+06
AGE_2*CAR_2 0.164914 2.173997e-01
AGE_2*CAR_3 0.572574 2.524014e-01
AGE_2*CAR_4 0.255636 3.876206e-01
AGE_3*CAR_2 0.204921 2.247971e-01
AGE_3*CAR_3 0.717246 2.574974e-01
AGE_3*CAR_4 0.239020 3.888001e-01
AGE_4*CAR_2 0.202070 1.730744e-01
AGE_4*CAR_3 0.485397 2.123608e-01
AGE_4*CAR_4 0.250490 3.398640e-01
AGE_2 * CAR_2 * DIST -1.184217 9.950476e-01
AGE_2 * CAR_3 * DIST -1.424807 1.051933e+00
AGE_2 * CAR_4 * DIST 28.334457 1.136949e+06
AGE_3 * CAR_2 * DIST -0.429822 9.944071e-01
AGE_3 * CAR_3 * DIST -0.883231 1.038655e+00
AGE_3 * CAR_4 * DIST 28.607983 1.136949e+06
AGE_4 * CAR_2 * DIST -0.804173 8.428485e-01
AGE_4 * CAR_3 * DIST -1.032332 8.881148e-01
AGE_4 * CAR_4 * DIST 28.365895 1.136949e+06

(c)

交互作用項を考慮せず, AGEやCARを連続変数とみなして推定する. つまり, 次のように定式化することを念頭においてモデリングを実施する.

\begin{aligned} \log(claim\_rate_i) &= \log (n_i) + \beta_0 + \beta^{age} AGE_i + \beta^{car} CAR_i + \beta^{dist} DIST_i \end{aligned}

係数推定に使用したコードは以下である.

# 使用する説明変数
expvar = ["CAR", "AGE", "DIST"]

# 係数の初期値
beta = np.zeros(len(expvar)+1)

# 反復回数
n_iters = 100

# Newton-Raphson法による推定
beta_c, J_c, std_error = NewtonRaphson(n_iters, beta, df_arrange, expvar)

# 係数推定値
df_coef = pd.DataFrame({"ポアソン回帰の係数推定値": beta_c, "標準誤差":std_error}, index = ["定数"] + expvar)

df_coefの出力結果は以下となった. CARの係数値が正, AGEの係数値が負であり, DISTが正である. これは, (a) の可視化結果と整合する.

係数 推定値 標準誤差
定数 -1.852532 0.079903
CAR 0.197769 0.020803
AGE -0.176740 0.018489
DIST 0.218646 0.058528

(追記:12/6) (b)における2次の交互作用までの場合

コメントで「3次の交互作用までいれるのは, 効果やメカニズムが複雑になりすぎるので, 2次の交互作用までの結果も確認すべきでは」というコメントをいただいていた. そこで, 2次までの交互作用項を用いて同様の推定を行う.

# 使用する説明変数
expvar = ['CAR_2',
 'CAR_3',
 'CAR_4',
 'AGE_2',
 'AGE_3',
 'AGE_4',
 'DIST',
 'AGE_2*DIST',
 'AGE_3*DIST',
 'AGE_4*DIST',
 'CAR_2*DIST',
 'CAR_3*DIST',
 'CAR_4*DIST',
 'AGE_2*CAR_2',
 'AGE_2*CAR_3',
 'AGE_2*CAR_4',
 'AGE_3*CAR_2',
 'AGE_3*CAR_3',
 'AGE_3*CAR_4',
 'AGE_4*CAR_2',
 'AGE_4*CAR_3',
 'AGE_4*CAR_4'
]

# 係数の初期値
beta = np.zeros(len(expvar)+1)

# 反復回数
n_iters = 30

# Newton-Raphson法による推定
beta_b2, J_b, std_error_b2 = NewtonRaphson(n_iters, beta, df_dummies, expvar)

# (b)と比較するため1つのテーブルにする.
df_coef_b2 = pd.DataFrame({"ポアソン回帰の係数推定値": beta_b2, "標準誤差":std_error_b2}, index = ["定数"] + expvar)
df_coef_b_table = pd.merge(df_coef_b, df_coef_b2, left_index = True, right_index = True, how = "outer", suffixes=('_with_3rd_inter_term', '_without_3rd_inter_term')).reindex(index = df_coef_b.index)

結果は以下のようになった. 係数の符号が反転したもののみ太字で記載した. (a) の結果では, CARの係数値が正, AGEの係数値が負が示唆されており, 今回の結果は (a) の結果に (b) よりも近づく結果となっている. ただ, 以前として CAR_3 の係数値が負のままなどの課題は残る.

係数名 推定値_with_3rd SE_with_3rd 推定値_without_3rd SE_without_3rd
定数 -1.584515 1.240347e-01 -1.608286 0.123205
CAR_2 -0.016727 1.599647e-01 0.016862 0.156663
CAR_3 -0.109085 1.994364e-01 -0.047719 0.191942
CAR_4 0.293530 3.260272e-01 0.220929 0.325960
AGE_2 -0.406516 1.754116e-01 -0.364350 0.171865
AGE_3 -0.650450 1.860521e-01 -0.642968 0.182222
AGE_4 -0.768089 1.364200e-01 -0.740581 0.134666
DIST -0.718071 7.179029e-01 -0.127023 0.302763
AGE_2*DIST 0.822032 8.548504e-01 -0.063619 0.339671
AGE_3*DIST 0.650450 8.857852e-01 0.266343 0.315739
AGE_4*DIST 0.898451 7.391808e-01 0.270904 0.285347
CAR_2*DIST 0.831235 8.175854e-01 0.081737 0.177277
CAR_3*DIST 1.130736 8.601017e-01 0.125001 0.189899
CAR_4*DIST -27.977383 1.136949e+06 0.426333 0.221617
AGE_2*CAR_2 0.164914 2.173997e-01 0.104071 0.211332
AGE_2*CAR_3 0.572574 2.524014e-01 0.485358 0.243014
AGE_2*CAR_4 0.255636 3.876206e-01 0.339813 0.380538
AGE_3*CAR_2 0.204921 2.247971e-01 0.199599 0.217782
AGE_3*CAR_3 0.717246 2.574974e-01 0.662997 0.247280
AGE_3*CAR_4 0.239020 3.888001e-01 0.327455 0.381310
AGE_4*CAR_2 0.202070 1.730744e-01 0.162836 0.168625
AGE_4*CAR_3 0.485397 2.123608e-01 0.421445 0.203750
AGE_4*CAR_4 0.250490 3.398640e-01 0.319137 0.337689
AGE_2CAR_2DIST -1.184217 9.950476e-01 NaN NaN
AGE_2CAR_3DIST -1.424807 1.051933e+00 NaN NaN
AGE_2CAR_4DIST 28.334457 1.136949e+06 NaN NaN
AGE_3CAR_2DIST -0.429822 9.944071e-01 NaN NaN
AGE_3CAR_3DIST -0.883231 1.038655e+00 NaN NaN
AGE_3CAR_4DIST 28.607983 1.136949e+06 NaN NaN
AGE_4CAR_2DIST -0.804173 8.428485e-01 NaN NaN
AGE_4CAR_3DIST -1.032332 8.881148e-01 NaN NaN
AGE_4CAR_4DIST 28.365895 1.136949e+06 NaN NaN
脚注

  1. AGEやCARの符号が, (a)で想定されている符号と逆になっており, 解釈がよくわからない. もしかしたら誤りかもしれない. https://www.academia.edu/5517676/Statistical_Modelling_in_GLIM が引用元となっているがアクセスできなかったため, もし誰か知っている人がいればコメントしていただけると助かります. ↩︎

Discussion

Kien Y. KnotKien Y. Knot

丁寧な演習大変すごいなあと思います。私も学生時代にDobsonを学んでいて、懐かしかったです。
脚注の部分について、私も再現できているわけではなく、リンク先もIndexのみのようなので、これはあくまで仮説なので、要検証なのですが……：

  • データ数nが16のデータに対して、切片、3次の交互作用を含む変数の数pは32なので、n<pとなっており、パラメータの解が一意には定まりません。結果の整合性がないことはこれに起因するもの(過剰適合・あるいは次元の呪い)かなと思います。

  • 交互作用の結果も標準誤差で雑に計算した信頼区間に0を含み、推定結果として「重要ではない」というcの導入に近いものを感じます。また、9章では二次の交互作用までを取っており、実務や実証研究でも3次の交互作用までを取る研究は少ない(メカニズムや効果の解釈が複雑になりすぎる)ので、2次の交互作用までにとどめた結果を確認することも必要かなと思いました。

ほんわかほんわか

コメントありがとうございます。
ご指摘いただいた点、いずれもその通りだと思います。

データ数nが16のデータに対して、切片、3次の交互作用を含む変数の数pは32なので、n<pとなっており、パラメータの解が一意には定まりません。結果の整合性がないことはこれに起因するもの(過剰適合・あるいは次元の呪い)かなと思います。

最小2乗推定量が一意に定まらないので、これはあると思います。何かしらの変数選択機構をモデルに組み込む（ex. 正則化項の導入）、相関の高い変数を除外（多重共線性の防止にもつながる）、など方法は色々あるのかなと思いました。

交互作用の結果も標準誤差で雑に計算した信頼区間に0を含み、推定結果として「重要ではない」というcの導入に近いものを感じます。また、9章では二次の交互作用までを取っており、実務や実証研究でも3次の交互作用までを取る研究は少ない(メカニズムや効果の解釈が複雑になりすぎる)ので、2次の交互作用までにとどめた結果を確認することも必要かなと思いました。

確かに本文では3次の項は取っていませんでしたね。解釈が複雑になり過ぎるのはそうだと思うので、今度2次の交互作用までの結果もそのうち出してみようと思います。信頼区間はちょっと面倒()で、計算していなかったのですが、0を含んでいるのであれば、(c)を検討する余地があったのだと思います。