2章ではとりあえず「ポアソン分布」を導入して、データがどういう分布に従っているのかを調べました。この章では、説明変数をどのようにして目的変数に取り入れるのかをポアソン回帰を用いて学習します。
3.2 観測されたデータの概要を調べる
- データの読み込み
import pandas as pd
d = pd.read_csv("data3a.csv")
- データの表示(jupyter notebook利用)
from IPython.display import display
display(d)
# dだけでも表示される
- 型を調べる
print(d.dtypes)
# y int64
# x float64
# f object
# dtype: object
- データフレームの概要
d.describe()
d.f.describe()
3.3 統計モデリングの前にデータを図示する
- 散布図
# 散布図の作成
import matplotlib.pyplot as plt
import matplotlib.font_manager
fontprop = matplotlib.font_manager.FontProperties(fname="/Library/Fonts/ヒラギノ丸ゴ ProN W4.ttc")
plt.rcParams['figure.figsize'] = (7.0, 5.0)
d1 = d[d.f=='C']
d2 = d[d.f=='T']
plt.plot(d1.x, d1.y, 'o')
plt.plot(d2.x, d2.y, '^')
# plt.scatter(data.x, data,y) #この書き方でも良い
plt.title("散布図", fontdict={"fontproperties":fontprop}, fontsize=12)
plt.xlabel("体サイズ", fontdict = {"fontproperties": fontprop}, fontsize=12)
plt.ylabel("植物の種子数", fontdict={"fontproperties":fontprop}, fontsize=12)
plt.legend(["C", "T"])
plt.show()
- 箱ひげ図
# 箱ひげ図の作成
fontprop = matplotlib.font_manager.FontProperties(fname="/Library/Fonts/ヒラギノ丸ゴ ProN W4.ttc")
data3 = (d1.y, d2.y)
fig = plt.figure()
ax = fig.add_subplot(111)
bp = ax.boxplot(data3,widths=0.6)
ax.set_xticklabels(['C', 'T'])
plt.grid()
plt.setp(bp['medians'], color='red', linewidth=2)
plt.title("箱ひげ図", fontdict={"fontproperties":fontprop}, fontsize=12)
plt.ylabel("植物の種子数", fontdict={"fontproperties":fontprop}, fontsize=12)
plt.show()
3.4 ポアソン回帰の統計モデル
- 線形予測子と対数リンク関数
ポアソン分布を特徴づけるパラメータである平均($\lambda$)が説明変数である植物の大きさ$(x_i)$に依存する(によって決まる)ようなモデルを考える。
\begin{align}
\lambda_i = \exp{(\beta_1 + \beta_2 x_i)}
\end{align}
\begin{align}
\log{\lambda_i} = \beta_1 + \beta_2 x_i
\end{align}
下の式の左辺の関数のことを「リンク関数」と呼び、今回の場合は$\log$がリンク関数に当たる。また右辺についてはパラメータ$\beta_{i}$が線形和で書かれているため「線形予測子」と呼ぶ。
fontprop = matplotlib.font_manager.FontProperties(fname="/Library/Fonts/ヒラギノ丸ゴ ProN W4.ttc")
x = np.arange(-5 ,6, 0.1)
y1 = [np.exp(-2 - _x * 0.8) for _x in x]
y2 = [np.exp(-1 + _x * 0.4) for _x in x]
y1 = pd.Series(y1, index = x)
y2 = pd.Series(y2, index = x)
plt.plot(y1, "-", label=r"{$\beta_{1}, \beta_{2}$}={-2, -0.8}")
plt.plot(y2, "--", label=r"{$\beta_{1}, \beta_{2}$}={-1, 0.4}")
plt.xlim(-4.2, 5.2)
plt.ylim(0,2.8)
plt.xlabel("個体$i$の体サイズ$x_{i}$", fontdict={"fontproperties":fontprop}, fontsize=12)
plt.ylabel("個体$i$の$\lambda_{i}$", fontdict={"fontproperties":fontprop}, fontsize=12)
plt.vlines([0], 0, 3.0, "black", linestyles='dashed')
plt.legend(loc='upper right')
plt.show()
- モデルを用いたfitting(ポアソン回帰)
# Poisson回帰
import statsmodels.api as sm
import statsmodels.formula.api as smf
model = smf.glm('y ~ x', data=d, family=sm.families.Poisson()) # ポアソン回帰
results = model.fit() # fitting
print(results.summary()) # fitting結果の表示
results.llf # 最大対数尤度の表示
# Generalized Linear Model Regression Results
# ==============================================================================
# Dep. Variable: y No. Observations: 100
# Model: GLM Df Residuals: 98
# Model Family: Poisson Df Model: 1
# Link Function: log Scale: 1.0
# Method: IRLS Log-Likelihood: -235.39
# Date: Sat, 30 Dec 2017 Deviance: 84.993
# Time: 14:50:06 Pearson chi2: 83.8
# No. Iterations: 4
# ==============================================================================
# coef std err z P>|z| [0.025 0.975]
# ------------------------------------------------------------------------------
# Intercept 1.2917 0.364 3.552 0.000 0.579 2.005
# x 0.0757 0.036 2.125 0.034 0.006 0.145
# ==============================================================================
# -235.38625076986077 #最大対数尤度
- ポアソン回帰モデルによる予測
# fittingの結果を図示
xx = np.arange(min(d.x), max(d.x), 0.1)
yy = [np.exp(1.29 + 0.0757 * _x) for _x in xx]
plt.plot(d1.x, d1.y, 'o')
plt.plot(d2.x, d2.y, '^')
plt.plot(xx, yy)
plt.xlabel("d.x", fontsize=12)
plt.ylabel("d.y", fontsize=12)
plt.show()
3.5 説明変数が因子型の統計モデル
上と同じことを因子型の説明変数について行う。
# 説明変数が因子型の統計モデル
model = smf.glm('y ~ f', data=data, family=sm.families.Poisson())
results2 = model.fit()
print(results2.summary()) # 結果の表示
results2.llf # 最大対数尤度を評価
# Generalized Linear Model Regression Results
# ==============================================================================
# Dep. Variable: y No. Observations: 100
# Model: GLM Df Residuals: 98
# Model Family: Poisson Df Model: 1
# Link Function: log Scale: 1.0
# Method: IRLS Log-Likelihood: -237.63
# Date: Mon, 01 Jan 2018 Deviance: 89.475
# Time: 21:10:43 Pearson chi2: 87.1
# No. Iterations: 4
# ==============================================================================
# coef std err z P>|z| [0.025 0.975]
# ------------------------------------------------------------------------------
# Intercept 2.0516 0.051 40.463 0.000 1.952 2.151
# f[T.T] 0.0128 0.071 0.179 0.858 -0.127 0.153
# ==============================================================================
# -237.62725696068688
- 対数リンク関数のわかりやすさ
# 恒等リンク関数を使ってPoisson回帰
model = smf.glm('y ~ x + f', data=data, family=sm.families.Poisson(link=sm.genmod.families.links.identity))
results4 = model.fit()
print(results4.summary())
results4.llf
# Generalized Linear Model Regression Results
# ==============================================================================
# Dep. Variable: y No. Observations: 100
# Model: GLM Df Residuals: 97
# Model Family: Poisson Df Model: 2
# Link Function: identity Scale: 1.0
# Method: IRLS Log-Likelihood: -235.16
# Date: Mon, 01 Jan 2018 Deviance: 84.538
# Time: 21:45:50 Pearson chi2: 83.6
# No. Iterations: 5
# ==============================================================================
# coef std err z P>|z| [0.025 0.975]
# ------------------------------------------------------------------------------
# ntercept 1.2671 2.843 0.446 0.656 -4.306 6.840
# f[T.T] -0.2048 0.582 -0.352 0.725 -1.346 0.936
# x 0.6606 0.290 2.281 0.023 0.093 1.228
# ==============================================================================
#
# /Users/knakamura/.pyenv/versions/anaconda3-4.4.0/envs/M3/lib/python3.6/site-packages/statsmodels/genmod/generalized_linear_model.py:244: DomainWarning: The identity link function does not respect the domain of the Poisson family.
# DomainWarning)
# -235.15897254420767
ポアソン回帰で恒等リンク関数を使おうとすると警告が表示されます。。。ちなみに、statsmodelsにおけるポアソン回帰で使うことができるリンク関数は以下のコマンドで調べることができる。
sm.families.family.Poisson.links
# [statsmodels.genmod.families.links.log,
# statsmodels.genmod.families.links.identity,
# statsmodels.genmod.families.links.sqrt]
- 図3.8
plt.rcParams['figure.figsize'] = (10.0, 5.0)
fontprop = matplotlib.font_manager.FontProperties(fname="/Library/Fonts/ヒラギノ丸ゴ ProN W4.ttc")
x = np.arange(5, 20.1, 0.1)
y1 = [np.exp(1.26 + 0.08 * _x - 0.032 * 3) for _x in x]
y2 = [np.exp(1.26 + 0.08 * _x ) for _x in x]
y1 = pd.DataFrame(y1, index = x)
y2 = pd.DataFrame(y2, index = x)
z1 = pd.DataFrame([(1.27 + 0.661 * _x - 0.205 * 3) for _x in x], index=x)
z2 = pd.DataFrame([(1.27 + 0.661 * _x ) for _x in x], index=x)
plt.subplot(1,2,1)
plt.plot(y1)
plt.plot(y2)
plt.xlim(4, 21)
plt.ylim(3, 19)
plt.xticks([5,10,15,20])
plt.yticks([5,10,15])
plt.xlabel("体サイズ$x_{i}$", fontdict={"fontProperties":fontprop}, fontsize=12)
plt.ylabel("平均種子数$\lambda_{i}$", fontdict={"fontProperties":fontprop}, fontsize=12)
plt.subplot(1,2,2)
plt.plot(z1)
plt.plot(z2)
plt.xlim(4, 21)
plt.ylim(3, 19)
plt.xticks([5,10,15,20])
plt.yticks([5,10,15])
plt.xlabel("体サイズ$x_{i}$", fontdict={"fontProperties":fontprop}, fontsize=12)
plt.ylabel("平均種子数$\lambda_{i}$", fontdict={"fontProperties":fontprop}, fontsize=12)
plt.show()
