「データ解析のための統計モデリング入門」をPythonで写経 Vol.5 ～ ３章「一般化線形モデル（GLM）」③ポアソン回帰と恒等リンク関数・直線回帰

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

書籍の著者　久保拓弥 先生

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

この記事は前回記事の流れを汲み、確率分布とリンク関数 に注目して、「統計モデル」がデータを生成した「現象」の表現に適する／適さない、といったことの意味合いを考えます。

具体的にはカウントデータに対して、① ポアソン分布＆恒等リンク関数の当てはめ、② 正規分布＆恒等リンク関数の当てはめ、の２つを試します。

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

はじめに

このブログシリーズは、書籍「データ解析のための統計モデリング入門 一般化線形モデル・階層ベイズモデル・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

# Rデータセットの読み込み
import rdata

# 可視化
import matplotlib.pyplot as plt
plt.rcParams['font.family'] = 'Meiryo'  # または import japanize_matplotlib

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

統計モデリングの概要です。

■ 統計モデル
ポアソン回帰と直線回帰（線形回帰）です。

■ モデリング手続き
今回はモデリング自体に注目しないため、モデリング手続きの深堀りを省略します。

その代わり、カウントデータに関する統計モデルの比較検討に注力します。

• ポアソン分布でリンク関数の異なる２モデルの比較

• ポアソン回帰モデルと直線回帰（線形回帰）モデルの比較

ポアソン分布が恒等リンク関数とLINKする…

テキスト p.59 ～ の「ポアソン回帰に恒等リンク関数を使ってみたら…」に取り組みます。

データの読み込み

前回記事と同じデータを利用します。
data3a.csv ファイルを pandas データフレームの data に読み込みます。

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

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

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

統計モデルをデータに当てはめ

ポアソン分布に関する２つのモデルをフィッティングします。

◆ ◆ ◆

１つ目は「対数リンク関数モデル」です。
前回記事と同じポアソン回帰の「体サイズ＋施肥有無モデル」です。
リンク関数は「対数リンク関数」です。

🔷 確率分布と確率質量関数
個体 $${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の当てはめ

# 対数リンク関数を用いた統計モデルの当てはめ ※statsmodelsのglmを利用

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

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

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

【実行結果】
当てはまりの良さの指標「最大対数尤度」は $${-235.29}$$ です。

平均種子数 $${\lambda_i}$$ の予測式は次のようになります。

◆ ◆ ◆

２つ目は「恒等リンク関数モデル」です。
その名の通り「恒等リンク関数」を使うモデルです。
確率分布と線形予測子は１つ前のモデルと同じです。

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

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

恒等リンク関数の場合、左辺のパラメータ $${\lambda_i}$$ は特段の関数形になっていないことに注目しましょう。
テキストの言葉をお借りすると「リンク関数がとくに何もない」です。

🔷 GLMの当てはめ

# 恒等リンク関数を用いた統計モデルの当てはめ ※statsmodelsのglmを利用

# 確率分布の設定（リンク関数は恒等関数Identity）
family_id = sm.families.Poisson(link=sm.families.links.Identity())

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

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

【実行結果】
当てはまりの良さの指標「最大対数尤度」は $${-235.16}$$ です。
実はこちらの恒等リンク関数モデルの方が最大対数尤度が大きいので、「当てはまりの良いモデル」になっています。

ワーニング「ポアソン分布は恒等リンク関数を尊重しない」的な英文が出ますが、気にしないことにします。

平均種子数 $${\lambda_i}$$ の予測式は次のようになります。

２つのリンク関数の違い

２つの統計モデルを比べます。
最初に最大対数尤度を比べます。

最大対数尤度が大きい恒等リンク関数の方が、当てはまりの良いモデルということになります。
しかしテキストは「どちらのリンク関数が妥当なモデルなのかは、当てはまりの良しあしだけで決まる問題ではありません」と警鐘を鳴らしています。
では何を手がかりにしてモデルを見たら良いのでしょう？？？

◆

ひとまず２つのモデルの予測値を可視化します。
テキスト p.60 図 3.8 に相当します。

（注意点）
テキストは図をわかりやすくする目的で「施肥処理の効果を３倍に設定」（$${\beta_3 \times 3}$$ の意味）しているとしています。
ところが、サポートサイトの R プログラムコードは 10 倍しています。
コードを尊重して、こちらのコードも 10 倍しました。

# 対数リンク関数・恒等リンク関数の予測 p.60 図3.8

## 設定と準備
# 対数リンク関数の統計モデルのパラメータ推定値
beta1_log, beta3_log, beta2_log = result_log.params
# 恒等リンク関数の統計モデルのパラメータ推定値
beta1_id, beta3_id, beta2_id = result_id.params
# xの値
x_val = np.linspace(5, 20, 101)

## 関数定義
# 対数リンク関数のλ算出ヘルパー関数 ※施肥処理の効果を10倍に設定
def log_link_func(x, f, multiplier=10):
    f = int(f == 'T')
    return np.exp(beta1_log + beta2_log * x + beta3_log * multiplier * f)
# 恒等リンク関数のλ算出ヘルパー関数 ※施肥処理の効果を10倍に設定
def id_link_func(x, f, multiplier=10):
    f = int(f == 'T')
    return beta1_id + beta2_id * x + beta3_id * multiplier * f

## 描画
# 描画領域の設定
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(8, 4), sharey=True,
                               tight_layout=True)
# 対数リンク関数の描画
ax1.plot(x_val, log_link_func(x_val, 'C'), label='無処理')       # 無処理の描画
ax1.plot(x_val, log_link_func(x_val, 'T'), label='施肥処理')     # 有処理の描画
ax1.set(xlabel='体サイズ $x_i$', ylabel='平均種子数 $\lambda_i$', # 修飾
        title='対数リンク関数') 
ax1.legend(loc='upper left')
# 恒等リンク関数の描画
ax2.plot(x_val, id_link_func(x_val, 'C'), label='無処理')        # 無処理の描画
ax2.plot(x_val, id_link_func(x_val, 'T'), label='施肥処理')      # 有処理の描画
ax2.set(xlabel='体サイズ $x_i$', title='恒等リンク関数')          # 修飾
ax2.legend(loc='upper left');

【実行結果】

【チャートの解釈】
２つのモデルは施肥処理をすると平均種子数が少なくなります。
ただし、少なくなる「幅」が異なっています。
対数リンク関数モデルは、平均種子数が多くなるほど、施肥処理による減少幅が大きくなっています。
具体的には施肥処理によって $${\exp(-0.032)}$$ 倍減少します。
一方で恒等リンク関数モデルは、２本の直線が平行になっていて、平均種子数に関係なく一定幅（$${0.205}$$ 個）が減少します。

◆

テキストは「対数リンク関数の方がましな統計モデルのような気もします」とコメントしています。

確かに、対数リンク関数モデルのような変動の方が「どちらかと言えば現実に合っていそうだ」という気がします。
恒等リンク関数モデルは、左側に直線を伸ばして「マイナスの平均種子数」という現実離れした予測を出しかねません。

現実離れ・・・、そうなんです！

統計モデルはデータを生成した現実・現象を理解・予測する目的で作られるので、現象を表現できる統計モデルを作りたいです！

テキストは「数式が現象をどのように表現しているかに注意して統計モデルを設計しよう」（p.60）と促しています。
恒等リンク関数モデルは、観測データを生み出した現象を表現できていないだろう、と考えられます。

◆

テキストは、統計モデルがデータを表現できているかという視点に加えて、対数リンク関数の良いところを次のように語っています。

ポアソン回帰と直線回帰を比べる

テキスト p.61 ～ の「カウントデータに直線回帰モデルを当てはめることに関する考察」に取り組みます。

データの読み込み

新しい例題データを利用します。
d0.RData ファイルを rdata ライブラリで変換して、pandas データフレームの data2 を作成します。

# データの読み込み

# RdataをPythonオブジェクトに変換
converted = rdata.read_rda('./data/ch02/d0.RData', default_encoding='ASCII')

# pandasデータフレーム化
data2 = pd.DataFrame(converted['d0']).astype({'y': 'int'})

# csvファイルの出力
# data2.to_csv('./data/ch02/data02_d0.csv')

# データフレームの表示
print('data2.shape: ', data2.shape)
data2.head()

【実行結果】
データの個数（標本サイズ）が 30 の仮想の観測データです。
変数 y は整数です。

【変数の説明】
テキストでは明らかにされていません。
x が説明変数、y が目的変数に用いられます。

観測データの基本統計量を確認しましょう。

# 基本統計量の表示
data2.describe().round(1)

【実行結果】

【データプロファイル】

0 以上の整数である目的変数 y のバラツキを表現する確率分布は、ポアソン分布が候補になりそうです。

統計モデルをデータに当てはめ

観測データに２つのモデルをフィッティングします。

◆ ◆ ◆

１つ目は「直線回帰モデル」です。
いわゆる線形回帰モデル・単回帰モデルです。

テキストは線形モデル（linear model. LM）あるいは一般線形モデル（general linear model）の中でも特によく使われるモデルの１つである「直線回帰（linear regression）」と紹介しています。

🔷 確率分布と確率密度関数
$${y_i}$$ はパラメータ $${\mu_i, \sigma^2}$$ の正規分布に従います。

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

🔷 線形回帰モデルの当てはめ
statsmodels の ols を利用します。

# 統計モデルの当てはめ ※statsmodelsのolsを利用(正規分布・恒等リンク関数のGLMと同等)

# 線形回帰モデルLMの当てはめ
result_norm = smf.ols(formula='y ~ x', data=data2).fit()

# LMの結果表示
result_norm.summary()

【実行結果】
当てはまりの良さの指標「最大対数尤度」は $${-53.497}$$ です。

$${y_i}$$ の予測式は次のようになります。

◆ ◆ ◆

２つ目は「ポアソン回帰モデル」です。

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

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

🔷 GLMの当てはめ

# 統計モデルの当てはめ ※statsmodelsのglmを利用 ポアソン分布・対数リンク関数

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

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

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

【実行結果】
当てはまりの良さの指標「最大対数尤度」は $${-29.093}$$ です。
このポアソン回帰モデルの方が「当てはまりが良い」ことになります。

$${\lambda_i}$$ の予測式は次のようになります。

２つのモデルを比べる

観測データと統計モデルによる予測値を可視化して、統計モデルのデータへの当てはまりを確認します。
テキスト p.61 図 3.9 に相当します。

### 回帰モデルと確率分布の関係 p.61 図3.9

### 設定
x_val = np.linspace(data2['x'].min(), data2['x'].max(), 101)  # x軸の値
sigma = result_norm.resid.std()                             # 正規分布の標準偏差
prob_vals = [0.5, 1.1, 1.7]                                 # 確率分布の垂直線

### 描画
# 描画領域の設定
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(10, 4))

### 左の領域
# 観測値の散布図の描画
ax1.scatter(data2['x'], data2['y'], s=70, ec='tab:blue', fc='white')
# 回帰直線の描画
ax1.plot(x_val, result_norm.predict(dict(x=x_val)), color='tab:red', ls='--',
         lw=2)
# y=0の水平線の描画
ax1.axhline(0, color='black', lw=0.5)
# 修飾
ax1.set(ylim=(-2, 8), xlabel='$x$', ylabel='$y$',
        title='正規分布・恒等リンク関数の統計モデル')

## 正規分布の確率密度関数の描画
# 確率変数の設定
ymin, ymax = -1.9, 7
yy = np.linspace(ymin - 3, ymax + 3, 100)
# 正規分布の確率密度関数の描画
for prob_val in prob_vals:
    # 予測値の平均値の算出
    mu = result_norm.predict(dict(x=prob_val))
    # 確率密度関数の算出
    xx = stats.norm.pdf(yy, loc=mu, scale=sigma)
    # 確率密度関数の描画
    ax1.plot(prob_val - xx * 0.5, yy, lw=0.7, color='tab:green', alpha=0.3)
    # 確率密度関数の0を示す垂直線の描画
    ax1.axvline(prob_val, lw=0.7, color='tab:green', alpha=0.3)
    # 確率密度関数内部の塗りつぶし
    ax1.fill_betweenx(yy, prob_val, prob_val - xx * 0.5, color='tab:green',
                      alpha=0.3)

### 右の領域
# 観測値の散布図の描画
ax2.scatter(data2['x'], data2['y'], s=70, ec='tab:blue', fc='white')
# λの予測値（指数関数）の描画
ax2.plot(x_val, result_pois.predict(dict(x=x_val)), color='tab:red', ls='--',
         lw=2)
# y=0の水平線の描画
ax2.axhline(0, color='black', lw=0.5)
# 修飾
ax2.set(ylim=(-2, 8), xlabel='$x$', ylabel='$y$',
        title='ポアソン分布・対数リンク関数の統計モデル')

## ポアソン分布の確率質量関数の描画
# 確率変数の設定
yy = list(range(0, 11))
## ポアソン分布の確率質量関数の描画
for prob_val in prob_vals:
    # λの予測値の算出
    lam = result_pois.predict(dict(x=prob_val))
    # 確率質量関数の算出
    xx = stats.poisson.pmf(yy, mu=lam)
    # 確率質量関数の棒グラフの描画
    ax2.barh(yy, xx * 0.5, left=prob_val - xx * 0.5 - 0.01, height=0.4, 
             color='tab:orange', ec='darkorange', alpha=0.5)
    # 確率質量関数の0を示す垂直線の描画
    ax2.axvline(prob_val, color='darkorange', lw=1.5, alpha=0.5)

plt.show()

【実行結果】

【チャートの見方】

• y の観測値は青丸です。
  y が離散値（整数）になっていることが分かります。

• y の平均の予測値は赤点線です。
  観測値との当てはまりをチェックしましょう。

• 緑とオレンジの分布は、x = $${[0.5, 1.1, 1.7]}$$ のときの y のばらつきのイメージです。観測値・赤点線の乖離具合と重ね合わせてご覧ください。

  • 緑の３つの正規分布は形状が同じであり、分散が同じことを示しています。

  • オレンジの３つのポアソン分布は x が大きくなるにつれて、形状が平たくなっていて、分散が大きくなっていることを示しています。

◆

【チャートの解釈】
y の観測値の特徴は次の３点と思われます。

• 離散値（整数）

• 非負値（０以上）

• ばらつきに傾向あり（値が大きいほどばらつきが大きい）

２つの統計モデルは y の観測値の特徴を適度に表しているか（当てはまっているか）の状況を表にまとめます。

ポアソン回帰の方が観測値の特徴に適しており、データを生成した現象をよりよく表せている感じがします。

直線回帰の適用に関するテキストのコメントを引用いたします。

テキストは「データの特徴を無視して、どんなデータにも直線回帰を当てはめる解析行為をやめよう」と提案しています。

◆

【モデル比較のまとめ】
２モデルに限定した話として、観測データのような離散値・非負値・ばらつき傾向がある場合には、ポアソン分布・対数リンク関数を仮定するポアソン回帰のほうが適していました。
正規分布・恒等リンク関数を仮定する直線回帰モデルは不適でした。

まとめ

今回は統計モデリングにおいて、データの特徴に適する確率分布やリンク関数を選ぶことの大切さを学びました。

テキストが多用するカウントデータ（整数値・非負値）の場合、正規分布や恒等リンク関数が適さない場合があります。
カウントデータの特徴をよく吟味して、ポアソン分布・対数リンク関数がデータをより良く表現できることも同時に学びました。

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

次回は、モデル選択規準のひとつ「AIC」を学びます。

シリーズの記事

次の記事

前の記事

目次

ブログの紹介

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

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