다중공선성과 영향점¶

import numpy as np
import pandas as pd
import statsmodels.api as sm
import matplotlib.pyplot as plt
from sklearn.datasets import make_regression

# 깨끗한 선형 자료를 먼저 만든다.
np.random.seed(0)
X, y = make_regression(n_samples=100, n_features=1, noise=10)
data = pd.DataFrame({'X': X.flatten(), 'y': y})

# 여기에 이상치 넷을 일부러 심는다. 위아래로 짝을 맞춰 넣었으므로
# 기울기보다는 잔차의 퍼짐이 크게 흔들린다.
data.loc[95, 'y'] += 80
data.loc[96, 'y'] -= 80
data.loc[97, 'y'] += 60
data.loc[98, 'y'] -= 60

# 이상치가 든 채로 적합한다.
X_with_const = sm.add_constant(data['X'])
model = sm.OLS(data['y'], X_with_const).fit()

# Cook 의 거리는 그 관측값 하나를 뺐을 때 적합값 전체가 얼마나 움직이는지를
# 잰다. 잔차가 크다고 다 영향점인 것은 아니다 — 지렛값도 함께 커야 한다.
influence = model.get_influence()
cooks_d, _ = influence.cooks_distance

# 4/n 은 널리 쓰이는 어림 기준일 뿐 검정이 아니다.
n = len(data)
threshold = 4 / n
outliers = np.where(cooks_d > threshold)[0]

# 이상치가 있을 때와 없을 때를 나란히 그린다
fig, axes = plt.subplots(2, 2, figsize=(14, 10))

# 윗줄: 이상치가 든 자료. 회귀선과 잔차 그림을 나란히 본다.
axes[0, 0].scatter(data['X'], data['y'], alpha=0.7, label='Data Points')
axes[0, 0].plot(data['X'], model.fittedvalues, color='orange', label='Regression Line')
axes[0, 0].set_title('Regression Plot (With Outliers)')
axes[0, 0].set_xlabel('Predictor (X)')
axes[0, 0].set_ylabel('Response (y)')
axes[0, 0].legend()

axes[0, 1].scatter(model.fittedvalues, model.resid, alpha=0.7)
axes[0, 1].axhline(0, color='red', linestyle='--')
axes[0, 1].set_title('Residual Plot (With Outliers)')
axes[0, 1].set_xlabel('Fitted Values')
axes[0, 1].set_ylabel('Residuals')

# 아랫줄: 문턱을 넘은 관측값을 뺀 뒤 다시 적합한 결과다.
# 두 줄을 견주는 것이 이 그림의 목적이지, 이상치를 지우라는 뜻이 아니다.
# 지울지 말지는 그 값이 왜 생겼는지를 알아본 뒤에 정할 일이다.
data_no_outliers = data.drop(index=outliers)
X_with_const_no_outliers = sm.add_constant(data_no_outliers['X'])
model_no_outliers = sm.OLS(data_no_outliers['y'], X_with_const_no_outliers).fit()

axes[1, 0].scatter(data_no_outliers['X'], data_no_outliers['y'], alpha=0.7, label='Data Points')
axes[1, 0].plot(data_no_outliers['X'], model_no_outliers.fittedvalues, color='orange', label='Regression Line')
axes[1, 0].set_title('Regression Plot (Without Outliers)')
axes[1, 0].set_xlabel('Predictor (X)')
axes[1, 0].set_ylabel('Response (y)')
axes[1, 0].legend()

axes[1, 1].scatter(model_no_outliers.fittedvalues, model_no_outliers.resid, alpha=0.7)
axes[1, 1].axhline(0, color='red', linestyle='--')
axes[1, 1].set_title('Residual Plot (Without Outliers)')
axes[1, 1].set_xlabel('Fitted Values')
axes[1, 1].set_ylabel('Residuals')

plt.tight_layout()
plt.show()

import numpy as np
import pandas as pd
import statsmodels.api as sm
from sklearn.datasets import make_regression

# 교란하기 전의 적합
np.random.seed(0)
X0, y0 = make_regression(n_samples=100, n_features=1, noise=10)
clean = sm.OLS(y0, sm.add_constant(X0)).fit()

x = data['X'].values
xbar, Sxx = x.mean(), ((x - x.mean()) ** 2).sum()
c = (x - xbar) / Sxx
delta = np.zeros(len(x))
delta[[95, 96, 97, 98]] = [80, -80, 60, -60]

print(f"S_xx = {Sxx:.4f}")
print(f"공식  d_beta1 = sum c_i delta_i = {(c * delta).sum():.12f}")
print(f"실제        = {model.params.iloc[1] - clean.params[1]:.12f}")
print(f"공식  d_beta0 = {delta.mean() - xbar * (c * delta).sum():.12f}")
print(f"실제        = {model.params.iloc[0] - clean.params[0]:.12f}")

# 짝을 이루면 x 의 차만 남는다
print(f"\n두 쌍의 기여:  80(x95 - x96)/S_xx = "
      f"{80 * (x[95] - x[96]) / Sxx:+.4f},  60(x97 - x98)/S_xx = "
      f"{60 * (x[97] - x[98]) / Sxx:+.4f}")

print(f"\n{'':>12}{'절편':>10}{'기울기':>10}{'sigma_hat':>12}{'R^2':>9}")
for name, m in [('교란 전', clean), ('교란 후', model), ('제거 후', model_no_outliers)]:
    b0, b1 = np.asarray(m.params)
    print(f"{name:>12}{b0:>10.4f}{b1:>10.4f}"
          f"{np.sqrt(m.mse_resid):>12.4f}{m.rsquared:>9.4f}")

print(f"\nCook 문턱 {threshold} 을 넘은 관측 = {list(outliers)}")
print(f"다섯째로 큰 Cook 거리 = {np.sort(cooks_d)[-5]:.4f}")

S_xx = 101.5827
공식  d_beta1 = sum c_i delta_i = 0.531814777033
실제        = 0.531814777033
공식  d_beta0 = -0.031806786446
실제        = -0.031806786446

두 쌍의 기여:  80(x95 - x96)/S_xx = -0.6319,  60(x97 - x98)/S_xx = +1.1637

                    절편       기울기   sigma_hat      R^2
        교란 전   -0.8142   42.6194     10.7936   0.9417
        교란 후   -0.8460   43.1512     16.5534   0.8757
        제거 후   -1.1852   43.0230     10.7420   0.9430

Cook 문턱 0.04 을 넘은 관측 = [95, 96, 97, 98]
다섯째로 큰 Cook 거리 = 0.0355