本シリーズのモチベーション
データ解析のための統計モデリング入門の内容について理解を深めることを目的に、本書で例示されていた統計モデリング・データ解析に関する処理をPythonでコーディングします（本書ではRでコーディングされていました）。
注：記事は簡潔に書こうと思いますので、本書をご購入・ご一読の上でないと、理解が難しいです。
２章：データ解析のための統計モデリング入門 2章 Pythonによるコーディング
３章（本ページ）：データ解析のための統計モデリング入門 3章 Pythonによるコーディング
６章：データ解析のための統計モデリング入門 6章 Pythonによるコーディング
７章：データ解析のための統計モデリング入門 7章 Pythonによるコーディング
3章を読むモチベーション
一般化線形モデル（Generalized Linear Model：GLM)の理解を目指す。GLMでは、確率分布の母数が定数とは限らず、説明変数の一次関数（線形予測子）に依存する。線形予測子における、一次関数の傾きや切片がGLMのパラメータとなる。母数の線形予測子に対する依存の仕方は、リンク関数で表現される。母数にリンク関数を作用させたものが、線形予測子と一致すると定式化する。GLMの線形予測子に含まれるパラメータは、最尤推定により求める。
繰り返しになるが、上記の説明をさらに簡潔に述べると、GLMは、確率分布・リンク関数・線形予測子を指定してモデルを定義する。確率分布のパラメータにリンク関数を作用させたものを線形予測子で表現する。線形予測子には、パラメータが含まれており、そのパラメータを最尤推定により決める。
確率分布として、ポアソン分布を例にとりGLMを説明する。
確率分布:
目的変数 $( y )$ はポアソン分布に従う。期待値（平均）は $( \lambda )$ で、ポアソン分布の確率質量関数（PMF）は次の通り：
P(y = k) = \frac{\lambda^k e^{-\lambda}}{k!}, \quad k = 0, 1, 2, \dots
線形予測子:
説明変数 $( x )$ とパラメータ $( \beta_0 )$（切片）、$( \beta_1 )$（傾き）を用いて次のように表される：
\eta = \beta_0 + \beta_1 x
リンク関数:
ポアソン分布の場合、リンク関数として対数リンク関数を用いる。リンク関数を使うことで、ポアソン分布の平均 $( \lambda )$ を線形予測子 $( \eta )$ で表現する：
\log(\lambda) = \eta = \beta_0 + \beta_1 x
よって、平均 $( \lambda )$ は次のように表される：
\lambda = \exp(\beta_0 + \beta_1 x)
このモデルをより深く理解するために、次の問題を考え、Pythonで統計モデリングとデータ解析を行う。
ある架空の植物について、同種の植物を複数用意して、各個体（一つ一つの植物）から種子がいくつ採取できるか、カウントする。この種子数がGLMに従うと仮定して、モデルのパラメータ（説明変数と確率分布の母数の関係性）を求める。種子数は、施肥処理をしたかどうか（f）、個体の大きさ（x）に依存する、つまり、これらが説明変数になると仮定する。また、種子数は、0以上なので確率分布はポアソン分布になると仮定する。このようなモデルを使う解析方法をポアソン回帰と呼ぶ。
前準備
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")
データの観察
データの要約を表示します。
df.describe(include='all')
データの散布図を描きます。
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()
GLMのパラメータの最尤推定
個体の大きさのみを説明変数にする場合
GLMを定義し、最尤推定を行います。
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)
beta0_hat: 1.2917267436823487
beta1_hat: 0.0756613313668656
最大対数尤度: -235.38625076999512
観測値と予測値を重ねたプロットを描きます。
# 推定されたパラメータを使って予測値を計算
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を定義し、最尤推定を行います。
# '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)
beta0_hat: 1.2631126290230017
beta1_hat: 0.08007183224693507
beta2_hat: -0.03199895642885918
最大対数尤度: -235.29371924270876
f施肥処理の有無のみ、説明変数に使う場合
これまでの例と同様に、最尤推定を行います。
# '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)
beta0_hat: 2.0515563353504285
beta1_hat: 0.012771548343621751
最大対数尤度: -237.62725696068685