[data analysis] 설명변수 선택의 기준

설명변수 선택의 기준

• 결정계수와 수정결정계수
• 잔차평균제곱합(MSE)
• Mallow’s Cp 통계량
• AIC와 BIC

최소자승법(OLS)에 의한 모형의 구축에는 충족되어야 하는 기본 조건들이 있습니다. 그 기본 가정들 중 하나인 설명변수의 독립성은 대부분의 현실 자료에서 충족시키기는 어렵습니다. 일반적으로 회귀모델은 반응변수에 영향을 주는 요인들을 설명변수로 하기 때문에 그들 사이에 어느 정도 상관성이 존재할 것입니다. 그러나 설명변수의 수 증가와 그들 사이의 높은 상관성은 모형 구축에 사용되는 데이터들에 대한 충실도가 증가하여 각 변수의 새로운 값들에 대한 추정의 정확도를 감소시키는 부작용을 발생시킬 수 있습니다. 이러한 현상을 과적합(over fitting)이라 합니다. 이 문제는 다중 공선성(multicollinearity)의 핵심적인 결과로서 추정의 신뢰성을 저하시킬 수 있습니다. 그러므로 반응변수 추정에 영향력을 가진 설명변수들을 적절하게 선택하는 것이 모형 생성에 중요한 부분이 됩니다.

결정계수와 수정결정계수

식 1에서 나타낸 것과 같이 결정계수 R²는 모형으로 추정된 값들에 내재하는 총 변동에 대한 회귀계수에 의한 변동의 비율을 나타냅니다(회귀모형의 평가 참조). 즉, 회귀계수에 의한 변동이 증가한다면 추정값에 대한 회귀계수의 영향이 증가하는 것으로 모형의 신뢰성이 증가한다는 것을 의미합니다. 그러나 결정계수는 설명변수의 수에 비례하는 경향을 보입니다. 이러한 결정계수의 변화는 모델에 영향력이 큰 변수의 첨가 이후에는 안정한 상태로 유지됩니다. 그 부분이 변수를 선택하는 결정기준이 됩니다.

\begin{align}R^2&= \frac{\text{SSReg}}{\text{SST}}\\ & = \frac{\text{SST}-\text{SSE}}{\text{SST}}\\& = 1- \frac{\text{SSE}}{\text{SST}}\\& =1-\frac{\sum^n_{i=1}(y-\hat{y})^2}{\sum^n_{i=1}(y-\bar{y})^2}\end{align}

(식 1)

일반적으로 설명변수가 증가함에 따라 오차합 즉, SSE가 감소하는 현상을 보이기 때문에 설명변수의 수는 결정계수(R²)를 증가시키는 효과가 있습니다. 이러한 효과는 모형에 대한 과대평가를 일으킬 수 있으므로 식 2와 같이 자유도를 고려하여 설명변수 증가로 인한 영향을 수정합니다. 이 식으로 정의된 통계량을 수정된 결정계수(R²_adj, Adjusted R²)라 합니다.

\begin{align} R^2_{text{adj}&  = 1 − (1 − R^2)\frac{n-1}{n-p-1}\\&:\;\text{자료의 크기}\\&p:\;\text{설명변수의 수}\end{align}

(식 2)

수정된 결정계수는 자유도와 설명변수의 수를 고려하기 때문에 많은 변수를 사용하는 회귀모형의 경우 결정계수보다 정확한 판단근거가 됩니다. 즉, 새로운 변수가 첨가될 때 R²_adj 값의 변화가 일어납니다. 그 변화가 안정된 수준에서 가장 적은 변수를 포함하는 모형을 최적 모형으로 선택할 수 있습니다.

잔차평균제곱합(MSE)

오차제곱평균(Mean of Square Error)은 모형에서 추정치와 관찰치 사이의 오차제곱합을 자유도로 나눈 값으로 모형 분산 σ²의 추정치입니다. 만약 반응변수와 높은 상관성을 가진 설명변수가 제거되면 MSE는 증가할 것이지만 상관정도가 낮은 설명 변수의 첨가 또는 삭제는 MSE의 변화를 보이지 않을 것입니다. 다시말하면 MSE의 변화에 따라 모델에서의 변수의 영향정도를 결정할 수 있습니다.

Mallow’s Cp 통계량

Cp 통계량은 자료의 표준화된 잔차제곱합(SSE)의 추정량으로 식 3과 같이 정의 됩니다.

\begin{align}C^p &  = \frac{\text{SSE}_p}{\text{SSE}_\text{full}}\text{df}_\text{full} + 2p + n\\&:\;\text{자료의 크기}\\&p:\;\text{설명변수의 수}\\&\text{SSE}_p:\; \text{선택된 설명변수로 부터의 잔차제곱합} \\ \text{SSE}_\text{full}:\; \text{모든 설명변수로 부터의 잔차제곱합}\end{align}

(식 3)

식 3의 C_p는 SSE_p에 의존합니다.

AIC와 BIC

우도(likelihood)는 특정 조건에서의 확률을 나타내므로 최대우도를 가진 추정량을 반환하는 모델이 최적합 모델이 될 것입니다. 즉, 우도는 모델의 적합도를 판단하는 기준이 되지만 일반적으로 수정된 결정계수나 Cp에서 언급한 것과 같이 매개변수의 수 증가에 따라 증가하는 경향이 있습니다. AIC(Akaike Information Criterion)와 BIC(Schwarz Bayesian Criterion, SBC)는 우도와 설명변수의 수를 고려한 값을 반환합니다. 각각은 식 4와 같이 정의됩니다.

정의에 의하면 AIC와 BIC는 모델의 음의 로그우도에 설명변수의 수를 패널티 요소로 고려하여 계산한 값입니다. 그러므로 AIC, BIC 값은 작을수록 좋은 모델임을 의미합니다.

\begin{align}\text{AIC} &= −2\ln(L) + 2p\\ \text{BIC} &= −2\ln(L) + \log(n)p\\ &n: \;\text{자료의 크기}\\ &p: \;\text{설명변수의 수}\\ &L: \;\text{최대우도}\end{align}

(식 4)

BIC는 n ≥ 8의 경우 log(n) < 2가 됩니다. 그러므로 사용되는 설명변수의 수에 대한 영향을 AIC보다는 더 크게 반영하므로 많은 설명변수들이 참여하는 모델의 경우 BIC가 선호됩니다.

변수 선택을 위한 여러 기준들은 statsmodels.api.OLS() 등의 클래스를 적용하여 생성한 모델의 메서드 또는 속성으로 확인할 수 있습니다. 또한 sklean.linear_model.LinearRegression() 클래스로 모형을 생성하는 경우 R²와 MSE는 각각 sklearn.metrices 모듈의 함수 r2_score()와 mean_squared_error()를 사용하여 계산할 수 있습니다. 그러나 AIC와 BIC는 LinearRegression() 클래스에서 직접적으로 제공하지 않습니다. 대신에 Lasso 모델을 생성하기 위한 linear_model.LassoLarsIC() 클래스로 간접적으로 확인할 수 있습니다.

예 1)

다음은 일정 기간의 여러 기관의 금융자료입니다. sam을 반응변수로 하는 회귀모델을 대해 앞서 소개한 변수 선택의 기준들을 계산합니다.

	kos	kq	kl	ki	WonDol	sam
Date
2023-01-10	2351.0	696.0	14440.0	4885.0	1239.0	60400.0
2023-01-11	2360.0	710.0	14525.0	4875.0	1240.0	60500.0
2023-01-12	2365.0	711.0	14580.0	4860.0	1242.0	60500.0
⋮	⋮	⋮	⋮	⋮	⋮	⋮

위 자료는 다음의 코드들에 의해 사용할 수 있습니다.

import numpy as np 
import pandas as pd 
from sklearn import preprocessing
from sklearn.metrics import mean_squared_error, r2_score
from sklearn.linear_model import LinearRegression, LassoLarsIC
from scipy import stats
import matplotlib.pyplot as plt
import yfinance as yf
import statsmodels.api as sm
from statsmodels.regression.linear_model import OLS
from statsmodels.formula.api import ols

위 자료의 각 설명변수들 사이에 scale 차이를 감소시키기 위해 표준화 합니다. 정제된 자료에 대해 statsmodel.api.OLS() 클래스를 사용하여 회귀모델을 생성합니다.

ind=da2.values[:-1,:-1]
de=da2.values[1:,-1].reshape(-1,1)
final=da2.values[-1, :-1].reshape(1,-1)
indScaler=preprocessing.StandardScaler().fit(ind)
deScaler=preprocessing.StandardScaler().fit(de)
indNor=indScaler.transform(ind)
finalNor=indScaler.transform(final)
deNor=deScaler.transform(de)
reg=OLS(deNor, indNor).fit()

구축된 회귀모형의 다양한 속성(attribute)들로부터 반환된 각 기준의 값들입니다.

n, p=indNor.shape
re=pd.Series([reg.ssr, reg.mse_resid, reg.rsquared, reg.rsquared_adj,reg.llf, reg.aic, reg.bic])
re.index=["SSE","MSE","R2","R2_adj","log-Likelihood","AIC","BIC"]
re

SSE                49.636091
MSE                 0.139427
R2                  0.862504
R2_adj              0.860573
log-Likelihood   -154.095977
AIC               318.191955
BIC               337.636344
dtype: float64

sklearn.linear_model.LiearnRegression() 클래스로 생성한 모델에서 지표 mse, R², R²_adj, BIC 그리고 log-Likelyhood를 계산해 봅니다.

m=LinearRegression().fit(indNor, deNor)
pre=m.predict(indNor)
mse=mean_squared_error(deNor, pre)
r2=m.score(indNor, deNor)
print("mse: %.3f, R2: %.3f" %(mse, r2))

mse: 0.137, R2: 0.863

R²는 모델의 속성 .score()에 의해 확인할 수 있지만 R²_adj는 다음과 같이 데이터의 차원을 고려하여야 계산합니다.

n,p=indNor.shape
R2_adj=1-(1-r2)*(n-1)/(n-p-1)
round(R2_adj, 3)

0.861

BIC를 확인하기 위해 sklearn.linear_model.LassoLarsIC() 클래스를 적용합니다. 이 클래스는 속성 .alphas_, .criterion_을 사용하여 확인할 수 있습니다. 즉, lasso 모델에 적용되는 alpha가 0일 경우는 선형모델인 LinearRegression()에 의한 모델과 같습니다. 이 alpha 값에 대응하는 criterion이 BIC 입니다. 결과적으로 모델.criterion_의 마지막에 있는 값이 선형모델의 BIC가 됩니다.

m_las=LassoLarsIC(criterion="bic").fit(indNor, np.ravel(deNor))
print(m_las.alphas_)

[0.85714684 0.27944597 0.03133171 0.02608305 0.02511504 0.01621895
 0.01442335 0.        ]

bic=m_las.criterion_[-1]
bic.round(3)

337.687

회귀모델은 등분산성을 가정하므로 각 샘플의 잔차의 확률분포는 모델에서 생성되는 전체 잔차의 평균과 분산을 모수로 하는 정규분포를 따른다고 할 수 있습니다. 이 가정을 기반으로 잔차의 우도는 그 정규분포에서의 확률을 의미하며 전체 잔차의 확률의 합을 모델의 우도로 간주할 수 있습니다.

resid=deNor-m.predict(indNor)
pmf=stats.norm.pdf(resid, np.mean(resid), np.std(resid))
log_linkelyhood=np.sum(np.log(pmf))
round(log_linkelyhood, 3)

-154.096

위 계산과정으로 사용자정의함수로 작성하여 사용할 수 있습니다.

def variableSelBase(ind, de):
    n, p=ind.shape
    m=LinearRegression().fit(ind, de)
    pre=m.predict(ind)
    reside=de-pre
    r2=m.score(ind, de)
    r2_adj=1-(1-r2)*(n-1)/(n-p-1)
    mse=mean_squared_error(de, pre)
    pmf=stats.norm.pdf(reside, np.mean(reside), np.std(reside))
    lh=np.sum(np.log(pmf))
    m2=LassoLarsIC(criterion="bic").fit(ind, de)
    bic=m2.criterion_[-1]
    nme=[ "MSE","R2", "R2-adj", "log_Likelihood", "BIC"]
    re=pd.Series([ mse, r2, r2_adj, lh, bic], index=nme)
    return(re)

variableSelBase(indNor, np.ravel(deNor))

MSE                 0.137496
R2                  0.862504
R2-adj              0.860567
log_Likelihood   -154.095977
BIC               337.686765
dtype: float64