선형회귀의 선형성 확인¶

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
import statsmodels.api as sm

# (1) 의 두 조각
var_signal = 1.5 ** 2 * (10 ** 2 / 12)
e_sigma2 = 0.25 + 0.35 * 5 + 0.35 ** 2 * (100 / 3)
print(f"Var(E[Y|X]) = {var_signal:.4f}")
print(f"E[sigma(X)^2] = {e_sigma2:.4f}   (= 73/12)")
print(f"R^2_pop = {var_signal / (var_signal + e_sigma2):.4f}")
print(f"표본 R^2 (n = 120) = {model.rsquared:.4f}")

# 같은 모형에서 40만 개를 뽑으면 어디로 가는가
big = np.random.default_rng(5)
Xb = big.uniform(0, 10, 400_000)
Yb = 2 + 1.5 * Xb + big.normal(0, 0.5 + 0.35 * Xb)
print(f"표본 R^2 (n = 400000) = {sm.OLS(Yb, sm.add_constant(Xb)).fit().rsquared:.4f}")

# 기울기가 참값에서 얼마나 떨어져 있는가 (참 SE 로 재면)
Sxx = ((X - X.mean()) ** 2).sum()
c = (X - X.mean()) / Sxx
se_true = np.sqrt((c ** 2 * (0.5 + 0.35 * X) ** 2).sum())
print(f"\nb1 = {model.params[1]:.4f},  참 SE = {se_true:.4f},  "
      f"(b1 - 1.5)/SE = {(model.params[1] - 1.5) / se_true:+.3f}")

Var(E[Y|X]) = 18.7500
E[sigma(X)^2] = 6.0833   (= 73/12)
R^2_pop = 0.7550
표본 R^2 (n = 120) = 0.7813
표본 R^2 (n = 400000) = 0.7542

b1 = 1.5317,  참 SE = 0.0825,  (b1 - 1.5)/SE = +0.384

import matplotlib.pyplot as plt
import numpy as np

# 설명변수가 하나뿐이면 산점도가 선형성을 보는 가장 곧은 방법이다.
# 변수가 여럿이면 다른 변수의 효과가 섞이므로 아래 잔차 그림으로 가야 한다.
plt.scatter(X, Y)
plt.xlabel('Independent Variable')
plt.ylabel('Dependent Variable')
plt.title('Scatterplot of Y vs X')
plt.show()

print(f"corr(X, Y) = {np.corrcoef(X, Y)[0, 1]:.4f}")
print(f"그 제곱    = {np.corrcoef(X, Y)[0, 1] ** 2:.4f}   (= 단순회귀의 R^2)")

corr(X, Y) = 0.8839
그 제곱    = 0.7813   (= 단순회귀의 R^2)

import numpy as np
import statsmodels.api as sm

print(f"{'a':>6}{'corr(X,Y)':>11}{'선형 R^2':>11}{'t(x^2)':>9}{'p':>11}{'dAIC':>9}")
for a in (0.00, 0.05, 0.10, 0.20):
    Yc = Y + a * (X - 5) ** 2
    m1 = sm.OLS(Yc, sm.add_constant(X)).fit()
    m2 = sm.OLS(Yc, sm.add_constant(np.column_stack((X, X ** 2)))).fit()
    print(f"{a:>6.2f}{np.corrcoef(X, Yc)[0, 1]:>11.4f}{m1.rsquared:>11.4f}"
          f"{m2.tvalues[2]:>9.3f}{m2.pvalues[2]:>11.2e}{m2.aic - m1.aic:>+9.2f}")

     a  corr(X,Y)     선형 R^2   t(x^2)          p     dAIC
  0.00     0.8839     0.7813    1.380   1.70e-01    +0.06
  0.05     0.8770     0.7691    3.161   2.00e-03    -7.83
  0.10     0.8653     0.7487    4.941   2.61e-06   -20.74
  0.20     0.8295     0.6881    8.502   7.04e-14   -55.73

import numpy as np
import statsmodels.api as sm

e = model.resid
predictions = model.fittedvalues

# 그림이 아니라 수로 확인한다. 직교 조건 두 개는 적합이 강제한다.
print(f"sum e_i      = {e.sum():+.2e}")
print(f"sum x_i e_i  = {(X * e).sum():+.2e}")
print(f"sum x_i^2 e_i = {(X ** 2 * e).sum():+.4f}   (이것만 자유롭다)")

# x^2 에서 1 과 x 로 설명되는 몫을 걷어 낸 나머지
x2_resid = sm.OLS(X ** 2, sm.add_constant(X)).fit().resid
fwl = sm.OLS(e, x2_resid).fit()
full = sm.OLS(Y, sm.add_constant(np.column_stack((X, X ** 2)))).fit()
print(f"\nFWL 로 구한 beta2   = {fwl.params[0]:.15f}")
print(f"완전모형의 beta2    = {full.params[2]:.15f}")
print(f"\nFWL 의 t  = {fwl.tvalues[0]:.4f}  (자유도 {int(fwl.df_resid)})")
print(f"완전모형의 t = {full.tvalues[2]:.4f}  (자유도 {int(full.df_resid)})")
print(f"두 t 의 비 = {fwl.tvalues[0] / full.tvalues[2]:.6f},  "
      f"sqrt(119/117) = {np.sqrt(119 / 117):.6f}")

sum e_i      = +4.23e-13
sum x_i e_i  = +2.79e-12
sum x_i^2 e_i = +272.5912   (이것만 자유롭다)

FWL 로 구한 beta2   = 0.038757372543221
완전모형의 beta2    = 0.038757372543221

FWL 의 t  = 1.3919  (자유도 119)
완전모형의 t = 1.3802  (자유도 117)
두 t 의 비 = 1.008511,  sqrt(119/117) = 1.008511

import statsmodels.api as sm
from statsmodels.graphics.regressionplots import plot_ccpr
import matplotlib.pyplot as plt

# 성분-잔차 그림(CCPR)은 다른 변수의 효과를 걷어 낸 뒤, 관심 있는 변수
# 하나와 반응의 관계만 남겨 보여 준다. 설명변수가 여럿일 때 선형성을
# 변수별로 확인하는 방법이다. exog_idx 로 볼 변수를 고른다.
results = sm.OLS(Y, sm.add_constant(X)).fit()

fig, ax = plt.subplots(figsize=(8, 6))
plot_ccpr(results, exog_idx=1, ax=ax)
plt.show()

import numpy as np
import statsmodels.api as sm

b0, b1 = model.params
partial = b1 * X + model.resid          # CCPR 의 세로축

print(f"max |부분잔차 - (Y - b0)| = {np.abs(partial - (Y - b0)).max():.2e}")
print(f"부분잔차에 그은 최소제곱 기울기 = "
      f"{sm.OLS(partial, sm.add_constant(X)).fit().params[1]:.12f}")
print(f"원래 모형의 b1                 = {b1:.12f}")

# 변수가 둘일 때 비로소 쓸모가 생긴다. X2 가 X1 과 섞여 있는 자료.
r = np.random.default_rng(11)
x1 = r.uniform(0, 10, 300)
x2 = x1 + r.normal(0, 1.0, 300)
y = 1 + 2.0 * x1 - 1.5 * x2 + r.normal(0, 2.0, 300)
mm = sm.OLS(y, sm.add_constant(np.column_stack((x1, x2)))).fit()

marg = sm.OLS(y, sm.add_constant(x2)).fit().params[1]
part = mm.params[2] * x2 + mm.resid
cc = sm.OLS(part, sm.add_constant(x2)).fit().params[1]
print(f"\n참 계수            = -1.5")
print(f"y 를 x2 에만 회귀   = {marg:+.4f}   (산점도가 보여 주는 기울기)")
print(f"완전모형의 beta_2   = {mm.params[2]:+.4f}")
print(f"부분잔차 그림의 기울기 = {cc:+.4f}")

max |부분잔차 - (Y - b0)| = 3.55e-15
부분잔차에 그은 최소제곱 기울기 = 1.531654179791
원래 모형의 b1                 = 1.531654179791

참 계수            = -1.5
y 를 x2 에만 회귀   = +0.2369   (산점도가 보여 주는 기울기)
완전모형의 beta_2   = -1.5705
부분잔차 그림의 기울기 = -1.5705

import numpy as np
import statsmodels.api as sm

# 곡선이 의심되면 이차항을 넣어 견준다. R^2 는 항을 더하면 반드시 오르므로
# 판단 근거가 되지 못한다. 벌점이 붙는 AIC 로 보아야 한다.
X_poly = np.column_stack((X, X**2))
model_quad = sm.OLS(Y, sm.add_constant(X_poly)).fit()

# summary()는 실행 날짜와 시각을 함께 찍으므로 계수 표만 인쇄한다.
print(model_quad.summary().tables[1])
print(f"R^2: 선형 {model.rsquared:.4f}  →  이차 {model_quad.rsquared:.4f}")
print(f"AIC: 선형 {model.aic:.2f}  →  이차 {model_quad.aic:.2f}")

==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
const          2.2252      0.628      3.545      0.001       0.982       3.468
x1             1.1466      0.289      3.971      0.000       0.575       1.718
x2             0.0388      0.028      1.380      0.170      -0.017       0.094
==============================================================================
R^2: 선형 0.7813  →  이차 0.7848
AIC: 선형 549.02  →  이차 549.08

import numpy as np

t = model_quad.tvalues[2]
nu = model_quad.df_resid          # n - p - 1 = 117
dAIC_t = n * np.log(1 - t ** 2 / (t ** 2 + nu)) + 2
dAIC_sse = n * np.log(model_quad.ssr / model.ssr) + 2

print(f"t(x^2) = {t:.4f},  t^2 = {t ** 2:.4f},  nu = {int(nu)}")
print(f"\ndAIC  (t 로)       = {dAIC_t:+.6f}")
print(f"dAIC  (SSE 비로)   = {dAIC_sse:+.6f}")
print(f"dAIC  (실제)       = {model_quad.aic - model.aic:+.6f}")
print(f"큰 n 근사 -t^2 + 2 = {-t ** 2 + 2:+.6f}")

print(f"\nAIC 문턱 sqrt(2)      = {np.sqrt(2):.4f}")
print(f"BIC 문턱 sqrt(log n)  = {np.sqrt(np.log(n)):.4f}")
print(f"dBIC (실제)           = {model_quad.bic - model.bic:+.6f}")

t(x^2) = 1.3802,  t^2 = 1.9048,  nu = 117

dAIC  (t 로)       = +0.062043
dAIC  (SSE 비로)   = +0.062043
dAIC  (실제)       = +0.062043
큰 n 근사 -t^2 + 2 = +0.095152

AIC 문턱 sqrt(2)      = 1.4142
BIC 문턱 sqrt(log n)  = 2.1880
dBIC (실제)           = +2.849535