「データ解析のための統計モデリング入門」をPythonで写経 Vol.11 ～ ６章「GLMの応用範囲を広げる」③ポアソン回帰とオフセット項

６章「GLMの応用範囲を広げる」

書籍の著者　久保拓弥 先生

書籍「データ解析のための統計モデリング入門」６章「GLMの応用範囲を広げる」の Python写経活動記録 です。 

この記事は 一般化線形モデル（GLM）の一種「ポアソン回帰」を実践 します。
今回は オフセット項 を線形予測子に取り入れます。

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

はじめに

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

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

準備と概要

準備

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

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

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

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

# インポート

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

# 統計計算
import statsmodels.api as sm
import statsmodels.formula.api as smf

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

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

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

■ 統計モデル
「オフセット項」を含むポアソン回帰です。

■ モデリング手続き

データの概要

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

### データの読み込み
# y:調査地の植物個体数、x:調査地の明るさ、A:調査地の面積
data4b = pd.read_csv('./data/ch06/data4b.csv')
print('data4b.shape: ', data4b.shape)
data4b.head()

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

【変数の説明】
調査地における植物の個体数 y、調査地の明るさ x、調査地の面積 A です。

解析の目的、ポイント

■ 解析の目的

■ ポイント
「人口密度」は面積 $${A_i}$$当たりの個体数 $${y_i}$$、つまり $${y_i / A_i}$$ です。
割り算で求めているので、テキストは「割算値」と呼んでいます。
割算値を目的変数にする統計モデルの是非（答えは非）と、対応策（オフセット項）を学びます。

■ 「人口密度」変数の追加
ポイントを受けて、人口密度に関するデータも見ていこうと思います。
data に人口密度 $${y_i / A_i}$$ を表す「y_per_A」列を追加します。

# 割算値「面積当たりの個体数y_per_A」を追加
data['y_per_A'] = data.y / data.A
print('data.shape: ', data.shape)
data.head()

【実行結果】

データの確認

■ 基本統計量の確認
基本的な統計量やチャートでデータを概観します。

① 要約統計量と標本分散の表示

# 要約統計量
data.describe().round(3)

【実行結果】
個体数 y の範囲が $${13}$$ から $${95}$$ とかなり広くなっています。

# 標本分散
data.var(ddof=1).rename('var').to_frame().T.round(3)

【実行結果】

【考察】
ポアソン分布は「平均＝分散」の分布です。
目的変数候補の標本平均と標本分散が近似しているか比較すると…

個体数 y の標本平均と標本分散は大きく乖離しています。
個体数は調査地の面積に比例すると考えられますが、面積が一定ではないので、個体数のバラツキが大きくなると思われます。
一方で、割算値の y_per_A は標本平均よりも標本分散が小さくなっています。
オフセット項を用いる場合の分散の評価を後で確認することにして、先に進みます。

■ 相関係数の表示

# 相関係数
data.corr().round(3)

【実行結果】

個体数 y と明るさ x の相関係数は $${0.530}$$ であり、中程度の正の相関があります。
人口密度 y_per_A と明るさ x の相関係数は $${0.820}$$ であり、強い正の相関があります。

■ ヒストグラムと散布図の描画
すべての変数のヒストグラムと散布図を seaborn の pairplot で描画します。

# ヒストグラムと散布図の描画
sns.pairplot(data=data, height=1.5,
             diag_kws={'edgecolor': 'white', 'alpha': 0.7});

【実行結果】
各変数のヒストグラムの峰は１つです。
目的変数候補と明るさ x との間には正の相関を確認できます。

■ 個体数と面積と明るさの散布図の描画
テキスト p.132 図 6.10 (A)の「x軸：面積 A、y軸：個体数 y、hue：明るさ x」の散布図を描画します。
予測時に使うチャートを先取りして、見慣れておきます。

# 散布図の描画

# 描画領域の設定
fig, ax = plt.subplots(figsize=(6, 4))

# 観測データの散布図の描画
sns.scatterplot(data=data, x='A', y='y', hue='x', s=80, ax=ax,
                edgecolor='lightblue', palette='Blues_r')
# 修飾
ax.set(xlim=(0, 18), ylim=(0, 100), title='観測データの散布図')
ax.set_xlabel('調査地の面積 $A_i$', fontsize=14)
ax.set_ylabel('植物の個体数 $y$', fontsize=14)
ax.legend(title='明るさ $x$');

【実行結果】
x 軸は面積 A、y 軸は個体数 y です。
面積の大きさに個体数が比例するのは直感的にそう思えます。
縦方向を確認すると個体数が大きくなるにつれて色が薄くなる、つまり明るいほど個体数が大きくなる傾向が分かります。

ポアソン回帰（オフセット項モデル）

テキスト p.130 ～ の「オフセット項の入った線形予測子」の統計モデルに取り組みます。

「割合の目的変数」を回避できるオフセット項

■ 割合を目的変数にすることの問題と回避策
今回の「人口密度」（面積当たり個体数：個体数÷面積）や、時間あたり回数（回数÷時間）、成功率（成功回数÷試行回数）など、「観測値どうしの割算値で作った割合や比率」を目的変数に使わないようにしよう、というのがテキストの大主張です。

観測値どうしの割合を目的変数に使うことの問題は：

• 個々の観測値が持つ確からしさなどの 情報が失われること
  例えば、「1 回の試行で 1 回成功」と「100 回の試行で 100 回成功」は、割合にすると両方１になります。
  ２つのケースは同じ確からしさではないです（後者のほうが成功の確からしさが大きいはず）が、割合に要約してしまうことで、確からしさの情報が失われます。

• 割合の確率分布を特定しにくい
  例えば、分子の個体数がポアソン分布に従い、分母の面積が定数や変数のとき、割合の「人口密度」が従う確率分布は何になるのでしょう？
  目的変数が従う確率分布が分からなかったり複雑になると、尤度関数が見つからずに最尤推定によるパラメータ推定が困難になる場合があります。

この割合を回避するカギがオフセット項なのです。

◆ ◆ ◆

■ オフセット項
GLM におけるオフセット項は、係数が１で固定されている変数です。
通常、線形予測子の各変数に対して係数（パラメータ）を推定しますが、オフセット項は「係数を推定しない変数」です。

オフセット項を使えば割合を目的変数にする必要が無くなることを、テキストは教えてくれます。

ポアソン分布の平均パラメータ $${\lambda_i}$$ は今回の例題データの場合、平均個体数になります。
「平均的な人口密度」は観測値 $${i}$$ の平均個体数 $${\lambda_i}$$ を面積 $${A_i}$$ で割った値です。
このことを対数リンク関数と線形予測子で表してみます。

$${\log \frac{M}{N} = \log M - \log N}$$ のように、割り算の対数は対数の引き算で表せることを利用して、上の式を変形していきます。

最後の式の $${\log A_i}$$ がオフセット項です。
割合の分母に置きたい面積の対数になっています。

対数リンク関数を使うGLMでは、割合の分母に使いたい変数の対数をオフセット項に加えることで、割り算の意味を正しくモデルに組み込むことができます。

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

オフセット項を使うポアソン回帰の GLM を観測データに当てはめます。
モデル名は「$${\mathtt{x + offset}}$$ モデル」（非公式）です。

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

🔷 確率分布と確率質量関数
調査地 $${i}$$ における個体数 $${y_i}$$ は平均個体数 $${\lambda_i}$$ のポアソン分布に従います。

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

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

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

🔷 今回の統計モデルの尤度関数 $${L}$$ と対数尤度関数 $${\log L}$$

◆ ◆ ◆

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

# オフセット項を指定してポアソン回帰の当てはめを行う p.133

# 設定
family = sm.families.Poisson()  # GLMの引数familyに与える確率分布＝ポアソン分布

# モデルの当てはめ
result = smf.glm(
    formula='y ~ x', offset=np.log(data.A), data=data, family=family).fit()
result.summary()

【実行結果】
最下２行の Intercept が $${\beta_1}$$、x が $${\beta_2}$$に対応しています。

【オフセット項の記述】
glm の引数 offset にオフセット項を記述します。
（formula には記述しません）

result = smf.glm(..., offset=np.log(data.A), ...)

◆ ◆ ◆

■ 当てはめ結果の分析
🔷 係数の推定値

切片 Intercept、体サイズ x の係数の推定値を線形予測子に当てはめてみます。

🔷 パラメータ推定値の評価
２つ係数の推定値の $${p}$$ 値（P>|z|）は $${0.000 < 0.05}$$、95% 信頼区間（[0.025  0.975]）はゼロを含んでいません。
統計的な視点では、パラメータ推定値は有意ということになります。
また、係数の推定値の標準誤差は、推定値と比べて小さな値であり、特段の問題は無いように思われます。

🔷 モデルの評価
AIC で予測の良さを確認します。

# Null Deviance, Residual Deviance, AICの表示
print(f'Null 逸脱度\t: {result.null_deviance:.1f}')
print(f'残差逸脱度\t: {result.deviance:.1f}')
print(f'AIC\t\t: {result.aic:.1f}')

【実行結果】
AIC は 650.3 です。
比較するモデルが無いので、ここでは参考値の扱いです。

平均個体数の予測

テキストにならって、平均個体数の予測を行って可視化します。
テキスト p.132 図 6.10 に相当します。
面白いのは、x 軸に面積 A を使い、明るさ x をカラーマッピングで表現するところです。

# オフセット項を利用するGLMを説明するための例題 p.132 図6.10

## 設定
xlabel, ylabel = '調査地の面積 $A_i$', '植物の個体数 $y$'  # x軸とy軸のラベル

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

## (A)観測データの描画
# 観測データの散布図の描画
sns.scatterplot(data=data, x='A', y='y', hue='x', s=80, ax=ax1,
                edgecolor='lightblue', palette='Blues_r')
# 修飾
ax1.set(xlim=(0, 18), ylim=(0, 100), xlabel=xlabel, ylabel=ylabel,
        title='(A) 観測データ')
ax1.legend(title='明るさ $x$')

## (B)推定されるモデルによる予測
# 観測データの散布図の描画
sns.scatterplot(data=data, x='A', y='y', hue='x', s=80, ax=ax2,
                edgecolor='lightblue', palette='Blues_r', legend=False)
# 予測に与えるAの値の設定
A_vals =np.linspace(0.01, 18, 100)
# xの値ごとに予測値の算出と直線の描画を繰り返し処理
for x_val in [0.1, 0.3, 0.5, 0.7, 0.9]:
    # yの予測値の算出
    y_vals = result.predict(
        dict(x=[x_val]*len(A_vals)), offset=np.log(A_vals)
    )
    # xの値に応じた直線の色の取得（値が小さい：暗い、値が大きい：明るい）
    color = cm.Blues_r(x_val)
    # 予測値の直線の描画
    sns.lineplot(x=A_vals, y=y_vals, color=color, lw=3, label=f'{x_val}', ax=ax2)
# 修飾
ax2.set(xlim=(0, 18), ylim=(0, 100), xlabel=xlabel, ylabel='', 
        title='(B) 推定されたモデルによる予測')
ax2.legend(title='明るさ $x$');

【実行結果】
(B) が平均個体数の予測値のプロットです。

まず、面積が大きくなるにつれて平均個体数が直線的に大きくなることが分かります。
そして、明るさごとに色を変えた５本の直線を比較することで、明るいほど（直線の色が薄いほど）平均個体数が大きくなることが分かります。

ちなみに x 軸を明るさ x にするバージョンはこちらです。

# オフセット項を利用するGLMを説明するための例題 ※x軸に明るさxを適用

# 描画領域の設定
fig, ax = plt.subplots(figsize=(6, 4))

# 観測データの散布図の描画
sns.scatterplot(data=data, x='x', y='y', hue='A', s=80, ax=ax,
                edgecolor='lightblue', palette='Greens_r', legend=False)
# 予測に与えるAの値の設定
x_vals =np.linspace(0.01, 18, 100)
# xの値ごとに予測値の算出と直線の描画を繰り返し処理
for A_val in [1, 5, 10, 15, 20]:
    # yの予測値の算出
    y_vals = result.predict(
        dict(x=x_vals), offset=np.log([A_val]*len(x_vals))
    )
    # 予測値の直線の描画
    sns.lineplot(x=x_vals, y=y_vals, color=cm.Greens_r(A_val*10), lw=3,
                 label=f'{A_val}', ax=ax)
# 修飾
ax.set(xlim=(0, 1), ylim=(0, 100), title='(C) 推定されたモデルによる予測')
ax.set_xlabel('調査地の明るさ $x$', fontsize=14)
ax.set_ylabel('植物の個体数 $y$', fontsize=14)
ax.legend(title='面積 $A$');

【実行結果】
明るさに比例して個体数が非線形に大きくなることが分かります。

テキストはオフセット項で割合を表すモデルについて、次のようにまとめています。

アディショナルタイム

テキストを飛び越えて、３つのテーマに取り組みます！

• 最適化ライブラリでパラメータ推定
  オフセット項を追加しても最適化ライブラリでパラメータ推定できることを実験します

• オフセット項を使う場合の分散の評価
  「データの確認」の際に先送りした分散の評価を深堀りします。

• 割合を目的変数にするGLMの適否
  テキストがNGと言っている「人口密度」を目的変数にして、線形回帰モデルに当てはめします。

（注意）
こちらは趣味的な深堀りとコードです。
ご興味ない方はスルーしてくださって大丈夫です。

最適化ライブラリでパラメータ推定

オフセット項を含むポアソン回帰のパラメータ推定を scipy の minimize で実践します。
最適化の目的関数内の線形予測子に $${\log(A)}$$ を含めるだけです。

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

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

# 目的関数の定義：負の対数尤度関数 ※負の対数尤度の最小化問題を解く
def llf(params, x, y, A):
    # パラメータの分解
    beta1, beta2 = params
    # λ = exp(線形予測子)
    lam = np.exp(beta1 + beta2 * x + np.log(A)) # ⬅️ オフセット項 log(A)
    # 戻り値：負の対数尤度
    return -np.sum(stats.poisson.logpmf(k=y, mu=lam))

# 最適化の実行
res_optim = minimize(fun=llf, x0=[0, 0], args=(data['x'], data['y'], data['A']))

# 結果の表示
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}')

【実行結果】
GLM のパラメータ推定値と一致しました。

オフセット項を使う場合の分散の評価

■ ChatGPTによると：
ポアソン回帰では「平均 = 分散」という仮定があるため、実データがこの仮定に反して過大分散（分散が平均より大きい）、過小分散（分散が平均より小さい）を持っていると、推定結果の信頼区間や仮説検定の正しさに影響します。
ただし、オフセット項を用いると単純な「目的変数の標本平均と標本分散の比較」では評価が不十分になるので、モデル残差の観点から過大・過小分散を評価する必要があります。

データの確認で見た標本平均と標本分散の表を再掲します。

個体数 y の標本平均と標本分散は、オフセット項による「実質的な割合化」を反映していないので、単純な比較では不十分です。
人口密度 y_per_A の場合は、割合（割算値）に変換済みなので、ポアソン分布に従っているとはいい難く、標本平均と標本分散の比較は無意味な気がします。

このような状況の場合の分散等の評価についてChatGPTに訊いたところ、残差を用いる次の方法があると教えてくれました。

やってみましょう！

◆ ◆ ◆

■ 分散パラメータ（dispersion parameter）
分散パラメータの計算過程を示しつつ、分散パラメータに迫ります。

①ピアソン残差 $${r_i}$$ の算出
目的変数の観測値 $${y_i}$$、モデルによる予測値 $${\hat{y}_i}$$、予測値の分散 $${V(\hat{y}_i)}$$ とするとき、ピアソン残差 $${r_i}$$ は次の式で計算されます。

今回のポアソン回帰の場合、「平均」個体数の予測値は  $${\widehat{\lambda}_i}$$ であり、また、ポアソン分布の分散は平均と等しいので、次のように書き換えできます。

② ピアソン $${\chi^2}$$ 統計量
ピアソン $${\chi^2}$$ 統計量はピアソン残差の２乗和です。

③ 分散パラメータ $${\phi}$$
分散パラメータ $${\phi}$$ はピアソン $${\chi^2}$$ 統計量をモデルの残差の自由度で割ったものです。

④ 分散パラメータの見方

⑤ Python で分散パラメータを計算
分散パラメータを計算します。
GLM の結果 resultの属性「pearson_chi2：ピアソン $${\chi^2}$$ 統計量」と「df_resid：残差の自由度」を利用します。

# 分散パラメータ（ピアソン残差ベース）
dispersion = result.pearson_chi2 / result.df_resid
print(f"分散パラメータ φ = {dispersion:.3f}")

【実行結果】

分散パラメータは $${0.832}$$ です。
ChatGPTによると：

⑥ Python で計算過程を追う
ピアソン残差 $${r_i}$$ とピアソン $${\chi^2}$$ 統計量の計算を Pythonで確認します。

・ピアソン残差 $${r_i}$$（先頭５行）
GLM の結果 resultの属性「resid_pearson：ピアソン残差」の取得と、ピアソン残差の計算式に基づく計算の２つを行います。

# ピアソン残差の計算

# GLMの結果からresid_pearson属性で取得
print('result.resid_pearson:')
print(result.resid_pearson.iloc[:5], '\n')

# ピアソン残差の公式で計算
print('ri = (yi - hat{yi}) / sqrt(hat{yi}):')
print((
    (result.model.endog - result.fittedvalues) / np.sqrt(result.fittedvalues)
).iloc[:5])

【実行結果】

・ピアソン $${\chi^2}$$ 統計量
GLM の結果 resultの属性「pearson_chi2：ピアソン $${\chi^2}$$ 統計量」の取得と、ピアソン残差の二乗和の計算の２つを行います。

# ピアソンΧ²統計量

# GLMの結果からpearson_chi2属性で取得
print('result.pearson_chi2\t:', result.pearson_chi2)

# ピアソン残差の二乗和で計算
print('ピアソン残差の二乗和\t:', sum(result.resid_pearson**2))

【実行結果】

◆ ◆ ◆

■ 残差プロット
x 軸：予測値、y 軸：ピアソン残差の散布図を描画します。

# 予測値とピアソン残差の可視化

# 予測値とピアソン残差の散布図の描画
plt.scatter(result.fittedvalues, result.resid_pearson)
# y=0, -2, 2 の水平点線の描画
plt.axhline(0, color='gray', linestyle='--')
plt.axhline(-2, color='tab:red', linestyle='--')
plt.axhline(2, color='tab:red', linestyle='--')
# 修飾
plt.xlabel('平均個体数の予測値', fontsize=12)
plt.ylabel('ピアソン残差', fontsize=12)
plt.title('ピアソン残差 vs. 予測値')
plt.grid(alpha=0.3)
plt.show()

【実行結果】

【残差プロットの評価ポイント】
ChatGPT に評価ポイントを聴きました。

残差プロットを確認しましょう。
① はランダムにばらついているのでＯＫです。
② は多くの残差が $${-2 ～ +2}$$ の範囲に収まっているのでＯＫです。
③ は予測値が大きくなるにつれて残差のばらつきが広がる傾向は見られないのでＯＫです。

割合を目的変数にする GLM の適否

割合である人口密度 $${y_i / A_i}$$ を目的変数にした線形回帰モデルを例題データに当てはめて、いろいろ眺めてみたいと思います。
モデル名は「人口密度モデル」（非公式）です。

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

🔷 確率分布
人口密度 $${y_i / A_i}$$ は正規分布に従うと仮定します。

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

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

# 目的変数を割合にして、線形回帰（正規分布・恒等リンク関数）の当てはめ

# 設定
family = sm.families.family.Gaussian()  # GLMの引数familyに与える確率分布＝正規分布

# モデルの当てはめ
res_lr = smf.glm(formula='y_per_A ~ x', data=data, family=family).fit()
res_lr.summary()

【実行結果】
最下２行の Intercept が $${\beta_1}$$、x が $${\beta_2}$$に対応しています。

■ 当てはめ結果の分析
🔷 係数の推定値

切片 Intercept、明るさ x の係数の推定値を線形予測子に当てはめてみます。

🔷 パラメータ推定値の評価
２つ係数の推定値の $${p}$$ 値（P>|z|）は $${0.000 < 0.05}$$、95% 信頼区間（[0.025  0.975]）はゼロを含んでいません。
統計的な視点では、パラメータ推定値は有意ということになります。
また、係数の推定値の標準誤差は、推定値と比べて小さな値であり、特段の問題は無いように思われます。

◆ ◆ ◆

■ モデルの選択
① AIC
AIC で予測の良さを確認します。

# Null Deviance, Residual Deviance, AICの表示
print(f'Null 逸脱度\t: {res_lr.null_deviance:5.1f}')
print(f'残差逸脱度\t: {res_lr.deviance:5.1f}')
print(f'AIC\t\t: {res_lr.aic:5.1f}')

【実行結果】
人口密度モデルの AIC は 196.7 です。

ポアソン回帰の$${\mathtt{x + offset}}$$ モデルと線形回帰の人口密度モデルの評価指標を比べてみます。

# ポアソン回帰と割算値の線形回帰の評価指標

# 追加インポート
from sklearn.metrics import r2_score as R2
from sklearn.metrics import mean_absolute_percentage_error as MAPE

# 疑似決定係数、目的変数の平均の予測値の算出
pseudo_r2_pois = result.pseudo_rsquared()
pseudo_r2_lr = res_lr.pseudo_rsquared()
y_pred_pois = result.fittedvalues
y_pred_lr = res_lr.fittedvalues

# データフレーム表示
pd.DataFrame({
    'AIC': [result.aic, res_lr.aic],
    'MAPE': [MAPE(data.y, y_pred_pois), MAPE(data.y_per_A, y_pred_lr)], 
    '（参考）疑似決定係数（Cox & Snell’s）': [pseudo_r2_pois, pseudo_r2_lr],
    '（参考）決定係数 R²': [R2(data.y, y_pred_pois), R2(data.y_per_A, y_pred_lr)],
    }, index=['x+offset：ポアソン回帰', '人口密度：線形回帰']
).round(4)

【実行結果】

まず AIC を見ます。
値としては人口密度モデルのほうが AIC が小さいのでドキッとします。
目的変数や確率分布が異なるモデルなので、単純に AIC を比較できないでしょう。
けれども、これだけ差が大きいと、AIC の小さい人口密度モデルに心惹かれる気持ちも理解できそうです。

続いて MAPE です。
MAPE は観測値に対する予測誤差の比率（絶対値）の平均です。
値が小さいほど予測誤差が小さい＝予測精度が高いことになります。
$${\mathtt{x+offset}}$$ モデルのほうが予測誤差が小さいモデルと言えます。
ただし、モデルの学習で使用した観測値で予測している点には注意が必要です。

参考値で疑似決定係数と決定係数を併記しました。
両指標は値が大きいほど当てはまりが良いことを示します。
線形回帰では決定係数を用い、線形回帰以外では疑似決定係数を用いるようです。
疑似決定係数は GLM 結果の 「Pseudo R-squ. (CS)」の値です。

◆

こうやって比較してみると、目的変数や確率分布が異なるモデルの比較は難しいなぁという印象です。
また AIC の単純な値比較によって、割合を用いる人口密度モデルを選択したい誘惑に駆られます…

そこで ChatGPT の冷静な意見を聴きましょう。

📌 理論的にはポアソン回帰が正当な選択です。

• 個体数（y）を目的変数とし、オフセット log(A) を使うことで、「面積あたりの個体生息率」という率の構造を適切にモデル化しています。

• モデルの確率分布の前提にも合っていて、GLMの基本にも忠実です。

📌 MAPEによる予測精度の差はわずかです。

• 線形回帰も悪くはありませんが、「y/A」のような割合の目的変数は、理論的な正しさに不安があります。

• 尤度や分布仮定を考えると、将来的に汎化性能や頑健性で劣る可能性もあります。

📌 したがって実務でも理論でも、「y~x, offset=log(A)」というGLMのほうが望ましいと言えるでしょう。

まとめ

今回はオフセット項を含むポアソン回帰を実践しました。

🔷 オフセット項
GLM におけるオフセット項は、係数が１で固定されている変数です。
オフセット項は、観測値どうしの割合・比率を目的変数にしないための回避策に活用できます。

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

🔷 リンク関数と線形予測子
オフセット項は線形予測子の中に、割合の分母の対数 $${\log A_i}$$ の形式で含められます。

🔷 statsmodels のポアソン回帰モデル構築と結果表示
オフセット項は glm の引数 offset に、割合の分母の対数 $${\mathtt{np.log(A)}}$$ の形式で設定されます。

family = sm.families.Poisson()  # GLMの引数familyに与える確率分布＝ポアソン分布
result = smf.glm(
    formula='y ~ x', offset=np.log(data.A), 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の教科書です。
よかったらぜひ、お試しくださいませ。

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