PandasとStatsModelsを使って重回帰分析をしてみた

Pythonで重回帰分析をやってみました

多変量解析法入門 (ライブラリ新数学大系)

この文献のP43～86を参考を参考にしています．

マンガでわかる統計学 回帰分析編

また，回帰分析の全体の流れは上の書籍を参考にしています．

今回はStatsModelsというライブラリを使いました．

データを可視化する

データファイルを読み込んでどんなデータかを観察します．
今回は３次元データなので３次元プロットをしてみます．

# データの読み込み
import pandas as pd

file = 'reg_test_data.csv'
df = pd.read_csv(file)
df = df.set_index('val')

import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

# 初期化
fig = plt.figure()
ax = Axes3D(fig)

# 軸ラベルを設定する
ax.set_xlabel("X-axis")
ax.set_ylabel("Y-axis")
ax.set_zlabel("Z-axis")

# データを無理やりarrayに変換する
xs = df.T.ix[:,0].values
ys = df.T.ix[:,1].values
zs = df.T.ix[:,2].values

# 描画
ax.scatter3D(xs, ys, zs)
plt.show()

３次元プロットをしてみるとこんな感じです．

横からみるとこんな感じです．x1とx2には相関はなさそうですが，x1とy，x2とyには相関がありそうです．

散布図行列を描く

描画をしただけだと感覚的にしかわからないので，定量的に評価します．
今回は散布図行列という可視化をしてみました．

from pandas.tools.plotting import scatter_matrix
scatter_matrix(df.T)
plt.show()

相関係数を計算する

同じく，各変数の相関を定量評価します．

# 相関係数を表示
df.T.corr()
print(df.T.corr())

val        x1        x2         y
val                              
x1   1.000000 -0.170384  0.675085
x2  -0.170384  1.000000  0.603907
y    0.675085  0.603907  1.000000

この通りx1とx2には相関はなさそうですが，x1とy，x2とyには相関がありそうというのがわかります．

重回帰分析の計算をする

ここからが実際の分析です．
今回使ったStatsModelsのドキュメントは下記のとおりです

• StatsModels: Statistics in Python — statsmodels 0.9.0 documentation

import statsmodels.formula.api as sm
reg = "y ~ x1 + x2"
model = sm.ols(formula=reg, data=df.T)

# 回帰分析を実行する
result = model.fit()

==============================================================================
Dep. Variable:                      y   R-squared:                       0.988
Model:                            OLS   Adj. R-squared:                  0.987
Method:                 Least Squares   F-statistic:                     703.1
Date:                Sun, 28 Jan 2018   Prob (F-statistic):           4.53e-17
Time:                        19:34:45   Log-Likelihood:                -25.297
No. Observations:                  20   AIC:                             56.59
Df Residuals:                      17   BIC:                             59.58
Df Model:                           2                                         
Covariance Type:            nonrobust                                         
==============================================================================
                 coef    std err          t      P>|t|      [95.0% Conf. Int.]
------------------------------------------------------------------------------
Intercept      0.7622      0.829      0.919      0.371        -0.988     2.512
x1             1.0248      0.034     29.785      0.000         0.952     1.097
x2             1.0218      0.037     27.524      0.000         0.943     1.100
==============================================================================
Omnibus:                        4.291   Durbin-Watson:                   1.543
Prob(Omnibus):                  0.117   Jarque-Bera (JB):                1.439
Skew:                           0.024   Prob(JB):                        0.487
Kurtosis:                       1.687   Cond. No.                         86.1
==============================================================================

回帰式を使ってもう一回描画する

# 先程の結果から値を獲得
b0,b1,b2= result.params

# もう一回描画する
import numpy as np

# 変数の区間の指定
x = np.arange(0, 21, 3)
y = np.arange(0, 21, 3)

# メッシュ表示
X, Y = np.meshgrid(x, y)

# 回帰式を代入する
Z = b0 + b1*X +b2*Y

# 初期化
fig = plt.figure()
ax = Axes3D(fig)
ax.plot_wireframe(X,Y,Z)

# 軸名を設定する
ax.set_xlabel("X-axis")
ax.set_ylabel("Y-axis")
ax.set_zlabel("Z-axis")

# 表示範囲を指定する
ax.set_xlim(-5, 25)
ax.set_ylim(0, 25)
ax.set_zlim(0, 35)

# データを無理やりarrayに変換する
xs = df.T.ix[:,0].values
ys = df.T.ix[:,1].values
zs = df.T.ix[:,2].values

結構上手くフィティングされているように見えます．

真横から見ると，ピタッと合っている気がします