잔차 분석¶

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import statsmodels.api as sm
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})

X_with_const = sm.add_constant(data['X'])
model = sm.OLS(data['y'], X_with_const).fit()

data['Fitted'] = model.fittedvalues
data['Residuals'] = model.resid

# 왼쪽은 회귀 그림, 오른쪽은 잔차 그림이다. 늘 짝으로 본다.
fig, (ax0, ax1) = plt.subplots(1, 2, figsize=(12, 3))

ax0.scatter(data['X'], data['y'], alpha=0.7, label='Data Points')
ax0.plot(data['X'], data['Fitted'], color='orange', label='Regression Line')
ax0.set_title('Regression Plot')
ax0.set_xlabel('Predictor (X)')
ax0.set_ylabel('Response (y)')
ax0.legend()

# 잔차가 0 선 둘레에 무늬 없이 흩어져 있으면 좋다. 이 자료가 그렇다.
ax1.scatter(data['Fitted'], data['Residuals'], alpha=0.7)
ax1.axhline(y=0, color='r', linestyle='--')
ax1.set_title('Residuals vs. Fitted Values Plot')
ax1.set_xlabel('Fitted Values')
ax1.set_ylabel('Residuals')

plt.tight_layout()
plt.show()

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

# make_regression 은 참 계수를 돌려줄 수 있다. 같은 씨앗이면 같은 자료다.
np.random.seed(0)
_, _, true_coef = make_regression(n_samples=100, n_features=1, noise=10, coef=True)
n, p, sigma = 100, 2, 10.0

print(f"참 기울기 = {true_coef:.4f},   적합값 = {model.params.iloc[1]:.4f}")
print(f"참 절편   = 0.0000,   적합값 = {model.params.iloc[0]:.4f}")
print(f"모집단 R^2 = {true_coef ** 2 / (true_coef ** 2 + sigma ** 2):.4f},"
      f"   표본 R^2 = {model.rsquared:.4f}")

e = model.resid.values
print(f"\n잔차제곱평균의 이론값 sigma^2 (n-p)/n = {sigma ** 2 * (n - p) / n:.3f}")
print(f"실측 = {(e ** 2).mean():.3f}")
print(f"sigma_hat = sqrt(SSE/(n-p)) = {np.sqrt(model.mse_resid):.4f}")
print(f"|e| > 2 sigma_hat 인 점의 개수 = {(np.abs(e) > 2 * np.sqrt(model.mse_resid)).sum()} / {n}")

# "무늬 없음" 을 눈 대신 수로 판정한다
f = model.fittedvalues.values
print(f"\ncorr(e, fitted)  = {np.corrcoef(e, f)[0, 1]:+.2e}   (항등적으로 0)")
print(f"corr(|e|, fitted) = {np.corrcoef(np.abs(e), f)[0, 1]:+.4f}   (등분산)")
aux = sm.OLS(e, sm.add_constant(np.column_stack([f, f ** 2]))).fit()
print(f"e 를 fitted, fitted^2 에 회귀한 2차항의 p = {aux.pvalues[2]:.4f}   (선형성)")

참 기울기 = 42.3855,   적합값 = 42.6194
참 절편   = 0.0000,   적합값 = -0.8142
모집단 R^2 = 0.9473,   표본 R^2 = 0.9417

잔차제곱평균의 이론값 sigma^2 (n-p)/n = 98.000
실측 = 114.171
sigma_hat = sqrt(SSE/(n-p)) = 10.7936
|e| > 2 sigma_hat 인 점의 개수 = 2 / 100

corr(e, fitted)  = -1.94e-15   (항등적으로 0)
corr(|e|, fitted) = +0.0223   (등분산)
e 를 fitted, fitted^2 에 회귀한 2차항의 p = 0.4335   (선형성)

import matplotlib.pyplot as plt
import numpy as np
from sklearn.linear_model import LinearRegression

def generate_data(n=50, noise_level=3.0, seed=0):
    """기울기 2 의 직선 자료를 만든다. noise_level 로 잡음 크기를 조절한다."""
    np.random.seed(seed)
    x = np.random.randn(n, 1)
    x.sort(axis=0)
    noise = np.random.normal(0, 1, size=x.shape)
    y = (1 + 2 * x + noise_level * noise).reshape((-1,))
    return x, y

def perform_regression(x, y):
    """최소제곱으로 적합하고 예측값까지 돌려준다."""
    model = LinearRegression()
    model.fit(x, y)
    y_pred = model.predict(x)
    return model, y_pred

def plot_regression_and_residuals(x, y, y_pred):
    """회귀 그림과 잔차 그림을 나란히 그린다. 아래에서 되풀이해 쓴다."""
    fig, (ax0, ax1) = plt.subplots(1, 2, figsize=(12, 3))

    ax0.plot(x, y, 'o', label="Data")
    ax0.plot(x, y_pred, '-b', label="Predicted")
    ax0.set_title('Regression Plot')
    ax0.legend()

    ax1.plot(x, y - y_pred, 'o', label="Residuals")
    ax1.set_title('Residual Plot')
    ax1.legend()

    for ax in (ax0, ax1):
        ax.spines['top'].set_visible(False)
        ax.spines['right'].set_visible(False)
        ax.spines['bottom'].set_position("zero")

    plt.tight_layout()
    plt.show()

x, y = generate_data()
model, y_pred = perform_regression(x, y)
plot_regression_and_residuals(x, y, y_pred)

import numpy as np
import statsmodels.api as sm

e = y - y_pred
xf = x.ravel()
Sxx = ((xf - xf.mean()) ** 2).sum()

print(f"참값 (절편, 기울기) = (1, 2)")
print(f"적합값              = ({model.intercept_:.4f}, {model.coef_[0]:.4f})")
print(f"참 SE(기울기) = sigma/sqrt(S_xx) = {3 / np.sqrt(Sxx):.4f}   (S_xx = {Sxx:.3f})")
print(f"(b1 - 2)/SE = {(model.coef_[0] - 2) / (3 / np.sqrt(Sxx)):+.3f}")

print(f"\n잔차제곱평균 이론 sigma^2 (n-p)/n = {9 * 48 / 50:.3f},  실측 {(e ** 2).mean():.3f}")

# 정말 깔때기인가. 적합값 사분위로 잘라 잔차의 폭을 잰다.
q = np.quantile(y_pred, [0, 0.25, 0.5, 0.75, 1.0])
print()
for k in range(4):
    sel = (y_pred >= q[k]) & (y_pred <= q[k + 1])
    print(f"  Q{k + 1}: n = {sel.sum()}, sd(e) = {e[sel].std(ddof=1):.4f}")
print(f"\ncorr(|e|, fitted) = {np.corrcoef(np.abs(e), y_pred)[0, 1]:+.4f}")
from statsmodels.stats.diagnostic import het_breuschpagan
lm, p = het_breuschpagan(e, sm.add_constant(xf))[:2]
print(f"Breusch-Pagan: LM = {lm:.4f}, p = {p:.4f}")

# 씨앗을 바꾸어 400번 되풀이하면 기울기는 2 둘레에 모인다.
slopes = []
for s in range(400):
    xs, ys = generate_data(seed=s)
    slopes.append(perform_regression(xs, ys)[0].coef_[0])
print(f"\n씨앗 400개의 기울기: 평균 {np.mean(slopes):.4f}, 표준편차 {np.std(slopes):.4f}")

참값 (절편, 기울기) = (1, 2)
적합값              = (0.8125, 2.8868)
참 SE(기울기) = sigma/sqrt(S_xx) = 0.3769   (S_xx = 63.340)
(b1 - 2)/SE = +2.352

잔차제곱평균 이론 sigma^2 (n-p)/n = 8.640,  실측 5.770

  Q1: n = 13, sd(e) = 1.6495
  Q2: n = 12, sd(e) = 2.7632
  Q3: n = 12, sd(e) = 3.2007
  Q4: n = 13, sd(e) = 2.1579

corr(|e|, fitted) = +0.1255
Breusch-Pagan: LM = 0.6560, p = 0.4180

씨앗 400개의 기울기: 평균 1.9988, 표준편차 0.4458

def generate_data(n=50, noise_level=3.0, d=1, seed=0):
    """차수 d 인 다항 자료를 만든다. d=1 이면 앞과 같은 직선이다."""
    np.random.seed(seed)
    x = np.random.randn(n, 1)
    x.sort(axis=0)
    noise = np.random.normal(0, 1, size=x.shape)
    y = (1 + np.sum([(k+1) * x**k for k in range(1, d+1)], axis=0) + noise_level * noise).reshape((-1,))
    return x, y

# 이차 자료를 직선으로 맞춘다. 회귀 그림만 보면 그럴듯해 보이지만
# 잔차 그림에는 굽은 무늬가 또렷하게 남는다. 잔차 그림을 보는 까닭이다.
x, y = generate_data(d=2)
model, y_pred = perform_regression(x, y)
plot_regression_and_residuals(x, y, y_pred)

import numpy as np
import statsmodels.api as sm

xf = x.ravel()
e = y - y_pred

print(f"직선 적합 (절편, 기울기) = ({model.intercept_:.4f}, {model.coef_[0]:.4f})")
print(f"모집단 극한값            = (4, 2)")
print(f"표본 적률: E[x] = {xf.mean():+.4f}, E[x^2] = {(xf ** 2).mean():.4f}, "
      f"E[x^3] = {(xf ** 3).mean():+.4f}")

# 잔차를 두 조각으로 정확히 가른다: 곡률의 몫과 잡음의 몫
curv = 3 * xf ** 2
curv_part = curv - sm.OLS(curv, sm.add_constant(xf)).fit().fittedvalues
noise_part = e - curv_part
print(f"\n두 조각의 합이 잔차와 같은가: 최대 오차 {np.abs(e - (curv_part + noise_part)).max():.1e}")
print(f"곡률 조각의 표준편차 = {curv_part.std():.4f}")
print(f"잡음 조각의 표준편차 = {noise_part.std():.4f}")
print(f"잔차의 분산 가운데 곡률이 차지하는 몫 = "
      f"{np.corrcoef(e, curv_part)[0, 1] ** 2:.4f}")
print(f"곡률 조각의 범위 = [{curv_part.min():.3f}, {curv_part.max():.3f}]")

직선 적합 (절편, 기울기) = (4.6161, 3.2857)
모집단 극한값            = (4, 2)
표본 적률: E[x] = +0.1406, E[x^2] = 1.2866, E[x^3] = +0.3493

두 조각의 합이 잔차와 같은가: 최대 오차 0.0e+00
곡률 조각의 표준편차 = 4.6361
잡음 조각의 표준편차 = 2.4021
잔차의 분산 가운데 곡률이 차지하는 몫 = 0.8136
곡률 조각의 범위 = [-3.816, 16.768]

import numpy as np
import matplotlib.pyplot as plt
import statsmodels.api as sm
from statsmodels.nonparametric.smoothers_lowess import lowess

# 같은 이야기를 평활선까지 얹어 더 또렷하게 본다.
np.random.seed(42)
x = np.random.uniform(-3, 3, 100)
y_true = 2 + 0.5 * x - 1.5 * x**2
y = y_true + np.random.normal(0, 1, len(x))

# 왼쪽에 쓸 선형 모형
X_linear = sm.add_constant(x)
model_linear = sm.OLS(y, X_linear).fit()
residuals_linear = model_linear.resid
y_pred_linear = model_linear.fittedvalues

# 오른쪽에 쓸 이차 모형
X_quad = sm.add_constant(np.column_stack([x, x**2]))
model_quad = sm.OLS(y, X_quad).fit()
residuals_quad = model_quad.resid
y_pred_quad = model_quad.fittedvalues

# 잔차에 평활선을 얹어 그린다
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 5))

# 선형 모형의 잔차
ax1.scatter(y_pred_linear, residuals_linear, alpha=0.6)
ax1.axhline(y=0, color='r', linestyle='--', linewidth=2)

# 평활선이 굽어 있으면 남은 구조가 있다는 신호다.
lowess_result = lowess(residuals_linear, y_pred_linear, frac=0.3)
ax1.plot(lowess_result[:, 0], lowess_result[:, 1], 'b-', linewidth=2.5,
         label='LOWESS Trend')

ax1.set_xlabel('Fitted Values')
ax1.set_ylabel('Residuals')
ax1.set_title('Linear Model: Clear Non-linearity Pattern')
ax1.legend()
ax1.grid(True, alpha=0.3)

# 이차항을 넣고 나면 평활선이 평평해진다. 무늬가 사라진 것이다.
ax2.scatter(y_pred_quad, residuals_quad, alpha=0.6)
ax2.axhline(y=0, color='r', linestyle='--', linewidth=2)

# 평활선을 얹는다
lowess_result_quad = lowess(residuals_quad, y_pred_quad, frac=0.3)
ax2.plot(lowess_result_quad[:, 0], lowess_result_quad[:, 1], 'b-', linewidth=2.5,
         label='LOWESS Trend')

ax2.set_xlabel('Fitted Values')
ax2.set_ylabel('Residuals')
ax2.set_title('Quadratic Model: Non-linearity Removed')
ax2.legend()
ax2.grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

# 요약 비교
print("Model Comparison:")
print(f"Linear Model R²:    {model_linear.rsquared:.4f}")
print(f"Quadratic Model R²: {model_quad.rsquared:.4f}")
print(f"Linear Model RSS:    {np.sum(residuals_linear**2):.2f}")
print(f"Quadratic Model RSS: {np.sum(residuals_quad**2):.2f}")

Model Comparison:
Linear Model R²:    0.0830
Quadratic Model R²: 0.9534
Linear Model RSS:    1530.56
Quadratic Model RSS: 77.72

import numpy as np
import statsmodels.api as sm

# 모집단 적률 (x ~ U(-3,3))
VarX = 36 / 12
EX4 = 3 ** 4 / 5
VarX2 = EX4 - VarX ** 2
VarY = 0.5 ** 2 * VarX + 1.5 ** 2 * VarX2 + 1.0
print(f"Var(x) = {VarX:.4f},  Var(x^2) = {VarX2:.4f},  Var(y) = {VarY:.4f}")
print(f"직선 적합의 극한 (절편, 기울기) = ({2 - 1.5 * VarX:.4f}, {0.5:.4f})")
print(f"실제 적합                      = ({model_linear.params[0]:.4f}, "
      f"{model_linear.params[1]:.4f})")
print(f"\n모집단 R^2:  직선 {0.5 ** 2 * VarX / VarY:.4f},"
      f"  이차 {(0.5 ** 2 * VarX + 1.5 ** 2 * VarX2) / VarY:.4f}")
print(f"표본   R^2:  직선 {model_linear.rsquared:.4f},  이차 {model_quad.rsquared:.4f}")

# 표본 R^2 이 0.083 인 것이 이상한가. 같은 모형에서 2000번 뽑아 본다.
r2s = []
for s in range(2000):
    rng = np.random.default_rng(s)
    xx = rng.uniform(-3, 3, 100)
    yy = 2 + 0.5 * xx - 1.5 * xx ** 2 + rng.normal(0, 1, 100)
    r2s.append(sm.OLS(yy, sm.add_constant(xx)).fit().rsquared)
r2s = np.array(r2s)
print(f"\n직선 R^2 의 모의분포: 평균 {r2s.mean():.4f}, 표준편차 {r2s.std():.4f}")
print(f"0.0830 은 그 분포의 {(r2s < model_linear.rsquared).mean():.1%} 분위")

Var(x) = 3.0000,  Var(x^2) = 7.2000,  Var(y) = 17.9500
직선 적합의 극한 (절편, 기울기) = (-2.5000, 0.5000)
실제 적합                      = (-2.7512, 0.6626)

모집단 R^2:  직선 0.0418,  이차 0.9443
표본   R^2:  직선 0.0830,  이차 0.9534

직선 R^2 의 모의분포: 평균 0.0553, 표준편차 0.0493
0.0830 은 그 분포의 75.9% 분위

def perform_regression(x, y, d=1):
    """차수 d 의 다항회귀. x, x^2, ... 를 열로 쌓아 넣기만 하면 된다.

    항이 x 의 거듭제곱일 뿐 계수에 대해서는 여전히 선형이므로,
    최소제곱을 그대로 쓸 수 있다. 다항회귀도 선형모형인 까닭이다.
    """
    x_poly = np.concatenate([x**k for k in range(1, d+1)], axis=1)
    model = LinearRegression()
    model.fit(x_poly, y)
    y_pred = model.predict(x_poly)
    return model, y_pred

x, y = generate_data(d=2)
model, y_pred = perform_regression(x, y, d=2)
plot_regression_and_residuals(x, y, y_pred)

import numpy as np
import statsmodels.api as sm

xf = x.ravel()
e2 = y - y_pred                                        # 이차 적합의 잔차
_, y_lin = perform_regression(x, y, d=1)
e1 = y - y_lin                                         # 보기 3 의 직선 적합 잔차

print(f"이차 적합 계수 = (절편 {model.intercept_:.4f}, "
      f"x {model.coef_[0]:.4f}, x^2 {model.coef_[1]:.4f})   참값 (1, 2, 3)")
print(f"\n{'':>10}{'직선':>12}{'이차':>12}")
print(f"{'RSS':>10}{(e1 ** 2).sum():>12.2f}{(e2 ** 2).sum():>12.2f}")
print(f"{'sigma_hat':>10}{np.sqrt((e1 ** 2).sum() / 48):>12.4f}"
      f"{np.sqrt((e2 ** 2).sum() / 47):>12.4f}   (참값 3)")

# 보기 3 의 곡률 조각이 정말 사라졌는가
curv = 3 * xf ** 2
print(f"\ncorr(e, x^2):  직선 {np.corrcoef(e1, curv)[0, 1]:+.4f},"
      f"  이차 {np.corrcoef(e2, curv)[0, 1]:+.2e}")
aux1 = sm.OLS(e1, sm.add_constant(np.column_stack([xf, xf ** 2]))).fit()
print(f"e 를 (1, x, x^2) 에 회귀한 R^2:  직선 {aux1.rsquared:.4f}")

이차 적합 계수 = (절편 0.5384, x 2.8580, x^2 3.2162)   참값 (1, 2, 3)

                    직선          이차
       RSS     1518.09      282.93
 sigma_hat      5.6238      2.4535   (참값 3)

corr(e, x^2):  직선 +0.8978,  이차 +4.02e-16
e 를 (1, x, x^2) 에 회귀한 R^2:  직선 0.8136

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import statsmodels.api as sm
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})

X_with_const = sm.add_constant(data['X'])
model = sm.OLS(data['y'], X_with_const).fit()

data['Fitted'] = model.fittedvalues
data['Residuals'] = model.resid
# 척도-위치 그림을 위해 잔차를 표준화하고 절댓값의 제곱근을 취한다.
data['Standardized Residuals'] = data['Residuals'] / np.std(data['Residuals'])
data['Sqrt Abs Standardized Residuals'] = np.sqrt(np.abs(data['Standardized Residuals']))

fig, (ax0, ax1, ax2) = plt.subplots(1, 3, figsize=(18, 5))

# 회귀 그림
ax0.scatter(data['X'], data['y'], alpha=0.7, label='Data Points')
ax0.plot(data['X'], data['Fitted'], color='orange', label='Regression Line')
ax0.set_title('Regression Plot')
ax0.set_xlabel('Predictor (X)')
ax0.set_ylabel('Response (y)')
ax0.legend()

# 잔차 대 적합값
ax1.scatter(data['Fitted'], data['Residuals'], alpha=0.7)
ax1.axhline(y=0, color='r', linestyle='--')
ax1.set_title('Residuals vs. Fitted Values Plot')
ax1.set_xlabel('Fitted Values')
ax1.set_ylabel('Residuals')

# 척도-위치 그림은 등분산성만 본다. 부호를 없앴으므로 점들이 이루는 띠의
# 높이가 일정한지만 보면 된다.
ax2.scatter(data['Fitted'], data['Sqrt Abs Standardized Residuals'], alpha=0.7)
ax2.axhline(y=0, color='r', linestyle='--')
ax2.set_title('Scale-Location Plot')
ax2.set_xlabel('Fitted Values')
ax2.set_ylabel(r'$\sqrt{|\text{Standardized Residuals}|}$')

plt.tight_layout()
plt.show()

import math
import numpy as np
from scipy.stats import skew
from statsmodels.nonparametric.smoothers_lowess import lowess

theory = 2 ** 0.25 * math.gamma(0.75) / math.sqrt(math.pi)
v = data['Sqrt Abs Standardized Residuals'].values
print(f"E[sqrt|Z|] = 2^(1/4) Gamma(3/4) / sqrt(pi) = {theory:.6f}")
print(f"그림의 세로축 평균                        = {v.mean():.6f}")
print(f"(n = {len(v)}, 표준오차 = {v.std(ddof=1) / math.sqrt(len(v)):.6f})")

# 띠가 평평한가: 평활선의 양 끝
f = data['Fitted'].values
sm_line = lowess(v, f, frac=0.6)
print(f"\n평활선 왼쪽 끝 {sm_line[0, 1]:.4f}  →  오른쪽 끝 {sm_line[-1, 1]:.4f}")
print(f"제곱한 비 (sd 의 비로 읽는다) = {(sm_line[-1, 1] / sm_line[0, 1]) ** 2:.4f}")

# 제곱근이 치우침을 없앴는가
print(f"\n왜도:  |r| = {skew(np.abs(data['Standardized Residuals'].values)):+.4f},"
      f"  sqrt|r| = {skew(v):+.4f}")

E[sqrt|Z|] = 2^(1/4) Gamma(3/4) / sqrt(pi) = 0.822179
그림의 세로축 평균                        = 0.843058
(n = 100, 표준오차 = 0.033396)

평활선 왼쪽 끝 0.7258  →  오른쪽 끝 0.6442
제곱한 비 (sd 의 비로 읽는다) = 0.7876

왜도:  |r| = +0.5204,  sqrt|r| = -0.0001