「データ解析のための統計モデリング入門」をPythonで写経 Vol.10 ～ ６章「GLMの応用範囲を広げる」②ロジスティック回帰と交互作用項

６章「GLMの応用範囲を広げる」

書籍の著者　久保拓弥 先生

書籍「データ解析のための統計モデリング入門」６章「GLMの応用範囲を広げる」の Python写経活動記録 です。 

この記事は 一般化線形モデル（GLM）の一種「ロジスティック回帰」を実践 します。
今回は 交互作用項 を線形予測子に取り入れます。

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

はじめに

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

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

準備と概要

準備

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

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

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

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

# インポート

# 数値計算
import numpy as np
import pandas as pd
from scipy.special import expit  # ロジスティック関数（シグモイド関数）

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

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

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

この記事で扱う統計モデリングの概要です。

■ 統計モデル
ロジスティック回帰と呼ばれる統計モデルです。

■ モデリング手続き

データの確認

■ データの読み込み
今回の統計モデルは前回記事と同じデータに当てはめします。
data4a.csv ファイルを pandas データフレームの data に読み込みます。

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

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

【変数の説明】
個体からの観察種子数 N と生存種子数 y、個体の体サイズ x、肥料を与えたかどうかの情報 f です。
今回の統計モデリングではすべての変数を利用します。
目的変数は生存種子数 y です。

生存種子数 y を観察種子数 N で割った $${q = y/N}$$ を生存確率 q と呼びます。

■ 体サイズ x と生存種子数 y の散布図を施肥処理別に描画
テキスト p.117 図 6.2 の散布図に相当します。

# 体サイズxと生存種子数yの散布図の描画 p.117 図6.2

# 散布図の描画
sns.scatterplot(
    data=data, x='x', y='y', hue='f', s=70, palette=['tab:blue', 'tab:red'],
    alpha=0.7)
# 修飾
plt.xlabel('植物の体サイズ $x_i$', fontsize=14)
plt.ylabel('生存種子数 $y_i$', fontsize=14)
plt.legend(title='施肥処理 f');

【実行結果】
次の２つの傾向を読み取れました。
・体サイズが大きいほど生存種子数が大きい
・体サイズが同じの場合に施肥ありの方が生存種子数が大きい

ロジスティック回帰（交互作用項モデル）

テキスト p.127 ～ の「交互作用項の入った線形予測子」の統計モデルに取り組みます。

交互作用項

交互作用項は 複数の説明変数の積で作られる新しい変数 であり、説明変数間の相互作用をモデルに組み込むために使用されます。
ChatGPTが「交互作用項の一般的な用途」をまとめてくれました。

今回の統計モデルは「体サイズ x と施肥処理 f の交互作用項」を含めます。
交互作用項とその係数を $${\beta_4 x_i f_i}$$ と表して、線形予測子は次のようになります。

ちなみに交互作用項との対比で、体サイズ x、施肥処理 f を「主効果項」と呼びます。

テキストは p.22 脚注 *22 で交互作用に関する２つの見方を教えてくれます。

• 生存種子数の体サイズ依存性が施肥処理の有無で変わると考える

• 生存種子数に対する施肥処理の効果が体サイズに依存すると考える

１点目の生存種子数の体サイズ依存性が施肥処理の有無で大きく変わる例を可視化します。
テキスト p.128 図 6.8 に相当します。
この図はあくまで架空の分析結果の描画であり、例題データに基づく分析結果を表していません。

# 交互作用項が大きいのでサイズ依存性が施肥処理によって大きく変わる例 p.128 図6.8

## モデルの当てはめ
# 交互作用項を含めるモデルの当てはめ
family = sm.families.Binomial()
res_plot = smf.glm(formula='y + I(N-y) ~ x * f', data=data, family=family).fit()
# 係数のβ1（切片）とβ2（体サイズ）の設定
beta1, beta2 = res_plot.params.iloc[0], res_plot.params.iloc[2]
# 作図用の係数のβ3(施肥処理)とβ4（体サイズと施肥処理の交互作用）の設定
beta3, beta4 = 25, -2.5

## 描画
# 設定
N = 8  # 種子数
x_val = np.linspace(data.x.min(), data.x.max(), 101)  # x軸の値
# 施肥なしCの描画
plt.plot(x_val, expit(beta1 + beta2 * x_val) * N, label='C:施肥処理なし')
# 施肥ありTの描画
plt.plot(x_val, expit(beta1 + beta2 * x_val + beta3 + beta4 * x_val) * N,
         color='tab:red', label='T:施肥処理あり')
# 修飾
plt.xlabel('植物の体サイズ $x$', fontsize=14)
plt.ylabel('生存種子数 $y$', fontsize=14)
plt.legend();

【実行結果】

【考察】
生存種子数の体サイズ依存性が施肥処理の有無で逆転しています。
施肥なしの場合（青い線）、体サイズが大きくなるにつれて、平均的な生存種子数が増加する傾向が見られます。
施肥ありの場合（赤い線）、体サイズが大きくなるにつれて、平均的な生存種子数が減少する傾向が見られます。

統計モデルをデータに当てはめ（モデルの理解）

ロジスティック回帰の GLM を観測データに当てはめます。
モデル名は「$${\mathtt{x * f}}$$ モデル」（非公式）です。

■ 確率分布、リンク関数、線形予測子
GLMの３要素である「確率分布」、「リンク関数」、「線形予測子」を定義します。

🔷 確率分布と確率質量関数
生存種子数 $${y_i}$$ は観察種子数 $${N_i}$$、生存確率 $${q_i}$$ の二項分布に従います。

$${\displaystyle \binom{N_i}{q_i}}$$ は二項係数です。組み合わせ $${{}_{N_i} \text{C}_{y_i}}$$ と等しいです。

🔷 リンク関数と線形予測子
リンク関数は「ロジット」、線形予測子は「$${\beta_1 + \beta_2 x_i + \beta_3 f_i + \beta_4 x_i f_i}$$」です。
施肥処理 $${f_i}$$ は説明の便宜上、$${\mathtt{C=0, T=1}}$$ と読み替えます。

統計モデルをデータに当てはめ（当てはめと評価）

■ 統計モデルをデータに当てはめ、の準備
GLMの３要素「確率分布」「リンク関数」「線形予測子」を用いる統計モデルをデータに当てはめします。
統計モデルの対数尤度 $${\log L}$$ が最大になるパラメータ $${\beta_1, \beta_2, \beta_3, \beta_4}$$（線形予測子のパラメータ）を推定します。

🔷 今回の統計モデルの尤度関数 $${L}$$ と対数尤度関数 $${\log L}$$

$${n}$$ は標本サイズです（観察種子数 $${N}$$ と被って…）

◆ ◆ ◆

■ Python ライブラリで統計モデルをデータに当てはめ
statsmodels の glm を利用して統計モデルを実装します。
引数 family には二項分布 sm.families.Binomial() を設定します。
引数 link は未設定とし、デフォルトのリンク関数 logit を使います。
当てはめ結果を変数 result に格納します。

# 生存種子数のデータの交互作用を推定 x * f p.127~128

# 設定
family = sm.families.Binomial()   # GLMの引数familyに与える確率分布＝二項分布

# モデルの当てはめ ※目的変数には生存種子数yと死滅種子数N-yを与える
result = smf.glm(formula='y + I(N-y) ~ x * f', data=data, family=family).fit()
result.summary()

【実行結果】
最下４行の Intercept が $${\beta_1}$$、x が $${\beta_2}$$、f[T.T] が $${\beta_3}$$、x:f[T.T] が $${\beta_4}$$ に対応しており、テキスト p.128  glm() の推定結果 Coefficients に相当します。

【交互作用項がある場合の formula】
$${\mathtt{x * f}}$$ と記述します。

formula='y + I(N-y) ~ x * f'

上の記述で以下の formula と同じことを表わせます。
$${\mathtt{x:f}}$$ は交互作用項単体の表現です。

formula='y + I(N-y) ~ x + f + x:f'

◆ ◆ ◆

■ 当てはめ結果の分析
🔷 係数の推定値
coef に注目します。

パラメータである係数の最尤推定値は 切片（Intercept）$${\beta_1 = -18.5233}$$、x の係数 $${\beta_2 = -0.0638}$$、f の係数 $${\beta_3 = 1.8525}$$、x と f の交互作用項の係数 $${\beta_4 = 0.2163}$$です。
この推定値を線形予測子に当てはめてみます。

施肥なしの場合の $${\text{logit}(q_i)}$$ を算出します。
$${f_i = \mathtt{C}= 0}$$ を代入します。

こちらは施肥ありの場合です。
$${f_i = \mathtt{T}= 1}$$ を代入します。

施肥処理のなし・ありで $${\text{logit}(q_i)}$$ はほとんど変わらないことが分かりました。

◆ ◆ ◆

🔷 パラメータ推定値の評価
係数の推定値の $${p}$$ 値（P>|z|）と 95% 信頼区間（[0.025  0.975]）に注目します。

体サイズ x の係数 $${\beta_2}$$ の $${p}$$ 値は $${0.000 < 0.05}$$ であり、有意水準 5% で統計的に有意です。
体サイズが１単位増加すると生存のオッズは 6.4 倍になります。

切片 $${\beta_1}$$ の $${p}$$ 値は $${0.000 < 0.05}$$ であり、有意水準 5% で統計的に有意です。
ただ、切片だけのときの生存確率 $${q_i}$$ はほぼ０です。

施肥処理 f、体サイズと施肥処理の交互作用項 の $${p}$$ 値は大きく、統計的有意性がありません。
施肥処理 f の主効果も体サイズ x との交互作用も無さそう、です。

係数の推定値の 95% 信頼区間を用いてフォレストプロットを描画して、統計的有意性を可視化してみます。
まずフォレストプロット描画関数を定義します。

# フォレストプロット（係数の推定値と95%信頼区間）の描画関数の定義
# 引数　result:GLMの結果、var_names: 変数名、x_adjust: xlimの両端を加減する調整値

def forest_plot(result, var_names=None, x_adjust=1):

    ## 設定
    # 係数・変数名
    VARS = result.model.exog_names if var_names is None else var_names
    # 係数の推定値
    params = result.params
    # 係数の95%信頼区間
    conf = result.conf_int()

    ## 描画設定
    # 信頼区間のエラーバーの値の算出
    errors = np.array([params - conf[0], conf[1] - params])
    # 変数の個数の取得
    y_len = len(params)
    # 説明変数のy軸上の位置の算出
    y_pos = np.arange(y_len)[::-1]

    ## 描画
    # 描画領域の設定
    plt.figure(figsize=(6, y_len))  # (横:6, 縦:変数の個数)
    # 係数の点と95%信頼区間バーの描画
    plt.errorbar(params, y_pos, xerr=errors, fmt='o')
    # 有意性を判断する偏回帰係数＝０の垂直線の描画
    plt.axvline(0, color='tab:red', ls='--')
    # 修飾
    plt.xlim(round(conf[0].min() - x_adjust), round(conf[1].max() + x_adjust))
    plt.ylim(-1, y_len)
    plt.yticks(y_pos, VARS, fontsize=12)
    plt.xlabel('係数', fontsize=14)
    plt.title('係数の推定値と95%信頼区間')
    plt.grid(lw=0.5, alpha=0.5, axis='x')

    plt.show()

【実行結果】なし

ではフォレストプロットを描画します。

# フォレストプロットの描画
var_names = ['切片', '施肥処理', '体サイズ', '体サイズ×施肥処理']
forest_plot(result, var_names)

【実行結果】
横バーで示す 95% 信頼区間にゼロを含まない体サイズと切片が統計的に有意であり、ゼロを含む施肥処理と交互作用が統計的に有意では無いことが分かります。

◆

（プチまとめ）
統計的な視点によれば、「交互作用がない」という仮説を棄却できません。
おそらく交互作用の効果はないのでしょう。

◆ ◆ ◆

🔷 モデルの評価
AIC で予測の良さを確認します。

# Null Deviance, Residual Deviance, AICの表示
print(f'Null 逸脱度\t: {result.null_deviance:.1f}')
print(f'残差逸脱度\t: {result.deviance:.1f}')
print(f'AIC\t\t: {result.aic:.1f}')

【実行結果】
AIC は 273.6 です。

前回記事の $${\mathtt{x+f}}$$ モデルの AIC は 272.2 でした。
交互作用を含めたことで、予測の良さは悪化しました。
パラメータ数が増えてモデルが複雑化したことが、悪化に繋がったようです。

（参考：前回記事のモデル評価表）

◆

（プチまとめ）
モデルの「予測の良さ」の観点では、交互作用項を追加しても予測能力は改善されず、AIC が大きくなりました。

生存種子数の予測

テキストにならって、生存種子数 y の予測を行って可視化します。
交互作用のないモデルと比べます。
テキスト p.129 図 6.9 に相当します。

■ $${\mathtt{x+f}}$$ モデルと $${\mathtt{x*f}}$$ モデルの予測
交互作用のない $${\mathtt{x+f}}$$ モデルと交互作用のある $${\mathtt{x*f}}$$ モデルの平均生存種子数の予測値を可視化します。

# 交互作用の有無を調べる図示 p.129

## 準備
# x + fモデルの当てはめ
family = sm.families.Binomial()
result_x_plus_f = smf.glm(
    formula='y + I(N-y) ~ x + f', data=data, family=family).fit()

## 描画
# 設定
N = 8  # 種子数
x_val = np.linspace(data.x.min(), data.x.max(), 101)  # x軸の値
xlabel, ylabel = '植物の体サイズ $x$', '生存種子数 $y$'

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

# (A)交互作用のないモデルの描画
ax1.plot(x_val, result_x_plus_f.predict(dict(x=x_val, f=['C']*len(x_val))) * N,
         label='C:施肥処理なし')
ax1.plot(x_val, result_x_plus_f.predict(dict(x=x_val, f=['T']*len(x_val))) * N,
         color='tab:red', label='T:施肥処理あり')
ax1.set(xlabel=xlabel, ylabel=ylabel, title='(A) 交互作用のないモデル')
ax1.legend()

# (B)交互作用のあるモデルの描画
ax2.plot(x_val, result.predict(dict(x=x_val, f=['C']*len(x_val))) * N,
         label='C:施肥処理なし')
ax2.plot(x_val, result.predict(dict(x=x_val, f=['T']*len(x_val))) * N,
         color='tab:red', label='T:施肥処理あり')
ax2.set(xlabel=xlabel, title='(B) 交互作用のあるモデル')
ax2.legend();

【実行結果】
両方のチャートはほぼ同じ曲線を描いています。
交互作用を追加してもモデルの予測はほとんど変化していないことが分かります。

◆

【交互作用項を追加した $${\mathtt{x*f}}$$ モデルの評価まとめ】
$${\mathtt{x*f}}$$ モデルは $${\mathtt{x+f}}$$ モデルに交互作用項を追加したモデルです。
交互作用項の追加によって、AIC は悪化し、モデルの予測値は追加前のモデルとほぼ同じです。
また交互作用項の係数は統計的に有意ではありませんでした。
以上から、交互作用項の追加による統計モデルへの効果（予測性能向上・統計的有意性）は無さそうです。

 アディショナルタイム【x:f モデル】 

■ $${\mathtt{x:f}}$$ モデル、降臨
テキスト p.129 脚注 *23 には神秘的なコメントが書かれています。

こ、こ、これは実践せねば！

（注意）
こちらは趣味のコードです。
ご興味ない方はスルーしてくださって大丈夫です。

■ 線形予測子
$${\mathtt{x:f}}$$ モデルは、線形予測子が切片と交互作用項で構成されます。
交互作用項に関する係数は、施肥ありのときの係数 $${\beta_2}$$ と施肥なしのときの係数 $${\beta_3}$$ に分かれます。

$${\mathbb{I}[\cdot]}$$ は指示関数です。
$${\mathbb{I}[f_i = \mathtt{C}]}$$ は施肥なしのとき $${1}$$、施肥ありのとき $${0}$$ になります。
$${\mathbb{I}[f_i = \mathtt{T}]}$$ は施肥ありのとき $${1}$$、施肥なしのとき $${0}$$ になります。

◆

■ Python ライブラリで統計モデルをデータに当てはめ

# AIC最良モデル x : f p.129脚注

# 設定
family = sm.families.Binomial()   # GLMの引数familyに与える確率分布＝二項分布

# モデルの当てはめ
result2 = smf.glm(formula='y + I(N-y) ~ x : f', data=data, family=family).fit()
result2.summary()

【実行結果】
$${p}$$ 値、95% 信頼区間は全ての係数が統計的に有意であることを示しています。

■ AIC

# AICの表示
print(f'AIC: {result2.aic:.1f}')

【実行結果】

今までの AIC 最小モデル $${\mathtt{x+f}}$$ モデルの AIC 272.2 よりも小さい値になりました。
$${\mathtt{x:f}}$$ モデルが AIC の観点で予測の良いモデルとして選択されます！

■ 対数オッズ $${\text{logit}(q_i)}$$
パラメータの推定値を対数オッズ $${\text{logit}(q_i)}$$ に当てはめます。

この対数オッズは以下の $${\mathtt{x+f}}$$ モデルの対数オッズととてもよく似ています。

■ オッズ

体サイズが１単位大きくなるときの生存のオッズの増加（倍）は次のようになります。

施肥あり $${f_i = \mathtt{T}}$$ の方が増加倍率が大きいです。

◆

■ 予測
$${\mathtt{x+f}}$$ モデルと $${\mathtt{x:f}}$$ モデルの平均的な生存種子数の予測値を可視化します。

# x+fモデルとx:fモデルの予測値を比べる図示

## 描画
# 設定
N = 8  # 種子数
x_val = np.linspace(data.x.min(), data.x.max(), 101)  # x軸の値
xlabel, ylabel = '植物の体サイズ $x$', '生存種子数 $y$'

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

# (A)x+fモデルの描画（切片を主効果項のモデル）
ax1.plot(x_val, result_x_plus_f.predict(dict(x=x_val, f=['C']*len(x_val))) * N,
         label='C:施肥処理なし')
ax1.plot(x_val, result_x_plus_f.predict(dict(x=x_val, f=['T']*len(x_val))) * N,
         color='tab:red', label='T:施肥処理あり')
ax1.set(xlabel=xlabel, ylabel=ylabel, title='(A) $\mathtt{x+f}$ モデル')
ax1.legend()

# (B)x:fモデルの描画（切片と交互作用項のモデル）
ax2.plot(x_val, result2.predict(dict(x=x_val, f=['C']*len(x_val))) * N,
         label='C:施肥処理なし')
ax2.plot(x_val, result2.predict(dict(x=x_val, f=['T']*len(x_val))) * N,
         color='tab:red', label='T:施肥処理あり')
ax2.set(xlabel=xlabel, title='(B) $\mathtt{x:f}$ モデル')
ax2.legend();

【実行結果】
両モデルの予測値はよく似ています。

結局、$${\mathtt{x+f}}$$ モデル、$${\mathtt{x:f}}$$ モデル、$${\mathtt{x*f}}$$ モデル の３モデルの予測値が類似していることになります。

◆

■ モデルの選択
前回記事の４モデルと今回記事の $${\mathtt{x*f}}$$ モデル、$${\mathtt{x:f}}$$ モデルの AIC 等の評価表を作成します。

# 種子の生存確率モデルのAICなど

## 設定
# モデル名のリスト（最後にフルモデルを表す'フル'を設定）
model_names = ['一定', 'f', 'x', 'x+f', 'x*f', 'x:f', 'フル']
# 各モデルのformulaのリスト（最後にフルモデルで用いる'一定モデル'を設定）
formulas = ['y + I(N-y) ~ 1', 'y + I(N-y) ~ f', 'y + I(N-y) ~ x',
            'y + I(N-y) ~ x + f', 'y + I(N-y) ~ x * f', 'y + I(N-y) ~ x:f',
            'y + I(N-y) ~ 1']
# GLMの引数familyに与える確率分布＝二項分布
family = sm.families.Binomial()
# 結果を格納するデータフレームの初期化
aic_df = pd.DataFrame()


## 関数定義
# モデルの結果からデータフレーム（１行）を作成する関数
def make_df(model_name, k, llf, deviance, resid_deviance, aic):
    return pd.DataFrame({'モデル': [model_name], '説明変数の数': k,
                         '最大対数尤度': llf, '逸脱度':  deviance,
                         '残差逸脱度': resid_deviance, 'AIC': aic})

# モデルのあてはめ～結果から１行のデータフレームを作成する関数
def model_fitting(data, model_name, formula, family):
    # GLM・ロジスティック回帰のあてはめ
    result = smf.glm(formula=formula, data=data, family=family).fit()
    # フルモデルの場合、一定モデルの結果から各種数値を算出する関数でデータフレーム化
    if model_name == 'フル':
        return full_model_result(result, data, model_name)
    # フルモデル以外の場合、あてはめ結果の各種数値をデータフレーム化
    else:
        return make_df(model_name, len(result.params), result.llf,
                       -2 * result.llf, result.deviance, result.aic)

# フルモデルの各種数値を算出する関数
def full_model_result(result, data, model_name):
    # 各種数値を算出
    k = len(data)                                      # 説明変数の数
    deviance = -2 * result.llf - result.null_deviance  # 逸脱度
    max_llf = deviance / -2                            # 最大対数尤度
    aic = -2 * (max_llf - k)                           # AIC
    # 戻り値：データフレーム化した各種数値
    return make_df(model_name, k, max_llf, deviance, 0, aic)


## モデル比較の実行
# 各モデルのあてはめと各種数値算出を繰り返し処理
for model_name, formula in zip(model_names, formulas):
    # モデルのあてはめと各種数値の算出
    tmp_df = model_fitting(data, model_name, formula, family)
    # 結果を格納するデータフレームに追加
    aic_df = pd.concat([aic_df, tmp_df], axis=0)

# データフレームの最終化と結果表示(最小AICをハイライト)
aic_df = aic_df.reset_index(drop=True)  # インデックスのリセット
aic_df.style.highlight_min(subset='AIC', color='lightpink').format(precision=1)

【実行結果】
AIC 最小は $${\mathtt{x:f}}$$ モデルの 271.6 でした。

テキストは p.129 脚注 *23 で「AICが真のモデル（この場合は $${\mathtt{x+f}}$$ モデル）を選ぶわけではない」と説明しています。

そういえば…
真のモデルに収束すると噂の情報量規準があったような…

◆

真のモデルを選ぶ、という観点で BIC さんが顔を出してくれました。
BIC（Bayesian information criterion：ベイズ情報量規準）は真のモデルに収束する性質を持つ情報量規準であり、次の式で算出されます。

$${\log L^*}$$ は最大対数尤度、$${k}$$ はパラメータ数、$${n}$$ は標本サイズです。
BIC 最小のモデルが良いモデル、ということになります。

BIC を含めて評価表を更新します。

# 種子の生存確率モデルの評価指標（BICを追加）

## 設定
# モデル名のリスト（最後にフルモデルを表す'フル'を設定）
model_names = ['一定', 'f', 'x', 'x+f', 'x*f', 'x:f', 'フル']
# 各モデルのformulaのリスト（最後にフルモデルで用いる'一定モデル'を設定）
formulas = ['y + I(N-y) ~ 1', 'y + I(N-y) ~ f', 'y + I(N-y) ~ x',
            'y + I(N-y) ~ x + f', 'y + I(N-y) ~ x * f', 'y + I(N-y) ~ x:f',
            'y + I(N-y) ~ 1']
# GLMの引数familyに与える確率分布＝二項分布
family = sm.families.Binomial()
# 結果を格納するデータフレームの初期化
aic_df = pd.DataFrame()


## 関数定義
# モデルの結果からデータフレーム（１行）を作成する関数
def make_df(model_name, k, llf, deviance, resid_deviance, aic, bic):
    return pd.DataFrame({'モデル': [model_name], '説明変数の数': k,
                         '最大対数尤度': llf, '逸脱度':  deviance,
                         '残差逸脱度': resid_deviance, 'AIC': aic,
                         'BIC': bic})

# モデルのあてはめ～結果から１行のデータフレームを作成する関数
def model_fitting(data, model_name, formula, family):
    # GLM・ロジスティック回帰のあてはめ
    result = smf.glm(formula=formula, data=data, family=family).fit()
    # フルモデルの場合、一定モデルの結果から各種数値を算出する関数でデータフレーム化
    if model_name == 'フル':
        return full_model_result(result, data, model_name)
    # フルモデル以外の場合、あてはめ結果の各種数値をデータフレーム化
    else:
        return make_df(model_name, len(result.params), result.llf,
                       -2 * result.llf, result.deviance, result.aic,
                       result.bic_llf)

# フルモデルの各種数値を算出する関数
def full_model_result(result, data, model_name):
    # 各種数値を算出
    k = n = len(data)                                  # 説明変数の数＝標本サイズ
    deviance = -2 * result.llf - result.null_deviance  # 逸脱度
    max_llf = deviance / -2                            # 最大対数尤度
    aic = -2 * (max_llf - k)                           # AIC
    bic = -2 * max_llf + k * np.log(n)                 # BIC
    # 戻り値：データフレーム化した各種数値
    return make_df(model_name, k, max_llf, deviance, 0, aic, bic)


## モデル比較の実行
# 各モデルのあてはめと各種数値算出を繰り返し処理
for model_name, formula in zip(model_names, formulas):
    # モデルのあてはめと各種数値の算出
    tmp_df = model_fitting(data, model_name, formula, family)
    # 結果を格納するデータフレームに追加
    aic_df = pd.concat([aic_df, tmp_df], axis=0)

# データフレームの最終化と結果表示(最小BICをハイライト)
aic_df = aic_df.reset_index(drop=True)  # インデックスのリセット
aic_df.style.highlight_min(subset='BIC', color='lightpink').format(precision=1)

【実行結果】
BIC 最小は $${\mathtt{x:f}}$$ モデルの 279.4 でした。

まとめ

今回は交互作用項を含むロジスティック回帰のモデリングを実践しました。

🔷 確率分布と確率質量関数

🔷 リンク関数と線形予測子

🔷 オッズ $${\cfrac{q_i}{1-q_i}}$$

🔷 statsmodels のロジスティック回帰モデル構築と結果表示

family = sm.families.Binomial()   # GLMの引数familyに与える確率分布＝二項分布
result = smf.glm(formula='y + I(N-y) ~ x * f', data=data, family=family).fit()
result.summary()

🔷 テキストからのメッセージ
テキストは「むやみに交互作用項をいれない」ことを勧めています。
交互作用項を多数含んだ統計モデルが 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の教科書です。
よかったらぜひ、お試しくださいませ。

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