説明変数の組み合わせは回帰モデルの適合に大きく関わります。
そこで組み合わせを求める「itertools.combinations」を使うと上手く組み合わせを抽出することができます。
スポンサーリンク
Pythonコード
データはsklearnのデフォルトで入っているボストン住宅価格データを使用します。
説明変数は13あるので考えられる組み合わせは以下の式から8178通りになります。
![\[ \sum_{r=2}^{13} {}_{13} \mathrm{C} _r = 8178\]](https://latotex.net/wp-content/plugins/lazy-load/images/1x1.trans.gif "Rendered by QuickLaTeX.com")
なお、単回帰の場合（r=1）はモデルの適合に大きく関わらないとみなし、r=2からスタートしています。
#必要なものを宣言
import re
import itertools
import numpy as np
import pandas as pd
from multiprocessing import Pool #並列化
from sklearn.datasets import load_boston
from sklearn import linear_model
from sklearn import model_selection
#csvの出力先
path = r'C:\Users\user\Desktop\comb_LR.csv'
#データの読み込み
boston = load_boston()
input = pd.DataFrame(boston.data, columns=boston.feature_names)
nf = len(input.columns)
def main_LR(itimes):
#結果をいれる空の辞書
dic = {}
#組み合わせ計算
comb = list(itertools.combinations(range(nf), itimes))
ncomb = len(comb)
data = [0]*ncomb
for i in range(ncomb):
element = comb[i]
element_list = np.array(element)
text = '{}'
#進行状況を確認したい場合、以下のコメントアウトを外す
#print('combinations = ', text.format(int_d), '_now_', i, '_in_', ncomb)
#組み合わせの順に説明変数を抽出
data[i] = input.iloc[:, element_list]
X = data[i]
Y = boston.target
#重回帰分析
clf = linear_model.LinearRegression()
clf.fit(X, Y)
Y_predicted = clf.predict(X)
#ピアソンの相関係数を計算
corr = np.corrcoef(Y, Y_predicted)[0,1]
#標準偏差と相関係数を空の辞書に追加
label = text.format(list(input.columns[element_list]))
dic.setdefault(label, []).append(Y_predicted.std())
dic.setdefault(label, []).append(corr)
#itimesごとにデータフレームを作成
output_se = pd.Series(data=dic)
output_df = pd.DataFrame.from_dict(dict(zip(output_se.index, output_se.values)),
orient='index', columns=['Standard Deviation', 'Correlation Coefficient'])
#csvに出力
if itimes == 2:
output_df.to_csv(path)
else:
output_df.to_csv(path, mode='a', header=False)
#目的変数の相関係数・標準偏差に近い順に並べる
def find_nearest(df,target):
arr = df[['Standard Deviation', 'Correlation Coefficient']]
newList = arr - target
sort = np.sum(np.power(newList, 2), axis=1).sort_values()
new_df = df.loc[sort.index]
new_df['difference'] = sort
return new_df
if __name__ == '__main__':
#並列化
p = Pool(4)
p.map(main_LR, range(2,nf+1))
#find_nearest用
df_output = pd.read_csv(path, index_col=0)
print(df_output)
target = np.array([boston.target.std(), 1])
print('Target = ', target)
new = find_nearest(df_output,target)
print(new)
print(new.index[:10])補足説明
print(output_df)をすると、以下のように表示されます。
Standard Deviation Correlation Coefficient
['CRIM', 'ZN'] 4.444457 0.483724
['CRIM', 'INDUS'] 4.844312 0.527243
['CRIM', 'CHAS'] 3.837402 0.417653
['CRIM', 'NOX'] 4.457800 0.485176
['CRIM', 'RM'] 6.764019 0.736179
... ... ...
['INDUS', 'CHAS', 'NOX', 'RM', 'DIS', 'RAD', 'T... 7.849168 0.854284
['INDUS', 'CHAS', 'NOX', 'AGE', 'DIS', 'RAD', '... 7.572242 0.824144
['INDUS', 'CHAS', 'RM', 'AGE', 'DIS', 'RAD', 'T... 7.789401 0.847779
['INDUS', 'NOX', 'RM', 'AGE', 'DIS', 'RAD', 'TA... 7.818189 0.850912
['CHAS', 'NOX', 'RM', 'AGE', 'DIS', 'RAD', 'TAX... 7.849451 0.854314
[8178 rows x 2 columns] しかし、このままではどの組み合わせが目的変数に近いのか分かりません。そこでoutput_dfを目的変数の標準偏差と相関係数に近い順に並べる作業をfind_nearestで行っています。
よって、print('Target = ', target)、print(new)、print(new.index[:10])の結果は以下の通りになります。
Target = [9.18801155 1. ]
Standard Deviation Correlation Coefficient difference
['CRIM', 'ZN', 'INDUS', 'CHAS', 'NOX', 'RM', 'A... 7.907258 0.860606 1.659761
['CRIM', 'ZN', 'INDUS', 'CHAS', 'NOX', 'RM', 'D... 7.907250 0.860605 1.659781
['CRIM', 'ZN', 'CHAS', 'NOX', 'RM', 'AGE', 'DIS... 7.906943 0.860572 1.660576
['CRIM', 'ZN', 'CHAS', 'NOX', 'RM', 'DIS', 'RAD... 7.906935 0.860571 1.660597
['CRIM', 'ZN', 'INDUS', 'NOX', 'RM', 'AGE', 'DI... 7.879846 0.857623 1.731568
... ... ... ...
['AGE', 'DIS'] 3.491647 0.380022 32.832943
['DIS', 'B'] 3.396905 0.369711 33.934177
['CHAS', 'B'] 3.395024 0.369506 33.956222
['ZN', 'DIS'] 3.314265 0.360716 34.909584
['CHAS', 'DIS'] 2.947057 0.320750 39.410894
[8178 rows x 3 columns]
Index(['['CRIM', 'ZN', 'INDUS', 'CHAS', 'NOX', 'RM', 'AGE', 'DIS', 'RAD', 'TAX', 'PTRATIO', 'B', 'LSTAT']',
'['CRIM', 'ZN', 'INDUS', 'CHAS', 'NOX', 'RM', 'DIS', 'RAD', 'TAX', 'PTRATIO', 'B', 'LSTAT']',
'['CRIM', 'ZN', 'CHAS', 'NOX', 'RM', 'AGE', 'DIS', 'RAD', 'TAX', 'PTRATIO', 'B', 'LSTAT']',
'['CRIM', 'ZN', 'CHAS', 'NOX', 'RM', 'DIS', 'RAD', 'TAX', 'PTRATIO', 'B', 'LSTAT']',
'['CRIM', 'ZN', 'INDUS', 'NOX', 'RM', 'AGE', 'DIS', 'RAD', 'TAX', 'PTRATIO', 'B', 'LSTAT']',
'['CRIM', 'ZN', 'INDUS', 'NOX', 'RM', 'DIS', 'RAD', 'TAX', 'PTRATIO', 'B', 'LSTAT']',
'['CRIM', 'ZN', 'NOX', 'RM', 'AGE', 'DIS', 'RAD', 'TAX', 'PTRATIO', 'B', 'LSTAT']',
'['CRIM', 'ZN', 'NOX', 'RM', 'DIS', 'RAD', 'TAX', 'PTRATIO', 'B', 'LSTAT']',
'['CRIM', 'ZN', 'INDUS', 'CHAS', 'NOX', 'RM', 'AGE', 'DIS', 'RAD', 'PTRATIO', 'B', 'LSTAT']',
'['CRIM', 'ZN', 'INDUS', 'CHAS', 'NOX', 'RM', 'DIS', 'RAD', 'PTRATIO', 'B', 'LSTAT']'],
dtype='object') 最初に説明したように、説明変数が多いほど目的変数の標準偏差と相関係数に近いことが読み取れます。nC_2の組み合わせなら、上位の組み合わせだけを散布図にして標準偏差と相関係数をその図にannotateしても良いですし、
今回のような重回帰分析ではTaylor Diagramを使って、出力したcsvファイルを１つの図に可視化することも可能です！
追記（2021/02/05）：Pool数について
今回のコードでは「p = Pool(4)」で全通りが重複なく出力されましたが、場合によっては重複した組み合わせが発生し、想定した組み合わせの数以上が出力される場合があります。その場合は、Pool数は説明変数の数と同じ数（今回の場合であれば「p = Pool(13)」）にすることで解決できます。
参考にしたサイト様
ある値に近い順にarrayを並び変える方法
https://stackoverflow.com/questions/22466273/sorting-a-2d-numpy-array-on-to-the-proximity-of-each-element-to-a-certain-point