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

# 本シリーズのモチベーション
[データ解析のための統計モデリング入門](https://amzn.asia/d/b9hqFc5)の内容について理解を深めることを目的に、本書で例示されていた統計モデリング・データ解析に関する処理をPythonでコーディングします（本書ではRでコーディングされていました）。
注：記事は簡潔に書こうと思いますので、本書をご購入・ご一読の上でないと、理解が難しいです。

２章：[データ解析のための統計モデリング入門　2章　Pythonによるコーディング](https://qiita.com/warper/items/82a4b7dd747690e7a530)
３章：[データ解析のための統計モデリング入門　3章　Pythonによるコーディング](https://qiita.com/warper/items/02a4436d6fa6919bc1eb)
６章（本ページ）：[データ解析のための統計モデリング入門　6章　Pythonによるコーディング](https://qiita.com/warper/items/f459bb214f5286869e89)
７章：[データ解析のための統計モデリング入門　7章　Pythonによるコーディング](https://qiita.com/warper/items/9b53db1b35f3b6c0bbef)

# 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は施肥処理なし）による項を用いる。説明変数の交互作用の取り扱いもこの問題で学ぶ。

## データの分析

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

https://kuboweb.github.io/-kubo/stat/iwanamibook/fig/binomial/data4a.csv

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

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/378787/1bfc237a-b1e0-c11f-17bb-01b051573bef.png)


```python
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()
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/378787/cc2d3477-980b-0e73-7497-8098ac9acaa8.png)

## 最尤推定

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

```python
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)
```

```bash
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の交互作用項を加えて、最尤推定を行う。

```python
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)
```

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