「データ解析のための統計モデリング入門」をPythonで写経 Vol.4 ～ ３章「一般化線形モデル（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

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

統計モデリングの概要です。
この記事は２つの統計モデルを取り扱います。

■ 統計モデル
ポアソン回帰と呼ばれる統計モデルです。

■ モデリング手続き

データの確認（簡略版）

■ データの読み込み
今回の２つの統計モデルは前回記事と同じデータに当てはめします。
data3a.csv ファイルを pandas データフレームの data に読み込みます。

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

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

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

今回の統計モデリングで活用する質的変数「施肥処理 f 」の概要を掴みます。
前回記事の振り返りにもなっています。

■ データ型の確認
最初にデータ型を確認します。
info メソッドでpandas データフレームの各変数の型を確認できます。

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

【実行結果】
f はオブジェクト型（文字など）です。数値ではありません。

■ 施肥処理の要素ごとの個数の確認
value_counts メソッドを利用して、施肥処理 f の要素「C, T」の個数を調べます。

# 変数 f のユニーク値と度数
data.f.value_counts().to_frame()

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

■ 可視化
施肥処理 f と目的変数である種子数 y の関係を箱ひげ図で確認します。

# 質的変数 f と量的変数 y の箱ひげ図の描画

# 箱ひげ図の描画
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 との間に、目に見える差は見られません。

ポアソン回帰（施肥有無モデル）

テキスト p.54 ～ の「説明変数が因子型の統計モデル」に取り組みます。

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

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

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

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

🔷 リンク関数と線形予測子
リンク関数は「対数」、線形予測子は「$${\beta_1 + \beta_3 d_i}$$」です。

変数 $${d_i}$$ は 施肥処理 $${f_i}$$ のカテゴリ値 C, T を数値 $${0, 1}$$ に変換したダミー変数です。
ダミー変数への変換は、機械学習の文脈では One-Hot エンコーディングで知られています。

なお、statsmodels による GLM の実装時（※）には、明示的にダミー変数化をしなくても大丈夫です。
※この記事で用いる Formula API の場合です。

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

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

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

◆ ◆ ◆

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

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

# 確率分布の設定（リンク関数はデフォルトの対数）
family=sm.families.Poisson()

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

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

【実行結果】
最下２行の Intercept が $${\beta_1}$$、f[T.T] が $${\beta_3}$$ に関連しており、テキスト p.56 上部の summary 関数の実行結果 Coefficients に相当します。

【コードの補足】
前回記事の統計モデルと異なるのは、formula を「y ~ f」とする点です。

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

formula='y ~ f'

【実行結果の補足】
「 f [T.T] 」は、変数 f をダミー変数化（Treatment Coding、最初のT）によって要素「T」（２番目のT）を値 $${1}$$ とし、暗黙で「C」を値 $${0}$$ にしたことを表現しています。

◆ ◆ ◆

■ 当てはめ結果の分析

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

施肥処理 f が C （施肥なし）のときは、$${d_i = 0}$$ となるため

となり、施肥処理 f が T （施肥あり）のときは、$${d_i = 1}$$ となるため

となります。
この統計モデルは、施肥ありの場合には施肥なしの場合よりも平均種子数 $${\lambda_i}$$ が $${0.1}$$ だけ増えることを示しています。

【補足情報】
「 f [T.T] 」の $${z}$$ 値の $${p}$$ 値が $${0.858}$$ であり、大きな値です（例えば $${5\%}$$ をはるかに超える）。
統計的検定の文脈では、施肥処理 f の係数は有意では無い、ということになります。

🔷 最大対数尤度
最大対数尤度 Log-Likelihood は $${-237.63}$$ です。
前回の体サイズモデルの最大対数尤度は $${-235.39}$$ でした。
最大対数尤度の大きいモデルが「データへの当てはまりが良い」ので、今回の施肥有無モデルの方が「当てはまりが悪い」モデルです。

🔷 AIC
AIC を計算します。第４章への布石です。

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

【実行結果】
AIC は $${479.25}$$ です。
前回の体サイズモデルの AIC は $${474.77}$$ でした。
AIC の小さいモデルが「予測が良い」ので、今回の施肥有無モデルの方が「予測が悪い」モデルです。

平均種子数の予測

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

# 平均種子数λの予測の可視化

## 設定
x_val = np.linspace(data.x.min(), data.x.max(), 100)  # x軸の値

## f=C,T別のλの予測値の算出
lam_hat1_C = result1.predict(dict(f=np.repeat('C', 100)))
lam_hat1_T = result1.predict(dict(f=np.repeat('T', 100)))

## 描画
# 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_hat1_C, color='tab:blue', label='Cの場合の$\lambda$の予測値')
plt.plot(x_val, lam_hat1_T, color='tab:red', label='Tの場合の$\lambda$の予測値')
# 凡例
plt.legend(bbox_to_anchor=(1, 1), title='施肥処理');

【実行結果】
y=8 付近の直線を凝視しましょう。
下側に施肥なし C の青い線、上側に施肥あり T の赤い線が見えるはずです！
施肥の影響はほとんど無いように見えます。

ポアソン回帰（体サイズ＋施肥有無モデル）

テキスト p.57 ～ の「説明変数が数量型＋因子型の統計モデル」に取り組みます。

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

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

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

🔷 確率分布と確率質量関数
施肥有無モデルと同じです。
個体 $${i}$$ の種子数 $${y_i}$$ はパラメータ $${\lambda_i}$$ のポアソン分布に従います。

🔷 リンク関数と線形予測子
リンク関数は「対数」、線形予測子は「$${\beta_1 + \beta_2 x_i + \beta_3 d_i}$$」です。

変数 $${x_i}$$ は体サイズです。
変数 $${d_i}$$ は 施肥処理 $${f_i}$$ のカテゴリ値 C, T を数値 $${0, 1}$$ に変換したダミー変数です。

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

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

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

◆ ◆ ◆

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

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

# 確率分布の設定（リンク関数はデフォルトの対数）
family=sm.families.Poisson()

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

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

【実行結果】
最下３行の Intercept が $${\beta_1}$$、f[T.T] が $${\beta_3}$$、x が $${\beta_2}$$ に関連しており、テキスト p.58 上部の summary 関数の実行結果 Coefficients に相当します。

【コードの補足】
施肥有無モデルと異なるのは、formula を「y ~ x + f」とする点です。

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

formula='y ~ x + f'

◆ ◆ ◆

■ 当てはめ結果の分析

🔷 係数の推定値
Intercept、x、 f [T.T] の coef に注目します。
パラメータである係数の最尤推定値は 切片（Intercept）$${\beta_1 = 1.2631}$$、体サイズ（ x ） $${\beta_2 = 0.0801}$$、施肥処理（ f [T.T] ）$${\beta_3 = -0.0320}$$ です。
施肥処理の係数がマイナスになってしまいました…

この推定値を線形予測子に当てはめてみます。

施肥処理 f が C （施肥なし）のときは、$${d_i = 0}$$ となるため

となり、施肥処理 f が T （施肥あり）のときは、$${d_i = 1}$$ となるため

となります。

【係数の推定値の効果】
この統計モデルは、施肥ありの場合には施肥なしの場合よりも平均種子数 $${\lambda_i}$$ が $${\exp(-0.0320) \approx 0.969}$$ 倍になることを示しています。
また、体サイズ $${x_i}$$ が $${1}$$ 単位大きくなると、$${\lambda_i}$$ が $${\exp(0.0801)  \approx 1.083}$$ 倍になることを示しています。

テキストは p.58 ～ で「対数リンク関数では要因の効果が積であらわされる」としています。

Python で効果（倍率）を計算してみます。
$${\lambda}$$ を予測するヘルパー関数を定義します。

# λ予測ヘルパー関数の定義
def pred_λ(x, f, result):
    # 戻り値：λの予測値
    return result.predict(dict(x=x, f=f))[0]

施肥ありの場合の効果を計算します。

# 説明変数fiが施肥ありのときλiは0.969倍になる、の検証 p.59
ans = pred_λ(x=1, f='T', result=result2) / pred_λ(x=1, f='C', result=result2)
print(f'説明変数fiが施肥ありのときλiは {ans:.3f} 倍になる')

【実行結果】

体サイズが１単位大きくなる場合の効果を計算します。

# 説明変数xiが1増加するとλiは1.8倍増える、の検証 p.59
ans = pred_λ(x=8, f='C', result=result2) / pred_λ(x=7, f='C', result=result2)
print(f'説明変数xiが1増加するとλiは {ans:.3f} 倍になる')

【実行結果】

【補足情報】
「 f [T.T] 」の $${z}$$ 値の $${p}$$ 値が $${0.667}$$ であり、大きな値です（例えば $${5\%}$$ をはるかに超える）。
統計的検定の文脈では、施肥処理 f の係数は有意では無い、ということになります。

🔷 最大対数尤度
最大対数尤度 Log-Likelihood は $${-235.29}$$ です。
施肥有無モデルは $${-237.63}$$、前回の体サイズモデルは $${-235.39}$$ でした。
最大対数尤度の大きいモデルが「データへの当てはまりが良い」ので、今回の体サイズ＋施肥有無モデルが「当てはまりが良い」モデルです。

🔷 AIC
AIC を計算します。第４章への布石です。

# AICの表示 p.58
print('AIC:', result2.aic)

【実行結果】
AIC は $${476.59}$$ です。
施肥有無モデルは $${479.25}$$、前回の体サイズモデルは $${474.77}$$ でした。
AIC の小さいモデルが「予測が良い」ので、前回の体サイズモデルが「予測が良い」モデルです。

平均種子数の予測

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

# 平均種子数λの予測の可視化

## 設定
x_val = np.linspace(data.x.min(), data.x.max(), 100)  # x軸の値

## f=C,T別のλの予測値の算出
lam_hat2_C = result2.predict(dict(x=x_val, f=np.repeat('C', 100)))
lam_hat2_T = result2.predict(dict(x=x_val, f=np.repeat('T', 100)))

## 描画
# 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_C, color='tab:blue', label='Cの場合の$\lambda$の予測値')
plt.plot(x_val, lam_hat2_T, color='tab:red', label='Tの場合の$\lambda$の予測値')
# 凡例
plt.legend(bbox_to_anchor=(1, 1), title='施肥処理');

【実行結果】
青・赤の$${\lambda}$$ の予測値の曲線は、なだらかに右上がりになっています。
体サイズ x が大きくなるにつれて、やや平均種子数が大きくなっています。
施肥なし C の青い線のほうが平均種子数が若干大きいことを示しています。

まとめ

今回は２つのポアソン回帰を実践しました。

1️⃣ 施肥有無モデル

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

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

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

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

family=sm.families.Poisson()
result1 = smf.glm(formula='y ~ f', data=data, family=family).fit()
result1.summary()

◆ ◆ ◆

2️⃣ 体サイズ＋施肥有無モデル

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

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

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

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

family=sm.families.Poisson()
result2 = smf.glm(formula='y ~ x + f', data=data, family=family).fit()
result2.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の教科書です。
よかったらぜひ、お試しくださいませ。

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