「データ解析のための統計モデリング入門」をPythonで写経 Vol.6 ～ ４章「GLMのモデル選択」①逸脱度・AIC

４章「GLMのモデル選択」

書籍の著者　久保拓弥 先生

書籍「データ解析のための統計モデリング入門」４章「GLMのモデル選択－AICとモデルの予測の良さ－」の Python写経活動記録 です。 

４章は AIC を学び、複数のモデルから「予測の良いモデル」を選択することを実践します。

この記事は 逸脱度 を学び、逸脱度の側面から見た AIC を考えます。

では書籍を開いて統計モデリングの旅に出かけましょう🚀

はじめに

このブログシリーズは、書籍「データ解析のための統計モデリング入門 一般化線形モデル・階層ベイズモデル・MCMC」（岩波書店、「テキスト」と呼びます）の Python 写経を通じて得た「統計モデリングの楽しさ」をご紹介します。

テキストの紹介と引用表記はリンク先の記事に掲載しています。

準備・サマリー

準備

■ 記事の範囲
この記事はテキスト４章の以下の節を取り扱います。

■ 利用データ
テキスト・サポートサイトのデータファイルを引用しています。
▶️ サポートサイト

Jupyter Notebook ファイルと同一フォルダ内に「data」フォルダを用意して、data フォルダ配下の章別フォルダにデータファイルを格納しています。

■ ライブラリのインポート
この記事で用いるライブラリをインポートします。

# インポート

# 数値計算
import numpy as np
import pandas as pd

# 統計計算
import scipy.stats as stats
import statsmodels.api as sm
import statsmodels.formula.api as smf

# 機械学習
from sklearn.metrics import r2_score

# Rデータセットの読み込み
import rdata

# 可視化
import matplotlib.pyplot as plt
import seaborn as sns
plt.rcParams['font.family'] = 'Meiryo'  # または import japanize_matplotlib

統計モデリング・サマリー

統計モデリングの概要です。
この記事は５つの統計モデルを取り扱います。

■ 統計モデル
線形予測子が異なる５つのポアソン回帰です。

■ モデリング手続き
今回はモデリング自体に注目しないため、モデリング手続きの深堀りを省略します。
その代わり、逸脱度とAICによるモデルの選択に注力して、次のテーマに取り組みます。

• 当てはまりの良さとモデルの複雑さの直感的理解

• ４つの統計モデルの逸脱度

• ４つの統計モデルの AIC

当てはまりの良さとモデルの複雑さの直感的理解

テキスト p.69 ～ の図 4.1「当てはまりの良さとモデルの複雑さ」を軸において、２つのモデルを比べます。

データの読み込み

前回記事と同じデータを利用します。
data3a.csv ファイルを pandas データフレームの data に読み込みます。

# データの読み込み
data = pd.read_csv('./data/ch03/data3a.csv')
print('data.shape: ', data.shape)
data.head()

【実行結果】
データの個数（標本サイズ）は 100 です。
100 個体の植物に関する仮想実験の観測データです。

【変数の説明】
データの各変数は、個体から採れた種子数、個体の体サイズ、肥料を与えたかどうかの情報です。

統計モデルをデータに当てはめ

２つのポアソン回帰モデルをフィッティングします。
推定するパラメータの数に着目します。

◆ ◆ ◆

１つ目は「切片モデル」（実は一定モデル）です。
線形予測子が切片のみのポアソン回帰です。

🔷 確率分布と確率質量関数
個体 $${i}$$ の種子数 $${y_i}$$ はパラメータ $${\lambda_i}$$（平均種子数）のポアソン分布に従います。

🔷 リンク関数と線形予測子
リンク関数は「対数」、線形予測子は「$${\beta_1}$$」です。
推定するパラメータは $${\beta_1}$$ の１つなので、パラメータ数 $${k=1}$$ です。

🔷 GLMの当てはめ
statsmodels の glm クラスを利用してポアソン回帰を実装します。

# ポアソン回帰モデルの当てはめ パラメータ数 k=1：切片のみモデル p.68

# 確率分布の設定（リンク関数はデフォルトの対数）
family=sm.families.Poisson()

# GLMの実行 ※formulaに線形予測子を指定
result0_k1 = smf.glm(formula='y ~ 1', data=data, family=family).fit()

# GLMの結果表示
result0_k1.summary()

【実行結果】
当てはまりの良さの指標「最大対数尤度」は $${-237.64}$$ です。

平均種子数 $${\lambda_i}$$ の予測式は次のようになります。

◆ ◆ ◆

２つ目は「6次式モデル」です。
線形予測子が 6 次の多項式のポアソン回帰です。

🔷 確率分布と確率質量関数
個体 $${i}$$ の種子数 $${y_i}$$ はパラメータ $${\lambda_i}$$ のポアソン分布に従います。

🔷 リンク関数と線形予測子
リンク関数は「対数」、線形予測子は「$${\beta_1 + \beta_2 x_i + \beta_3 x_i^2 + \beta_4 x_i^3 + \beta_5 x_i^4 + \beta_6 x_i^5 + \beta_7 x_i^6}$$」です。
変数 $${x_i}$$ は体サイズです。
推定するパラメータは $${\beta_1, \cdots, \beta_7}$$ の７つなので、パラメータ数 $${k=7}$$ です。

🔷 GLMの当てはめ
formula に変数の演算を記述する場合、その演算を $${\mathtt{I(\cdot)}}$$ で囲います。
例えば $${x^2}$$ は $${\mathtt{I(x**2)}}$$ と記述します。

# ポアソン回帰モデルの当てはめ パラメータ数 k=7：6次の多項式 p.68

# 確率分布の設定（リンク関数はデフォルトの対数）
family=sm.families.Poisson()

# GLMの実行 ※formulaに線形予測子を指定
result0_k7 = smf.glm(
    formula='y ~ x + I(x**2) + I(x**3) + I(x**4) + I(x**5) + I(x**6)',
    data=data, family=family).fit()

# GLMの結果表示
result0_k7.summary()

【実行結果】
当てはまりの良さの指標「最大対数尤度」は $${-232.58}$$ です。
実はこちらのモデルの方が最大対数尤度が大きいので、「当てはまりの良いモデル」になっています。

平均種子数 $${\lambda_i}$$ の予測式は次のようになります。

２つの統計モデルの比較

２つの統計モデルを比べましょう。
最初に最大対数尤度を比べます。

「最大対数尤度の大きさ」＝「データへの当てはまりの良さ」で評価するならば「6次式モデル」の方が良いモデルということになります。

◆ ◆ ◆

続いて２つのモデルの予測値を可視化します。
テキスト p.69 図 4.1 に相当します。

# 当てはまりの良さとモデルの複雑さ p.69 図4.1

## 設定
# x軸の値の設定
x_val = np.linspace(data.x.min(), data.x.max(), 100)
# 描画共通設定
params = dict(color='tab:red', lw=3, alpha=0.7)


## 描画処理
# 描画領域の設定
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(10, 4))

# 左の描画
# 観測値の散布図の描画
sns.scatterplot(data=data, x='x', y='y', s=70, alpha=0.7, ax=ax1)
# k=1モデルによる予測値の描画
ax1.plot(data.x, result0_k1.predict(), **params)
# 修飾：x軸ラベル、y軸ラベル、タイトル
ax1.set(xlabel='体サイズ $x$', ylabel='種子数 $y$', title='(A)パラメータ数 $k=1$')

# 左の描画
# 観測値の散布図の描画
sns.scatterplot(data=data, x='x', y='y', s=70, alpha=0.7, ax=ax2)
# k=7モデルによる予測値の描画
ax2.plot(x_val, result0_k7.predict(dict(x=x_val)), **params)
# 修飾：x軸ラベル、y軸ラベル、タイトル
ax2.set(xlabel='体サイズ $x$', ylabel='種子数 $y$', title='(B)パラメータ数 $k=7$');

【実行結果】
左が切片モデル、右が 6 次式モデルです。

【チャートの解釈】
左の切片モデルはシンプルすぎる「水平線」。
改善の余地がありそうです。
右の 6 次式モデルは複雑な曲線を描いています。

この複雑さは良いモデルの構築に必要なのでしょうか？
テキストは「線形予測子を複雑化するにつれて、データへの当てはまりは改善するが、果たして望ましい統計モデルなのでしょうか？」と疑問を投げかけています。

テキストが考える良いモデルとは…
「データへの当てはまりが良いモデルは、必ずしも良いモデルではない」がテキストの主張です。
複数モデルの中から１モデルを選択する際、テキストは「良い予測をするモデルが良いモデルである」ことを最重要視します。
AIC はこの「良い予測をするモデルが良いモデルである」という考えに基づいて設計された「モデル選択規準」である、とテキストは説明しています。

もしかすると…
機械学習の文脈には、モデルが学習データに過度に適合した過学習に陥ると予測性能（汎化性能）が悪化する、があります。
過学習を回避して汎化性能の高いモデルを選択することは、テキストの「当てはまりの良さではなく予測の良さで選択すること」と似ているかもです。

◆ ◆ ◆

テキストの p.68 脚注に気になる一文が…
「何でも直線回帰な人たちが、$${R^2}$$ という指標をモデルの説明力と信じているとか…」とあります。

ひとまず２つのモデルの決定係数 $${R^2}$$ を計算して、次に進みましょう。
scikit-learn の r2_score(y_true, y_pred) で$${R^2}$$ を計算します。

# 学習データによる予測値の決定係数 R² の算出
print(f'k=1 モデル R² = {r2_score(data.y, result0_k1.fittedvalues):8.5f}')
print(f'k=7 モデル R² = {r2_score(data.y, result0_k7.fittedvalues):8.5f}')

【実行結果】
両モデルともに当てはまりが良いとは言えません…

４つの統計モデルの逸脱度

例題データを引き続き利用して、テキスト p.69 ～ の４つのモデルを構築し、逸脱度を計算します。

統計モデルをデータに当てはめ

観測データに４つのモデルをフィッティングします。

◆ ◆ ◆

１つ目は「一定モデル」です。
先ほどの「切片モデル」と同じモデルです。

🔷 確率分布と確率質量関数
個体 $${i}$$ の種子数 $${y_i}$$ はパラメータ $${\lambda_i}$$（平均種子数）のポアソン分布に従います。

🔷 リンク関数と線形予測子
リンク関数は「対数」、線形予測子は「$${\beta_1}$$」です。
推定するパラメータは $${\beta_1}$$ の１つなので、パラメータ数 $${k=1}$$ です。

🔷 GLMの当てはめ
statsmodels の glm クラスを利用してポアソン回帰を実装します。

# ポアソン回帰モデル (A)一定モデル（k=1） p.70

# 確率分布の設定（リンク関数はデフォルトの対数）
family=sm.families.Poisson()

# GLMの実行 ※formulaに線形予測子を指定
result1_a = smf.glm(formula='y ~ 1', data=data, family=family).fit()

# GLMの結果表示
result1_a.summary()

【実行結果】
当てはまりの良さの指標「最大対数尤度」は $${-237.64}$$ です。

平均種子数 $${\lambda_i}$$ の予測式は次のようになります。

◆ ◆ ◆

２つ目は「 $${\mathtt{f}}$$ モデル」です。
施肥効果 $${f_i}$$ の影響を仮定するモデルです。
施肥効果 $${f_i}$$ は $${\mathtt{C}}$$：施肥なし、$${\mathtt{T}}$$：施肥ありの二値の変数です。
説明の便宜上、$${\mathtt{C=0,\ T=1}}$$ と読み替えます。

🔷 確率分布と確率質量関数
個体 $${i}$$ の種子数 $${y_i}$$ はパラメータ $${\lambda_i}$$（平均種子数）のポアソン分布に従います。

🔷 リンク関数と線形予測子
リンク関数は「対数」、線形予測子は「$${\beta_1 + \beta_3 f_i}$$」です。
推定するパラメータは $${\beta_1, \beta_3}$$ の２つなので、パラメータ数 $${k=2}$$ です。

🔷 GLMの当てはめ
statsmodels の glm クラスを利用してポアソン回帰を実装します。

### ポアソン回帰モデル (B)fモデル（k=2） p.70

# 確率分布の設定（リンク関数はデフォルトの対数）
family=sm.families.Poisson()

# GLMの実行 ※formulaに線形予測子を指定
result1_b = smf.glm(formula='y ~ f', data=data, family=family).fit()

# GLMの結果表示
result1_b.summary()

【実行結果】
当てはまりの良さの指標「最大対数尤度」は $${-237.63}$$ です。

平均種子数 $${\lambda_i}$$ の予測式は次のようになります。

◆ ◆ ◆

３つ目は「 $${\mathtt{x}}$$ モデル」です。
体サイズ $${x_i}$$ の影響を仮定するモデルです。
体サイズ $${x_i}$$ は連続値の変数です。

🔷 確率分布と確率質量関数
個体 $${i}$$ の種子数 $${y_i}$$ はパラメータ $${\lambda_i}$$（平均種子数）のポアソン分布に従います。

🔷 リンク関数と線形予測子
リンク関数は「対数」、線形予測子は「$${\beta_1 + \beta_2 x_i}$$」です。
推定するパラメータは $${\beta_1, \beta_2}$$ の２つなので、パラメータ数 $${k=2}$$ です。

🔷 GLMの当てはめ
statsmodels の glm クラスを利用してポアソン回帰を実装します。

# ポアソン回帰モデル (C)xモデル（k=2） p.70

# 確率分布の設定（リンク関数はデフォルトの対数）
family=sm.families.Poisson()

# GLMの実行 ※formulaに線形予測子を指定
result1_c = smf.glm(formula='y ~ x', data=data, family=family).fit()

# GLMの結果表示
result1_c.summary()

【実行結果】
当てはまりの良さの指標「最大対数尤度」は $${-235.39}$$ です。

平均種子数 $${\lambda_i}$$ の予測式は次のようになります。

◆ ◆ ◆

４つ目は「 $${\mathtt{x+f}}$$ モデル」です。
体サイズ $${x_i}$$ と施肥効果 $${f_i}$$ 影響を仮定するモデルです。

🔷 確率分布と確率質量関数
個体 $${i}$$ の種子数 $${y_i}$$ はパラメータ $${\lambda_i}$$（平均種子数）のポアソン分布に従います。

🔷 リンク関数と線形予測子
リンク関数は「対数」、線形予測子は「$${\beta_1 + \beta_2 x_i + \beta_3 f_i}$$」です。
推定するパラメータは $${\beta_1, \beta_2, \beta_3}$$ の３つなので、パラメータ数 $${k=3}$$ です。

🔷 GLMの当てはめ
statsmodels の glm クラスを利用してポアソン回帰を実装します。

# ポアソン回帰モデル (D)x+fモデル（k=3） p.70

# 確率分布の設定（リンク関数はデフォルトの対数）
family=sm.families.Poisson()

# GLMの実行 ※formulaに線形予測子を指定
result1_d = smf.glm(formula='y ~ x + f', data=data, family=family).fit()

# GLMの結果表示
result1_d.summary()

【実行結果】
当てはまりの良さの指標「最大対数尤度」は $${-235.29}$$ です。

平均種子数 $${\lambda_i}$$ の予測式は次のようになります。

４つの統計モデルの可視化

４つの統計モデルによる予測値を可視化します。
テキスト p.70 図 4.2 に相当します。

# 第３章の例題データを説明する４種類のポアソン回帰モデル p.70 図4.2

# x軸の値の設定
x_val = np.linspace(data.x.min(), data.x.max(), 100)

# 描画領域の設定
fig, ax = plt.subplots(2, 2, figsize=(8, 8), tight_layout=True)

# (A)一定モデル
# 予測値の描画
ax[0, 0].plot(data.x, result1_a.predict(), color='tab:blue')
# 修飾：x軸ラベル、y軸ラベル、タイトル
ax[0, 0].set(title='(A) 一定モデル ($k=1$)', ylim=(6, 10))

# (B)fモデル(k=2)
# 処理無しの予測値の描画
ax[0, 1].plot(x_val, result1_b.predict(dict(f=['C']*100)), color='tab:blue',
              ls='--', label='処理なし')
# 施肥処理の予測値の描画
ax[0, 1].plot(x_val, result1_b.predict(dict(f=['T']*100)), color='tab:green',
              label='施肥処理')
# 修飾：x軸ラベル、y軸ラベル、タイトル
ax[0, 1].set(title='(B) f モデル ($k=2$)', ylim=(6, 10))
ax[0, 1].legend()

# (C)xモデル(k=2)
# 処理無しの予測値の描画
ax[1, 0].plot(x_val, result1_c.predict(dict(x=x_val)), color='tab:blue')
# 修飾：x軸ラベル、y軸ラベル、タイトル
ax[1, 0].set(title='(C) x モデル ($k=2$)', ylim=(6, 10))

# (D)fモデル(k=3)
# 処理無しの予測値の描画
ax[1, 1].plot(x_val, result1_d.predict(dict(x=x_val, f=['C']*100)),
              color='tab:blue', ls='--', label='処理なし')
# 施肥処理の予測値の描画
ax[1, 1].plot(x_val, result1_d.predict(dict(x=x_val, f=['T']*100)),
              color='tab:green', label='施肥処理')
# 修飾：x軸ラベル、y軸ラベル、タイトル
ax[1, 1].set(title='(D) x+f モデル ($k=3$)', ylim=(6, 10))
ax[1, 1].legend()

fig.suptitle('種子数データを説明する４種類のポアソン回帰モデル', fontsize=16)
plt.show()

【実行結果】
施肥効果を含むモデルは、施肥処理なし／ありに分けて描画しています。
x 軸は 体サイズ、y 軸は種子数です。

【チャートの解釈】
どの統計モデルが良いのか、この図だけでは判断しかねます…

◆ ◆ ◆

【テキストの提案】
最大対数尤度で示される「データへの当てはまりの良さ」はモデルの良さではない、というのが主張です。
モデルの良さに向かう第一歩は、「当てはまりの悪さ」を示す 逸脱度 からはじめます。

逸脱度

逸脱度（deviance）$${D}$$ は、最大対数尤度 $${\log L^*}$$ に $${\times -2}$$ して求めます。
「当てはまりの良さ」を表す最大対数尤度 $${\log L^*}$$ にマイナスを掛けているので、逸脱度は「当てはまりの悪さ」を表します。

📊 逸脱度 $${D}$$ の公式

「$${-2}$$ を掛ける理由は $${\chi^2}$$ 分布との対応関係が良くなるから」と、テキスト p.71 の脚注 *6 に掲載されています。

statsmodels の GLM の結果を使って、逸脱度を計算します。
最大対数尤度 $${\log L^*}$$ を「llf 属性」で取り出して $${\times -2}$$ します。

# 逸脱度の計算

# GLM結果のリスト
results = [result1_a, result1_b, result1_c, result1_d]
# モデル名のリスト
model_names = ['(a) 一定', '(b) f', '(c) x', '(d) x+f']
# 各モデルの逸脱度の算出
for res, nam in zip(results, model_names):
    print(f'{nam} モデル\t: {-2 * res.llf:6.2f}') 

【実行結果】

以上で逸脱度の計算は完了です！
逸脱度を使えば、AIC も簡単に計算できます！

ところが・・・

逸脱度は statsmodels のGLM の結果に格納されていません。
GLM のサマリーにも表示されていません。
それどころか、逸脱度 deviance には別の値が表示されています！

テキストも「R の glm 関数の実行結果に逸脱度が表示されない」ことに言及しています。
そこで、テキスト p.72 ～に則って、さまざまな逸脱度を確認していきます。

次の図は (c) $${\mathtt{x}}$$ モデルのサマリー表示です。
Deviance 欄には「残差逸脱度」（residual deviance）が表示されています。

statsmodels のGLM の結果から「deviance 属性」で取り出せます。

# 残差逸脱度の取得

# 各モデルの残差逸脱度の取得
for res, nam in zip(results, model_names):
    print(f'{nam} モデル\t: {res.deviance:6.2f}')

【実行結果】

また、テキストの glm 関数の結果に表示されている Null deviance（Null 逸脱度）は、statsmodels のGLM の結果から「null_deviance 属性」で取り出せます。

# Null逸脱度の取得

# 各モデルのNull逸脱度の取得
for res, nam in zip(results, model_names):
    print(f'{nam} モデル\t: {res.null_deviance:6.2f}')

【実行結果】
全てのモデルで同じ値が入っています。
実はNull 逸脱度は、一定モデルの残差逸脱度なのです。

その他にも「最小の逸脱度」「最大の逸脱度」が出てきます！

(c) x モデルで逸脱度を深堀り

$${\mathtt{x}}$$ モデルを例にとって、さまざまな逸脱度の内容を確かめます。

■ 逸脱度のサマリー
テキスト p.72 表 4.1 の逸脱度の表を引用いたします。

◆ ◆ ◆

■ 最小の逸脱度
最小の逸脱度はフルモデルの逸脱度です。
フルモデルは、例題データの場合「データ数 $${100}$$ と等しいパラメータ $${100}$$ 個で当てはめたモデル」です。
テキストによると「フルモデルとは全データを読みあげているようなもの」です。

フルモデルの逸脱度の計算に必要な最大対数尤度は、ポアソン分布の平均パラメータ $${\lambda_i}$$ に 目的変数 $${y_i}$$ をそのまま当てはめて得られる「対数確率」の和で求められます。

フルモデルの最大対数尤度を計算します。
scipy.stats の poisson.logpmf で対数確率を算出します。

# フルモデルの最大対数尤度 p.74
max_lf_full = sum(stats.poisson.logpmf(k=data.y, mu=data.y))
print(f'フルモデルの最大対数尤度: {max_lf_full:.2f}')

【実行結果】

フルモデルの逸脱度＝最小の逸脱度を計算します。

# フルモデルの逸脱度＝最小の逸脱度 p.74
D_full = -2 * max_lf_full
print(f'フルモデルの逸脱度: {D_full:.2f}')

【実行結果】

フルモデルの逸脱度は、４つのモデルで共通の統計量になります。

◆ ◆ ◆

■ 残差逸脱度
残差逸脱度はフルモデルの逸脱度から各モデルの逸脱度を引いたものです。
テキストによると「フルモデルの逸脱度を基準にした、当てはまりの悪さの相対値」です。

フルモデルの逸脱度と $${\mathtt{x}}$$ モデルの逸脱度を使って、$${\mathtt{x}}$$ モデル の残差逸脱度を計算します。

# xモデルの残差逸脱度 p.74
D = -2 * result1_c.llf  # 逸脱度 D
D_resid = D - D_full    # 残差逸脱度 = D - フルモデルの逸脱度
print(f'xモデルの残差逸脱度: {D_resid:.2f}')

【実行結果】
残差逸脱度 $${84.99}$$ は先ほど GLM の結果から取得した残差逸脱度の値と一致しています。

◆ ◆ ◆

■ 最大の逸脱度
最大の逸脱度は Null モデルの逸脱度です。
今回の例題データに当てはめたポアソン回帰の場合、Null モデルは最もパラメータ数が少ないモデル ＝ (a) 一定モデル（切片モデル）です。

Null モデル＝一定モデルの最大対数尤度を取得します。

# nullモデルの最大対数尤度 ※(A)一定モデル p.75
max_lf_null = result1_a.llf
print(f'nullモデルの最大対数尤度: {max_lf_null:.4f}')

【実行結果】

Null モデルの逸脱度＝最大の逸脱度を計算します。

# nullモデルの逸脱度 ※(A)一定モデル 
D_null = -2 * max_lf_null
print(f'nullモデルの逸脱度: {D_null:.2f}\n')

【実行結果】

Null モデルの逸脱度は、４つのモデルで共通の統計量になります。

◆ ◆ ◆

■ Null 逸脱度
Null 逸脱度は最大の逸脱度と最小の逸脱度の差です。
Null モデルの残差逸脱度と一致します。

# Null逸脱度＝最大の逸脱度－最小の逸脱度 ※(A)一定モデル p.74, 75

# null逸脱度、残差逸脱度、最大逸脱度－最小逸脱度の確認　３つの逸脱度は一致する
print(f'xモデルのNull逸脱度\t\t: {result1_c.null_deviance:5.2f}')
print(f'一定モデルの残差逸脱度\t\t: {result1_a.deviance:5.2f}')
print(f'最大の逸脱度－最小の逸脱度\t: {D_null - D_full:.2f}')

【実行結果】

Null 逸脱度は、４つのモデルで共通の統計量になります。

◆ ◆ ◆

ここまで計算した逸脱度を一表にまとめましょう。
テキスト p.76 表 4.2 に相当します。

# 種子数モデルの最大対数尤度と逸脱度 p.75, 76 表4.2

## 各列の値の設定
# モデル
model_names = ['一定', 'f', 'x', 'x+f', 'フル']
# パラメータ数
ks = [len(result1_a.params), len(result1_b.params), len(result1_c.params),
     len(result1_d.params), 100]
# 最大対数尤度
max_lfs = np.array([result1_a.llf, result1_b.llf, result1_c.llf, result1_d.llf,
                    max_lf_full])
# 逸脱度
deviances = -2 * max_lfs
# 残差逸脱度
resid_deviances = [result1_a.deviance, result1_b.deviance, result1_c.deviance,
                   result1_d.deviance, 0]

## データフレーム化
stats1_df = pd.DataFrame({
    'モデル': model_names,
    'k パラメータ数': ks,
    '最大対数尤度': max_lfs,
    '逸脱度': deviances,
    '残差逸脱度': resid_deviances
})
stats1_df.round(1)

【実行結果】

当てはまりの悪さを表す逸脱度 $${D}$$ で評価できると仮定すると（フルモデルを除いて）、逸脱度 $${D}$$ が最も小さい $${\mathtt{x+f}}$$ モデルが良い、ということになります。
ただし、$${\mathtt{x+f}}$$ モデルはパラメータ数が最も多いモデルです。

テキストは、パラメータ数が多いことについて、次のように警鐘を鳴らしています。

いよいよ予測の良さでモデルを選択する規準 AIC に進みます。

４つの統計モデルのAIC

AIC とは

AIC （Akaike's Information Criterion）はテキストによると「よく使われているモデル選択規準のひとつであり、予測の良さを重視するモデル選択規準」です。
モデル選択は、複数のモデルの中から最も適切なモデルを選択するプロセスであり、テキストは「何らかの規準で良いモデルを選択すること」としています。
なお、テキスト p.76 脚注 *14 では「他にもさまざまなモデル選択基準があります」と添えられています。

AIC は、最尤推定したパラメータの個数 $${k}$$、最大対数尤度 $${\log L^*}$$、逸脱度 $${D}$$ を用いて、以下のように計算されます。

📊 AIC の公式

AIC が最も小さなモデルが「予測の良いモデル」となります。
ですので、AIC は値が大きいほど「予測の悪さ」を示します。

データの当てはまりの悪さを表す逸脱度に、パラメータ数の２倍を足すことで予測の悪さを表せるとは、とても不思議です。
次回の記事で、AIC と予測の良さの関係を深堀りします。

４つのモデルの AIC

AIC の値は、statsmodels のGLM の結果から「aic 属性」で取り出せます。

# 各モデルのAICの取得
for res, nam in zip(results, model_names):
    print(f'{nam} モデル\t: {res.aic:6.2f}')

【実行結果】

逸脱度の表に AIC 列を追加します。
$${逸脱度 + 2 \times パラメータ数}$$ で計算しましょう。
テキスト p.77 表 4.3 に相当します。

# AICを追加した表 p.76, 77 表4.3
stats1_df['AIC'] = stats1_df['逸脱度'] + 2 * stats1_df['k パラメータ数']
stats1_df.round(1)

【実行結果】

【モデル選択】
AIC の最も小さなモデルは「$${\mathtt{x}}$$ モデル」です。
４つのモデルの中では $${\mathtt{x}}$$ モデルが最も予測の良いモデルです。

最大対数尤度や逸脱度による「データの当てはまり」の観点では $${\mathtt{x+f}}$$ モデル（パラメータ数 $${2}$$）が良いモデルですが、パラメータ数を罰則のように反映した AIC ではパラメータ数 $${1}$$ の　$${\mathtt{x}}$$ モデルが予測の観点で良いモデルとなりました。

◆ ◆ ◆

■ おまけ
なんとなく…決定係数 $${R^2}$$ も計算しておきます。

# 学習データによる予測値の決定係数 R² の算出
for res, nam in zip(results, model_names):
    print(f'{nam} モデル\t: {r2_score(data.y, res.fittedvalues):8.5f}')
print(f'(e) フル モデル\t: {r2_score(data.y, data.y):8.5f}')

【実行結果】
フルモデルを除くと、どのモデルも当てはまりは良くないようです。

まとめ

今回はモデル選択の観点で、良いモデルを選択する指標を学びました。

逸脱度 $${D}$$は「データへの当てはまりの悪さ」を表します。
最大対数尤度 $${\log L^*}$$を $${\times -2}$$ して求めます。

statsmodels や R のGLM の結果には $${D}$$ とは異なる逸脱度が表示されていることを起点にして、さまざまな逸脱度を学びました。

AIC は「モデルの予測の良さ」を表すモデル選択規準です。
AIC の小さなモデルが予測の良いモデルです。
逸脱度 $${D}$$ からパラメータ数 $${k}$$ の２倍を差し引いて計算できます。

今回のブログは以上です。

次回は、AIC が予測の良さを表す理由を確認します。

シリーズの記事

次の記事

前の記事

目次

ブログの紹介

note で７つのシリーズ記事を書いています。
ぜひ覗いていってくださいね！

１．のんびり統計

統計検定２級の問題集を手がかりにして、確率・統計をざっくり掘り下げるブログです。
雑談感覚で大丈夫です。ぜひ覗いていってくださいね。
統計検定２級公式問題集CBT対応版に対応しています。
Python、EXCELのサンプルコードの配布もあります。

２．実験！たのしいベイズモデリング１＆２をPyMC Ver.5で

書籍「たのしいベイズモデリング」・「たのしいベイズモデリング２」の心理学研究に用いられたベイズモデルを PyMC Ver.5で描いて分析します。
この書籍をはじめ、多くのベイズモデルは Ｒ言語＋Stanで書かれています。
PyMCの可能性を探り出し、手軽にベイズモデリングを実践できるように努めます。
身近なテーマ、イメージしやすいテーマですので、ぜひぜひPyMCで動かして、一緒に楽しみましょう！

３．実験！岩波データサイエンス1のベイズモデリングをPyMC Ver.5で

書籍「実験！岩波データサイエンスvol.1」の４人のベイジアンによるベイズモデルを PyMC Ver.5で描いて分析します。
この書籍はベイズプログラミングのイロハをざっくりと学ぶことができる良書です。
楽しくPyMCモデルを動かして、ベイズと仲良しになれた気がします。
みなさんもぜひぜひPyMCで動かして、一緒に遊んで学びましょう！

４．楽しい写経 ベイズ・Python等

ベイズ、Python、その他の「書籍の写経活動」の成果をブログにします。
主にPythonへの翻訳に取り組んでいます。
写経に取り組むお仲間さんのサンプルコードになれば幸いです🍀

５．RとStanではじめる心理学のための時系列分析入門 を PythonとPyMC Ver.5 で

書籍「RとStanではじめる心理学のための時系列分析入門」の時系列分析をPythonとPyMC Ver.5 で実践します。
この書籍には時系列分析のテーマが盛りだくさん！
時系列分析の懐の深さを実感いたしました。
大好きなPythonで楽しく時系列分析を学びます。

６．データサイエンスっぽいことを綴る

統計、データ分析、AI、機械学習、Pythonのコラムを不定期に綴っています。
統計・データサイエンス書籍にまつわる記事が多いです。
「統計」「Python」「数学とPython」「R」のシリーズが生まれています。

７．Python機械学習プログラミング実践記

書籍「Python機械学習プログラミング PyTorch & scikit-learn編」を学んだときのさまざまな思いを記事にしました。
この書籍は、scikit-learnとPyTorchの教科書です。
よかったらぜひ、お試しくださいませ。

最後までお読みいただきまして、ありがとうございました。