[통계] 다중공선성(Multicollinearity)의 개념부터 검증까지

| 다중공선성 Multicollinearity 

통계학의 회귀분석에서 독립변수들 간에 강한 상관관계가 나타나는 문제.

왜 문제일까? 

회귀분석은 기본적으로 피처 간의 '독립성'을 전제하기 때문에, 각 피처간의 상관관계가 높으면 분석에 부정적인 영향을 미친다. 독립변수가 서로 의존하게 되면 보통 over-fitting문제가 발생한다. 

어떻게 찾고 해결할까? 

상관관계가 높은 변수를 찾아거나, VIF 검정을 통해 다중공선성을 유발하는 변수를 찾고,

1.제거하거나, 파생변수로 만들거나

2.PCA 기법 사용하여 서로 독립인 새로운 변수들을 생성 🔽 

[ML] 주성분 분석 PCA와 차원의 저주

 

 

Python으로 같이 찾고, 해당 변수를 제거해보도록 하자 

해당 포스팅에서 PCA를 통한 변수축소는 따로 진행하지 않는다.

| 데이터셋 불러오기

git에 올라와있는 보스턴 집값 예측 데이터를 받아온다. 우리의 목표는 MEDV (주택 가격)을 예측하는 것이다.

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns

df=pd.read_csv("https://raw.githubusercontent.com/signature95/tistory/main/dataset/boston.csv")

 

 

 

| 다중공선성 유발 변수 찾기 (1) 변수 간 상관관계 분석 

df.corr()

ax=sns.heatmap(df.corr(),cmap='Blues',annot=True)

 

상관분석 결과, 다중공선성을 유발할 수 있는 컬럼들이 (정말) 많이 보인다.

(0.3이상이면 어느정도 뚜렷한 상관관계가 있다고 보지만, 해당 포스팅에서는 0.6 이상을 기준으로 삼았다.) 

 

변수1	변수2	상관관계
RAD (방사형 도로까지의 접근성 지수)	CRIM (범죄율)	0.63
DIS (직업센터의 거리)	ZN (25,000ft² 초과 거주지역 비율)	0.66
NOX (일산화질소 농도)	INDUS (비소매상업지역 면적 비율)	0.76
AGE (1940년 이전 건축 주택 비율)	INDUS (비소매상업지역 면적 비율)	0.64
RAD (방사형 도로까지의 접근성 지수)	INDUS (비소매상업지역 면적 비율)	0.6
TAX  (재산세율)	INDUS (비소매상업지역 면적 비율)	0.72
LSTAT (인구 중 하위 계층 비율)	INDUS (비소매상업지역 면적 비율)	0.6
AGE  (1940년 이전 건축 주택 비율)	NOX (일산화질소 농도)	0.73
RAD (방사형 도로까지의 접근성 지수)	NOX (일산화질소 농도)	0.61
TAX (재산세율)	NOX (일산화질소 농도)	0.67
MEDV (본인 소유 주택 가격, 중앙값 )	RM (주택당 방 수)	0.7
CAT.MEDV (MEDV 30,000$가  넘는지 여부 1,0)	RM (주택당 방 수)	0.64
LSTAT (인구 중 하위 계층 비율)	AGE  (1940년 이전 건축 주택 비율)	0.6
TAX (재산세율)	RAD (방사형 도로까지의 접근성 지수)	0.91
CAT.MEDV (MEDV 30,000$가  넘는지 여부 1,0)	MEDV  (본인 소유 주택 가격,중앙값)	0.79

 

 

 

| 다중공선성 유발 변수 찾기 (2) VIF (Variance Inflation Fector, 분산확장요인) 

VIF는 다중 회귀 모델에서 독립 변수가 상관 관계가 있는지 측정하는 척도이다. 

VIF가 10이 넘으면 다중공선성이 있다고 판단하고, 5가 넘으면 주의할 필요가 있는 것으로 본다. 

#https://aliencoder.tistory.com/17

def feature_engineering_XbyVIF(df):
    vif = pd.DataFrame()
    vif['VIF_Factor'] = [variance_inflation_factor(df.values,i)
                        for i in range(df.shape[1])]
    vif['Feature'] = df.columns
    return vif 
    
df_vif = df.drop(columns='MEDV',axis=1)
vif=feature_engineering_XbyVIF(df_vif)
vif.sort_values('VIF_Factor',ascending=False).reset_index(drop=True)

 

다중공선성을 유발하는 변수를 모두 지우기 보다는, 변수를 하나하나 제거해가면서 VIF가 어떻게 변하는지 하나하나 씩 확인해보자. 순서대로 RM제거 -> PTRATIO -> NOX 제거 -> TAX 제거 -> MEDV 제거 

 

df_not_rm = df.drop(columns=['MEDV','RM'],axis=1)
vif_not_rm=feature_engineering_XbyVIF(df_not_rm)
vif_not_rm.sort_values('VIF_Factor',ascending=False).reset_index(drop=True)

df_not_PTRATIO = df.drop(columns=['MEDV','RM','PTRATIO'],axis=1)
vif_not_PTRATIO=feature_engineering_XbyVIF(df_not_PTRATIO)
vif_not_PTRATIO.sort_values('VIF_Factor',ascending=False).reset_index(drop=True)

df_not_NOX = df.drop(columns=['MEDV','RM','PTRATIO','NOX'],axis=1)
vif_not_NOX=feature_engineering_XbyVIF(df_not_NOX)
vif_not_NOX.sort_values('VIF_Factor',ascending=False).reset_index(drop=True)

df_not_TAX = df.drop(columns=['MEDV','RM','PTRATIO','NOX','TAX'],axis=1)
vif_not_TAX=feature_engineering_XbyVIF(df_not_TAX)
vif_not_TAX.sort_values('VIF_Factor',ascending=False).reset_index(drop=True)

 

 

 

| 그렇다면 뭘 지우지? 

이렇게 지우다보니 끝이 없다.. 도대체 어디까지 지워야할까? 다중공선성이 높은것만 지워야할까? 

하면 그건 아니다. 다중공선성을 가장 많이 유발하는 RM은 MEDV와 상관이 높기 때문에, 이걸 지우면 예측력이 떨어질 것이다. 따라서 VIF가 가장 높은 변수부터 삭제하기 보다는, 이렇게 종속변수에 영향을 미치는 변수와의 관계 등을 파악해서 선택해야한다. PCA 또한 좋은 선택지가 될 것이다. 

 

이번 포스팅에서 종속변수와 높은 상관을 보이는 RM, LSTAT,PTRATIO 중 한가지는 살린 채로 다중공선성을 최대로 제거하며 예측 전후를 비교해보도록 하겠다.

 

- 여러 비교 결과 TAX,PTRATIO,NOX,AGE,B,DIS,RAD,INDU를 제거해서 선형회귀분석을 실시하도록 한다. (물론 완벽하게 제거하지는 못했다)

PTRATIO,NOX,AGE,B,DIS,RAD,UNUDS 제거

TAX,NOX,AGE,B,DIS,RAD,INDUS 제거

 

RM,NOX,AGE,B,DIS,RAD,INDUS 제거

 

최적

 

| 모델 평가 (1) 모든 변수 사용

모든 변수를 사용했을 때 모델은 어떤 성능을 보일까.

from sklearn.linear_model import LinearRegression 
from sklearn.model_selection import train_test_split
from sklearn.metrics import mean_squared_error
from sklearn.metrics import r2_score

df_basic_train,df_basic_test=train_test_split(df,test_size=0.3)

df_basic_ml_x_train  = df_basic_train.drop(columns='MEDV',axis=1)
df_basic_ml_y_train  = df_basic_train['MEDV']

df_basic_ml_x_test  = df_basic_test.drop(columns='MEDV',axis=1)
df_basic_ml_y_test  = df_basic_test['MEDV']

#test-trian 분리

model=LinearRegression()
model.fit(X=df_basic_ml_x_train,y=df_basic_ml_y_train)

model.coef_
model.intercept_

y_pred=model.predict(X=df_basic_ml_x_test)
mse = mean_squared_error(y_true=df_basic_ml_y_test,y_pred=y_pred)
print('MSE=',mse)
r2_score_baisc= r2_score(y_true=df_basic_ml_y_test,y_pred=y_pred)
print('r2_score=',r2_score_baisc)

MSE= 14.744115761410495
r2_score= 0.8448388442095156

import matplotlib.pyplot as plt
from sklearn.model_selection import learning_curve
from sklearn.model_selection import ShuffleSplit
from sklearn.linear_model import LinearRegression
import numpy as np

def plot_learning_curve(estimator, title, X, y, ylim=None, cv=None, n_jobs=1, train_sizes=np.linspace(.1, 1.0, 5)):
    plt.figure()
    plt.title(title)
    if ylim is not None:
        plt.ylim(*ylim)
    plt.xlabel("Training examples")
    plt.ylabel("Score")
    train_sizes, train_scores, test_scores = learning_curve(
        estimator, X, y, cv=cv, n_jobs=n_jobs, train_sizes=train_sizes)
    train_scores_mean = np.mean(train_scores, axis=1)
    train_scores_std = np.std(train_scores, axis=1)
    test_scores_mean = np.mean(test_scores, axis=1)
    test_scores_std = np.std(test_scores, axis=1)
    plt.grid()

    plt.fill_between(train_sizes, train_scores_mean - train_scores_std,
                     train_scores_mean + train_scores_std, alpha=0.1,
                     color="r")
    plt.fill_between(train_sizes, test_scores_mean - test_scores_std,
                     test_scores_mean + test_scores_std, alpha=0.1, color="g")
    plt.plot(train_sizes, train_scores_mean, 'o-', color="r",
             label="Training score")
    plt.plot(train_sizes, test_scores_mean, 'o-', color="g",
             label="Cross-validation score")

    plt.legend(loc="best")
    return plt

 

 

 

 

| 모델 평가 (2) 다중공선성을 없애도록 일부 컬럼만 추출한 데이터 사용 

df_preprocess = df.drop(columns=['TAX','PTRATIO','NOX','AGE','B','DIS','RAD','INDUS'],axis=1)

df_preprocess_train,df_preprocess_test=train_test_split(df_preprocess,test_size=0.3)

df_preprocess_ml_x_train  = df_preprocess_train.drop(columns='MEDV',axis=1)
df_preprocess_ml_y_train  = df_preprocess_train['MEDV']

df_preprocess_ml_x_test  = df_preprocess_test.drop(columns='MEDV',axis=1)
df_preprocess_ml_y_test  = df_preprocess_test['MEDV']

#test-trian 분리

model=LinearRegression()
model.fit(X=df_preprocess_ml_x_train,y=df_preprocess_ml_y_train)

model.coef_
model.intercept_

y_pred=model.predict(X=df_preprocess_ml_x_test)
mse = mean_squared_error(y_true=df_preprocess_ml_y_test,y_pred=y_pred)
print('MSE=',mse)
# 수정: 변수의 이름을 함수 이름과 다르게 설정합니다.
r2_score_value = r2_score(y_true=df_preprocess_ml_y_test, y_pred=y_pred)
print('r2_score=', r2_score_value)

MSE= 13.25548556560361
r2_score= 0.786896261830193

# 예측 결과와 실제 결과를 산점도로 나타내어 모델의 성능을 시각적으로 확인합니다
plt.scatter(df_preprocess_ml_y_train, train_predictions, label='Training Data')
plt.scatter(df_preprocess_ml_y_test, test_predictions, label='Testing Data')
plt.xlabel('True Values')
plt.ylabel('Predictions')
plt.legend()
plt.show()

plot_learning_curve(model, 'Linear Regression Learning Curve (Housing Prices Raw)', df_preprocess_ml_x_train, df_preprocess_ml_y_train, ylim=(0, 1.01), cv=5, n_jobs=4)

 

 

 

 

| 평과 결과 해석 

1.정확도 비교 

	MSE	R^2
raw data	14.744115761410495	0.8448388442095156
가공 데이터	13.25548556560361	0.786896261830193

 

*MSE = 모델의 예측값과 실제값 간의 평균 제곱 오차 

*R-square = 모델이 데이터를 얼마나 잘 설명하는지 나타내는 지표 

 

comment 

MSE의 경우 가공 데이터가 더 낮지만, R^2의 경우 raw data에서 좀 더 좋은수치가 나왔다. 

 

2.예측값 그래프 확인

comment 

predict 값에 대한 산점도에서 유의미한 차이를 발견하기는 힘들다. 

 

 

3.학습 곡선 확인

 

*Training 곡선 : 훈련 데이터셋에 대한 성능이 표시, 훈련 데이터셋이 커질수록 성능이 어떻게 변하는지 보여준다.

 -> 초기에는 훈련데이터가 적기 때문에 모델이 데이터에 더 적응하기 쉽다. 따라서 처음에는 훈련 성능이 높게 나타날 수 있으나 데이터셋이 커질수록 성능이 안정화되거나 감소할 수 있다.

 

*Cross-Validation Score 곡선 : 검증 데이터셋에 대한 성능이 표시된 곡선, 모델이 얼마나 일반화 되었는지 

-> 초기에는 검증 성능이 개선되지만, 훈련 데이터셋이 특정 규모 이상으로 커지면 더 이상 성능이 향상되지 않거나 감소할 수 있다. 

 

*과적합 여부 확인 : 훈련 성능과 검증 성능 간의 간격 , 간격이 크게 벌어지면 과적합이 발생할 가능성이 잇다. 

 

*데이터 크기에 따른 성능 변화 : x축은 훈련 데이터셋의 크기. 훈련 데이터셋의 크기가 증가함에 따라 성능이 어떻게 변하는지 주목 

 

comment 

raw 데이터로 training 한 경우, 가공한 데이터보다 overfitting될 가능성이 더 높으며 최소한 80~90개 이상의 표본을 확보한 후 모델링을 진행하는 것이 좋다. 오른쪽의 경우 매우 우수한 그래프 모양을 보여주고 있다. 90개 이상의 표본 부터 성능이 안정화 되는 것을 확인할 수 있다. 

 

| 결론 

단순 R^2의 경우, raw데이터로 생성한 모델이 높은 성능을 보여주지만, 그 외의 지표를 보면 다가 아님을 알 수 있다. 

학습 곡선에서 다중공선성 변수를 제거한 모델이 훨씬 안정적이고 overfitting에서 벗어난 성능을 보여주고 있음을 확인할 수 있었다. 그럼 다중공선성을 제거한게 좋은 것이네? 하면 또 그건아니다. 내가 이 모델을 사용할 상황, 선형 회귀가 아닌 다른 방법으로의 모델링 등 여러 방면으로 검토 후 최종 모델을 선택하면 되는 것이다. PCA를 통해 주성분을 파악해서 모델링을 돌리면 더 나은 결과가 나올 수도 있고, raw 데이터에 L2 규제 등을 넣으면 훨씬 안정적인 모델로 변할 수도 있는 것이다. 하지만 회귀분석에서 다중공선성은 굉장히 부정적인 결과를 초래할 수 있기 때문에 인지하여 모델링 시 처리할 수 있도록 해야한다. 

 

 

참고 자료

https://gradient-descent.tistory.com/6

https://dnai-deny.tistory.com/22