선형회귀의 등분산성 확인¶

import numpy as np
import pandas as pd
import statsmodels.api as sm

rng = np.random.default_rng(7)
n = 120

# X는 균등, Y는 X에 선형으로 의존하되 오차의 분산이 X와 함께 커진다.
# 이렇게 두면 선형성은 성립하고 등분산성만 깨져, 각 진단이 무엇을
# 잡아내고 무엇을 놓치는지 구분해 볼 수 있다.
X = rng.uniform(0, 10, n)
Y = 2.0 + 1.5 * X + rng.normal(0, 0.5 + 0.35 * X, n)

df = pd.DataFrame({"X": X, "Y": Y})
model = sm.OLS(Y, sm.add_constant(X)).fit()
residuals = model.resid
fitted = model.fittedvalues

print(f"beta_hat = {model.params.round(4)}")
print(f"R^2 = {model.rsquared:.4f}")

beta_hat = [1.5933 1.5317]
R^2 = 0.7813

import numpy as np

# 모자행렬의 대각 h_ii. 멱등성 때문에 합이 계수의 개수 p = 2 다.
Xd = sm.add_constant(X)
H = Xd @ np.linalg.inv(Xd.T @ Xd) @ Xd.T
h = np.diag(H)
print(f"sum h_ii = {h.sum():.6f}   (p = 2 여야 한다)")
print(f"h_ii 범위 = [{h.min():.4f}, {h.max():.4f}],  평균 {h.mean():.4f} = p/n")

# (1) 의 식. 이 자료는 sigma_i 를 우리가 안다.
sigma = 0.5 + 0.35 * X
Es2 = ((1 - h) * sigma ** 2).sum() / (n - 2)
print(f"\nE[s^2] = sum (1-h_ii) sigma_i^2 / (n-2) = {Es2:.4f}")
print(f"참고: sigma_i^2 의 단순평균            = {(sigma ** 2).mean():.4f}")
print(f"이번 표본에서 실제로 나온 s^2         = {model.mse_resid:.4f}")

# 같은 X 를 고정한 채 오차만 20000번 다시 뽑아 E[s^2] 을 직접 재 본다.
sim = np.random.default_rng(99)
E = sim.normal(0, 1, (20000, n)) * sigma
R = E - E @ H.T
print(f"모의실험이 준 E[s^2]                  = {(R ** 2).sum(axis=1).mean() / (n - 2):.4f}")

# OLS 와 WLS 의 기울기 분산
Sxx = ((X - X.mean()) ** 2).sum()
c = (X - X.mean()) / Sxx
w = 1 / sigma ** 2
xbar_w = (w * X).sum() / w.sum()
var_ols = (c ** 2 * sigma ** 2).sum()
var_wls = 1 / (w * (X - xbar_w) ** 2).sum()
print(f"\nSE(b1) OLS = {np.sqrt(var_ols):.6f}")
print(f"SE(b1) WLS = {np.sqrt(var_wls):.6f}")
print(f"OLS 의 분산이 WLS 의 {var_ols / var_wls:.3f} 배  (효율 {var_wls / var_ols:.1%})")

sum h_ii = 2.000000   (p = 2 여야 한다)
h_ii 범위 = [0.0083, 0.0329],  평균 0.0167 = p/n

E[s^2] = sum (1-h_ii) sigma_i^2 / (n-2) = 6.0853
참고: sigma_i^2 의 단순평균            = 6.0916
이번 표본에서 실제로 나온 s^2         = 5.5889
모의실험이 준 E[s^2]                  = 6.0743

SE(b1) OLS = 0.082525
SE(b1) WLS = 0.056809
OLS 의 분산이 WLS 의 2.110 배  (효율 47.4%)

import matplotlib.pyplot as plt
import numpy as np

# 앞에서 적합한 model 을 그대로 쓴다.
residuals = model.resid
fitted = model.fittedvalues

# 등분산성 진단의 기본 그림이다. 점들의 세로 폭이 왼쪽부터 오른쪽까지
# 일정해야 한다. 깔때기처럼 벌어지면 분산이 적합값에 따라 커진다는 뜻이다.
plt.scatter(fitted, residuals, alpha=0.5)
plt.axhline(y=0, color='red', linestyle='--')
plt.xlabel('Fitted Values')
plt.ylabel('Residuals')
plt.title('Residuals vs. Fitted Values')
plt.show()

# 눈으로 본 깔때기를 수치로 바꾼다. 적합값의 사분위로 잘라 잔차의 폭을 잰다.
q = np.quantile(fitted, [0, 0.25, 0.5, 0.75, 1.0])
print(f"{'구간':>6}{'n':>5}{'sd(e)':>9}{'평균 sigma_i':>14}")
for k in range(4):
    sel = (fitted >= q[k]) & (fitted <= q[k + 1])
    print(f"{'Q' + str(k + 1):>6}{sel.sum():>5}{residuals[sel].std(ddof=1):>9.4f}"
          f"{(0.5 + 0.35 * X[sel]).mean():>14.4f}")

# 세로축의 0 선은 자료가 아니라 적합이 만든 것이다.
print(f"\n잔차의 합          = {residuals.sum():+.2e}")
print(f"잔차와 적합값의 상관 = {np.corrcoef(residuals, fitted)[0, 1]:+.2e}")

    구간    n    sd(e)    평균 sigma_i
    Q1   30   0.7828        0.9259
    Q2   30   1.6717        1.8428
    Q3   30   2.6705        2.6569
    Q4   30   3.4243        3.5777

잔차의 합          = +4.23e-13
잔차와 적합값의 상관 = +7.97e-16

from statsmodels.stats.diagnostic import het_breuschpagan

# Breusch-Pagan 검정: 잔차의 제곱을 설명변수에 회귀해, 설명되는 몫이
# 있는지를 본다. p-값이 작으면 등분산을 기각한다.
# 분산이 설명변수의 선형함수로 커지는 경우를 잘 잡아낸다.
bp_test = het_breuschpagan(model.resid, model.model.exog)
labels = ['LM Statistic', 'LM p-value', 'F-Statistic', 'F p-value']
for label, value in zip(labels, bp_test):
    print(f'{label}: {value:.4f}')

LM Statistic: 25.4256
LM p-value: 0.0000
F-Statistic: 31.7235
F p-value: 0.0000

import numpy as np
import statsmodels.api as sm
from scipy.stats import chi2

e = model.resid
Xd = model.model.exog

# 보조회귀: e^2 을 상수와 X 에 회귀시킨다.
aux = sm.OLS(e ** 2, Xd).fit()
print(f"보조회귀 R^2      = {aux.rsquared:.10f}")
print(f"단순회귀이므로 corr(e^2, X)^2 = {np.corrcoef(e ** 2, Xd[:, 1])[0, 1] ** 2:.10f}")
print(f"LM = n R^2        = {len(e) * aux.rsquared:.10f}")
print(f"statsmodels 의 LM = {het_breuschpagan(e, Xd)[0]:.10f}")
print(f"\np = chi2(1).sf(LM) = {chi2.sf(len(e) * aux.rsquared, 1):.4e}")
print(f"보조회귀의 F       = {aux.fvalue:.4f},  p = {aux.f_pvalue:.4e}")

보조회귀 R^2      = 0.2118803813
단순회귀이므로 corr(e^2, X)^2 = 0.2118803813
LM = n R^2        = 25.4256457519
statsmodels 의 LM = 25.4256457519

p = chi2(1).sf(LM) = 4.5977e-07
보조회귀의 F       = 31.7235,  p = 1.2258e-07

from statsmodels.stats.diagnostic import het_white

# White 검정: 설명변수의 제곱과 교차항까지 넣어 회귀한다. 그래서 선형이
# 아닌 형태의 이분산도 잡지만, 항이 많아지는 만큼 검정력이 흩어진다.
white_test = het_white(model.resid, model.model.exog)
labels = ['LM Statistic', 'LM p-value', 'F-Statistic', 'F p-value']
for label, value in zip(labels, white_test):
    print(f'{label}: {value:.4f}')

LM Statistic: 26.2253
LM p-value: 0.0000
F-Statistic: 16.3603
F p-value: 0.0000

import numpy as np
import statsmodels.api as sm
from scipy.stats import chi2
from statsmodels.stats.diagnostic import het_breuschpagan

e = model.resid
Xd = model.model.exog

# White 의 보조회귀는 상수, X, X^2 이다. 설명변수가 하나뿐이라 교차항은 없다.
Z = np.column_stack([np.ones(len(e)), Xd[:, 1], Xd[:, 1] ** 2])
auxw = sm.OLS(e ** 2, Z).fit()
print(f"White 보조회귀 R^2 = {auxw.rsquared:.10f}")
print(f"LM = n R^2         = {len(e) * auxw.rsquared:.10f}")
print(f"statsmodels 의 LM  = {het_white(e, Xd)[0]:.10f}")

bp = het_breuschpagan(e, Xd)[0]
wh = het_white(e, Xd)[0]
print(f"\n{'검정':>14}{'LM':>10}{'df':>5}{'p':>13}")
print(f"{'Breusch-Pagan':>14}{bp:>10.4f}{1:>5}{chi2.sf(bp, 1):>13.3e}")
print(f"{'White':>14}{wh:>10.4f}{2:>5}{chi2.sf(wh, 2):>13.3e}")
print(f"\nLM 은 {wh - bp:.4f} 커졌지만 p-값은 {chi2.sf(wh, 2) / chi2.sf(bp, 1):.1f} 배 커졌다")

White 보조회귀 R^2 = 0.2185443843
LM = n R^2         = 26.2253261126
statsmodels 의 LM  = 26.2253261126

            검정        LM   df            p
 Breusch-Pagan   25.4256    1    4.598e-07
         White   26.2253    2    2.019e-06

LM 은 0.7997 커졌지만 p-값은 4.4 배 커졌다

import numpy as np
import matplotlib.pyplot as plt

residuals = model.resid
fitted = model.fittedvalues

# 척도-위치 그림은 부호를 없애고 크기만 본다. 제곱근을 씌워 큰 값이
# 그림을 독차지하지 않게 한다.
standardized_residuals = (residuals - np.mean(residuals)) / np.std(residuals)

plt.scatter(fitted, np.sqrt(np.abs(standardized_residuals)), alpha=0.5)
plt.xlabel('Fitted Values')
plt.ylabel('√|Standardized Residuals|')
plt.title('Scale-Location Plot')

# 눈으로만 판단하기 어려우므로 평활선을 얹는다. 이 선이 평평하면
# 등분산이고, 올라가거나 내려가면 분산이 적합값에 따라 달라진다는 뜻이다.
from statsmodels.nonparametric.smoothers_lowess import lowess
smooth = lowess(np.sqrt(np.abs(standardized_residuals)), fitted, frac=0.6)
plt.plot(smooth[:, 0], smooth[:, 1], color='red', linewidth=2)
plt.show()

# (1) 제곱근이 왜 필요한가. 표준정규 Z 로 세 변환의 왜도를 잰다.
z = np.random.default_rng(0).normal(size=2_000_000)
from scipy.stats import skew
print(f"왜도:  |Z| = {skew(np.abs(z)):.4f},  "
      f"sqrt|Z| = {skew(np.sqrt(np.abs(z))):.4f},  Z^2 = {skew(z ** 2):.4f}")

# (2) 평활선이 얼마나 올라갔는가. 세로축이 sqrt 척도이므로 제곱하면 sd 비가 된다.
y0, y1 = smooth[0, 1], smooth[-1, 1]
print(f"\n평활선 왼쪽 끝 {y0:.4f}  →  오른쪽 끝 {y1:.4f}")
print(f"비 {y1 / y0:.4f},  제곱하면 {(y1 / y0) ** 2:.4f}  (sd 의 비로 읽는다)")

b0, b1 = model.params
x0, x1 = (smooth[0, 0] - b0) / b1, (smooth[-1, 0] - b0) / b1
print(f"같은 자리의 참 sigma = {0.5 + 0.35 * x0:.4f}, {0.5 + 0.35 * x1:.4f}  "
      f"→ 비 {(0.5 + 0.35 * x1) / (0.5 + 0.35 * x0):.4f}")

# 간단 표준화와 내부 스튜던트화 잔차의 차이
ri = model.get_influence().resid_studentized_internal
print(f"\n두 표준화의 최대 차이 = {np.abs(standardized_residuals - ri).max():.4f}")

왜도:  |Z| = 0.9971,  sqrt|Z| = 0.0852,  Z^2 = 2.8321

평활선 왼쪽 끝 0.4066  →  오른쪽 끝 1.1431
비 2.8111,  제곱하면 7.9021  (sd 의 비로 읽는다)
같은 자리의 참 sigma = 0.5131, 3.9843  → 비 7.7655

두 표준화의 최대 차이 = 0.0206