データ解析のための統計モデリング入門　3章　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)

# 3章を読むモチベーション
一般化線形モデル（Generalized Linear Model：GLM)の理解を目指す。GLMでは、確率分布の母数が定数とは限らず、説明変数の一次関数（線形予測子）に依存する。線形予測子における、一次関数の傾きや切片がGLMのパラメータとなる。母数の線形予測子に対する依存の仕方は、リンク関数で表現される。母数にリンク関数を作用させたものが、線形予測子と一致すると定式化する。GLMの線形予測子に含まれるパラメータは、最尤推定により求める。

繰り返しになるが、上記の説明をさらに簡潔に述べると、GLMは、確率分布・リンク関数・線形予測子を指定してモデルを定義する。確率分布のパラメータにリンク関数を作用させたものを線形予測子で表現する。線形予測子には、パラメータが含まれており、そのパラメータを最尤推定により決める。

:::note info
確率分布として、ポアソン分布を例にとりGLMを説明する。

確率分布: 
目的変数 $\( y \)$ はポアソン分布に従う。期待値（平均）は $\( \lambda \)$ で、ポアソン分布の確率質量関数（PMF）は次の通り：
```math
P(y = k) = \frac{\lambda^k e^{-\lambda}}{k!}, \quad k = 0, 1, 2, \dots
```

線形予測子:
説明変数 $\( x \)$ とパラメータ $\( \beta_0 \)$（切片）、$\( \beta_1 \)$（傾き）を用いて次のように表される：
```math
\eta = \beta_0 + \beta_1 x
```

リンク関数:
ポアソン分布の場合、リンク関数として対数リンク関数を用いる。リンク関数を使うことで、ポアソン分布の平均 $\( \lambda \)$ を線形予測子 $\( \eta \)$ で表現する：
```math
\log(\lambda) = \eta = \beta_0 + \beta_1 x
```
よって、平均 $\( \lambda \)$ は次のように表される：
```math
\lambda = \exp(\beta_0 + \beta_1 x)
```
:::

このモデルをより深く理解するために、次の問題を考え、Pythonで統計モデリングとデータ解析を行う。

ある架空の植物について、同種の植物を複数用意して、各個体（一つ一つの植物）から種子がいくつ採取できるか、カウントする。この種子数がGLMに従うと仮定して、モデルのパラメータ（説明変数と確率分布の母数の関係性）を求める。種子数は、施肥処理をしたかどうか（f）、個体の大きさ（x）に依存する、つまり、これらが説明変数になると仮定する。また、種子数は、0以上なので確率分布はポアソン分布になると仮定する。このようなモデルを使う解析方法をポアソン回帰と呼ぶ。

## 前準備

```python
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
import numpy as np
from scipy.optimize import minimize
from scipy.stats import poisson

def download_csv_to_dataframe(url):
  """
  CSVがダウンロードできるリンクからCSVをダウンロードし、データフレームに格納する関数

  Args:
    url: CSVがダウンロードできるリンク (httpsで指定)

  Returns:
    pandas DataFrame: ダウンロードしたCSVデータを含むデータフレーム
  """
  try:
    df = pd.read_csv(url)
    return df
  except Exception as e:
    print(f"Error downloading or reading CSV: {e}")
    return None
    
df = download_csv_to_dataframe("https://kuboweb.github.io/-kubo/stat/iwanamibook/fig/poisson/data3a.csv")
```

## データの観察
データの要約を表示します。
```python
df.describe(include='all')
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/378787/6bb53d94-f198-8a18-bc10-b25e839d07b9.png)


データの散布図を描きます。

```python
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/148f6ca4-120a-0173-24a8-f9a1401c68bc.png)


## GLMのパラメータの最尤推定

### 個体の大きさのみを説明変数にする場合

GLMを定義し、最尤推定を行います。
```python
def log_likelihood(params, x, y):
  """
  対数尤度関数
  params: [beta0, beta1] - 線形モデルのパラメータ
  x: 説明変数
  y: 観測値 (ポアソン分布に従う)
  """
  beta0, beta1 = params
  lambda_ = np.exp(beta0 + beta1 * x)
  return -np.sum(poisson.logpmf(y, lambda_))

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

# 最尤推定
initial_params = [0, 0]  # 初期値
result = minimize(log_likelihood, initial_params, args=(x, y))

# 推定されたパラメータ
beta0_hat, beta1_hat = result.x
print("beta0_hat:", beta0_hat)
print("beta1_hat:", beta1_hat)
max_log_likelihood = -result.fun
print("最大対数尤度:", max_log_likelihood)
```

```bash
beta0_hat: 1.2917267436823487
beta1_hat: 0.0756613313668656
最大対数尤度: -235.38625076999512
```

観測値と予測値を重ねたプロットを描きます。
```python
# 推定されたパラメータを使って予測値を計算
lambda_pred = np.exp(beta0_hat + beta1_hat * x)

# 散布図と予測値のプロット
plt.scatter(x, y, label='Observed Data')
plt.plot(x, lambda_pred, color='red', label='Predicted Values')
plt.xlabel('x')
plt.ylabel('y')
plt.title('Scatter Plot with Predicted Values from Poisson GLM')
plt.legend()
plt.show()
```

### 個体の大きさと施肥処理の有無を説明変数にする場合
説明変数を変えて、GLMを定義し、最尤推定を行います。
```python
# 'f' 列をダミー変数に変換
df['f_dummy'] = df['f'].map({'C': 0, 'T': 1})

def log_likelihood_with_dummy(params, x, f, y):
  """
  ダミー変数を含む対数尤度関数
  params: [beta0, beta1, beta2] - 線形モデルのパラメータ
  x: 説明変数
  f: ダミー変数
  y: 観測値 (ポアソン分布に従う)
  """
  beta0, beta1, beta2 = params
  lambda_ = np.exp(beta0 + beta1 * x + beta2 * f)
  return -np.sum(poisson.logpmf(y, lambda_))

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

# 最尤推定
initial_params = [0, 0, 0]  # 初期値
result = minimize(log_likelihood_with_dummy, initial_params, args=(x, f, y))

# 推定されたパラメータ
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("最大対数尤度:", max_log_likelihood)
```

```bash
beta0_hat: 1.2631126290230017
beta1_hat: 0.08007183224693507
beta2_hat: -0.03199895642885918
最大対数尤度: -235.29371924270876
```

### f施肥処理の有無のみ、説明変数に使う場合
これまでの例と同様に、最尤推定を行います。
```python
# 'f' 列をダミー変数に変換
df['f_dummy'] = df['f'].map({'C': 0, 'T': 1})

def log_likelihood_with_dummy(params, f, y):
  """
  ダミー変数のみを含む対数尤度関数
  params: [beta0, beta1] - 線形モデルのパラメータ
  f: ダミー変数
  y: 観測値 (ポアソン分布に従う)
  """
  beta0, beta1 = params
  lambda_ = np.exp(beta0 + beta1 * f)
  return -np.sum(poisson.logpmf(y, lambda_))

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

# 最尤推定
initial_params = [0, 0]  # 初期値
result = minimize(log_likelihood_with_dummy, initial_params, args=(f, y))

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

max_log_likelihood = -result.fun
print("最大対数尤度:", max_log_likelihood)
```

```bash
beta0_hat: 2.0515563353504285
beta1_hat: 0.012771548343621751
最大対数尤度: -237.62725696068685
```