다중공선성¶

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

# 캘리포니아 주택 자료. 위도와 경도처럼 서로 얽힌 변수가 들어 있다.
from sklearn.datasets import fetch_california_housing

housing = fetch_california_housing()
df = pd.DataFrame(housing.data, columns=housing.feature_names)

# 상관행렬은 다중공선성을 훑어보는 첫걸음이다. 다만 쌍끼리의 상관만
# 보므로, 세 변수 이상이 얽힌 경우는 놓친다. 그래서 VIF 가 필요하다.
corr_matrix = df.corr()

# 열지도로 한눈에 본다.
plt.figure(figsize=(10, 8))
plt.imshow(corr_matrix, cmap='coolwarm', vmin=-1, vmax=1)
plt.colorbar(label='Correlation')
plt.xticks(range(len(housing.feature_names)), housing.feature_names, rotation=45)
plt.yticks(range(len(housing.feature_names)), housing.feature_names)
plt.title('Correlation Matrix: Housing Features')
plt.tight_layout()
plt.show()

# 눈으로 놓치기 쉬우므로 0.7 을 넘는 쌍만 따로 찍는다.
print("High Correlations (|r| > 0.7):")
for i in range(len(corr_matrix.columns)):
    for j in range(i+1, len(corr_matrix.columns)):
        if abs(corr_matrix.iloc[i, j]) > 0.7:
            print(f"  {corr_matrix.columns[i]} <-> {corr_matrix.columns[j]}: {corr_matrix.iloc[i, j]:.3f}")

High Correlations (|r| > 0.7):
  AveRooms <-> AveBedrms: 0.848
  Latitude <-> Longitude: -0.925

import numpy as np
from statsmodels.stats.outliers_influence import variance_inflation_factor

# 쌍별 상관은 모두 1/sqrt(m) 인데 완전공선인 자료를 만든다.
rng = np.random.default_rng(0)
m, N = 10, 2000
Z = rng.normal(size=(N, m))
Z = np.column_stack([Z, Z.sum(axis=1) + rng.normal(0, 0.01, N)])   # 마지막 열 = 앞의 합

C = np.corrcoef(Z, rowvar=False)
off = C[~np.eye(m + 1, dtype=bool)]
print(f"변수 {m + 1}개, 쌍별 상관의 최대 절대값 = {np.abs(off).max():.4f}")
print(f"이론값 1/sqrt(m) = {1 / np.sqrt(m):.4f}")
print(f"0.7 을 넘는 쌍의 개수 = {(np.abs(off) > 0.7).sum()}")

vifs = [variance_inflation_factor(Z, j) for j in range(m + 1)]
print(f"\nVIF 최대 = {max(vifs):.1f}  (마지막 열)")
print(f"VIF 중앙값 = {np.median(vifs):.1f}")

변수 11개, 쌍별 상관의 최대 절대값 = 0.3569
이론값 1/sqrt(m) = 0.3162
0.7 을 넘는 쌍의 개수 = 0

VIF 최대 = 98068.6  (마지막 열)
VIF 중앙값 = 9979.2

import statsmodels.api as sm
from statsmodels.stats.outliers_influence import variance_inflation_factor

# 위도와 경도를 함께 넣으면 VIF 가 치솟는다. 캘리포니아가 북서에서 남동으로
# 비스듬히 뻗어 있어 두 좌표가 강하게 얽히기 때문이다.
X = sm.add_constant(df[['MedInc', 'AveRooms', 'AveOccup', 'Latitude', 'Longitude']])

vif_data = pd.DataFrame()
vif_data['Feature'] = X.columns[1:]  # Skip constant
vif_data['VIF'] = [variance_inflation_factor(X.values, i+1) for i in range(X.shape[1]-1)]

print(vif_data)

     Feature       VIF
0     MedInc  1.269059
1   AveRooms  1.248489
2   AveOccup  1.000990
3   Latitude  8.184505
4  Longitude  7.977739

import numpy as np
import statsmodels.api as sm

features = ['MedInc', 'AveRooms', 'AveOccup', 'Latitude', 'Longitude']
y = housing.target
full = sm.OLS(y, sm.add_constant(df[features])).fit()
s2 = full.mse_resid

print(f"{'변수':>10}{'R_j^2':>10}{'VIF':>9}{'공식 SE':>13}{'보고 SE':>13}")
for f in features:
    others = [g for g in features if g != f]
    r2 = sm.OLS(df[f], sm.add_constant(df[others])).fit().rsquared
    Sjj = ((df[f] - df[f].mean()) ** 2).sum()
    se_formula = np.sqrt(s2 / (Sjj * (1 - r2)))
    print(f"{f:>10}{r2:>10.6f}{1 / (1 - r2):>9.4f}"
          f"{se_formula:>13.8f}{full.bse[f]:>13.8f}")

        변수     R_j^2      VIF        공식 SE        보고 SE
    MedInc  0.212015   1.2691   0.00306483   0.00306483
  AveRooms  0.199032   1.2485   0.00233421   0.00233421
  AveOccup  0.000989   1.0010   0.00049790   0.00049790
  Latitude  0.877818   8.1845   0.00692281   0.00692281
 Longitude  0.874651   7.9777   0.00728654   0.00728654

import numpy as np
import pandas as pd
from sklearn.linear_model import LinearRegression
from sklearn.datasets import fetch_california_housing

# 자료 읽기
housing = fetch_california_housing()
df = pd.DataFrame(housing.data, columns=housing.feature_names)

features = ['MedInc', 'AveRooms', 'AveOccup', 'Latitude', 'Longitude']
X = df[features]

# VIF 를 정의대로 직접 구해 본다. 변수 하나를 반응으로 두고 나머지로
# 회귀한 뒤, 그 R^2 에서 1/(1-R^2) 을 계산하는 것이 전부다.
print("Manual VIF Calculation:")
print("=" * 60)

for j, target_feature in enumerate(features):
    # 관심 변수를 나머지 변수들로 회귀한다.
    other_features = [f for f in features if f != target_feature]
    X_j = X[target_feature].values.reshape(-1, 1)
    X_others = X[other_features].values

    # 관심 변수를 나머지 변수들로 회귀한다
    model = LinearRegression()
    model.fit(X_others, X_j.ravel())

    # 그 회귀의 R^2. 1 에 가까울수록 그 변수가 나머지로 거의 설명된다는 뜻이다.
    y_pred_j = model.predict(X_others)
    ss_res = np.sum((X_j.ravel() - y_pred_j) ** 2)
    ss_tot = np.sum((X_j.ravel() - X_j.mean()) ** 2)
    r2_j = 1 - (ss_res / ss_tot)

    # R^2 가 0.9 면 VIF 가 10, 0.99 면 100 이 된다. 계수의 표준오차가
    # 그 제곱근만큼 부풀려진다는 뜻이다.
    vif_j = 1 / (1 - r2_j)

    print(f"{target_feature:12s}:  R² = {r2_j:.4f},  VIF = {vif_j:7.2f}")

print("=" * 60)

Manual VIF Calculation:
============================================================
MedInc      :  R² = 0.2120,  VIF =    1.27
AveRooms    :  R² = 0.1990,  VIF =    1.25
AveOccup    :  R² = 0.0010,  VIF =    1.00
Latitude    :  R² = 0.8778,  VIF =    8.18
Longitude   :  R² = 0.8747,  VIF =    7.98
============================================================

# 대책 1: 얽힌 변수를 빼 버린다. 가장 간단하지만, 뺀 변수가 실제로
# 반응을 설명하고 있었다면 남은 계수에 누락변수 편향이 생긴다.
y = housing.target                                     # 주택가격 중앙값
features_reduced = ['MedInc', 'AveRooms', 'AveOccup']  # 위도·경도를 뺀다
X_reduced = sm.add_constant(df[features_reduced])
model_reduced = sm.OLS(y, X_reduced).fit()

import numpy as np
import statsmodels.api as sm

features_full = ['MedInc', 'AveRooms', 'AveOccup', 'Latitude', 'Longitude']
dropped = ['Latitude', 'Longitude']
model_full = sm.OLS(y, sm.add_constant(df[features_full])).fit()

# 뺀 변수들을 남긴 변수들에 회귀시킨 계수행렬 Delta
Delta = np.column_stack([
    sm.OLS(df[d], X_reduced).fit().params.values for d in dropped
])
predicted = (model_full.params[['const'] + features_reduced].values
             + Delta @ model_full.params[dropped].values)

print(f"{'':>10}{'완전모형':>12}{'축소모형':>12}{'항등식 예측':>14}")
for k, name in enumerate(['const'] + features_reduced):
    print(f"{name:>10}{model_full.params[name]:>12.6f}"
          f"{model_reduced.params[name]:>12.6f}{predicted[k]:>14.6f}")
print(f"\n최대 오차 = {np.abs(model_reduced.params.values - predicted).max():.2e}")
print(f"\nR^2: 완전 {model_full.rsquared:.4f}  →  축소 {model_reduced.rsquared:.4f}")

                  완전모형        축소모형        항등식 예측
     const  -42.829801    0.606932      0.606932
    MedInc    0.359016    0.434683      0.434683
  AveRooms    0.015013   -0.038326     -0.038326
  AveOccup   -0.003366   -0.004174     -0.004174

최대 오차 = 1.58e-14

R^2: 완전 0.5860  →  축소 0.4808

# 대책 2: 얽힌 변수를 하나로 합친다. 여기서는 위도와 경도를 평균 내
# "위치" 지표 하나로 만들었다. 뜻이 통하는 합성이어야 쓸모가 있다.
df['Location'] = (df['Latitude'] + df['Longitude']) / 2

import numpy as np
import statsmodels.api as sm

lat, lon = df['Latitude'], df['Longitude']
r = lat.corr(lon)
print(f"sd(Latitude) = {lat.std():.4f},  sd(Longitude) = {lon.std():.4f},  r = {r:.4f}")
print(f"sd(합/2) 공식 = {np.sqrt((lat.var() + lon.var() + 2 * lat.cov(lon)) / 4):.4f}"
      f"   실측 = {((lat + lon) / 2).std():.4f}")
print(f"sd(차/2) 공식 = {np.sqrt((lat.var() + lon.var() - 2 * lat.cov(lon)) / 4):.4f}"
      f"   실측 = {((lat - lon) / 2).std():.4f}")

df['Diff'] = (lat - lon) / 2
print(f"\n{'모형':>26}{'R^2':>9}{'합성변수의 계수':>16}")
print(f"{'MedInc+AveRooms+AveOccup':>26}{model_reduced.rsquared:>9.4f}{'-':>16}")
for name in ('Location', 'Diff'):
    mm = sm.OLS(y, sm.add_constant(df[features_reduced + [name]])).fit()
    print(f"{'  + ' + name:>26}{mm.rsquared:>9.4f}{mm.params[name]:>16.4f}")
print(f"{'  + Latitude, Longitude':>26}{model_full.rsquared:>9.4f}"
      f"{model_full.params['Latitude'] + model_full.params['Longitude']:>16.4f}")

sd(Latitude) = 2.1360,  sd(Longitude) = 2.0035,  r = -0.9247
sd(합/2) 공식 = 0.4069   실측 = 0.4069
sd(차/2) 공식 = 2.0304   실측 = 2.0304

                        모형      R^2        합성변수의 계수
  MedInc+AveRooms+AveOccup   0.4808               -
                + Location   0.5855         -0.9991
                    + Diff   0.4813         -0.0130
     + Latitude, Longitude   0.5860         -1.0097

from sklearn.linear_model import Ridge, Lasso

# 대책 3: 벌점회귀. 능형회귀는 계수를 0 쪽으로 줄여 분산을 낮춘다.
# 편향이 조금 생기는 대신 분산이 크게 줄어드는 맞바꿈이다.
ridge = Ridge(alpha=1.0)
ridge.fit(X, y)
print(ridge.coef_)

# 라쏘는 일부 계수를 정확히 0 으로 만들어 변수 선택까지 해 준다.
lasso = Lasso(alpha=0.1)
lasso.fit(X, y)
print(lasso.coef_)

[ 0.35902368  0.01500534 -0.00336581 -0.4979378  -0.51160538]
[ 0.37150935 -0.         -0.00281434 -0.18079981 -0.1744082 ]

import numpy as np
from sklearn.linear_model import LinearRegression

# 능형의 닫힌 꼴. sklearn 은 절편을 벌하지 않으므로 중심화한 뒤 푼다.
Xc = X.values - X.values.mean(axis=0)
yc = y - y.mean()
lam = 1.0
b_closed = np.linalg.solve(Xc.T @ Xc + lam * np.eye(X.shape[1]), Xc.T @ yc)
print("닫힌 꼴      :", np.round(b_closed, 8))
print("sklearn Ridge:", np.round(ridge.coef_, 8))
print(f"최대 오차 = {np.abs(b_closed - ridge.coef_).max():.2e}")

# 수축이 실제로 얼마나 일어났는가
b_ols = LinearRegression().fit(X, y).coef_
print(f"\n||b_OLS||   = {np.linalg.norm(b_ols):.6f}")
print(f"||b_Ridge|| = {np.linalg.norm(ridge.coef_):.6f}")
print(f"||b_Lasso|| = {np.linalg.norm(lasso.coef_):.6f}")

# 왜 능형이 거의 아무 일도 하지 않았는가: X'X 의 고윳값과 견주어 lambda 가 너무 작다
ev = np.linalg.eigvalsh(Xc.T @ Xc)[::-1]
print(f"\nX'X 의 고윳값 = {np.array2string(ev, precision=1)}")
print(f"수축비 d/(d+lambda) = {np.array2string(ev / (ev + lam), precision=6)}")

# 벌점은 척도에 의존한다. 표준화하면 결과가 완전히 달라진다.
Xs = (X - X.mean()) / X.std()
print(f"\n표준화 뒤 OLS   : {np.round(LinearRegression().fit(Xs, y).coef_, 4)}")
print(f"표준화 뒤 Ridge : {np.round(Ridge(alpha=1.0).fit(Xs, y).coef_, 4)}")
print(f"표준화 뒤 Lasso : {np.round(Lasso(alpha=0.1).fit(Xs, y).coef_, 4)}")

닫힌 꼴      : [ 0.35902368  0.01500534 -0.00336581 -0.4979378  -0.51160538]
sklearn Ridge: [ 0.35902368  0.01500534 -0.00336581 -0.4979378  -0.51160538]
최대 오차 = 1.11e-16

||b_OLS||   = 0.799364
||b_Ridge|| = 0.799260
||b_Lasso|| = 0.448480

X'X 의 고윳값 = [2226360.5  172828.2  139581.8   59862.2    5540.9]
수축비 d/(d+lambda) = [1.       0.999994 0.999993 0.999983 0.99982 ]

표준화 뒤 OLS   : [ 0.6821  0.0371 -0.035  -1.0637 -1.0252]
표준화 뒤 Ridge : [ 0.6821  0.0371 -0.035  -1.063  -1.0244]
표준화 뒤 Lasso : [ 0.6931 -0.     -0.     -0.011  -0.    ]

from sklearn.decomposition import PCA

# 대책 4: 주성분회귀. 서로 직교하는 성분으로 바꾸므로 다중공선성이
# 원리적으로 사라진다. 대신 성분이 원래 변수의 섞임이라 해석이 어려워진다.
pca = PCA(n_components=3)
X_pca = pca.fit_transform(X)

# 주성분으로 모형을 적합한다
model_pca = LinearRegression()
model_pca.fit(X_pca, y)
print(f"Explained variance ratio: {pca.explained_variance_ratio_}")

Explained variance ratio: [0.85492016 0.06636584 0.05359929]

import numpy as np
from sklearn.decomposition import PCA

# explained_variance 는 공분산행렬의 고윳값이다.
pca_all = PCA().fit(X)
eig = np.linalg.eigvalsh(np.cov(X.values, rowvar=False))[::-1]
print("sklearn explained_variance :", np.round(pca_all.explained_variance_, 6))
print("공분산행렬의 고윳값        :", np.round(eig, 6))
print("비율                       :", np.round(pca_all.explained_variance_ratio_, 6))
print(f"조건수 sqrt(l_1/l_p)       = {np.sqrt(eig[0] / eig[-1]):.3f}")

# 첫 성분이 무엇인지 적재로 확인한다.
print("\n첫 성분의 적재 :", dict(zip(features, np.round(pca_all.components_[0], 4))))
print(f"AveOccup 의 분산 = {X['AveOccup'].var():.4f},  첫 고윳값 = {eig[0]:.4f}")

# 표준화하면 완전히 다른 그림이 된다.
Xs = (X - X.mean()) / X.std()
pca_s = PCA().fit(Xs)
print("\n표준화 뒤 비율 :", np.round(pca_s.explained_variance_ratio_, 6))
print("표준화 뒤 조건수 =",
      f"{np.sqrt(pca_s.explained_variance_[0] / pca_s.explained_variance_[-1]):.3f}")
print("가장 작은 성분의 적재 :",
      dict(zip(features, np.round(pca_s.components_[-1], 4))))

sklearn explained_variance : [107.871529   8.373863   6.763014   2.900441   0.268466]
공분산행렬의 고윳값        : [107.871529   8.373863   6.763014   2.900441   0.268466]
비율                       : [0.85492  0.066366 0.053599 0.022987 0.002128]
조건수 sqrt(l_1/l_p)       = 20.045

첫 성분의 적재 : {'MedInc': 0.0035, 'AveRooms': -0.0012, 'AveOccup': 1.0, 'Latitude': 0.0005, 'Longitude': 0.0005}
AveOccup 의 분산 = 107.8700,  첫 고윳값 = 107.8715

표준화 뒤 비율 : [0.386895 0.265062 0.200103 0.135386 0.012553]
표준화 뒤 조건수 = 5.552
가장 작은 성분의 적재 : {'MedInc': -0.105, 'AveRooms': 0.0963, 'AveOccup': 0.0062, 'Latitude': -0.7051, 'Longitude': -0.6946}