6.4 다중공선성과 변수 선택¶
다중공선성(multicollinearity)란 독립 변수의 일부가 다른 독립 변수의 조합으로 표현될 수 있는 경우이다. 독립 변수들이 서로 독립이 아니라 상호상관관계가 강한 경우에 발생한다. 이는 독립 변수의 공분산 행렬이 full rank 이어야 한다는 조건을 침해한다.
다음 데이터는 미국의 거시경제지표를 나타낸 것이다.
TOTEMP - Total Employment
GNPDEFL - GNP deflator
GNP - GNP
UNEMP - Number of unemployed
ARMED - Size of armed forces
POP - Population
YEAR - Year (1947 - 1962)
스캐터 플롯에서 보듯이 독립변수간의 상관관계가 강하다.
from statsmodels.datasets.longley import load_pandas
dfy = load_pandas().endog
dfX = load_pandas().exog
df = pd.concat([dfy, dfX], axis=1)
sns.pairplot(dfX)
plt.show()
상관관계는 상관계수 행렬로도 살펴볼 수 있다.
dfX.corr()
| GNPDEFL | GNP | UNEMP | ARMED | POP | YEAR | |
|---|---|---|---|---|---|---|
| GNPDEFL | 1.000000 | 0.991589 | 0.620633 | 0.464744 | 0.979163 | 0.991149 |
| GNP | 0.991589 | 1.000000 | 0.604261 | 0.446437 | 0.991090 | 0.995273 |
| UNEMP | 0.620633 | 0.604261 | 1.000000 | -0.177421 | 0.686552 | 0.668257 |
| ARMED | 0.464744 | 0.446437 | -0.177421 | 1.000000 | 0.364416 | 0.417245 |
| POP | 0.979163 | 0.991090 | 0.686552 | 0.364416 | 1.000000 | 0.993953 |
| YEAR | 0.991149 | 0.995273 | 0.668257 | 0.417245 | 0.993953 | 1.000000 |
cmap = sns.light_palette("darkgray", as_cmap=True)
sns.heatmap(dfX.corr(), annot=True, cmap=cmap)
plt.show()
다중 공선성이 있으면 독립변수의 공분산 행렬의 조건수(conditional number)가 증가한다.
from sklearn.model_selection import train_test_split
def get_model1(seed):
df_train, df_test = train_test_split(df, test_size=0.5, random_state=seed)
model = sm.OLS.from_formula("TOTEMP ~ GNPDEFL + POP + GNP + YEAR + ARMED + UNEMP", data=df_train)
return df_train, df_test, model.fit()
df_train, df_test, result1 = get_model1(3)
print(result1.summary())
OLS Regression Results
==============================================================================
Dep. Variable: TOTEMP R-squared: 1.000
Model: OLS Adj. R-squared: 0.997
Method: Least Squares F-statistic: 437.5
Date: Sun, 23 Jun 2019 Prob (F-statistic): 0.0366
Time: 18:13:38 Log-Likelihood: -44.199
No. Observations: 8 AIC: 102.4
Df Residuals: 1 BIC: 103.0
Df Model: 6
Covariance Type: nonrobust
==============================================================================
coef std err t P>|t| [0.025 0.975]
------------------------------------------------------------------------------
Intercept -1.235e+07 2.97e+06 -4.165 0.150 -5e+07 2.53e+07
GNPDEFL 106.2620 75.709 1.404 0.394 -855.708 1068.232
POP 2.2959 0.725 3.167 0.195 -6.915 11.506
GNP -0.3997 0.120 -3.339 0.185 -1.920 1.121
YEAR 6300.6231 1498.900 4.203 0.149 -1.27e+04 2.53e+04
ARMED -0.2450 0.402 -0.609 0.652 -5.354 4.864
UNEMP -6.3311 1.324 -4.782 0.131 -23.153 10.491
==============================================================================
Omnibus: 0.258 Durbin-Watson: 1.713
Prob(Omnibus): 0.879 Jarque-Bera (JB): 0.304
Skew: 0.300 Prob(JB): 0.859
Kurtosis: 2.258 Cond. No. 2.01e+10
==============================================================================
Warnings:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
[2] The condition number is large, 2.01e+10. This might indicate that there are
strong multicollinearity or other numerical problems.
또한 다음처럼 학습용 데이터와 검증용 데이터로 나누어 회귀분석 성능을 비교하면 과최적화가 발생하였음을 알 수 있다.
def calc_r2(df_test, result):
target = df.loc[df_test.index].TOTEMP
predict_test = result.predict(df_test)
RSS = ((predict_test - target)**2).sum()
TSS = ((target - target.mean())**2).sum()
return 1 - RSS / TSS
test1 = []
for i in range(10):
df_train, df_test, result = get_model1(i)
test1.append(calc_r2(df_test, result))
test1
[0.9815050656837723,
0.9738497543069347,
0.9879366369871746,
0.7588861967897188,
0.980720608930437,
0.8937889315168234,
0.8798563810651999,
0.9314665778963799,
0.8608525682180641,
0.9677198735170137]
독립변수가 서로 의존하게 되면 이렇게 과최적화(over-fitting) 문제가 발생하여 회귀 결과의 안정성을 해칠 가능성이 높아진다. 이를 방지하는 방법들은 다음과 같다.
변수 선택법으로 의존적인 변수 삭제
PCA(principal component analysis) 방법으로 의존적인 성분 삭제
정규화(regularized) 방법 사용
VIF¶
다중 공선성을 없애는 가장 기본적인 방법은 다른 독립변수에 의존하는 변수를 없애는 것이다. 가장 의존적인 독립변수를 선택하는 방법으로는 VIF(Variance Inflation Factor)를 사용할 수 있다. VIF는 독립변수를 다른 독립변수로 선형회귀한 성능을 나타낸 것이다. \(i\)번째 변수의 VIF는 다음과 같이 계산한다.
\[ \text{VIF}_i = \frac{\sigma^2}{(n-1)\text{Var}[X_i]}\cdot \frac{1}{1-R_i^2} \]
여기에서 \(R^2_i\)는 다른 변수로 \(i\)번째 변수를 선형회귀한 성능(결정 계수)이다. 다른 변수에 의존적일 수록 VIF가 커진다.
StatsModels에서는 variance_inflation_factor 명령으로 VIF를 계산한다.
from statsmodels.stats.outliers_influence import variance_inflation_factor
vif = pd.DataFrame()
vif["VIF Factor"] = [variance_inflation_factor(
dfX.values, i) for i in range(dfX.shape[1])]
vif["features"] = dfX.columns
vif
| VIF Factor | features | |
|---|---|---|
| 0 | 12425.514335 | GNPDEFL |
| 1 | 10290.435437 | GNP |
| 2 | 136.224354 | UNEMP |
| 3 | 39.983386 | ARMED |
| 4 | 101193.161993 | POP |
| 5 | 84709.950443 | YEAR |
상관계수와 VIF를 사용하여 독립 변수를 선택하면 GNP, ARMED, UNEMP 세가지 변수만으로도 비슷한 수준의 성능이 나온다는 것을 알 수 있다.
def get_model2(seed):
df_train, df_test = train_test_split(df, test_size=0.5, random_state=seed)
model = sm.OLS.from_formula("TOTEMP ~ scale(GNP) + scale(ARMED) + scale(UNEMP)", data=df_train)
return df_train, df_test, model.fit()
df_train, df_test, result2 = get_model2(3)
print(result2.summary())
OLS Regression Results
==============================================================================
Dep. Variable: TOTEMP R-squared: 0.989
Model: OLS Adj. R-squared: 0.981
Method: Least Squares F-statistic: 118.6
Date: Sun, 23 Jun 2019 Prob (F-statistic): 0.000231
Time: 18:13:39 Log-Likelihood: -57.695
No. Observations: 8 AIC: 123.4
Df Residuals: 4 BIC: 123.7
Df Model: 3
Covariance Type: nonrobust
================================================================================
coef std err t P>|t| [0.025 0.975]
--------------------------------------------------------------------------------
Intercept 6.538e+04 163.988 398.686 0.000 6.49e+04 6.58e+04
scale(GNP) 4338.7051 406.683 10.669 0.000 3209.571 5467.839
scale(ARMED) -812.1407 315.538 -2.574 0.062 -1688.215 63.933
scale(UNEMP) -1373.0426 349.316 -3.931 0.017 -2342.898 -403.187
==============================================================================
Omnibus: 0.628 Durbin-Watson: 2.032
Prob(Omnibus): 0.731 Jarque-Bera (JB): 0.565
Skew: 0.390 Prob(JB): 0.754
Kurtosis: 1.958 Cond. No. 4.77
==============================================================================
Warnings:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
다중공선성을 제거한 경우에는 학습 성능과 검증 성능간의 차이가 줄어들었음을 확인할 수 있다. 즉, 과최적화가 발생하지 않는다.
test2 = []
for i in range(10):
df_train, df_test, result = get_model2(i)
test2.append(calc_r2(df_test, result))
test2
[0.9763608388904907,
0.9841984331185702,
0.9687069366140136,
0.9397304053201785,
0.9773357061188462,
0.9561262155732316,
0.980385249669863,
0.9917361722470804,
0.9837134067639467,
0.9789512977093212]
plt.subplot(121)
plt.plot(test1, 'ro', label="검증 성능")
plt.hlines(result1.rsquared, 0, 9, label="학습 성능")
plt.legend()
plt.xlabel("시드값")
plt.ylabel("성능(결정계수)")
plt.title("다중공선성 제거 전")
plt.ylim(0.5, 1.2)
plt.subplot(122)
plt.plot(test2, 'ro', label="검증 성능")
plt.hlines(result2.rsquared, 0, 9, label="학습 성능")
plt.legend()
plt.xlabel("시드값")
plt.ylabel("성능(결정계수)")
plt.title("다중공선성 제거 후")
plt.ylim(0.5, 1.2)
plt.suptitle("다중공선성 제거 전과 제거 후의 성능 비교", y=1.04)
plt.tight_layout()
plt.show()
보스턴 집값 예측 문제에 응용¶
from sklearn.datasets import load_boston
boston = load_boston()
dfX0 = pd.DataFrame(boston.data, columns=boston.feature_names)
from patsy import dmatrix
formula = "scale(CRIM) + scale(I(CRIM ** 2)) + " + \
"scale(ZN) + scale(I(ZN ** 2)) + scale(INDUS) + " + \
"scale(NOX) + scale(RM) + scale(AGE) + " + \
"scale(np.log(DIS)) + scale(RAD) + scale(TAX) + " + \
"scale(np.log(PTRATIO)) + scale(B) + scale(np.log(LSTAT)) + CHAS"
dfX = dmatrix(formula, dfX0, return_type="dataframe")
dfy = pd.DataFrame(boston.target, columns=["MEDV"])
idx_outlier = \
np.array([7, 54, 148, 152, 160, 214, 253, 267, 364, 365, 367, 368, 369,
371, 372, 374, 380, 385, 397, 398, 399, 400, 401, 405, 409, 410,
412, 413, 414, 415, 416, 418, 419, 426, 445, 489, 490, 492, 505,
161, 162, 163, 166, 186, 195, 204, 225, 257, 267, 283, 368, 369,
370, 371, 372])
idx = list(set(range(len(dfX))).difference(idx_outlier))
dfX = dfX.iloc[idx, :].reset_index(drop=True)
dfy = dfy.iloc[idx, :].reset_index(drop=True)
cmap = sns.light_palette("black", as_cmap=True)
sns.heatmap(dfX.corr(), annot=True, fmt='3.1f', cmap=cmap)
plt.show()
vif = pd.DataFrame()
vif["VIF Factor"] = [variance_inflation_factor(
dfX.values, i) for i in range(dfX.shape[1])]
vif["features"] = dfX.columns
vif = vif.sort_values("VIF Factor").reset_index(drop=True)
vif
| VIF Factor | features | |
|---|---|---|
| 0 | 1.061624 | CHAS |
| 1 | 1.338325 | scale(B) |
| 2 | 1.478553 | Intercept |
| 3 | 1.780320 | scale(np.log(PTRATIO)) |
| 4 | 2.596496 | scale(RM) |
| 5 | 3.748931 | scale(AGE) |
| 6 | 3.807459 | scale(INDUS) |
| 7 | 4.682812 | scale(np.log(LSTAT)) |
| 8 | 5.071802 | scale(NOX) |
| 9 | 5.215025 | scale(np.log(DIS)) |
| 10 | 9.107858 | scale(TAX) |
| 11 | 10.218588 | scale(I(CRIM ** 2)) |
| 12 | 11.254736 | scale(RAD) |
| 13 | 11.751869 | scale(I(ZN ** 2)) |
| 14 | 14.646056 | scale(ZN) |
| 15 | 21.260182 | scale(CRIM) |
model_boston1 = sm.OLS(np.log(dfy), dfX)
result_boston1 = model_boston1.fit()
print(result_boston1.summary())
OLS Regression Results
==============================================================================
Dep. Variable: MEDV R-squared: 0.872
Model: OLS Adj. R-squared: 0.868
Method: Least Squares F-statistic: 199.9
Date: Sun, 23 Jun 2019 Prob (F-statistic): 1.56e-185
Time: 18:13:41 Log-Likelihood: 317.45
No. Observations: 456 AIC: -602.9
Df Residuals: 440 BIC: -536.9
Df Model: 15
Covariance Type: nonrobust
==========================================================================================
coef std err t P>|t| [0.025 0.975]
------------------------------------------------------------------------------------------
Intercept 3.0338 0.007 433.880 0.000 3.020 3.048
scale(CRIM) -0.3471 0.044 -7.976 0.000 -0.433 -0.262
scale(I(CRIM ** 2)) 0.3075 0.071 4.331 0.000 0.168 0.447
scale(ZN) -0.0465 0.022 -2.110 0.035 -0.090 -0.003
scale(I(ZN ** 2)) 0.0440 0.020 2.206 0.028 0.005 0.083
scale(INDUS) 0.0037 0.012 0.323 0.747 -0.019 0.026
scale(NOX) -0.0652 0.013 -5.001 0.000 -0.091 -0.040
scale(RM) 0.0999 0.011 9.195 0.000 0.079 0.121
scale(AGE) -0.0273 0.011 -2.438 0.015 -0.049 -0.005
scale(np.log(DIS)) -0.1008 0.014 -7.368 0.000 -0.128 -0.074
scale(RAD) 0.1634 0.020 8.106 0.000 0.124 0.203
scale(TAX) -0.0934 0.018 -5.153 0.000 -0.129 -0.058
scale(np.log(PTRATIO)) -0.0699 0.008 -8.872 0.000 -0.085 -0.054
scale(B) 0.0492 0.007 6.699 0.000 0.035 0.064
scale(np.log(LSTAT)) -0.1487 0.013 -11.074 0.000 -0.175 -0.122
CHAS 0.0659 0.026 2.580 0.010 0.016 0.116
==============================================================================
Omnibus: 28.653 Durbin-Watson: 1.309
Prob(Omnibus): 0.000 Jarque-Bera (JB): 43.266
Skew: 0.465 Prob(JB): 4.03e-10
Kurtosis: 4.188 Cond. No. 35.2
==============================================================================
Warnings:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
cols = ["Intercept", "CHAS", "scale(B)", "scale(CRIM)",
"scale(np.log(PTRATIO))", "scale(RM)", "scale(np.log(LSTAT))"]
model_boston2 = sm.OLS(np.log(dfy), dfX[cols])
result_boston2 = model_boston2.fit()
print(result_boston2.summary())
OLS Regression Results
==============================================================================
Dep. Variable: MEDV R-squared: 0.836
Model: OLS Adj. R-squared: 0.834
Method: Least Squares F-statistic: 380.7
Date: Sun, 23 Jun 2019 Prob (F-statistic): 1.42e-172
Time: 18:13:41 Log-Likelihood: 260.52
No. Observations: 456 AIC: -507.0
Df Residuals: 449 BIC: -478.2
Df Model: 6
Covariance Type: nonrobust
==========================================================================================
coef std err t P>|t| [0.025 0.975]
------------------------------------------------------------------------------------------
Intercept 3.0192 0.007 445.252 0.000 3.006 3.033
CHAS 0.0884 0.028 3.141 0.002 0.033 0.144
scale(B) 0.0558 0.008 6.989 0.000 0.040 0.072
scale(CRIM) -0.1179 0.013 -9.120 0.000 -0.143 -0.092
scale(np.log(PTRATIO)) -0.0508 0.007 -6.936 0.000 -0.065 -0.036
scale(RM) 0.1153 0.011 10.828 0.000 0.094 0.136
scale(np.log(LSTAT)) -0.1570 0.011 -14.179 0.000 -0.179 -0.135
==============================================================================
Omnibus: 29.141 Durbin-Watson: 1.113
Prob(Omnibus): 0.000 Jarque-Bera (JB): 42.637
Skew: 0.483 Prob(JB): 5.51e-10
Kurtosis: 4.145 Cond. No. 5.91
==============================================================================
Warnings:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.