[data analysis] 일변량 변수의 선택

일변량 변수의 선택

2개 이상의 설명변수에서 각 변수가 반응변수에 미치는 영향을 고려하여 선택합니다. 즉, 설명변수 한개와 반응변수의 단순선형모델을 생성하여 각 모델의 f 검정 결과를 비교하여 변수를 선태합니다. 이 과정은 sklearn.feature_selection.SelectKBest(function, k=10) 클래스를 적용하여 계산할 수 있습니다. 이 클래스의 인수 function은 서브모듈인 feature_selection에서 제공하는 f_regression() 와 r_regression() 함수입니다. 이 함수는 각 모델의 f-통계량과 pearson's r(상관계수)를 반환하는 함수입니다. 인수 k는 선택되는 최대변수의 수로서 10이 기본값입니다.

예 1)

코스피지수(kos), 코스탁지수(kq), kodex 레버리지(kl), kodex 인버스(ki), 그리고 원달러환율(WonDol)의 일일 시가, 고가, 저가, 종가(o,h,p,c)들을 설명변수로 사용하여 삼성전자(sam)의 일일 종가를 추정하는 회귀모델을 위해 설명변수를 선택합니다.

이 모델에 적합한 설명변수들을 선택하기 위해 sklearn.feature_selection 모듈의 여러 클래스를 적용하여봅니다.

import numpy as np
import pandas as pd
from scipy import stats
from sklearn import linear_model
from sklearn.preprocessing import StandardScaler, MinMaxScaler
import statsmodels.api as sm
from sklearn.metrics import mean_squared_error
from sklearn.model_selection import train_test_split
import yfinance as yf
from sklearn import feature_selection
import matplotlib.pyplot as plt
plt.rcParams['font.family'] ='NanumGothic'
plt.rcParams['axes.unicode_minus'] =False

st=pd.Timestamp(2023,1, 10)
et=pd.Timestamp(2024, 9, 13)
code=["^KS11", "^KQ11", "122630.KS", "114800.KS","KRW=X","005930.KS"]
nme=["kos","kq","kl", "ki", "WonDol","sam" ]
da=pd.DataFrame()
for i, j in zip(nme,code):
    d=yf.download(j,st, et)[["Open","High","Low","Close"]]
    d.columns=[i+"_"+k for k in ["o","h","l","c"]]
    da=pd.concat([da, d], axis=1)
da=da.ffill()
da.columns

Index(['kos_o', 'kos_h', 'kos_l', 'kos_c', 'kq_o', 'kq_h', 'kq_l', 'kq_c',
       'kl_o', 'kl_h', 'kl_l', 'kl_c', 'ki_o', 'ki_h', 'ki_l', 'ki_c',
       'WonDol_o', 'WonDol_h', 'WonDol_l', 'WonDol_c', 'sam_o', 'sam_h',
       'sam_l', 'sam_c'],
      dtype='object')

ind=da.values[:-1,:-1]
de=da.values[1:,-1].reshape(-1,1)
final=da.values[-1, :-1].reshape(1,-1)
[i.shape for i in [ind, de, final]]

[(437, 23), (437, 1), (1, 23)]

위 결과와 같이 자료의 설명변수는 23개로 반응변수와 분리된 구조입니다. 회귀모델은 sklearn.linear_model.LinearRegression() 클래스를 사용하여 생성합니다. 모델 생성전에 설명변수들을 표준화하였습니다. 또한 설명변수의 증가로 인한 다중공선성의 영향은 모델 구축에 사용된 자료(train data)와 사용되지 않은 자료(test data)들 사이에 결과를 비교하는 것으로 판단할 수 있습니다. 그러므로 train_test_split() 함수에 의해 train과 test를 분리합니다.

indScaler=StandardScaler().fit(ind)
indNor=indScaler.transform(ind)
finalNor=indScaler.transform(final)
Xtr, Xte, ytr, yte=train_test_split(indNor, de, test_size=0.3, random_state=3)
[i.shape for i in [Xtr, Xte]]

[(305, 23), (132, 23)]

모든 설명변수를 사용하는 회귀 모델(m_F)와 5개와 10개의 선택된 설명변수들을 적용한 회귀모델(m_5, m_10)을 비교하기 위한 항목으로 $R^2$와 rmse(=$\sqrt{\text{mse}}$) 그리고 BIC를 사용합니다. $R^2$와 rmse는 모델.score() 속성, sklearn.metrics.mean_squared_error() 함수에 의해 계산하고 BIC는 sklearn.linear_model.LassoLarsIC() 클래스를 사용하여 계산합니다.

m_F=linear_model.LinearRegression().fit(Xtr, ytr)
R2_trF=m_F.score(Xtr, ytr)
R2_teF=m_F.score(Xte, yte)
rmse_trF=mean_squared_error(ytr, m_F.predict(Xtr), squared=False)
rmse_teF=mean_squared_error(yte, m_F.predict(Xte), squared=False)
print("train_R2: %.4f, test_R2:%.4f" %(R2_trF, R2_teF))
print("train_rmse: %.2f, test_rmse:%.2f" %(rmse_trF, rmse_teF))

train_R2: 0.9695, test_R2:0.9680
train_rmse: 1161.16, test_rmse:1081.66

las_F=linear_model.LassoLarsIC(criterion="bic").fit(Xtr, np.ravel(ytr))
bic_F=las_F.criterion_[-1]
bic_F.round(2)

5297.27

변수 선택후 $R^2$, rmse, bic를 비교할 것입니다. 모든 계산을 편리하게 하기 위해 이 결과를 계산하기 위한 함수를 작성합니다.

def r_rmse_bic(model, X, y, bicT=True):
    r2=model.score(X,y)
    rmse=mean_squared_error(y, model.predict(X), squared=False)
    if bicT==True:
        las=linear_model.LassoLarsIC(criterion="bic").fit(X, y)
        bic=las.criterion_[-1]
    else:
        bic=np.nan
    return(pd.DataFrame([r2, rmse, bic], index=["R2","RMSE","BIC"]))

위 함수를 적용하여 회귀모델 m_F의 결과를 확인합니다.

tr_re=r_rmse_bic(m_F, Xtr, np.ravel(ytr))
te_re=r_rmse_bic(m_F, Xte, np.ravel(yte), bicT=False)
reF=pd.concat([tr_re, te_re], axis=1).round(3)
reF.columns=["Train","Test"]
reF

	Train	Test
R2	0.970	0.968
RMSE	1161.162	1081.662
BIC	5297.274	NaN

먼저 각 설명변수에 대한 F-통계량을 기준으로 5개의 설명변수를 선택합니다. 각 단일 변수에 대한 F-통계량과 p-value(유의확률)은 sklearn.feature_selection.f_regression() 함수로 확인할 수 있습니다. sklearn의 회귀모델 등에 관련된 모든 클래스와 함수에 인수로 전달되는 자료는2차원 행렬구조이어야 합니다.

변수 선택 과정을 소개하기 위해 각 변수의 이름을 명시할 필요가 있습니다. 이 측면에서 데이터의 구조는 numpy 배열형보다 pandas의 DataFrame 또는 Series 형이 유리합니다. 후자의 두 형태는 numpy의 배열형을 기본으로 하므로 상호 호환됩니다.

inddf=pd.DataFrame(indNor, columns=da.columns[:-1])
dedf=pd.Series(np.ravel(de), name=da.columns[-1])
fstat, pval=feature_selection.f_regression(inddf, dedf)
re=pd.DataFrame([fstat, pval], index=['fstat','p-value']).T
re.head(3).round(3)

	fstat	p-value
0	1050.308	0.0
1	1112.190	0.0
2	1039.036	0.0

위 결과는 각 설명변수와 반응변수의 분산비를 기준으로 F 검정을 실시한 것으로 모든 변수에서 유의한 차이를 나타냅니다. 즉, 반응변수를 설명하기 위한 모든 설명변수들이 효과가 있음을 의미합니다. 이 결과를 내림차순으로 정리하여 지정한 변수의 수만큼 선택하기 위해 SelectKBest() 클래스를 적용합니다. skleran의 변수선택을 위한 클래스의 다양한 속성이나 메소드를 제공합니다.

ksel5_f=feature_selection.SelectKBest(feature_selection.f_regression, k=5).fit(inddf, dedf)
X_new5_f=ksel5_f.fit_transform(inddf, dedf)    #새로운 변수의 생성
print(ksel5_f.pvalues_.round(3))    #유의확률
print(ksel5_f.scores_.round(2))   #통계량

[0. 0. 0. 0. 0. 0. …]
[ 1050.31  1112.19  1039.04  1070.56    46.31    …]

print(ksel_f.get_feature_names_out()) #선택된 변수의 열이름

['kl_h' 'kl_c' 'sam_o' 'sam_h' 'sam_l']

X_new_f.shape #변형된 변수의 모양

(437, 5)

선택된 5개의 설명변수와 반응변수를 훈련(train)세트와 검정(test)세트로 분리하여 훈련세트를 사용하여 회귀모델을 생성합니다. 작성한 함수 r_rmse_bic()를 사용하여 결과를 계산합니다.

Xtr5, Xte5, ytr, yte= train_test_split(X_new5_f, de, test_size=0.3, random_state=3)m_5=linear_model.LinearRegression().fit(Xtr5, ytr)
tr5_re=r_rmse_bic(m_5, Xtr5, np.ravel(ytr))
te5_re=r_rmse_bic(m_5, Xte5, np.ravel(yte), bicT=False)
re_5=pd.concat([tr5_re, te5_re], axis=1).round(3)
re_5.columns=["Train","Test"]
re_5

	Train	Test
R2	0.967	0.966
RMSE	1203.845	1115.328
BIC	5221.112	NaN

위와 같은 방법으로 10개의 설명변수를 선택하고 회귀모델을 생성합니다. 비교항목들에 대한 그 모델의 결과는 다음과 같습니다.

ksel10_f=feature_selection.SelectKBest(feature_selection.f_regression, k=10).fit(inddf, dedf)
X_new10_f=ksel10_f.fit_transform(inddf, dedf)
Xtr10, Xte10, ytr, yte= train_test_split(X_new10_f, de, test_size=0.3, random_state=3)
m_10=linear_model.LinearRegression().fit(Xtr10, ytr)
tr10_re=r_rmse_bic(m_10, Xtr10, np.ravel(ytr))
te10_re=r_rmse_bic(m_10, Xte10, np.ravel(yte), bicT=False)
re_10=pd.concat([tr10_re, te10_re], axis=1).round(3)
re_10.columns=["Train","Test"]
re_10

	Train	Test
R2	0.968	0.966
RMSE	1189.601	1109.225
BIC	5242.597	NaN

그림 1은 각 모델의 $R^2$와 BIC를 비교를 시각화 한 것입니다.

그림 1. 5개와 10개, 그리고 설명변수를 가진 회귀모델

plt.figure(figsize=(7,2))
plt.subplots_adjust(wspace=0.3)
plt.subplot(1,2,1)
plt.scatter([5, 10, 23],[re_5.iloc[0,0], re_10.iloc[0,0], reF.iloc[0,0]], s=30, color="b",label="Train")
plt.scatter([5, 10, 23],[re_5.iloc[0,1], re_10.iloc[0,1], reF.iloc[0,1]], s=30, color="r", label="Test")
plt.xlabel("설명변수의 수")
plt.xticks([5, 10, 23])
plt.grid(True)
plt.legend(loc="best")
plt.title("R2의 변화")
plt.subplot(1,2,2)
plt.scatter([5, 10, 23],[re_5.iloc[2,0], re_10.iloc[2,0], reF.iloc[2,0]], s=30, color="b",label="Train")
plt.xlabel("설명변수의 수")
plt.xticks([5, 10, 23])
plt.grid(True)
plt.legend(loc="best")
plt.title("bic의 변화")
plt.show()

이 결과에 의하면 모든 변수를 사용한 Full-Model에 비해 모델의 정확도를 나타내는 결정계수(R²)는 약간 감소했지만 설명변수의 수가 감소한 경우 훈련과 검정세트 사이의 차이가 감소한 것으로 과적합(다중공선성)이 개선되었음을 알 수 있습니다. 이것은 BIC의 개선으로 이어집니다. 즉, 모델의 적합도에 대해 영향력이 낮은 변수들을 제거함으로서 모델을 단순화 할 수 있습니다. 이 결과는 설명변수의 증가로 인해 나타나는 과적합문제를 감소시킬 수 있습니다. 그러나 개선된 모델의 $R^2$는 감소한 상태로 사용자는 모델의 목적등을 고려하여 어떤 모델을 선택할지에 대한 결정을 하여야 합니다.

다음은 변수들 사이의 상관계수를 기준으로 위와 동일한 과정으로 변수를 선택해 봅니다. sklearn.feature_selection.r_regression(ind,de)를 분류기준 함수로 적용합니다. 이 함수는 설명변수와 반응변수의 Pearson's r(상관계수)를 반환합니다.

print(feature_selection.r_regression(inddf, dedf).round(3))

[0.841  0.848  0.84   0.843  0.31  …]

SelectKBest(--, k) 함수는 위 결과의 반응변수에 대한 각 설명변수간의 상관계수를 기준으로 상위 k개를 선택합니다.

ksel_r=feature_selection.SelectKBest(feature_selection.r_regression, k=5).fit(inddf, dedf)
print(ksel_r.get_feature_names_out())

['kl_h' 'kl_c' 'sam_o' 'sam_h' 'sam_l']

r_regression() 함수로 확인할 수 있는 상관계수는 객체 ksel_r의 속성 .scores를 사용하여 확인할 수 있습니다.

print(ksel_r.scores_.round(3))

[ 0.841  0.848  0.84   0.843  0.31  …]

위 r_regression()를 적용한 결과는 f_regression() 함수를 사용한 결과와 같습니다. 그러므로 선택한 설명변수들로 구성된 모델의 결과 역시 같습니다.

Xtr5r, Xte5r, ytr5, yte5= train_test_split(X_new5_r, de, test_size=0.3, random_state=3)
m_5r=linear_model.LinearRegression().fit(Xtr5r, ytr5)
tr5r_re=r_rmse_bic(m_5r, Xtr5r, np.ravel(ytr5))
te5r_re=r_rmse_bic(m_5r, Xte5r, np.ravel(yte5), bicT=False)
re_5r=pd.concat([tr5r_re, te5r_re], axis=1).round(3)
re_5r.columns=["Train","Test"]
re_5r

	Train	Test
R2	0.967	0.966
RMSE	1203.845	1115.328
BIC	5221.112	NaN