投稿者: Yutaro

以前の記事で、statsmodelsを使って線形回帰をやったので、今回はポアソン回帰をやってみます。
参考:statsmodelsで重回帰分析

データは久保拓弥先生の、データ解析のための統計モデリング入門 (通称緑本)の第３章から拝借し、
本に載っているのと同じ結果を得ることを目指します。
ちなみに本のコードはRで実装されています。
Rは勉強中なので写経はしましたがそれはそれとして、
今回はいつも使っているpythonのstatsmodelsでやってみます。

データはこちらのサポートページからダウンロードできます。
生態学のデータ解析 – 本/データ解析のための統計モデリング入門
３章の data3a.csv を保存しておきましょう。

このデータは、架空の植物の種子の数に関するもので、
種子の数が$y_i$, 植物のサイズが$x_i$, 施肥処理を行ったかどうかが$f_i$列に格納されています。

そして、種子の個数$y_i$が期待値 $\lambda_i=\exp(\beta_1+\beta_2x_i+\beta_3f)$のポアソン分布に従うと仮定した場合の尤度を最大化するように係数を決定します。
($f$は施肥処理を行ったら1, 行ってない場合は0)

早速やってみましょう。
ドキュメントはこのあたりが参考になります。
statsmodels.genmod.generalized_linear_model.GLM
Model Family に　Poisson を指定すると、リンク関数は自動的にlogが選ばれます。

import pandas as pd
import statsmodels.api as sm

# データの読み込み
# ファイルの取得元
# http://hosho.ees.hokudai.ac.jp/~kubo/ce/IwanamiBook.html
df = pd.read_csv("./data3a.csv")

# f 列の値を Tならば 1 , その他(Cのみ)ならば0に符号化
df["fT"] = (df.f == "T").astype(int)

y = df.y
X = df[["x", "fT"]]
# 定数列(1)を作成
X = sm.add_constant(X)

# モデル生成と学習
model = sm.GLM(y, X, family=sm.families.Poisson())
result = model.fit()

# 結果出力
print(result.summary())

# 以下出力結果
                 Generalized Linear Model Regression Results                  
==============================================================================
Dep. Variable:                      y   No. Observations:                  100
Model:                            GLM   Df Residuals:                       97
Model Family:                 Poisson   Df Model:                            2
Link Function:                    log   Scale:                          1.0000
Method:                          IRLS   Log-Likelihood:                -235.29
Date:                Sun, 21 Apr 2019   Deviance:                       84.808
Time:                        16:36:24   Pearson chi2:                     83.8
No. Iterations:                     4   Covariance Type:             nonrobust
==============================================================================
                 coef    std err          z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------
const          1.2631      0.370      3.417      0.001       0.539       1.988
x              0.0801      0.037      2.162      0.031       0.007       0.153
fT            -0.0320      0.074     -0.430      0.667      -0.178       0.114
==============================================================================

本のP58に載っているのと全く同じ係数を得ることができました。

ついでに実際のデータとこのモデルから予測される種子数をプロットして可視化してみます。

import matplotlib.pyplot as plt
import numpy as np

# プロット用のデータ作成
xx = np.linspace(7, 13, 101)
yy0 = np.exp(result.params["const"] + result.params["x"]*xx)
yy1 = np.exp(result.params["const"] + result.params["x"]*xx + result.params["fT"])

# 可視化
fig = plt.figure(figsize=(8, 7))
ax = fig.add_subplot(1, 1, 1)
ax.set_xlabel("個体のサイズ")
ax.set_ylabel("種子の個数")
ax.scatter(X[X.fT == 0]["x"], y[X.fT == 0].values, marker="x", label="肥料無し")
ax.scatter(X[X.fT == 1]["x"], y[X.fT == 1].values, alpha=0.7, label="肥料有り")
ax.plot(xx, yy0, linestyle="--",  label="肥料無し")
ax.plot(xx, yy1, label="肥料有り")
plt.legend()
plt.show()

出力結果がこちら。

前回に引き続きpythonで重回帰分析を行い、回帰係数と切片を求める方法の紹介です。
今回はstatsmodelsを使います。

ドキュメントはこちら。
statsmodels.regression.linear_model.OLS

データは同じものを使い、結果が一致することを確認したいので
保存してたものを読み込みます。

import numpy as np
import statsmodels.api as sm

# データの読み込み
npzfile = np.load("sample_data.npz")
X = npzfile["X"]
y = npzfile["y"]

つぎに、回帰分析に入りますが、statsmodelsで回帰分析する場合には、一点注意があります。
それは、上で読み込んだ Xとyをそのままモデルに食わせると、切片(定数項)を0として結果を返してくることです。
それを回避するために、 X に対して、全ての値が1の列を追加する必要があります。
それは、次の関数を使って実行しています。
statsmodels.tools.tools.add_constant

# 配列Xに1だけからなる列を追加する。
X0 = sm.add_constant(X)
# 回帰分析実行
st_model = sm.OLS(y, X0)
results = st_model.fit()

# 切片と回帰係数
print(results.params)
# [ 4.94346432  2.00044153 -2.99255801  3.98231315  0.01708309]

無事にscikit-learnで行った時と全く同じ結果を得ることができましたね。
これだけだと、add_constantをやらなくてよかった分scikit-learnの圧勝に見えますが、
statsmodels を使うと一つメリットがあります。
それは、各係数が0でないことの検定ができることです。
resultsの sumary()って関数を実行してみましょう。

print(results.summary())

# 以下出力です。
                            OLS Regression Results                            
==============================================================================
Dep. Variable:                      y   R-squared:                       0.999
Model:                            OLS   Adj. R-squared:                  0.999
Method:                 Least Squares   F-statistic:                 2.136e+04
Date:                Mon, 15 Apr 2019   Prob (F-statistic):          2.24e-139
Time:                        00:59:22   Log-Likelihood:                -143.94
No. Observations:                 100   AIC:                             297.9
Df Residuals:                      95   BIC:                             310.9
Df Model:                           4                                         
Covariance Type:            nonrobust                                         
==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
const          4.9435      0.108     45.726      0.000       4.729       5.158
x1             2.0004      0.020    101.534      0.000       1.961       2.040
x2            -2.9926      0.018   -166.158      0.000      -3.028      -2.957
x3             3.9823      0.017    231.575      0.000       3.948       4.016
x4             0.0171      0.019      0.914      0.363      -0.020       0.054
==============================================================================
Omnibus:                        2.170   Durbin-Watson:                   1.994
Prob(Omnibus):                  0.338   Jarque-Bera (JB):                1.862
Skew:                          -0.208   Prob(JB):                        0.394
Kurtosis:                       2.477   Cond. No.                         7.12
==============================================================================

Warnings:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.

出力の真ん中にある通り、定数項および回帰係数のt値と、p値が出力されています。
そして、x4 は P値が 0.05より大きいので、x4の係数が0であるという帰無仮説を棄却できていないことがわかります。
もともと、x4の回帰係数は0としてサンプルデータを作っていたので、これは良い結果です。
そのほか、AICなどの値も出してくれています。

scikit-learnでこれらの情報を収集するのは手間なので、
目的によってはstatsmodelsが便利な場面もあります。