「データ解析のための統計モデリング入門」をPythonで写経 Vol.3 ～ ３章「一般化線形モデル（GLM）」①ポアソン回帰（体サイズモデル）

３章「一般化線形モデル（GLM）」

書籍の著者　久保拓弥 先生

書籍「データ解析のための統計モデリング入門」３章「一般化線形モデル（GLM）－ポアソン回帰－」の Python写経活動記録 です。 

この記事は 一般化線形モデル（GLM）の一種「ポアソン回帰」を実践 します。
１つの統計モデルに集中して、モデリング手続きを丁寧に辿っていきます。

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

はじめに

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

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

ポアソン回帰（体サイズモデル）

準備

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

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

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

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

# インポート

# 数値計算
import numpy as np
import pandas as pd

# 統計計算
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

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

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

■ 統計モデル
ポアソン回帰と呼ばれる統計モデルです。
平均パラメータ $${\lambda_i}$$ は個体 $${i}$$ ごとに異なります。

■ モデリング手続き

データの確認

データを読み込み、データの外観を眺めてから、統計モデルを選択するためのデータの特徴確認を行います。

■ CSVファイルの読み込み
data3a.csv ファイルを pandas データフレームの data に読み込みます。

# データの読み込み p.42
data = pd.read_csv('./data/ch03/data3a.csv')
print('data.shape: ', data.shape)
data.head()

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

【変数の説明】
データの各変数は、個体から採れた種子数、個体の体サイズ、肥料を与えたかどうかの情報です。

■ データの確認
テキスト p.42～45 のデータ要約・可視化を通じて、データの特徴を確認します。
今回もデータ操作に関する説明が厚い（熱い）です。

① data の列ごとにデータを表示
values 属性で numpy 配列に変換して表示します。

列 x（体サイズ）

# dataの列ごとにデータを表示 p.42, 43

# x列（体サイズ）
print(data['x'].values)

【実行結果】
10 前後の値が多いようです！

列 y（種子数）

# y列（種子数）
print(data['y'].values)

【実行結果】
正の整数のみです。１桁の数が多い印象でしょうか。

列 f（施肥処理）

# f列（施肥処理）　※因子型ではありません
print(data['f'].values) 

【実行結果】
文字列 C と T のみです。因子型ではありません。
前半が C（施肥なし）、後半が T（施肥あり）のようです。

② データの型の表示
info メソッドでpandas データフレームの各変数の型を確認できます。

# 列ごとの型を確認 p.43, 44
data.info()

【実行結果】
y は64ビット整数型、x は64ビット浮動小数点型、f はオブジェクト型（文字など）です。

③ 要約統計量の表示
describe メソッドで要約統計量を表示します。
数値変数のみ表示可能です。

# データフレームの概要を調べる 数値型：要約統計量の表示 p.44
data.describe().round(3)

【実行結果】
最小値 min と最大値 max を確認すると、種子数 y は 2 ～ 15、体サイズ x は およそ 1 ～ 12.5 であることが分かります。

標本分散を調べます。
var メソッドを利用します。

# 分散を調べる
data.var(ddof=1, numeric_only=True).to_frame(name='分散').round(3)

【実行結果】
種子数 y の標本分散 $${6.890}$$ は標本平均 $${7.830}$$ と近いです。

オブジェクト型の f （施肥処理）は変数の要素、つまり C と T の度数を計算します。
value_counts メソッドを利用します。

# データフレームの概要を調べる 文字列型：要素ごとの個数のカウント p.44
data['f'].value_counts().to_frame()

【実行結果】
C と T は 50 個づつ含まれています。

④ 散布図の描画
可視化ライブラリ seaborn を利用して散布図を描画します。
x 軸：体サイズ x、y 軸：種子数 y、データ点の色分け（hue）：施肥処理 f を指定します。

# 散布図の描画 p.45 図3.2

# 散布図の描画（hueで施肥処理有無を設定）
sns.scatterplot(data=data, x='x', y='y', hue='f', 
                palette=['tab:blue', 'tab:red'], s=70, alpha=0.7)
# 修飾：x軸・y軸ラベル、凡例
plt.gca().set(xlabel='x: 体サイズ', ylabel='y: 種子数')
plt.legend(loc='upper left', title='施肥処理');

【実行結果】
体サイズと種子数の相関は・・・無い？うっすら有る？

体サイズと種子数の相関係数を計算します。

# 種子数と体サイズの相関係数
data.corr(numeric_only=True).round(3)

【実行結果】
弱い正の相関があります。

⑤ 箱ひげ図の描画
seaborn で箱ひげ図を描画します。
分布が見えるようにスウォームプロットを重ねます。
x 軸：施肥処理 f、y 軸：種子数 y を指定します。

# 箱ひげ図の描画 p.45 図3.3

# 箱ひげ図の描画
sns.boxplot(data=data, x='f', y='y', hue='f', fill=False)
# スウォームプロットの重ね描き
sns.swarmplot(data=data, x='f', y='y', hue='f')
# 修飾：x軸・y軸ラベル
plt.gca().set(xlabel='f: 施肥処理の有無', ylabel='y: 種子数');

【実行結果】
施肥なし C と施肥あり T の種子数はほとんど変わらない感じです。

⑥ 念の為、種子数のヒストグラムを seaborn で描画します。

# 種子数のヒストグラム
sns.histplot(data=data, x='y', edgecolor='white');

【実行結果】
峰が１つです。左半分に偏っている風に見えます。

◆ ◆ ◆

■ データの特徴まとめ
２章と同じようにデータの特徴を整理します。

種子数のばらつきを「確率分布」で表現します。
種子数の特徴を表現できる確率分布は「ポアソン分布」です。

今回の統計モデルには、説明変数として体サイズ x を使っていきます。
説明変数を使うことで、個体ごとに平均種子数 $${\lambda_i}$$ が異なる統計モデルを構築することになります！

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

今回はポアソン回帰の GLM を観測データに当てはめます。

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

🔷 確率分布と確率質量関数
種子数 $${y_i}$$ はパラメータ $${\lambda_i}$$ のポアソン分布に従います。
パラメータ $${\lambda_i}$$ は個体 $${i}$$ ごとに異なります。

🔷 リンク関数と線形予測子
リンク関数は「対数」、線形予測子は「$${\beta_1 + \beta_2 x_i}$$」です。
パラメータ $${\lambda_i}$$ は指数関数「$${\exp(線形予測子)}$$」で表現されます。

$${\log  \lambda_i = \beta_1 + \beta_2 x_i}$$ の形式に注目します。
（$${\lambda_i}$$ の関数）＝（線形の式）になっています。
左辺の関数を「リンク関数」、右辺の線形の式を「線形予測子」と呼びます。

◆ ◆ ◆

■ 対数リンク関数とポアソン分布
リンク関数が対数 $${\log}$$ の場合「対数リンク関数」と呼びます。
確率分布がポアソン分布のときはたいてい、対数リンク関数を用いるようです。
テキストによると、ポアソン回帰の GLM で対数リンク関数を用いる理由は、$${\lambda_i = \exp(線形予測子)}$$ が必ず $${0}$$ 以上になるので、ポアソン分布のパラメータ $${\lambda_i}$$ を非負とする条件を満たせるからです。
GLM を扱うライブラリでは、確率分布にポアソン分布を指定すると、デフォルトで対数リンク関数を使うようになっています。

◆ ◆ ◆

■ 線形予測子と対数リンク関数の可視化
次の２つの $${\lambda_i}$$ を可視化します。
テキスト p.48 図 3.4 に相当するチャートを描画します。

# 個体iの平均種子数λiと体サイズx_iの関係 p.48 図3.4

## 設定
beta1s = [-2, -1]                 # 切片β1
beta2s = [-0.8, 0.4]              # 傾きβ2
colors = ['tab:blue', 'tab:red']  # グラフの線の色
linestyles = ['--', '-']          # グラフの線種
x_val = np.linspace(-4, 5, 101)   # x_iの値

## 描画
# ２組のパラメータごとにlam算出とグラフ描画を繰り返し処理
for beta1, beta2, color, linestyle in zip(beta1s, beta2s, colors, linestyles):
    # lam算出
    lam_val = np.exp(beta1 + beta2 * x_val)
    # 折れ線グラフの描画
    plt.plot(x_val, lam_val, c=color, ls=linestyle,
             label=rf'{{$\beta_1, \beta2$}}={{{beta1}, {beta2}}}')
# x_i=0の垂直線の描画
plt.axvline(0, color='black', ls=':')
# 修飾：y軸の範囲、x軸ラベル、y軸ラベル、凡例
plt.ylim(0, 2.8)
plt.xlabel('個体 $i$ の体サイズ $x_i$', fontsize=14)
plt.ylabel('個体 $i$ の平均種子数 $\lambda_i$', fontsize=14)
plt.legend(loc='upper center');

【実行結果】

式１は x が大きくなるにつれて線形予測子の値が小さくなるので$${\lambda_i}$$ は「指数関数的減衰（減少）」（青点線）になります。
式２は  x が大きくなるにつれて線形予測子の値が大きくなるので$${\lambda_i}$$ は「指数関数的成長（増加）」（赤実線）になります。

$${\exp(線形予測子)}$$ は $${e^{\ 線形予測子}}$$ です。
$${e = 2.71828 \cdots}$$ はネイピア数、オイラー数、自然対数の底などと呼ばれる定数です。
線形予測子の値によって、ネイピア数 $${e}$$ の「線形予測子」乗が変化します。
このチャートは「指数関数的な変化」を示しています。

線形予測子のパラメータ $${\beta_1, \beta_2}$$ の推定値が気になりますね！

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

まずテキスト 2.2 節の「データとポアソン分布の対応」を実践します。
次にテキスト 2.4 節の「最尤推定」を実践します。

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

🔷 今回の統計モデルの対数尤度

◆ ◆ ◆

■ Python ライブラリで統計モデルをデータに当てはめ
statsmodels の glm を利用して統計モデルを実装します。
引数 family で確率分布、link でリンク関数、formula で線形予測子を指定します。
当てはめ結果を変数 result に格納します。

# 統計モデルの当てはめ ※statsmodelsのglmを利用 p.49, 50

# 確率分布とリンク関数の設定
# ※ポアソン分布の場合、リンク関数のデフォルトは対数 Log（省略可能）
family = sm.families.Poisson(link=sm.families.links.Log())

# GLMの実行 ※formulaに線形予測子を指定
result = smf.glm(formula='y ~ x', data=data, family=family).fit()

# GLMの結果表示
result.summary()

【実行結果】
最下２行の Intercept が $${\beta_1}$$、x が $${\beta_2}$$ に対応しており、テキスト p.50 下部の summary 関数の実行結果 Coefficients に相当します。

【コード：３要素の設定】
🔷モデル表現
「smf.glm」でGLMモデルを指定します

smf.glm(formula=formula, family=family, ...)

🔷 確率分布
smf.glm の引数 family に「sm.families.確率分布名」と指定します。

family = sm.families.Poisson(...)

🔷 リンク関数
sm.families.Poisson() の引数 link で与えます。
ポアソン分布の場合、デフォルトで対数リンク関数を用いるので、引数 link を省略できます。

link=sm.families.links.Log()

🔷 線形予測子
smf.glm の引数 formula に、目的変数とともに、線形予測子を指定します。
formulaの「y ~ x」は「目的変数 y $${\sim}$$ 切片 $${+}$$ 係数 $${\times}$$ 説明変数 x 」を指定しています。

formula='y ~ x'

【コード：当てはめ】
glm() の後ろに fit メソッドを付けることでモデルをデータに当てはめ（フィット）します。

result = smf.glm(formula='y ~ x', data=data, family=family).fit()

上述のコードではモデルと当てはめを一文で書きましたが、モデルと当てはめを次のように別々に書くこともできます。

model = smf.glm(formula='y ~ x', data=data, family=family)
result = model.fit()

【コード：当てはめ結果】
当てはめ結果を格納した変数 result 対して、属性 summary() とすることで、「実行結果」の当てはめ結果表を表示できます。

result.summary()

◆ ◆ ◆

■ 当てはめ結果の分析
当てはめ結果の出力を再掲します。
テキスト p.51 ～ p.53 に現れる推定値・統計量を確認します。
上段の対数尤度と下段の係数の推定値関連を探ります。

🔷 係数の推定値
Intercept と x の coef に注目します。
パラメータである係数の最尤推定値は 切片（Intercept）$${\beta_1 = 1.2917}$$、傾き（x）$${\beta_2 = 0.0757}$$ です。
この推定値を線形予測子に当てはめてみます。

当てはめ結果 result から取り出すこともできます。

# 係数の推定値
print('beta_1:', result.params.Intercept)
print('beta_2:', result.params.x)

【実行結果】

🔷 標準誤差の推定値
続いて標準誤差 std err です。
パラメータ $${\beta_1, \beta_2}$$ の推定値の標準偏差であり、それぞれ $${0.364, 0.036}$$ です。
パラメータの標準誤差の推定値に関するテキストの説明を引用いたします。

当てはめ結果 result から取り出しましょう。

# 標準誤差の推定値
print('beta_1:', result.bse.Intercept)
print('beta_2:', result.bse.x)

【実行結果】

⏰️ アディショナルタイム：標準誤差を計算してみる
ChatGPTに教えてもらった標準誤差の計算方法をPythonで実装してみます。
「観測情報行列 $${J}$$」に基づくものであり、ポアソン回帰のケースにフォーカスしています。

# 標準誤差の算出（観測情報行列 J で） ※ChatGPTに訊いた計算方法に基づき実装

# 説明変数の行列を作成：1の列を追加
vars = np.column_stack([np.ones(len(data)), data.x.values])

# 係数ベクトルの推定値の取得
betas = result.params

# パラメータλiの推定値の算出
lams = np.exp(vars @ betas)

# 観測情報行列 J の算出
J = sum([lams[i] * np.outer(vars[i], vars[i]) for i in range(len(vars))])

# 係数ベクトルの推定値の分散共分散行列の算出
V_hat = np.linalg.inv(J)

# 係数ベクトルの標準誤差の推定値の算出
beta1_se, beta2_se = np.sqrt(np.diag(V_hat))

# 結果の表示
print('beta_1の標準誤差:', beta1_se)
print('beta_2の標準誤差:', beta2_se)

【実行結果】
小数点第８位くらいまでは statsmodels の標準誤差と同じになりました。

観測情報行列のことや、文脈上に現れるフィッシャーの情報行列 $${I}$$ については、ぜひご自身に合った情報を探してみて下さい🍀
（やさしい解説記事を見つけられませんでした…）

◆ ◆ ◆

🔷 $${z}$$ 値
$${z}$$ 値はパラメータの推定値を標準誤差で割って求めた統計量です。
パラメータ $${\beta_1, \beta_2}$$ の $${z}$$ 値はそれぞれ $${3.552, 2.125}$$ です。
テキストによると「$${z}$$ 値はパラメータ推定値がゼロから十分に離れているかどうかの粗い目安になります」とのこと。
当てはめ結果 result から取り出しましょう。

# z値 ※属性名はtvaluesですが、summary()のzを取得できます
print('beta_1:', result.tvalues.Intercept)
print('beta_2:', result.tvalues.x)

【実行結果】

⏰️ アディショナルタイム：Wald 統計量を調べる
テキストは「$${z}$$ 値は Wald統計量 とも呼ばれる」としています。
Web の情報や ChatGPT との会話で調べてみると「$${z}$$ 値と Wald 統計量は別物」との文脈がありました。
ちょいと Wald 統計量を調べたので、メモ的に記録を残します（調べたことはちょっぴりです）。

$${j}$$ 番目の係数の Wald 統計量 $${W_j}$$ は、係数の推定値 $${\hat{\beta_j}}$$ を係数の標準誤差の推定値 $${\text{SE}(\hat{\beta})}$$ で割った値の二乗値です。
つまり、$${z}$$ 値の二乗値です。
$${W_j}$$ は自由度 $${1}$$ のカイ二乗分布に従います。

説明変数 x の係数 $${\beta_2}$$ について、Wald 統計量とその $${p}$$ 値を計算してみます。

# xの係数の推定値 beta2 のWald統計量とそのp値

# Wald 統計量 W の算出
wald_stats = (result.params.x / result.bse.x)**2
print('Wald統計量           :', wald_stats)
print('Wald統計量の正の平方根:', np.sqrt(wald_stats))

# Wald統計量 W のp値の算出
p_value = stats.chi2.sf(x=wald_stats, df=1)
print('Wald統計量のp値      :', p_value)

【実行結果】

２行目の値は $${\beta_2}$$ の $${z}$$ 値と一致しています。
３行目の $${p}$$ 値は次に出現する $${z}$$ 値の $${p}$$ 値と一致します。

◆ ◆ ◆

🔷 P>|z|（$${p}$$ 値）
P>|z| は、statsmodels の公式サイトによると「$${z}$$ 値の $${p}$$ 値」です。
テキストはこの項目について、今回のケースに限定して言えば「平均が $${z}$$ 値の絶対値であり標準偏差が $${1}$$ の正規分布におけるマイナス無限大からゼロまでの値をとる確率の２倍です」と説明しています。

なお、テキストは $${p}$$ 値と呼んでいません（テキストはアンチ $${p}$$ 値の立場です）。
この記事は $${p}$$ 値と呼ぶことにします。

ひとまず、当てはめ結果 result から取り出しましょう。

# p値 
print('beta_1:', result.pvalues.Intercept)
print('beta_2:', result.pvalues.x)

【実行結果】

$${\beta_2}$$ の $${p}$$ 値は、先の Wald 統計量の $${p}$$ 値と一致しています（ほぼほぼ）。

◆ ◆ ◆

🔷 「パラメータの推定値のばらつき」の可視化
テキスト p.52 図 3.6 に相当するチャートを描画します。
このチャートは、平均がパラメータの推定値、標準偏差がパラメータの標準誤差の推定値の「正規分布の確率密度関数」のプロットです。
なおテキストは「パラメータの最尤推定値のばらつきが正規分布で近似できると仮定すると…」と、前置きしています。

### パラメータ推定値のばらつきの評価 p.52 図3.6

## 設定
x_val = np.linspace(-0.1, 1.6, 1001)  # x軸の値
x_val_0 = np.linspace(-0.1, 0, 101)   # x軸の0未満の値（塗りつぶし用）

## 正規分布の確率密度関数の算出 ※bse: パラメータ推定値の標準誤差
# β1
beta1_y_val = stats.norm.pdf(
    x_val, loc=result.params.Intercept, scale=result.bse.Intercept)
# β2
beta2_y_val = stats.norm.pdf(x_val, loc=result.params.x, scale=result.bse.x)
# β2（塗りつぶし用）
beta2_y_val_0 = stats.norm.pdf(x_val_0, loc=result.params.x, scale=result.bse.x)

## 描画
# β1の確率密度関数の描画
plt.plot(x_val, beta1_y_val, color='tab:blue', label=r'$\hat{\beta}_1$')
# β2の確率密度関数の描画
plt.plot(x_val, beta2_y_val, color='tab:red', label=r'$\hat{\beta}_2$')
# β2のx<=0の領域の塗りつぶし描画
plt.fill_between(x_val_0, beta2_y_val_0, 0, color='tab:red', alpha=0.4)
# x=0の垂直線の描画
plt.axvline(0, color='black', lw=0.8)
# 修飾
plt.xlabel('パラメータの推定値')
plt.ylabel('確率密度')
plt.legend();

【実行結果】

テキストに掲載の「$${p}$$ 値が大きいほど $${z}$$ 値がゼロに近くなり、パラメータ推定値がゼロに近いことを表現できる」を頭に置いてチャートを眺めます。

$${\beta_2}$$ の推定値（赤い線）の「ゼロ以下の面積」（赤塗り領域）の２倍が  $${p}$$ 値 $${0.03358}$$ に該当します。
$${p}$$ 値が「そこそこ小さいので $${z}$$ 値がゼロから遠ざかり、パラメータ推定値がゼロではなさそうだ」をイメージします。
$${\beta_1}$$ の推定値（青い線）の「ゼロ以下の面積」の２倍＝ $${p}$$ 値は $${0.00038}$$ です。
$${p}$$ 値が「めっちゃ小さいので $${z}$$ 値がゼロからすごく遠ざかり、パラメータ推定値がゼロではなさそうだ」をイメージします。

テキストが推さない $${p}$$ 値系の文脈だと、こんな感じでしょうか。

最後に$${p}$$ 値を正規分布の確率計算で求めてみます。
$${\beta_1}$$ から。

# β1の推定値のp値 ※上記の青い正規分布でx=0のときの下側確率 × 2
stats.norm.cdf(x=0, loc=result.params.Intercept, scale=result.bse.Intercept) * 2

【実行結果】

続いて $${\beta_2}$$ です。

# β2の推定値のp値 ※上記の赤い正規分布でx=0のときの下側確率 × 2 ※赤塗り部分の確率
stats.norm.cdf(x=0, loc=result.params.x, scale=result.bse.x) * 2

【実行結果】

◆ ◆ ◆

🔷 パラメータ推定値の 95% 信頼区間
statsmodels の GLM の結果にはパラメータ推定値の 95% 信頼区間が表示されています。
95% 信頼区間は、区間が 95% の確率でパラメータの真値を含む、そんな区間です。
パラメータ $${\beta_1, \beta_2}$$ の 95% 信頼区間はそれぞれ $${[0.579\ 2.005],\ [0.006,\ 0.145]}$$ です。

$${p}$$ 値系の文脈だとこんな感じでしょうか。

当てはめ結果 result から 95% 信頼区間を取り出しましょう。

# 係数の95%信頼区間の計算

# GLMの結果から取得
result.conf_int().set_axis(['2.5%', '97.5%'], axis=1)

【実行結果】

最後に、標準正規分布の上側 $${2.5\%}$$ 点（$${\approx 1.96}$$）を用いてパラメータ推定値の 95% 信頼区間を計算してみます。
$${\beta_1}$$ から。

# 標準正規分布の上側2.5%点の取得
ppf = stats.norm.isf(q=0.025, loc=0, scale=1)

# パラメータβ1の95%信頼区間の算出
(result.params.Intercept - ppf * result.bse.Intercept,
 result.params.Intercept + ppf * result.bse.Intercept)

【実行結果】

続いて $${\beta_2}$$ です。

# パラメータβ2の95%信頼区間の算出
(result.params.x - ppf * result.bse.x,
 result.params.x + ppf * result.bse.x)

【実行結果】

◆ ◆ ◆

🔷 最大対数尤度
大トリは、対数尤度 Log-Likelihood です。値は $${-235.39}$$ です。
最尤推定＝対数尤度が最大となるパラメータを推定していますので、この対数尤度は最大対数尤度です。
ひとまず、当てはめ結果 result から取り出しましょう。

# 最大対数尤度の表示
print('最大対数尤度:', result.llf)

【実行結果】

R の logLik 関数で併記される自由度も result から取り出します。
実際には推定したパラメータ数を表示しています。

# 自由度（パラメータ数）
print('パラメータ数:', result.params.size)

【実行結果】

テキストは最大対数尤度について

としています。

モデルのデータに対する当てはまりの文脈で語っているのだと思います。
テキストは、観測データへの「当てはまりの良さ」を重視せず、代わりに未知データへの「予測の良さ」を重視しています。
テキストは予測の良さの指標に AIC を用います。
ということで、result からこのモデルの AIC を取り出しましょう。

# AIC ※-2 * (最大対数尤度 - パラメータ数)
print('AIC:', result.aic)

【実行結果】

パラメータ数の表示は AIC への布石なのかも…、ですね。

◆ ◆ ◆

🔷 statsmodels の GLM の結果
「result から取り出す」流れでさまざまな推定値・統計量を見てきました。
GLMの結果から、もっともっと多くの情報を引き出すことができます。
引き出せる情報の詳細は statsmodels のGLM結果クラス GLMResults の公式サイトをご覧ください！

🔷 statmodels の GLM あるあるメモ
GLM の結果サマリー上部に掲載された「Method: IRLS」は、最尤推定を数値に解く手法として IRLS を使用していることを示しています。
ChatGPT によると、IRLS（Iteratively Reweighted Least Squares：反復再重み付け最小二乗法）とは、

とのこと📝

⏰️ 当てはめのアディショナルタイム

ここまでは statsmodels の GLM クラスを利用してポアソン回帰のモデルを構築しました。
GLM クラスはポアソン回帰に加えて、さまざまなGLMのモデルを構築できます。
汎用性の高い GLM に加えて、statsmodels にはポアソン回帰専用の Poisson クラスが準備されています。

🔷 statsmodels の Poisson クラス

# （参考）統計モデルの当てはめ ※statsmodelsのpoisson（ポアソン回帰）を利用 p.49, 50

# ポアソン回帰の実行
result2 = smf.poisson(formula='y ~ x', data=data).fit()
# ポアソン回帰の結果表示
result2.summary()

【実行結果】
対数尤度およびパラメータ関連の推定値・統計量は GLM と同じ値です。

🔷 最適化ライブラリで最尤推定
statsmodels を使わないで最尤推定を行ってみます。
scipy の最小化問題ソルバー mininize で、負の対数尤度が最小となるパラメータを推定します。

# 最適化ライブラリで最尤推定

# 追加インポート
from scipy.optimize import minimize

# 対数尤度関数（目的関数）※負の対数尤度の最小化問題を解く
def llf(params, x, y):
    beta1, beta2 = params
    lam = np.exp(beta1 + beta2 * x)
    return -np.sum(stats.poisson.logpmf(k=y, mu=lam))

# β1, β2の最尤推定値を算出
res_optim = minimize(fun=llf, x0=[0.1, 0.1], args=(data['x'], data['y']))

# 結果の表示
print(f'β1の最尤推定値 = {res_optim.x[0]:9.4f}')
print(f'β2の最尤推定値 = {res_optim.x[1]:9.4f}')
print(f'最大対数尤度   = {-res_optim.fun:9.4f}')

【実行結果】
statsmodels の当てはめ結果と同じになっています。

平均種子数の予測

テキストにならい、平均種子数 $${\lambda}$$ の予測を行って可視化します。
p.54 図 3.7 に相当します。

🔷 $${\lambda}$$ の関数に代入

まずは、テキストの方法「$${\lambda}$$ の関数にパラメータ推定値を代入する」（次式参照）で予測します。

# 平均種子数λiの予測の可視化 p.53, 54 図3.7 λの数式で計算

## 設定
beta1_hat, beta2_hat = result.params                  # β1、β2の推定値
x_val = np.linspace(data.x.min(), data.x.max(), 100)  # x軸の値

## λiの予測値の算出 ※ λi = exp(β1 + β2 * x) <= テキストの数式
lam_hat1 = np.exp(beta1_hat + beta2_hat * x_val)

## 描画
# x,yの観測値の散布図の描画
sns.scatterplot(data=data, x='x', y='y', hue='f', 
                palette=['tab:blue', 'tab:red'], s=70, alpha=0.7)
# λiの予測値の描画
plt.plot(x_val, lam_hat1, color='black', lw=2, label='$\lambda_i$の予測値')
# 凡例
plt.legend(loc='upper left', title='施肥処理');

【実行結果】
黒い線が $${\lambda}$$ の予測値です。
体サイズ x が大きくなるにつれて、緩やかに平均種子数が増えることが分かります。
青い点は施肥なしのデータ点、赤い点は施肥ありのデータ点です。

🔷 GLMの結果から predict する
GLMの結果を格納した変数 result に対して predict メソッドを実行することで、$${\lambda}$$ の予測値を得られます。
新しい説明変数 x_val を与えて予測します。

モデルの当てはめの際、変数名（列名）が設定された pandas データフレームを利用しました。
予測の際も変数名が必要です（変数名が無いとstatsmodelsに怒られます）。
辞書の key に変数名 x を設定することで正常に動くようになります。

# 平均種子数λの予測の可視化 p.53, 54 図3.7 predictで計算

## 設定
beta1_hat, beta2_hat = result.params                  # β1、β2の推定値
x_val = np.linspace(data.x.min(), data.x.max(), 100)  # x軸の値

## λの予測値の算出 ※ result に対して predictメソッドで予測
lam_hat2 = result.predict(dict(x=x_val)).values  # 変数名xを与える目的で辞書化

## 描画
# x,yの観測値の散布図の描画
sns.scatterplot(data=data, x='x', y='y', hue='f', 
                palette=['tab:blue', 'tab:red'], s=70, alpha=0.7)
# λの予測値の描画
plt.plot(x_val, lam_hat2, color='black', lw=2, label='$\lambda$の予測値')
# 凡例
plt.legend(loc='upper left', title='施肥処理');

【実行結果】
先に描画したチャートと同じ結果になっています。

最初の予測値を格納した lam_hat1 と次の予測値を格納した lam_hat2 を比べていましょう。

# λの計算値とget_predictメソッドの予測値は等しい
lam_hat1 - lam_hat2

【実行結果】
両者が一致していることを確認できました。

🔷 予測区間も推定してみる
調子に乗って種子数の 95% 予測区間も可視化します。
ChatGPT と相談しながら「正規分布近似前提の線形化した近似」を目指します。
予測に用いるメソッドは get_prediction です。

summary_frame() で 予測区間算出に必要な追加情報を取得しています。

# 95% 予測区間（の近似値）をあわせて描画

## 予測の実行
# 予測の実行：summaryに実行結果を格納
pred = result.get_prediction(dict(x=x_val))
summary = pred.summary_frame()
# 平均 λ の予測値・平均の予測値の標準誤差（パラメータの不確実性）の取得
mean_pred = summary['mean'].values
se_mean = summary['mean_se'].values

## 平均 λ の予測値の標準誤差の算出
# 平均の予測値の分散の算出 ※正規分布近似前提の線形化した近似です
# (平均の予測値の標準誤差)² + 目的変数の分散（ポアソン分布の分散は予測平均に等しい）
var_pred = se_mean**2 + mean_pred
# 平均の予測値の標準誤差の算出
se_pred = np.sqrt(var_pred)

## 目的変数 y の95%予測区間の算出
ppt = stats.norm.isf(q=0.025)      # 標準正規分布の上側 2.5%点 ≈ 1.96
lower = mean_pred - ppt * se_pred  # 95% 信頼区間の下端
upper = mean_pred + ppt * se_pred  # 95% 信頼区間の上端

## 描画
# x,yの観測値の散布図の描画
sns.scatterplot(data=data, x='x', y='y', hue='f', 
                palette=['tab:blue', 'tab:red'], s=70, alpha=0.7)
# λの予測値の描画
plt.plot(x_val, mean_pred, color='black', lw=2, label='$\lambda$の予測値')
# yの95%予測区間の塗りつぶし描画
plt.fill_between(x_val, lower, upper, color='lightpink', alpha=0.3,
                 label='95%予測区間')
# 凡例
plt.legend(bbox_to_anchor=(1, 1), title='施肥処理');

【実行結果】
薄赤色の塗りつぶし部分が種子数 y の 95% 予測区間です。

summary_frame() の結果を見てみましょう。

# get_prediction()の実行結果の表示：すべて平均に対する推定値
summary

【実行結果】
左から平均、平均の標準誤差、平均の95%信頼区間の下端・上端です。
今回のモデルの場合、平均は $${\lambda}$$ の予測値に相当します。

まとめ

今回はポアソン回帰と呼ばれる GLM を実践しました。

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

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

🔷 平均パラメータ $${\lambda}$$ の予測

🔷 statsmodels のポアソン回帰モデル構築と結果表示

family = sm.families.Poisson(link=sm.families.links.Log())
result = smf.glm(formula='y ~ x', data=data, family=family).fit()
result.summary()

今回のブログは以上です。

次回は、今回の線形予測子を変形したポアソン回帰を実践します。

シリーズの記事

次の記事

前の記事

目次

ブログの紹介

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の教科書です。
よかったらぜひ、お試しくださいませ。

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