重回帰分析をしてみよう-2
ここでは重回帰分析がうまく行かない例を示してみる．しかしながらこのデータはニューラルネットワークによる回帰分析を行えば予測が可能になるデータです．
データの読み込み
まずは，必要なモジュールをインポートする．
モジュールのインポートimport pandas as pd
import numpy as np
import seaborn as sns
import statsmodels.api as sm
import matplotlib.pyplot as plt
# 高解像度ディスプレイ用
from IPython.display import set_matplotlib_formats
# from matplotlib_inline.backend_inline import set_matplotlib_formats # バージョンによってはこちらを有効に
set_matplotlib_formats('retina')
次に，CSVファイルを読み込んで表示してみる．ここで利用するデータは GitHub のリポジトリで公開しているので，プログラムの中で読み込むことができる．このデータは，\(x_1\), \(x_2\), \(x_3\) の値が決まれば何らかの法則で \(y\) の値がおおよそ決まるような500件のデータである．
ファイルを読み込んで表示するurl = "https://github.com/rinsaka/sample-data-sets/blob/master/mra-02.csv?raw=true"
df = pd.read_csv(url)
print(df)
no x1 x2 x3 y 0 0 0.412 3.968 3.083 -0.764 1 1 2.562 4.609 0.683 -0.648 2 2 0.874 4.577 2.704 1.005 3 3 0.869 3.275 4.463 0.200 4 4 1.326 4.764 4.729 0.212 .. ... ... ... ... ... 495 495 3.552 4.549 1.786 -0.806 496 496 4.754 1.685 1.406 1.019 497 497 0.826 4.907 3.840 0.641 498 498 1.948 3.455 1.806 1.088 499 499 2.671 2.491 4.405 -0.398 [500 rows x 5 columns]
相関係数と散布図行列
重回帰分析を行う前に，変数間の相関係数を計算し，散布図行列も描いて確認してみよう．まずは，ここと同じ方法で相関係数を計算する．これには Pandas のデータフレームを NumPy 配列に変換し，これを転置して np.corrcoef 関数に与えると良い．
相関係数# x1, x2, x3, y 列だけとりだして NumPy 配列に変換
xy = df.loc[:,['x1','x2','x3','y']].values
# NumPy配列を転置して相関係数を求める
print(np.corrcoef(xy.T))
[[ 1. 0.02061493 -0.0424868 0.03623328] [ 0.02061493 1. 0.00665035 0.00484668] [-0.0424868 0.00665035 1. 0.06805707] [ 0.03623328 0.00484668 0.06805707 1. ]]
上の結果から，どの変数間にも相関関係がなさそうであることがわかる．
次に，変数ごとの散布図を描いて変数間の関連を確認する．複数の散布図をまとめた散布図行列を描くには seaborn パッケージを使うとよい．このとき Pandas のデータフレームには no 列が含まれているので，この列を削除してから散布図を描く．なお，df.drop でデータフレームの行や列を削除できるが，axis=1とすると列を削除し，axis=0や省略すると行を削除する．
データフレームから散布図行列を描く# no列を削除したデータフレームについて， x1, x2, x3, y それぞれの組み合わせで散布図行列を描く
sns.pairplot(df.drop('no', axis=1))
plt.show()
上の散布図を確認すると，やはり変数間の相関関係はなさそうである．
重回帰分析（失敗例）
ここでは上のようなデータに対して重回帰分析が失敗する例を確認しよう．まずは，データフレームから必要な列を取り出して NumPy 配列に格納する．
x1, x2, x3 列を NumPy 配列 x_data に格納x_data = df.loc[:,['x1', 'x2', 'x3']].values
print(x_data)
[[0.412 3.968 3.083] [2.562 4.609 0.683] [0.874 4.577 2.704] ... [0.826 4.907 3.84 ] [1.948 3.455 1.806] [2.671 2.491 4.405]]
y 列を NumPy 配列 y_data に格納y_data = df.loc[:, 'y'].values
print(y_data)
[-0.764 -0.648 1.005 ．．．（中略）．．． 0.641 1.088 -0.398]
y 切片も含めて重回帰分析をしたいので，定数項を変数xに加える．加えなければ当てはめた直線（平面，超平面）は原点を通ることになる．（つまり，\(x_1 = 0\)，\(x_2 = 0\)，\(x_3 = 0\) のとき，\(y = 0\) となる．）
定数項を加えるx = sm.add_constant(x_data)
モデルを定義する．
モデルの定義model = sm.OLS(y_data, x)
モデルの当てはめ（パラメータ推定）を行う．
パラメータ推定results = model.fit()
結果を表示する．推定はできている．
print(results.params)
[-0.01425276 0.01948354 0.00178876 0.03489884]
詳細な推定結果を表示する．
print(results.summary2())
Results: Ordinary least squares
==================================================================
Model: OLS Adj. R-squared: 0.000
Dependent Variable: y AIC: 1088.6260
Date: 20yy-mm-dd hh:mm BIC: 1105.4844
No. Observations: 500 Log-Likelihood: -540.31
Df Model: 3 F-statistic: 1.028
Df Residuals: 496 Prob (F-statistic): 0.380
R-squared: 0.006 Scale: 0.51243
--------------------------------------------------------------------
Coef. Std.Err. t P>|t| [0.025 0.975]
--------------------------------------------------------------------
const -0.0143 0.1033 -0.1380 0.8903 -0.2172 0.1887
x1 0.0195 0.0223 0.8730 0.3831 -0.0244 0.0633
x2 0.0018 0.0224 0.0799 0.9364 -0.0422 0.0458
x3 0.0349 0.0224 1.5556 0.1205 -0.0092 0.0790
------------------------------------------------------------------
Omnibus: 4428.886 Durbin-Watson: 2.107
Prob(Omnibus): 0.000 Jarque-Bera (JB): 46.960
Skew: -0.048 Prob(JB): 0.000
Kurtosis: 1.502 Condition No.: 15
==================================================================
上の結果を確認すると，決定係数 (R-squared : \(R^2\)) が 0.006，自由度補正済み決定係数 (Adj. R-squared) が 0.000 となっており，全く当てはまっていないことがわかる．また各変数のP値 (P>|t|) も非常に大きな値になっている．したがって，このモデルで予測をしたとしても結果は全く当てにならはいはずです．確認のため，yのデータと重回帰モデルによるyの予測値を散布図としてプロットすると，やはり全く当てにならない予測モデルであることがわかります．
散布図を作成するfig, ax = plt.subplots(1, 1, figsize=(6, 4))
ax.scatter(y_data, results.fittedvalues)
ax.set_xlabel('y_data')
ax.set_ylabel('y_pred')
ax.set_xlim(-1.1,1.3)
ax.set_ylim(-0.02,0.27)
# plt.savefig('mra2.png', dpi=300, facecolor='white')
plt.show()
データの詳細
ここで使用したデータは，重回帰分析に失敗するように恣意的に作成したダミーデータです．このデータに1つの列を加えたデータファイル (mra-02b.csv) を GitHub のリポジトリから読み込み，中身を確認してください．
データの確認url = "https://github.com/rinsaka/sample-data-sets/blob/master/mra-02b.csv?raw=true"
df = pd.read_csv(url)
print(df)
no x1 x2 x3 x123 y 0 0 0.412 3.968 3.083 5.347 -0.764 1 1 2.562 4.609 0.683 11.609 -0.648 2 2 0.874 4.577 2.704 7.499 1.005 3 3 0.869 3.275 4.463 3.130 0.200 4 4 1.326 4.764 4.729 6.390 0.212 .. ... ... ... ... ... ... 495 495 3.552 4.549 1.786 11.574 -0.806 496 496 4.754 1.685 1.406 7.669 1.019 497 497 0.826 4.907 3.840 6.965 0.641 498 498 1.948 3.455 1.806 7.442 1.088 499 499 2.671 2.491 4.405 3.782 -0.398 [500 rows x 6 columns]
上のデータは元のデータに列x123が追加されています．具体的には，\(x_1\)，\(x_2\)，\(x_3\) は最小値0，最大値5の一様乱数です．次に，\(x_{123}\) は次の式で求まるような構造になっています． \begin{eqnarray} x_{123} &=& 1.2 * x_1 + 2 * x_2 - x_3 \end{eqnarray} 例えば，\( \rm{no} = 0 \) のデータは，\( x_1 = 0.412\)，\( x_2 = 3.968 \)，\( x_3 = 3.083 \) であるので， \begin{eqnarray} x_{123} &=& 1.2 * 0.412 + 2 * 3.968 - 3.083 \\ &=& 5.3474 \end{eqnarray} となる．さらに \(y\) は \(y = \sin(x_{123}) \) で決まるので， \begin{eqnarray} y &=& \sin(5.3474) \\ &=& -0.8051 \end{eqnarray} に僅かな誤差を加えたような構造のデータになっています．
ここで，\(x_{123} \) も含めたデータで相関分析を行ってみる．相関係数を確認すると，\(x_{123} \) と \(y \) は無相関（相関係数は -0.00567841）である．
データフレームから相関分析# x1, x2, x3, y 列だけとりだして NumPy 配列に変換
xy = df.loc[:,['x1','x2','x3','x123','y']].values
# NumPy配列を転置して相関係数を求める
print(np.corrcoef(xy.T))
[[ 1. 0.02061493 -0.0424868 0.50067719 0.03623328] [ 0.02061493 1. 0.00665035 0.78412477 0.00484668] [-0.0424868 0.00665035 1. -0.40276532 0.06805707] [ 0.50067719 0.78412477 -0.40276532 1. -0.00567841] [ 0.03623328 0.00484668 0.06805707 -0.00567841 1. ]]
しかしながら，散布図行列を描いてみると，\(x_{123} \) と\(y \) の関連（Sinカーブ）が確認できた．
データフレームから散布図行列を描くsns.pairplot(df.drop('no', axis=1))
plt.show()
散布図から，\( x_{123} \) の値を観測できれば予測ができそうである．しかしながら，単回帰分析や重回帰分析では線形の予測はできるが，非線形のデータでは予測ができません．ニューラルネットワークを使えば，非線形の回帰分析も可能になります．たとえば，\( x_{123} \) の値から \( y \) を推定することもでき，さらに，\( x_{123} \) が観測できなくても，\( x_{1} \)，\( x_{2} \)，\( x_{3} \) から \( y \) を予測できる可能性があります．