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

다중공선성

다중공선성(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 [1]:
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 [2]:
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 [3]:
cmap = sns.light_palette("darkgray", as_cmap=True)
sns.heatmap(dfX.corr(), annot=True, cmap=cmap)
plt.show()
In [4]:
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:                Tue, 06 Nov 2018   Prob (F-statistic):           4.98e-10
Time:                        09:43:32   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.
/Users/joelkim/anaconda3/lib/python3.6/site-packages/scipy/stats/stats.py:1394: UserWarning: kurtosistest only valid for n>=20 ... continuing anyway, n=16
  "anyway, n=%i" % int(n))

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

In [5]:
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:                Tue, 06 Nov 2018   Prob (F-statistic):           4.98e-10
Time:                        09:43:32   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.
/Users/joelkim/anaconda3/lib/python3.6/site-packages/scipy/stats/stats.py:1394: UserWarning: kurtosistest only valid for n>=20 ... continuing anyway, n=16
  "anyway, n=%i" % int(n))

독립 변수가 서로 의존하게 되면 이른바 과최적화(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 [6]:
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 [7]:
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:                Tue, 06 Nov 2018   Prob (F-statistic):           3.19e-11
Time:                        09:43:32   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.
/Users/joelkim/anaconda3/lib/python3.6/site-packages/scipy/stats/stats.py:1394: UserWarning: kurtosistest only valid for n>=20 ... continuing anyway, n=16
  "anyway, n=%i" % int(n))

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

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