前回に引き続き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が便利な場面もあります。