重回帰における相関係数と特徴量選択

重回帰を扱えるようになりたかったので以下のことを調べました。

• 多項式特徴量の自動生成
• 特徴量選択の自動化
• 偏相関係数・偏回帰係数の計算

備忘録も兼ねてサンプルコードを投稿します。これから重回帰分析を使うという方の参考になれば幸いです。また、もっと良いアプローチがあれば教えていただけるとありがたいです。

目的

今回は、以下のような回帰式で表すことができる目的変数yを用意しました。

y = 2.5 * {X1の二乗} - 1.5 * {X2とX3の積} + 3.0 * {X4の三乗} + 10

複数の特徴量の中から、目的変数と線形関係にある特徴量の組み合わせを見つけることを目標とします。

また、今回は元の特徴量ではなく、特徴量の累乗や特徴量同士の積が目的変数と線形関係にある、という設定にしました。このような特徴量は多項式特徴量（Polynomial Features）と呼ばれ、sklearnで自動生成が可能らしいので、それも試したいです。

最後に、各特徴量の偏回帰係数と偏相関係数を求めます。重回帰では、特定の特徴量と目的変数の関係を評価するとき、他の特徴量の影響を除外して考える必要があるため、その具体的な方法を調べました。

データセットの用意

まず、正規分布に従う4つの特徴量（X1~X4）を生成します。

import numpy as np
import pandas as pd

num_samples = 100
feature_names = ['X1', 'X2', 'X3', 'X4']

data = np.random.randn(num_samples, len(feature_names))
df = pd.DataFrame(data, columns=feature_names)

ヒストグラムで可視化して、分布を確認しました。

import seaborn as sns
import matplotlib.pyplot as plt

sns.histplot(df['X1'], bins=20)

次に、これらの特徴量から導出できる従属変数yを生成します。

先述の通り、今回は以下の回帰式に従うものとします。

y = 2.5 * {X1の二乗} - 1.5 * {X2とX3の積} + 3.0 * {X4の三乗} + 10

random_error = np.random.randn(num_samples)
df['y'] = 2.5 * (df['X1'] ** 2) - 1.5 * (df['X2'] * df['X3']) + 3.0 * (df['X4'] ** 3) + random_error + 10

生成されたデータセットはこちらです。

index	X1	X2	X3	X4	y
0	-0.5986539369229861	-1.1158969859603944	0.7666631816450861	0.3562928174722889	11.210313760853275
1	-1.7685384506770307	0.35548179274376907	0.8145198224878664	0.05892558918162996	17.354677730171204
2	-0.18505367100934153	-0.8076484876163557	-1.4465346995633879	0.8002979493400275	11.449756831170482
3	-0.3091144447717088	-0.23346666154369272	1.7327211869191332	0.6845011068591904	11.01232981706235
4	0.3708250012811021	0.14206180518723566	1.5199948607657727	1.7195893074161945	24.70784995383167

多項式特徴量の生成

sklearnを使用して多項式特徴量を自動生成します。

多項式特徴量とは、X1の二乗やX2とX3の積のようにオリジナルの特徴量と次元の異なる特徴量のことです。このような特徴量を生成することで、重回帰でも特徴量と目的変数の非線形な関係を捉えることができるようになります。

from sklearn.preprocessing import PolynomialFeatures

poly = PolynomialFeatures(degree=3, include_bias=False)

poly_features = poly.fit_transform(df[feature_names])
poly_feature_names = poly.get_feature_names_out(feature_names)
df_poly = pd.DataFrame(poly_features, columns=poly_feature_names)

degreeで生成される多項式特徴量の最大の次元数を指定できます。今回はdegreeを3に設定したので、X4の三乗も生成されました。

生成されたdf_polyは以下のようになります。（一部省略）

index	X1	X2	X3	X4	X1^2	X1 X2	X1 X3	X1 X4	X2^2	X2 X3	X2 X4	X3^2	X3 X4	X4^2	X1^3	X1^2 X2	X1^2 X3	X1^2 X4	X1 X2^2	X1 X2 X3
0	-0.5986539369229861	-1.1158969859603944	0.7666631816450861	0.3562928174722889	0.3583865361933906	0.6680361238456842	-0.45896593198573316	-0.21329609787716863	1.2452260832754927	-0.855517133644558	-0.39758608113666416	0.5877724340901662	0.2731565850405969	0.12694457178234178	-0.21454951083236554	-0.39992245554699035	0.2747617620967866	0.12769054872447758	-0.7454594971120638	0.5121587001613831
1	-1.7685384506770307	0.35548179274376907	0.8145198224878664	0.05892558918162996	3.1277282515231124	-0.6286832189829586	-1.4405096249084213	-0.10421217019651105	0.12636730497232398	0.2895469667233233	0.020946974080768663	0.6634425412256654	0.04799606044021418	0.003472225060402226	-5.531507676087464	1.1118504460667702	2.54759666022089	0.1843032300210287	-0.2234854377519857	-0.5120749439271
2	-0.18505367100934153	-0.8076484876163557	-1.4465346995633879	0.8002979493400275	0.03424486115403361	0.14945831751854932	0.2676865563965999	-0.14809807342662012	0.6522960795489866	1.1682915623869496	-0.646359428426944	2.0924626370409407	-1.157658753709772	0.6404768077178532	-0.006337137269759115	-0.027657810319687334	-0.049536379941039946	0.027406092157007067	-0.12070978410554144	-0.21619664242894418
3	-0.3091144447717088	-0.23346666154369272	1.7327211869191332	0.6845011068591904	0.0955517399665218	0.07216791745578303	-0.5356091476386842	-0.21158917959239876	0.05450668205235717	-0.4045326308960348	-0.15980818824137766	3.00232271159845	1.1860495703245169	0.4685417652914568	-0.029536423046722084	-0.022308145734674885	0.16556452428698007	0.06540527176940572	-0.01684880275896245	0.1250468795914664
4	0.3708250012811021	0.14206180518723566	1.5199948607657727	1.7195893074161945	0.13751118157512934	0.05268006909055233	0.5636520961907362	0.6376667071255797	0.020181556493056096	0.2159332137957066	0.2442879611922129	2.310384376754361	2.61376690990039	2.9569873861801077	0.0509925840837632	0.019535086687992614	0.20901628929202562	0.23646275748675924	0.00748382571239216	0.08007343428242539

重回帰モデルの適用と特徴量選択

次に、重回帰モデルを適用し、最も決定係数が大きくなる特徴量の組み合わせを見つけます。

基本的な方針としては、特徴量の可能な組み合わせを全て試したいのですが、それだと数が多すぎるので特徴量を最大1〜3つに制限しています。最も決定係数が高い特徴量の組み合わせが、こちらが想定している正解X1^2,X2 X3,X4^3と一致するかを確認します。

特徴量の組み合わせを生成する際は、itertoolsが便利そうです。

from itertools import combinations

subset_list = []
for L in range(1, 4):
    for subset in combinations(df_poly.columns, L):
        subset_list.append(subset)

上記のスクリプトは、最大3つの特徴量を使った全組み合わせパターンを生成します。さきほども書きましたが、今回は特徴量を最大3つに制限していますが、本来はX1~X4すべてを含むパターンも試したほうが良いと思います。

次に、それぞれの特徴量の組み合わせで決定係数を計算し、上位5件を出力してみます。

from sklearn.linear_model import LinearRegression
from sklearn.metrics import r2_score
from sklearn.model_selection import train_test_split

X = df_poly.values
y = df['y'].values

X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=0)

results = []
for subset in subset_list:
    subset_indices = [df_poly.columns.get_loc(col) for col in subset]
    X_train_sub = X_train[:, subset_indices]
    X_test_sub = X_test[:, subset_indices]
    model = LinearRegression().fit(X_train_sub, y_train)
    score = model.score(X_test_sub, y_test)
    results.append((score, subset))

result_df = pd.DataFrame(results, columns=['決定係数', '特徴量']).sort_values(by='決定係数', ascending=False)
result_df.head(5)

出力結果は以下のとおりです。

index	決定係数	特徴量
2648	0.9895386159984231	X1^2,X2 X3,X4^3
2903	0.9750920182539076	X1^2,X2^2 X4,X4^3
2624	0.9745906963871167	X1^2,X2^2,X4^3
2914	0.9744648543587847	X1^2,X2 X3 X4,X4^3
2833	0.9742858631948628	X1^2,X1 X2 X3,X4^3

最も決定係数が高かった特徴量の組み合わせはX1^2,X2 X3,X4^3であることがわかります。これは期待していた正しい特徴量の組み合わせなので、上手くいきました。

しかし、他の特徴量の組み合わせでも高い決定係数が出てしまうことがわかりました。決定係数が上位の組み合わせをよく見ると、どれもX4^3が含まれていることがわかります。このように特定の特徴量が過度に目的変数に影響を与えている場合、決定係数だけを見て特徴量選択をするのはやめたほうが良さそうです。

偏相関係数の計算

特徴量と目的変数の線形関係の強さの指標といえば相関係数ですが、目的変数に影響する特徴量が複数ある場合は解釈が複雑になります。単純に相関係数を計算すると、以下のようになります。

df_poly['y'] = df['y']
correlation_matrix = df_poly.corr()
correlation_matrix = correlation_matrix[['y']].sort_values(by='y', ascending=False)
correlation_matrix.head()

影響の大きい特徴量X4^3は相関係数が大きいので検出できますが、線形関係にあるにも関わらずX1^2とX2 X3は上位にありません。

index	y
y	1.0
X4^3	0.9745094564614605
X4	0.7973728041847046
X3^2 X4	0.4553346579325712
X2^2 X4	0.3632763040276862

このようなときは、偏相関係数を用いると各特徴量と目的変数の関係がより理解しやすくなります。
偏相関係数は、他の特徴量の値を固定されているという条件の下で、特定の特徴量と目的変数の相関係数を求めたものです。偏相関係数はpingouinというライブラリを使用して以下のように計算できるようです。

import pandas as pd
import numpy as np
import pingouin as pg

df_poly['y'] = df['y']
p_corr_results_list = []

for feature in df_poly.columns[:-1]:
    covar_list = df_poly.columns.difference([feature, 'y']).tolist()
    p_corr = pg.partial_corr(data=df_poly, x=feature, y='y', covar=covar_list)
    p_corr['feature'] = feature
    p_corr_results_list.append(p_corr[['feature', 'r', 'CI95%']])

p_corr_df = pd.concat(p_corr_results_list, ignore_index=True)
p_corr_df['abs_r'] = p_corr_df['r'].abs()
p_corr_df = p_corr_df.sort_values(by='abs_r', ascending=False).drop(columns=['abs_r'])

p_corr_df.head()

実行結果は次のようになりました。偏相関係数の絶対値上位5件のみ表示しています。

index	feature	r	CI95%
33	X4^3	0.9878938180354511	[0.98 0.99]
4	X1^2	0.8665786424516077	[0.79 0.92]
9	X2 X3	-0.7770194477726097	[-0.86 -0.66]
2	X3	0.25644568237111043	[0.02 0.47]
18	X1 X2^2	0.2504446978022831	[0.01 0.46]

rが偏相関係数、CI95%は偏相関係数の95%信頼区間です。

これを見ると、目的変数と線形関係にある特徴量X1^2,X2 X3,X4^3がいずれも高い偏相関係数となっていることがわかります。重回帰の特徴量選択では、この偏相関係数も利用すると良さそうです。

偏回帰係数

回帰係数というのは、特定の特徴量が一単位変化したときに目的変数がどの程度変化するのかを表す係数です。特徴量が複数ある重回帰においては、他の特徴量の影響を固定した上での特定の特徴量の影響を示す偏回帰係数が使われます。

y = 2.5 * {X1の二乗} - 1.5 * {X2とX3の積} + 3.0 * {X4の三乗} + 10

上記の回帰式では2.5、-1.5、3.0が偏回帰係数となります。

偏回帰係数は、訓練済みの回帰モデルのcoef_属性で調べることができます。

import pandas as pd
from sklearn.linear_model import LinearRegression

best_features = ['X1^2', 'X2 X3', 'X4^3']
best_feature_indices = [df_poly.columns.get_loc(col) for col in best_features]

X_train_best = X_train[:, best_feature_indices]
X_test_best = X_test[:, best_feature_indices]

model_best = LinearRegression().fit(X_train_best, y_train)

coef_df = pd.DataFrame({
    '特徴量': best_features,
    '偏回帰係数': model_best.coef_
})

coef_df

index	特徴量	偏回帰係数
0	X1^2	2.5855396044907195
1	X2 X3	-1.5411878447743166
2	X4^3	2.969040448316783

このように、期待している偏回帰係数2.5、-1.5、3.0に近い値になることがわかりました。

まとめ

重回帰では、特定の特徴量が目的変数に過度な影響を与えている場合、不適切な特徴量が含まれていても決定係数が大きくなってしまうようです。よって、決定係数だけを見て特徴量選択をするのは危険そうです。

一方で、偏相関係数は特徴量選択に有用そうでした。