データ解析のための統計モデリング入門 6章 Pythonによるコーディング

本シリーズのモチベーション

データ解析のための統計モデリング入門の内容について理解を深めることを目的に、本書で例示されていた統計モデリング・データ解析に関する処理をPythonでコーディングします（本書ではRでコーディングされていました）。
注：記事は簡潔に書こうと思いますので、本書をご購入・ご一読の上でないと、理解が難しいです。

２章：データ解析のための統計モデリング入門　2章　Pythonによるコーディング
３章：データ解析のための統計モデリング入門　3章　Pythonによるコーディング
６章（本ページ）：データ解析のための統計モデリング入門　6章　Pythonによるコーディング
７章：データ解析のための統計モデリング入門　7章　Pythonによるコーディング

6章を読むモチベーション

3章に続き、一般化線形モデル（Generalized Linear Model：GLM)の理解を深める。

具体的には、以下の問題を考える。

架空植物の個体$i$において、観察した$N_i$個の種子の内、生きている種子（発芽能力がある種子）と死んでいる種子の個数を数え、それぞれ$y_i$個、$N_{i} - y_{i}$個とする。$y_i$は 0 ~ 8までの整数値を取るとする。
上限があるカウントデータには、確率分布として二項分布を利用するのが適している。また、確率分布の母数と線形予測子を結びつけるリンク関数には、二項分布の母数が確率にあたるので、ロジット関数を用いる。そして、モデルの線形予測子には、個体の大きさxi、施肥処理の有無を表す$f_i$（$f_i$=Tは施肥処理あり、$f_i$=Cは施肥処理なし）による項を用いる。説明変数の交互作用の取り扱いもこの問題で学ぶ。

データの分析

以下からデータを取得し、データフレームに格納する。

# データの概要を表示
df.describe(include='all')

import matplotlib.pyplot as plt
import seaborn as sns

# 'f' 列の値に基づいて色分け
sns.scatterplot(x='x', y='y', hue='f', data=df)
plt.xlabel('x')
plt.ylabel('y')
plt.title('Scatter Plot of x and y, colored by f')
plt.show()

最尤推定

線形予測子の形を変えながら、最尤推定を行う。

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import poisson, norm, binom
from scipy.stats import poisson
from scipy.optimize import minimize_scalar
import pandas as pd
import seaborn as sns
from scipy.optimize import minimize
from scipy.stats import binom

# 'f' 列をダミー変数に変換
df['f_dummy'] = df['f'].map({'C': 0, 'T': 1})

# scipy.stats モジュールから binom オブジェクトをインポート
from scipy.stats import binom

def log_likelihood_binom(params, x, f, y, N):
  """
  ロジスティック回帰モデルの対数尤度関数
  params: [beta0, beta1, beta2] - 線形モデルのパラメータ
  x: 説明変数 (個体の大きさ)
  f: ダミー変数 (施肥処理の有無)
  y: 生きている種子の数
  N: 全種子の数
  """
  beta0, beta1, beta2 = params
  logit_p = beta0 + beta1 * x + beta2 * f
  p = 1 / (1 + np.exp(-logit_p))  # ロジスティック関数で確率に変換
  return -np.sum(binom.logpmf(y, N, p))

# dfからx, f_dummy, y, Nのデータを取り出す
x = df['x'].values
f = df['f_dummy'].values
y = df['y'].values
N = df['N'].values

# 最尤推定 (xとfの両方を使う場合)
initial_params = [0, 0, 0]  # 初期値
result = minimize(log_likelihood_binom, initial_params, args=(x, f, y, N))

# 推定されたパラメータ
beta0_hat, beta1_hat, beta2_hat = result.x
print("beta0_hat:", beta0_hat)
print("beta1_hat:", beta1_hat)
print("beta2_hat:", beta2_hat)

max_log_likelihood = -result.fun
print("最大対数尤度 (xとfを使う):", max_log_likelihood)

# 最尤推定 (xのみを使う場合)
initial_params = [0, 0, 0]  # 初期値、ただし3番目のパラメータ (fに対応) は固定
result_x = minimize(log_likelihood_binom, initial_params, args=(x, np.zeros_like(x), y, N), bounds=[(None, None), (None, None), (0, 0)])  # fを0に固定

# 推定されたパラメータ
beta0_hat_x, beta1_hat_x, _ = result_x.x  # 3番目のパラメータは無視
print("beta0_hat (xのみ):", beta0_hat_x)
print("beta1_hat (xのみ):", beta1_hat_x)

max_log_likelihood_x = -result_x.fun
print("最大対数尤度 (xのみ):", max_log_likelihood_x)

# 最尤推定 (fのみを使う場合)
initial_params = [0, 0, 0]  # 初期値、ただし2番目のパラメータ (xに対応) は固定
result_f = minimize(log_likelihood_binom, initial_params, args=(np.zeros_like(f), f, y, N), bounds=[(None, None), (0, 0), (None, None)])  # xを0に固定

# 推定されたパラメータ
beta0_hat_f, _, beta1_hat_f = result_f.x  # 2番目のパラメータは無視
print("beta0_hat (fのみ):", beta0_hat_f)
print("beta1_hat (fのみ):", beta1_hat_f)

max_log_likelihood_f = -result_f.fun
print("最大対数尤度 (fのみ):", max_log_likelihood_f)

beta0_hat: -19.536051473419874
beta1_hat: 1.9524049563994703
beta2_hat: 2.0215042843386195
最大対数尤度 (xとfを使う): -133.1055646426629
beta0_hat (xのみ): -13.778457416129898
beta1_hat (xのみ): 1.4625956218496687
最大対数尤度 (xのみ): -180.17272164191053
beta0_hat (fのみ): 0.34333330747765917
beta1_hat (fのみ): 0.43351283034217425
最大対数尤度 (fのみ): -316.87987672833407

さらに、線形予測子にxとfの交互作用項を加えて、最尤推定を行う。

def log_likelihood_binom_interaction(params, x, f, y, N):
  """
  ロジスティック回帰モデルの対数尤度関数（交互作用項を含む）
  params: [beta0, beta1, beta2, beta3] - 線形モデルのパラメータ
  x: 説明変数 (個体の大きさ)
  f: ダミー変数 (施肥処理の有無)
  y: 生きている種子の数
  N: 全種子の数
  """
  beta0, beta1, beta2, beta3 = params
  logit_p = beta0 + beta1 * x + beta2 * f + beta3 * x * f  # 交互作用項を追加
  p = 1 / (1 + np.exp(-logit_p))  # ロジスティック関数で確率に変換
  return -np.sum(binom.logpmf(y, N, p))

# dfからx, f_dummy, y, Nのデータを取り出す
x = df['x'].values
f = df['f_dummy'].values
y = df['y'].values
N = df['N'].values

# 最尤推定 (交互作用項を含む)
initial_params = [0, 0, 0, 0]  # 初期値
result_interaction = minimize(log_likelihood_binom_interaction, initial_params, args=(x, f, y, N))

# 推定されたパラメータ
beta0_hat_interaction, beta1_hat_interaction, beta2_hat_interaction, beta3_hat_interaction = result_interaction.x
print("beta0_hat (交互作用項を含む):", beta0_hat_interaction)
print("beta1_hat (交互作用項を含む):", beta1_hat_interaction)
print("beta2_hat (交互作用項を含む):", beta2_hat_interaction)
print("beta3_hat (交互作用項を含む):", beta3_hat_interaction)  # 交互作用項のパラメータ

max_log_likelihood_interaction = -result_interaction.fun
print("最大対数尤度 (交互作用項を含む):", max_log_likelihood_interaction)

beta0_hat (交互作用項を含む): -18.523310270871317
beta1_hat (交互作用項を含む): 1.8525078464928062
beta2_hat (交互作用項を含む): -0.06376486498652352
beta3_hat (交互作用項を含む): 0.2163374404018216
最大対数尤度 (交互作用項を含む): -132.80529836302244