잔차 분석¶

import numpy as np
import pandas as pd
from statsmodels.formula.api import ols

# 이 페이지의 진단은 모두 아래 모형 하나를 놓고 수행한다.
# 집단마다 표준편차를 1.0, 1.3, 1.6으로 다르게 주었고,
# 집단 C에 이상점을 하나 심어 두었다.
rng = np.random.default_rng(42)
n = 20
response = np.concatenate([
    rng.normal(10.0, 1.0, n),
    rng.normal(10.8, 1.3, n),
    rng.normal(12.0, 1.6, n),
])
response[-1] = 20.0                     # 마지막 관측값을 이상점으로 만든다
data = pd.DataFrame({
    "group": np.repeat(["A", "B", "C"], n),
    "response": response,
})
group1 = data.loc[data["group"] == "A", "response"]
group2 = data.loc[data["group"] == "B", "response"]
group3 = data.loc[data["group"] == "C", "response"]

model = ols("response ~ C(group)", data=data).fit()

print(data.groupby("group").response.agg(["count", "mean", "std"]).round(3))
print(f"\nF = {model.fvalue:.4f}, p = {model.f_pvalue:.4f}")

       count    mean    std
group                      
A         20   9.967  0.870
B         20  10.942  1.034
C         20  12.513  2.077

F = 16.1314, p = 0.0000

import numpy as np

# 집단별 요약만으로 MSE, MSB, F 를 모두 재구성한다.
n, k = 20, 3
means = np.array([9.967068, 10.942301, 12.513460])
sds = np.array([0.870178, 1.033538, 2.077021])
N = n * k

MSE = (sds ** 2).sum() / k                       # 균형설계에서는 단순평균
grand = means.mean()                             # 균형설계에서는 집단평균의 평균
MSB = n * ((means - grand) ** 2).sum() / (k - 1)
print(f"집단분산 s_i^2 = {np.round(sds ** 2, 4)}")
print(f"MSE  = {MSE:.6f},   sigma_hat = {np.sqrt(MSE):.6f}")
print(f"전체평균 = {grand:.6f}")
print(f"MSB  = {MSB:.6f}")
print(f"F    = MSB/MSE = {MSB / MSE:.4f}      (쪽의 적합 결과 16.1314)")

print(f"\n합동분산에서 각 집단이 차지하는 몫")
for g, v in zip("ABC", sds ** 2):
    print(f"  {g}: {v:>8.4f} / {(sds ** 2).sum():.4f} = {v / (sds ** 2).sum():>7.2%}")

# 이상점을 빼면 합동분산이 어떻게 되는가
sd_c_clean = 1.145019
MSE_clean = (sds[0] ** 2 + sds[1] ** 2 + sd_c_clean ** 2) / k
print(f"\n집단 C 의 s 를 이상점 없는 값 {sd_c_clean} 로 바꾸면")
print(f"  MSE = {MSE_clean:.6f},  sigma_hat = {np.sqrt(MSE_clean):.6f}  "
      f"({np.sqrt(MSE) / np.sqrt(MSE_clean):.3f} 배 작다)")

집단분산 s_i^2 = [0.7572 1.0682 4.314 ]
MSE  = 2.046476,   sigma_hat = 1.430551
전체평균 = 11.140943
MSB  = 33.012441
F    = MSB/MSE = 16.1314      (쪽의 적합 결과 16.1314)

합동분산에서 각 집단이 차지하는 몫
  A:   0.7572 / 6.1394 =  12.33%
  B:   1.0682 / 6.1394 =  17.40%
  C:   4.3140 / 6.1394 =  70.27%

집단 C 의 s 를 이상점 없는 값 1.145019 로 바꾸면
  MSE = 1.045493,  sigma_hat = 1.022494  (1.399 배 작다)

import matplotlib.pyplot as plt

# 잔차 대 적합값 그림은 진단의 출발점이다. 점들이 0 선 둘레에 폭을 일정하게
# 유지하며 흩어져 있으면 좋다. 깔때기 모양이면 등분산이 깨진 것이고, 굽은
# 모양이면 모형이 놓친 구조가 남아 있다는 뜻이다.
plt.scatter(model.fittedvalues, model.resid, alpha=0.6)
plt.axhline(y=0, color='r', linestyle='--')
plt.xlabel("Fitted Values")
plt.ylabel("Residuals")
plt.title("Residuals vs. Fitted Values")
plt.show()

import numpy as np
import pandas as pd
from statsmodels.formula.api import ols

rng = np.random.default_rng(42)
n = 20
response = np.concatenate([
    rng.normal(10.0, 1.0, n), rng.normal(10.8, 1.3, n), rng.normal(12.0, 1.6, n)])
response[-1] = 20.0
data = pd.DataFrame({"group": np.repeat(["A", "B", "C"], n), "response": response})
model = ols("response ~ C(group)", data=data).fit()
e = model.resid.values
fit = model.fittedvalues.values

print(f"서로 다른 적합값: {np.unique(np.round(fit, 6))}")
print(f"\n{'띠':>3}{'가로 위치':>11}{'합':>12}{'최소':>9}{'최대':>9}"
      f"{'범위':>8}{'잔차 sd':>9}{'s_i*sqrt(19/20)':>17}")
for i, g in enumerate("ABC"):
    ei = e[i * n:(i + 1) * n]
    si = response[i * n:(i + 1) * n].std(ddof=1)
    print(f"{g:>3}{fit[i * n]:>11.4f}{ei.sum():>12.2e}{ei.min():>9.4f}{ei.max():>9.4f}"
          f"{ei.max() - ei.min():>8.4f}{ei.std(ddof=0):>9.4f}{si * np.sqrt((n - 1) / n):>17.4f}")

c = np.sort(e[2 * n:])[::-1]
print(f"\n집단 C 의 큰 잔차 넷: {np.round(c[:4], 4)}")
print(f"이상점은 둘째로 큰 잔차의 {c[0] / c[1]:.2f} 배")
print(f"이상점을 뺀 집단 C 잔차의 범위: "
      f"[{e[2 * n:-1].min():.4f}, {e[2 * n:-1].max():.4f}]")

서로 다른 적합값: [ 9.967068 10.942301 12.51346 ]

  띠      가로 위치           합       최소       최대      범위    잔차 sd  s_i*sqrt(19/20)
  A     9.9671   -7.11e-14  -1.9181   1.1602  3.0783   0.8481           0.8481
  B    10.9423   -5.33e-15  -1.2345   2.6418  3.8763   1.0074           1.0074
  C    12.5135   -1.81e-13  -2.8449   7.4865 10.3314   2.0244           2.0244

집단 C 의 큰 잔차 넷: [7.4865 1.8784 0.8808 0.6757]
이상점은 둘째로 큰 잔차의 3.99 배
이상점을 뺀 집단 C 잔차의 범위: [-2.8449, 1.8784]

import numpy as np

# 표준화 잔차는 잔차를 그 표준오차로 나눈 것이다. 단위가 사라지므로
# 어느 자료에서든 ±2 를 같은 뜻으로 읽을 수 있다.
influence = model.get_influence()
standardized_resid = influence.resid_studentized_internal

plt.scatter(model.fittedvalues, standardized_resid, alpha=0.6)
plt.axhline(y=0, color='r', linestyle='--')
plt.axhline(y=2, color='gray', linestyle=':', alpha=0.5)
plt.axhline(y=-2, color='gray', linestyle=':', alpha=0.5)
plt.xlabel("Fitted Values")
plt.ylabel("Standardized Residuals")
plt.title("Standardized Residuals vs. Fitted Values")
plt.show()

import numpy as np
import pandas as pd
from scipy import stats
from statsmodels.formula.api import ols

rng = np.random.default_rng(42)
n, k = 20, 3
response = np.concatenate([
    rng.normal(10.0, 1.0, n), rng.normal(10.8, 1.3, n), rng.normal(12.0, 1.6, n)])
response[-1] = 20.0
data = pd.DataFrame({"group": np.repeat(["A", "B", "C"], n), "response": response})
model = ols("response ~ C(group)", data=data).fit()
N = 3 * n

inf = model.get_influence()
h = inf.hat_matrix_diag
r = inf.resid_studentized_internal
e = model.resid.values
sigma = np.sqrt(model.mse_resid)

print(f"지렛값이 갖는 서로 다른 값: {np.unique(np.round(h, 12))}   (1/n = {1 / n})")
print(f"지렛값의 합 = {h.sum():.6f}   (= 모수 개수 k = {k})")

const = sigma * np.sqrt(1 - 1 / n)
print(f"\nsigma_hat = {sigma:.6f},  sqrt(1-1/n) = {np.sqrt(1 - 1 / n):.6f}")
print(f"나누는 상수 = {const:.6f}")
print(f"\nsum r^2 = {(r ** 2).sum():.6f}   (이론 (N-k)/(1-1/n) = {(N - k) / (1 - 1 / n):.6f} = N)")

print(f"\n이상점: e = {e[-1]:.6f}")
print(f"  손계산 r = e/상수 = {e[-1] / const:.6f}")
print(f"  statsmodels  r = {r[-1]:.6f}")

t_ext = inf.resid_studentized_external
t_hand = r[-1] * np.sqrt((N - k - 1) / (N - k - r[-1] ** 2))
print(f"\n외부 스튜던트화 t = r*sqrt((N-k-1)/(N-k-r^2))")
print(f"  손계산 = {t_hand:.6f},  statsmodels = {t_ext[-1]:.6f}")
p_raw = 2 * stats.t.sf(abs(t_ext[-1]), N - k - 1)
print(f"  t(56) 양측 p = {p_raw:.4e},  본페로니 보정 (x{N}) = {min(1, N * p_raw):.4e}")

print(f"\n|r| > 2 인 관측: {np.where(np.abs(r) > 2)[0]} -> r = {np.round(r[np.abs(r) > 2], 4)}")
print(f"|r| > 3 인 관측 수 = {(np.abs(r) > 3).sum()}")

지렛값이 갖는 서로 다른 값: [0.05]   (1/n = 0.05)
지렛값의 합 = 3.000000   (= 모수 개수 k = 3)

sigma_hat = 1.430551,  sqrt(1-1/n) = 0.974679
나누는 상수 = 1.394329

sum r^2 = 60.000000   (이론 (N-k)/(1-1/n) = 60.000000 = N)

이상점: e = 7.486540
  손계산 r = e/상수 = 5.369280
  statsmodels  r = 5.369280

외부 스튜던트화 t = r*sqrt((N-k-1)/(N-k-r^2))
  손계산 = 7.570250,  statsmodels = 7.570250
  t(56) 양측 p = 3.9517e-10,  본페로니 보정 (x60) = 2.3710e-08

|r| > 2 인 관측: [52 59] -> r = [-2.0403  5.3693]
|r| > 3 인 관측 수 = 1

# 척도-위치 그림은 부호를 없애고 크기만 본다. 제곱근을 씌우는 것은 큰 값이
# 그림을 독차지하지 않게 하려는 것이다. 추세선이 평평하면 등분산이다.
plt.scatter(model.fittedvalues, np.sqrt(np.abs(standardized_resid)), alpha=0.6)
plt.xlabel("Fitted Values")
plt.ylabel(r"$\sqrt{|\mathrm{Standardized\ Residuals}|}$")
plt.title("Scale-Location Plot")
plt.show()

import numpy as np
import pandas as pd
from math import gamma, pi
from statsmodels.formula.api import ols

rng = np.random.default_rng(42)
n = 20
response = np.concatenate([
    rng.normal(10.0, 1.0, n), rng.normal(10.8, 1.3, n), rng.normal(12.0, 1.6, n)])
response[-1] = 20.0
data = pd.DataFrame({"group": np.repeat(["A", "B", "C"], n), "response": response})
model = ols("response ~ C(group)", data=data).fit()
r = model.get_influence().resid_studentized_internal
sq = np.sqrt(np.abs(r))
sigma = np.sqrt(model.mse_resid)

ref = 2 ** 0.25 * gamma(0.75) / np.sqrt(pi)
print(f"기준 높이 E sqrt|Z| = 2^(1/4) Gamma(3/4) / sqrt(pi) = {ref:.6f}")
print(f"몬테카를로 확인: {np.sqrt(np.abs(np.random.default_rng(0).normal(size=2_000_000))).mean():.6f}")

print(f"\n{'띠':>3}{'s_i':>9}{'s_i/sigma':>11}{'예측 평균':>11}{'실제 평균':>11}{'실제 중앙값':>12}")
for i, g in enumerate("ABC"):
    si = response[i * n:(i + 1) * n].std(ddof=1)
    v = sq[i * n:(i + 1) * n]
    print(f"{g:>3}{si:>9.4f}{si / sigma:>11.4f}{ref * np.sqrt(si / sigma):>11.4f}"
          f"{v.mean():>11.4f}{np.median(v):>12.4f}")

# 집단 C 에서 이상점을 빼면
si_c = response[2 * n:-1].std(ddof=1)
v = sq[2 * n:-1]
print(f"\n집단 C 에서 이상점을 뺀 19 개:")
print(f"  s = {si_c:.4f},  예측 평균 = {ref * np.sqrt(si_c / sigma):.4f},  실제 평균 = {v.mean():.4f}")
print(f"\n이상점의 sqrt|r| = {sq[-1]:.4f}  (기준 높이의 {sq[-1] / ref:.2f} 배)")

기준 높이 E sqrt|Z| = 2^(1/4) Gamma(3/4) / sqrt(pi) = 0.822179
몬테카를로 확인: 0.822008

  띠      s_i  s_i/sigma      예측 평균      실제 평균      실제 중앙값
  A   0.8702     0.6083     0.6412     0.6554      0.7647
  B   1.0335     0.7225     0.6988     0.7459      0.7008
  C   2.0770     1.4519     0.9907     0.8134      0.7455

집단 C 에서 이상점을 뺀 19 개:
  s = 1.1296,  예측 평균 = 0.7306,  실제 평균 = 0.7342

이상점의 sqrt|r| = 2.3172  (기준 높이의 2.82 배)

import warnings
warnings.filterwarnings("ignore")

import numpy as np
import pandas as pd
from statsmodels.formula.api import ols

rng = np.random.default_rng(9001)
B = 30_000
NS = [5, 15, 40]
acc = {n: [] for n in NS}
for _ in range(B):
    rows = [pd.DataFrame({"y": rng.normal(0, 1.0, n), "g": f"G{i}"})
            for i, n in enumerate(NS)]
    df = pd.concat(rows, ignore_index=True)
    r = ols("y ~ C(g)", data=df).fit().resid.values
    i = 0
    for n in NS:
        acc[n].append(r[i:i + n])
        i += n

print("모든 집단의 참 σ = 1.0 (집단 크기만 다름)")
print(f"{'집단 크기 n':>11s} {'잔차의 SD':>10s} {'이론 σ√(1-1/n)':>15s} {'h=1/n':>8s}")
for n in NS:
    a = np.concatenate(acc[n])
    print(f"{n:11d} {a.std(ddof=0):10.4f} {np.sqrt(1 - 1 / n):15.4f} {1 / n:8.4f}")

모든 집단의 참 σ = 1.0 (집단 크기만 다름)
    집단 크기 n     잔차의 SD    이론 σ√(1-1/n)    h=1/n
          5     0.8960          0.8944   0.2000
         15     0.9658          0.9661   0.0667
         40     0.9880          0.9874   0.0250

import warnings
warnings.filterwarnings("ignore")

import numpy as np
import pandas as pd
from scipy import stats
from statsmodels.formula.api import ols

rng = np.random.default_rng(9002)
cases = {}

rows = []                                   # (1) 숨은 이분 잠재변수
for g in range(3):
    h = rng.binomial(1, 0.5, 60)
    rows.append(pd.DataFrame({"y": 10 + 2 * g + 4 * h + rng.normal(0, 1, 60),
                              "g": f"G{g}"}))
cases["숨은 하위집단 (이분 잠재변수)"] = pd.concat(rows, ignore_index=True)

rows = [pd.DataFrame({"y": 10 + 2 * g + rng.normal(0, 1, 60), "g": f"G{g}"})
        for g in range(3)]                  # (2) 문제 없음
cases["문제 없음 (비교용)"] = pd.concat(rows, ignore_index=True)

rows = [pd.DataFrame({"y": np.clip(10 + 2 * g + rng.normal(0, 3, 60), 8, 16),
                      "g": f"G{g}"}) for g in range(3)]    # (3) 절단
cases["천장·바닥 효과 (절단)"] = pd.concat(rows, ignore_index=True)

rows = [pd.DataFrame({"y": rng.binomial(1, 0.3 + 0.15 * g, 60).astype(float),
                      "g": f"G{g}"}) for g in range(3)]    # (4) 이분 반응
cases["사실은 이분 반응"] = pd.concat(rows, ignore_index=True)

def center_mass(r):
    """표준화 잔차가 0 근처에 얼마나 몰려 있는가. 정규면 약 0.383."""
    z = (r - r.mean()) / r.std(ddof=1)
    return np.mean(np.abs(z) < 0.5)

print(f"{'상황':>26s} {'치우침':>7s} {'초과첨도':>8s} {'중심밀도':>9s} {'샤피로 p':>9s}")
for lab, df in cases.items():
    r = ols("y ~ C(g)", data=df).fit().resid.values
    print(f"{lab:>26s} {stats.skew(r):7.3f} {stats.kurtosis(r):8.3f} "
          f"{center_mass(r):9.3f} {stats.shapiro(r).pvalue:9.4f}")
print(f"\n  (정규 잔차라면 |z|<0.5 의 비율이 약 {2 * stats.norm.cdf(0.5) - 1:.3f})")

                        상황     치우침     초과첨도      중심밀도     샤피로 p
         숨은 하위집단 (이분 잠재변수)  -0.289   -1.116     0.228    0.0000
               문제 없음 (비교용)  -0.199   -0.027     0.400    0.1374
             천장·바닥 효과 (절단)  -0.126   -0.925     0.278    0.0008
                 사실은 이분 반응   0.172   -1.778     0.000    0.0000

  (정규 잔차라면 |z|<0.5 의 비율이 약 0.383)

import numpy as np

rng = np.random.default_rng(9003)
B = 20_000
print("모든 집단의 참 σ 가 같을 때, 표본 SD 의 최대/최소 비")
print(f"{'k':>3s} {'집단당 n':>8s} {'비의 중앙값':>11s} {'90 백분위':>10s} "
      f"{'비>2 확률':>10s} {'비>3 확률':>10s}")
for k, n in [(3, 5), (3, 10), (3, 20), (3, 50), (5, 10), (8, 10)]:
    r = []
    for _ in range(B // 4):
        s = np.array([rng.normal(0, 1, n).std(ddof=1) for _ in range(k)])
        r.append(s.max() / s.min())
    r = np.array(r)
    print(f"{k:3d} {n:8d} {np.median(r):11.3f} {np.quantile(r, 0.9):10.3f} "
          f"{(r > 2).mean():10.4f} {(r > 3).mean():10.4f}")

모든 집단의 참 σ 가 같을 때, 표본 SD 의 최대/최소 비
  k    집단당 n      비의 중앙값     90 백분위     비>2 확률     비>3 확률
  3        5       1.831      3.245     0.4092     0.1294
  3       10       1.475      2.067     0.1184     0.0078
  3       20       1.305      1.629     0.0112     0.0000
  3       50       1.176      1.347     0.0000     0.0000
  5       10       1.743      2.392     0.2756     0.0212
  8       10       1.991      2.693     0.4888     0.0504

import warnings
warnings.filterwarnings("ignore")

import numpy as np
import pandas as pd
from statsmodels.formula.api import ols

rng = np.random.default_rng(9001)
B = 30_000
NS = [5, 15, 40]
acc = {n: [] for n in NS}
for _ in range(B):
    rows = [pd.DataFrame({"y": rng.normal(0, 1.0, n), "g": f"G{i}"})
            for i, n in enumerate(NS)]
    df = pd.concat(rows, ignore_index=True)
    r = ols("y ~ C(g)", data=df).fit().resid.values
    i = 0
    for n in NS:
        acc[n].append(r[i:i + n])
        i += n

a5 = np.array(acc[5])
a40 = np.array(acc[40])
print(f"집단 크기 5 의 e1, e2 상관  = {np.corrcoef(a5[:, 0], a5[:, 1])[0, 1]:.4f}"
      f"   (이론 -1/(n-1) = {-1 / 4:.4f})")
print(f"집단 크기 40 의 e1, e2 상관 = "
      f"{np.corrcoef(a40[:, 0], a40[:, 1])[0, 1]:.4f}   (이론 {-1 / 39:.4f})")
print(f"집단 크기 5 의 잔차합: 평균 {a5.sum(1).mean():.2e}, "
      f"SD {a5.sum(1).std():.2e}   (항등적으로 0)")

집단 크기 5 의 e1, e2 상관  = -0.2447   (이론 -1/(n-1) = -0.2500)
집단 크기 40 의 e1, e2 상관 = -0.0230   (이론 -0.0256)
집단 크기 5 의 잔차합: 평균 -1.53e-17, SD 7.67e-15   (항등적으로 0)

import warnings
warnings.filterwarnings("ignore")

import numpy as np
import pandas as pd
from scipy import stats
from statsmodels.formula.api import ols

rng = np.random.default_rng(9004)
B = 8_000
print("(가) 정규성 검정에는 차이가 거의 없다 (모두 정규, 명목 0.05)")
print(f"{'집단 크기':>18s} {'원 잔차':>9s} {'표준화 잔차':>11s}")
for ns in [[20, 20, 20], [5, 20, 60], [3, 10, 100]]:
    a = b = 0
    for _ in range(B):
        rows = [pd.DataFrame({"y": rng.normal(0, 1, n), "g": f"G{i}"})
                for i, n in enumerate(ns)]
        df = pd.concat(rows, ignore_index=True)
        f = ols("y ~ C(g)", data=df).fit()
        a += stats.shapiro(f.resid).pvalue < 0.05
        b += stats.shapiro(
            f.get_influence().resid_studentized_internal).pvalue < 0.05
    print(f"{str(ns):>18s} {a / B:9.4f} {b / B:11.4f}")

rng = np.random.default_rng(9005)
B = 20_000
NS = [4, 12, 40]
cnt_raw = np.zeros(3)
cnt_std = np.zeros(3)
for _ in range(B):
    rows = [pd.DataFrame({"y": rng.normal(0, 1, n), "g": f"G{i}"})
            for i, n in enumerate(NS)]
    df = pd.concat(rows, ignore_index=True)
    f = ols("y ~ C(g)", data=df).fit()
    gi = np.repeat([0, 1, 2], NS)
    cnt_raw[gi[np.argmax(np.abs(f.resid.values))]] += 1
    cnt_std[gi[np.argmax(np.abs(
        f.get_influence().resid_studentized_internal))]] += 1

print("\n(나) 가장 큰 |잔차| 가 어느 집단에서 나오는가")
print(f"{'집단':>8s} {'n':>4s} {'기대 비율':>9s} {'원 잔차 최대':>12s} {'표준화 최대':>11s}")
for i, n in enumerate(NS):
    print(f"{'G' + str(i):>8s} {n:4d} {n / sum(NS):9.4f} "
          f"{cnt_raw[i] / B:12.4f} {cnt_std[i] / B:11.4f}")

(가) 정규성 검정에는 차이가 거의 없다 (모두 정규, 명목 0.05)
             집단 크기      원 잔차      표준화 잔차
      [20, 20, 20]    0.0512      0.0512
       [5, 20, 60]    0.0534      0.0516
      [3, 10, 100]    0.0511      0.0500

(나) 가장 큰 |잔차| 가 어느 집단에서 나오는가
      집단    n     기대 비율      원 잔차 최대      표준화 최대
      G0    4    0.0714       0.0329      0.0664
      G1   12    0.2143       0.1900      0.2114
      G2   40    0.7143       0.7772      0.7221

fit.resid                                   원 잔차
fit.get_influence().resid_studentized_internal    표준화
fit.get_influence().resid_studentized_external    삭제(이상점 탐지에 최선)

1. 집단별 n, 평균, s 를 표로 먼저 본다
     ↓
2. 집단별 상자그림 — 산포·이상점·치우침
     ↓  (SD 비를 n 과 함께 해석 — 연습문제 6)
3. 잔차의 Q-Q — 벗어남의 모양
     ↓
4. 표준화 잔차의 최댓값 — 본페로니 임계값과 비교
     ↓
5. 수집 순서가 있으면 순서에 대한 잔차
     ↓
6. 이상한 것이 있으면 그 관측을 조사하고,
   빼고 다시 적합해 결론이 바뀌는지 확인