書記が数学やるだけ#602 ポアソン回帰

ポアソン回帰についての実装例を示す。

問題

説明

ポアソン回帰では，対数関数を連結関数としてポアソン分布で推測する。よく似たものとして対数正規モデルがあるが，連結関数が対数関数である点は同じだが，正規分布により推測する点が異なる。

ポアソン分布ではなさそうな場合のうち，過分散が疑われる場合には負の2項分布モデルを用いることがある。

使い分けについては以下を参照：

「使い分け」ではなく「妥当かどうか」が大事：重回帰分析＆一般化線形モデル選択まわりの再まとめ - 渋谷駅前で働くデータサイエンティストのブログ

解答

import math
import numpy as np
import scipy.stats as stats
import statsmodels.api as sm
import statsmodels.formula.api as smf
import seaborn as sns
import matplotlib.pyplot as plt
import pandas as pd

# チップのデータセット
df = sns.load_dataset("tips")
df

可視化をするにはSeabornのjointplotが便利である。

sns.set(style="darkgrid")
sns.jointplot(x="total_bill", y="tip", data=df,
             kind="scatter",
             xlim=(0, 60), ylim=(0, 12),
             color="b",
             height=7);

ここで標本のパラメータを求めておく。

mean = df['tip'].mean()
var = df['tip'].var()
scale = math.sqrt(var)

print("平均:", mean, "分散:", var, "標準偏差:", scale)

平均: 2.9982786885245902 分散: 1.9144546380624725 標準偏差: 1.3836381890011826

説明変数の分布を見ると，裾が長い分正規分布には見えない気もする。正規分布とポアソン分布でのサンプリングを行なってみると，その違いが出てくる。

norm = np.random.normal(mean,scale,244)
poisson = np.random.poisson(mean,244)
bins = np.linspace(-4, 12, 10)

plt.hist(norm, bins=bins, alpha=0.6, label="noem")
plt.hist(poisson, bins=bins,alpha=0.6, label="poisson")
plt.legend()

正規分布でのフィッテイングが妥当かどうか，いくつか試験してみる。

#Q-Qプロット，正規分布ではなさそう
stats.probplot(df['tip'], dist="norm", plot=plt)

#シャピロ・ウィルク検定，帰無仮説「正規分布である」が棄却される
stats.shapiro(df["tip"])

ShapiroResult(statistic=0.897811233997345, pvalue=8.20057563521992e-12)

これらより，単純に正規分布でフィッティングするのは危ういかもしれない。

では，実際にいくつかモデルを立ててみる。まずは正規線形モデル，ここで定数項なしに設定していることに注意（会計を払わない人にチップはないだろう）。

#正規線形モデル，定数項なし
y = df['tip']
x = df['total_bill']
link = sm.families.links.identity()
family = sm.families.Gaussian(link)

model_norm = sm.GLM(y, x, family=family)
results_norm = model_norm.fit()
results_norm.summary()

この予測値を可視化してみると，割といい当てはまりに見える。

y_hat_norm = results_norm.predict(x)

plt.plot(x, y, "o")
plt.plot(x, y_hat_norm, "*", color="r")
plt.xlabel('x (total_bill)'), plt.ylabel('y (tips)')
plt.xlim(0, 60), plt.ylim(0, 12)
plt.show()

次にポアソンモデルを示す。

#ポアソンモデル
link = sm.families.links.log()
family = sm.families.Poisson(link)

model_poisson = sm.GLM(y, x, family=family)
results_poisson = model_poisson.fit()
results_poisson.summary()

y_hat_poisson = results_poisson.predict(x)

plt.plot(x, y, "o")
plt.plot(x, y_hat_poisson, "*", color="r")
plt.xlabel('x (total_bill)'), plt.ylabel('y (tips)')
plt.xlim(0, 60), plt.ylim(0, 12)
plt.show()

対数正規モデルも組んでみる。

#対数正規モデル
link = sm.families.links.log()
family = sm.families.Gaussian(link)

model_lognorm = sm.GLM(y, x, family=family)
results_lognorm = model_lognorm.fit()
results_lognorm.summary()

y_hat_lognorm = results_lognorm.predict(x)

plt.plot(x, y, "o")
plt.plot(x, y_hat_lognorm, "*", color="r")
plt.xlabel('x (total_bill)'), plt.ylabel('y (tips)')
plt.xlim(0, 60), plt.ylim(0, 12)
plt.show()

本記事のもくじはこちら：