다운로드
작성자: admin 작성일시: 2017-10-16 20:39:19 조회수: 2513 다운로드: 155
카테고리: 머신 러닝 태그목록:

다중공선성

다중공선성(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)

스캐터 플롯에서 보듯이 독립변수간의 상관관계가 강하다.

In [2]:
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()

상관관계는 상관계수 행렬로도 살펴볼 수 있다.

In [3]:
dfX.corr()
Out:
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
In [4]:
cmap = sns.light_palette("darkgray", as_cmap=True)
sns.heatmap(dfX.corr(), annot=True, cmap=cmap)
plt.show()
In [5]:
model = sm.OLS.from_formula("TOTEMP ~ "
                            "GNPDEFL + POP + GNP + "
                            "YEAR + ARMED + UNEMP", data=df)
result = model.fit()
print(result.summary())
                            OLS Regression Results                            
==============================================================================
Dep. Variable:                 TOTEMP   R-squared:                       0.995
Model:                            OLS   Adj. R-squared:                  0.992
Method:                 Least Squares   F-statistic:                     330.3
Date:                Sat, 30 Jun 2018   Prob (F-statistic):           4.98e-10
Time:                        17:17:07   Log-Likelihood:                -109.62
No. Observations:                  16   AIC:                             233.2
Df Residuals:                       9   BIC:                             238.6
Df Model:                           6                                         
Covariance Type:            nonrobust                                         
==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
Intercept  -3.482e+06    8.9e+05     -3.911      0.004    -5.5e+06   -1.47e+06
GNPDEFL       15.0619     84.915      0.177      0.863    -177.029     207.153
POP           -0.0511      0.226     -0.226      0.826      -0.563       0.460
GNP           -0.0358      0.033     -1.070      0.313      -0.112       0.040
YEAR        1829.1515    455.478      4.016      0.003     798.788    2859.515
ARMED         -1.0332      0.214     -4.822      0.001      -1.518      -0.549
UNEMP         -2.0202      0.488     -4.136      0.003      -3.125      -0.915
==============================================================================
Omnibus:                        0.749   Durbin-Watson:                   2.559
Prob(Omnibus):                  0.688   Jarque-Bera (JB):                0.684
Skew:                           0.420   Prob(JB):                        0.710
Kurtosis:                       2.434   Cond. No.                     4.86e+09
==============================================================================

Warnings:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
[2] The condition number is large, 4.86e+09. This might indicate that there are
strong multicollinearity or other numerical problems.

다중 공선성이 있으면 독립변수의 공분산 행렬의 조건수(conditional number)가 증가한다. 조건수는 독립변수가 스케일링이 되어 있지 않아도 증가하므로 일단 스케일링을 통해 조건수를 감소시킨다.

In [6]:
model = sm.OLS.from_formula("TOTEMP ~ "
                            "scale(GNPDEFL) + scale(POP) + scale(GNP) + "
                            "scale(YEAR) + scale(ARMED) + scale(UNEMP)", data=df)
result = model.fit()
print(result.summary())
                            OLS Regression Results                            
==============================================================================
Dep. Variable:                 TOTEMP   R-squared:                       0.995
Model:                            OLS   Adj. R-squared:                  0.992
Method:                 Least Squares   F-statistic:                     330.3
Date:                Sat, 30 Jun 2018   Prob (F-statistic):           4.98e-10
Time:                        17:17:07   Log-Likelihood:                -109.62
No. Observations:                  16   AIC:                             233.2
Df Residuals:                       9   BIC:                             238.6
Df Model:                           6                                         
Covariance Type:            nonrobust                                         
==================================================================================
                     coef    std err          t      P>|t|      [0.025      0.975]
----------------------------------------------------------------------------------
Intercept       6.532e+04     76.214    857.026      0.000    6.51e+04    6.55e+04
scale(GNPDEFL)   157.3796    887.266      0.177      0.863   -1849.755    2164.514
scale(POP)      -344.1972   1522.652     -0.226      0.826   -3788.675    3100.281
scale(GNP)     -3447.1925   3223.132     -1.070      0.313   -1.07e+04    3844.039
scale(YEAR)     8431.9716   2099.652      4.016      0.003    3682.229    1.32e+04
scale(ARMED)    -696.2102    144.382     -4.822      0.001   -1022.826    -369.594
scale(UNEMP)   -1827.8860    441.900     -4.136      0.003   -2827.533    -828.239
==============================================================================
Omnibus:                        0.749   Durbin-Watson:                   2.559
Prob(Omnibus):                  0.688   Jarque-Bera (JB):                0.684
Skew:                           0.420   Prob(JB):                        0.710
Kurtosis:                       2.434   Cond. No.                         111.
==============================================================================

Warnings:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.

독립 변수가 서로 의존하게 되면 이른바 과최적화(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를 계산한다.

In [7]:
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
Out:
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 세가지 변수만으로도 비슷한 수준의 성능이 나온다는 것을 알 수 있다.

In [8]:
model2 = sm.OLS.from_formula("TOTEMP ~ scale(GNP) + scale(ARMED) + scale(UNEMP)", data=df)
result2 = model2.fit()
print(result2.summary())
                            OLS Regression Results                            
==============================================================================
Dep. Variable:                 TOTEMP   R-squared:                       0.985
Model:                            OLS   Adj. R-squared:                  0.981
Method:                 Least Squares   F-statistic:                     264.4
Date:                Sat, 30 Jun 2018   Prob (F-statistic):           3.19e-11
Time:                        17:17:07   Log-Likelihood:                -119.16
No. Observations:                  16   AIC:                             246.3
Df Residuals:                      12   BIC:                             249.4
Df Model:                           3                                         
Covariance Type:            nonrobust                                         
================================================================================
                   coef    std err          t      P>|t|      [0.025      0.975]
--------------------------------------------------------------------------------
Intercept     6.532e+04    119.824    545.106      0.000    6.51e+04    6.56e+04
scale(GNP)    3925.3853    212.359     18.485      0.000    3462.696    4388.075
scale(ARMED)  -325.2979    171.932     -1.892      0.083    -699.906      49.310
scale(UNEMP)  -720.9526    193.085     -3.734      0.003   -1141.649    -300.257
==============================================================================
Omnibus:                        4.210   Durbin-Watson:                   0.904
Prob(Omnibus):                  0.122   Jarque-Bera (JB):                1.788
Skew:                           0.532   Prob(JB):                        0.409
Kurtosis:                       4.245   Cond. No.                         3.26
==============================================================================

Warnings:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.

보스턴 집값 예측 문제에 응용

In [9]:
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"])
In [10]:
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)
In [11]:
cmap = sns.light_palette("black", as_cmap=True)
sns.heatmap(dfX.corr(), annot=True, fmt='3.1f', cmap=cmap)
plt.show()