重回帰分析をしてみよう-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 \) を予測できる可能性があります．