다중회귀 진단¶

import numpy as np
import pandas as pd
import statsmodels.api as sm
from statsmodels.api import OLS, add_constant
from statsmodels.stats.outliers_influence import variance_inflation_factor
from sklearn.datasets import fetch_california_housing

housing = fetch_california_housing()
df = pd.DataFrame(housing.data, columns=housing.feature_names)
df['PRICE'] = housing.target

# 설명변수 셋으로 시작한다. 뒤에서 변수를 늘려 가며 견줄 것이다.
features = ['MedInc', 'AveRooms', 'AveOccup']
X = add_constant(df[features])
y = df['PRICE']

model = OLS(y, X).fit()

# 분산팽창인자(VIF)는 그 변수를 나머지 변수들로 회귀했을 때의 R^2 로
# 정해진다. VIF = 1/(1-R^2) 이므로, 다른 변수들로 잘 설명될수록 커진다.
# 보통 5 나 10 을 넘으면 다중공선성을 의심한다. 상수항은 셈에서 뺀다.
vif_data = pd.DataFrame()
vif_data['Feature'] = X.columns[1:]
vif_data['VIF'] = [variance_inflation_factor(X.values, i)
                   for i in range(1, X.shape[1])]
print(vif_data)

    Feature       VIF
0    MedInc  1.120166
1  AveRooms  1.119797
2  AveOccup  1.000488

import numpy as np

R = df[features].corr().values
print("상관행렬 R:")
print(np.round(R, 4))
print("\nR^-1 의 대각 =", np.round(np.diag(np.linalg.inv(R)), 6))
print("variance_inflation_factor =",
      np.round([variance_inflation_factor(X.values, i) for i in range(1, X.shape[1])], 6))
print(f"\n이 모형의 R^2 = {model.rsquared:.4f}")

상관행렬 R:
[[ 1.      0.3269  0.0188]
 [ 0.3269  1.     -0.0049]
 [ 0.0188 -0.0049  1.    ]]

R^-1 의 대각 = [1.120166 1.119797 1.000488]
variance_inflation_factor = [1.120166 1.119797 1.000488]

이 모형의 R^2 = 0.4808

y_pred = model.predict(X)
residuals = y - y_pred

# 왼쪽은 등분산성, 오른쪽은 정규성을 본다. 이 자료는 집값이 50만 달러에서
# 잘려 있어 잔차 그림 오른쪽에 뚜렷한 사선이 생긴다.
import matplotlib.pyplot as plt
fig, axes = plt.subplots(1, 2, figsize=(14, 5))
axes[0].scatter(y_pred, residuals, alpha=0.3, s=10)
axes[0].axhline(y=0, color='red', linestyle='--')
axes[0].set_xlabel('Fitted Values')
axes[0].set_ylabel('Residuals')

# line='45' 는 기울기 1 의 기준선이다. 표준화하지 않은 잔차에 쓰면
# 점들이 그 선에서 벗어나 보이므로, 척도까지 맞추려면 line='s' 를 쓴다.
sm.qqplot(residuals, line='45', ax=axes[1])
plt.tight_layout()
plt.show()

import numpy as np
from scipy.stats import skew, kurtosis

cap = y.max()
capped = (y == cap)
print(f"반응변수의 최댓값 = {cap}")
print(f"그 값에 붙어 있는 관측 = {capped.sum()}건 ({capped.mean():.2%})")

# 그 점들의 잔차는 적합값의 1차함수여야 한다
band = residuals[capped]
print(f"\n절단된 점들에서  잔차 + 적합값 의 범위 = "
      f"[{(band + y_pred[capped]).min():.5f}, {(band + y_pred[capped]).max():.5f}]")
print(f"그 점들의 corr(적합값, 잔차) = {np.corrcoef(y_pred[capped], band)[0, 1]:+.10f}")
slope = np.polyfit(y_pred[capped], band, 1)
print(f"그 점들에 맞춘 직선: 기울기 {slope[0]:+.6f},  절편 {slope[1]:.5f}")

print(f"\n잔차의 왜도 = {skew(residuals):.4f},  첨도 = {kurtosis(residuals, fisher=False):.4f}")
print(f"잔차의 표준편차 = {residuals.std():.4f}  (qqplot 의 line='45' 가 가정하는 값은 1)")

반응변수의 최댓값 = 5.00001
그 값에 붙어 있는 관측 = 965건 (4.68%)

절단된 점들에서  잔차 + 적합값 의 범위 = [5.00001, 5.00001]
그 점들의 corr(적합값, 잔차) = -1.0000000000
그 점들에 맞춘 직선: 기울기 -1.000000,  절편 5.00001

잔차의 왜도 = 1.2564,  첨도 = 5.9654
잔차의 표준편차 = 0.8315  (qqplot 의 line='45' 가 가정하는 값은 1)

# 변수를 늘려 가며 네 모형을 견준다. R^2 는 변수를 더하면 반드시 오르므로
# 모형 고르기에 쓸 수 없다. 벌점이 붙는 AIC·BIC 를 본다.
# BIC 의 벌점이 더 무거우므로 대개 더 작은 모형을 고른다.
feature_sets = {
    'Model 1': ['MedInc'],
    'Model 2': ['MedInc', 'AveRooms'],
    'Model 3': ['MedInc', 'AveRooms', 'AveOccup'],
    'Model 4': list(housing.feature_names),
}

for name, feats in feature_sets.items():
    X_temp = add_constant(df[feats])
    m = OLS(y, X_temp).fit()
    print(f"{name}: AIC={m.aic:.1f}, BIC={m.bic:.1f}, "
          f"R2={m.rsquared:.4f}")

Model 1: AIC=51249.3, BIC=51265.2, R2=0.4734
Model 2: AIC=51016.2, BIC=51040.1, R2=0.4794
Model 3: AIC=50962.2, BIC=50994.0, R2=0.4808
Model 4: AIC=45265.5, BIC=45337.0, R2=0.6062

import numpy as np
from scipy.stats import f as fdist

n = len(y)
fits = {name: OLS(y, add_constant(df[feats])).fit()
        for name, feats in feature_sets.items()}

# 간이 공식과 statsmodels 의 차이가 상수인가
print(f"{'모형':>9}{'statsmodels AIC':>18}{'간이 공식':>14}{'차':>14}")
for name, m in fits.items():
    simple = n * np.log(m.ssr / n) + 2 * len(m.params)
    print(f"{name:>9}{m.aic:>18.1f}{simple:>14.1f}{m.aic - simple:>14.4f}")
print(f"n ln(2 pi) + n = {n * np.log(2 * np.pi) + n:.4f}")

# Model 3 -> Model 4 를 F 로 읽는다
m3, m4 = fits['Model 3'], fits['Model 4']
dk = len(m4.params) - len(m3.params)
F = ((m3.ssr - m4.ssr) / dk) / (m4.ssr / m4.df_resid)
print(f"\n추가한 변수 {dk}개,  F = {F:.2f},  p = {fdist.sf(F, dk, m4.df_resid):.3g}")
print(f"dAIC: 공식 {n * np.log(m4.ssr / m3.ssr) + 2 * dk:.4f},"
      f"  실제 {m4.aic - m3.aic:.4f}")
print(f"dBIC: 공식 {n * np.log(m4.ssr / m3.ssr) + dk * np.log(n):.4f},"
      f"  실제 {m4.bic - m3.bic:.4f}")
print(f"\n문턱:  AIC 는 F > 2,  BIC 는 F > ln n = {np.log(n):.4f}")

       모형   statsmodels AIC         간이 공식             차
  Model 1           51249.3       -7324.4    58573.7827
  Model 2           51016.2       -7557.5    58573.7827
  Model 3           50962.2       -7611.5    58573.7827
  Model 4           45265.5      -13308.2    58573.7827
n ln(2 pi) + n = 58573.7827

추가한 변수 5개,  F = 1314.16,  p = 0
dAIC: 공식 -5696.7080,  실제 -5696.7080
dBIC: 공식 -5657.0331,  실제 -5657.0331

문턱:  AIC 는 F > 2,  BIC 는 F > ln n = 9.9350

X_full = add_constant(df[list(housing.feature_names)])
for i, col in enumerate(X_full.columns):
    if col == 'const':
        continue          # VIF of the intercept is meaningless
    vif = variance_inflation_factor(X_full.values, i)
    print(f"{col}: VIF = {vif:.2f}")

MedInc: VIF = 2.50
HouseAge: VIF = 1.24
AveRooms: VIF = 8.34
AveBedrms: VIF = 6.99
Population: VIF = 1.14
AveOccup: VIF = 1.01
Latitude: VIF = 9.30
Longitude: VIF = 8.96

y_log = np.log(df['PRICE'])
model_log = OLS(y_log, X).fit()
sm.qqplot(model_log.resid, line='45')

influence = model.get_influence()
cooks_d = influence.cooks_distance[0]
mask = cooks_d < 4 / len(y)
X_clean = X[mask]
y_clean = y[mask]
model_clean = OLS(y_clean, X_clean).fit()
print(f"Original R2: {model.rsquared:.4f}")
print(f"Cleaned R2:  {model_clean.rsquared:.4f}")

Original R2: 0.4808
Cleaned R2:  0.5784

df['MedInc_x_AveRooms'] = df['MedInc'] * df['AveRooms']
df['MedInc_sq'] = df['MedInc'] ** 2
X_ext = add_constant(df[features + ['MedInc_x_AveRooms', 'MedInc_sq']])
model_ext = OLS(y, X_ext).fit()
print(f"Base AIC: {model.aic:.1f}")
print(f"Extended AIC: {model_ext.aic:.1f}")

Base AIC: 50962.2
Extended AIC: 50767.5