영향점¶

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})

full = ols("response ~ C(group)", data=data).fit()
drop = ols("response ~ C(group)", data=data.drop(index=59)).fit()
e59 = full.resid.values[-1]

print(f"e_59 = {e59:.6f}")
print(f"공식  집단 C 평균의 변화 = -e_59/(n-1) = {-e59 / (n - 1):.6f}")
print(f"실제  {response[40:59].mean():.6f} - {response[40:].mean():.6f} = "
      f"{response[40:59].mean() - response[40:].mean():.6f}")

print(f"\n{'':>8}{'F':>10}{'p':>12}{'sigma_hat':>12}{'eta^2':>9}{'평균 C':>11}")
print(f"{'전체 60':>8}{full.fvalue:>10.4f}{full.f_pvalue:>12.3e}"
      f"{np.sqrt(full.mse_resid):>12.4f}{full.rsquared:>9.4f}{response[40:].mean():>11.4f}")
print(f"{'59 빼고':>8}{drop.fvalue:>10.4f}{drop.f_pvalue:>12.3e}"
      f"{np.sqrt(drop.mse_resid):>12.4f}{drop.rsquared:>9.4f}{response[40:59].mean():>11.4f}")

e_59 = 7.486540
공식  집단 C 평균의 변화 = -e_59/(n-1) = -0.394028
실제  12.119432 - 12.513460 = -0.394028

                 F           p   sigma_hat    eta^2       평균 C
   전체 60   16.1314   2.808e-06      1.4306   0.3614    12.5135
   59 빼고   21.9577   9.111e-08      1.0146   0.4395    12.1194

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
import matplotlib.pyplot as plt

# 관측값마다 Cook 의 거리를 막대로 세운다. 유독 솟은 막대가 있는지를 본다.
influence = model.get_influence()
cooks_d = influence.cooks_distance[0]

plt.stem(range(len(cooks_d)), cooks_d, markerfmt=",")
plt.xlabel("Observation Index")
plt.ylabel("Cook's Distance")
plt.title("Cook's Distance")
plt.axhline(y=4/len(cooks_d), color='r', linestyle='--', label=f'Threshold = {4/len(cooks_d):.3f}')
plt.legend()
plt.show()

print(f"threshold = {4/len(cooks_d):.4f}")
print(f"max Cook's D = {cooks_d.max():.4f} at obs {cooks_d.argmax()}")
print(f"flagged = {np.where(cooks_d > 4/len(cooks_d))[0]}")

threshold = 0.0667
max Cook's D = 0.5058 at obs 59
flagged = [52 59]

import numpy as np
import pandas as pd
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()
e = model.resid.values
fit = model.fittedvalues.values
s2 = model.mse_resid
D = model.get_influence().cooks_distance[0]

def cook_by_deletion(i):
    """정의 그대로: i 를 빼고 다시 적합해 60 개 적합값이 얼마나 움직였는지 잰다."""
    refit = ols("response ~ C(group)", data=data.drop(index=i)).fit()
    moved = fit - refit.predict(data).values
    return (moved ** 2).sum() / (k * s2)

print(f"{'관측':>6}{'e_i':>11}{'삭제 정의':>13}{'공식':>13}{'cooks_distance':>16}")
for i in [59, 52, 30, 0]:
    formula = n * e[i] ** 2 / (k * s2 * (n - 1) ** 2)
    print(f"{i:>6}{e[i]:>11.4f}{cook_by_deletion(i):>13.8f}{formula:>13.8f}{D[i]:>16.8f}")

# 59 번을 뺐을 때 실제로 움직인 적합값
refit = ols("response ~ C(group)", data=data.drop(index=59)).fit()
moved = fit - refit.predict(data).values
print(f"\n59 번을 뺐을 때 적합값이 움직인 집단: "
      f"{sorted(set(data['group'][np.abs(moved) > 1e-12]))}")
print(f"  움직인 크기 = {moved[-1]:.6f}  (= e_59/(n-1) = {e[-1] / (n - 1):.6f})")
print(f"  움직인 관측 수 = {(np.abs(moved) > 1e-12).sum()}")
print(f"  제곱합 = {(moved ** 2).sum():.6f} = n*(e/(n-1))^2 = {n * (e[-1] / (n - 1)) ** 2:.6f}")
print(f"  p*sigma^2 = {k * s2:.6f}")
print(f"  나누면 {(moved ** 2).sum() / (k * s2):.8f}")

    관측        e_i        삭제 정의           공식  cooks_distance
    59     7.4865   0.50577484   0.50577484      0.50577484
    52    -2.8449   0.07303512   0.07303512      0.07303512
    30     2.6418   0.06298082   0.06298082      0.06298082
     0     0.3376   0.00102879   0.00102879      0.00102879

59 번을 뺐을 때 적합값이 움직인 집단: ['C']
  움직인 크기 = 0.394028  (= e_59/(n-1) = 0.394028)
  움직인 관측 수 = 20
  제곱합 = 3.105168 = n*(e/(n-1))^2 = 3.105168
  p*sigma^2 = 6.139427
  나누면 0.50577484

# 지렛값은 설명변수 쪽에서 그 점이 얼마나 외따로 있는지를 잰다.
# 지렛값이 크고 잔차도 큰 점이 가장 위험하다. 그 둘을 곱해 놓은 것이
# 앞의 Cook 거리라고 보면 된다.
leverage = influence.hat_matrix_diag

plt.scatter(leverage, influence.resid_studentized_internal, alpha=0.6)
plt.xlabel("Leverage")
plt.ylabel("Studentized Residuals")
plt.title("Leverage vs. Studentized Residuals")
plt.axhline(y=0, color='r', linestyle='--')
plt.show()

print(f"leverage: min = {leverage.min():.4f}, max = {leverage.max():.4f}")

leverage: min = 0.0500, max = 0.0500

import numpy as np
import pandas as pd
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})
h = ols("response ~ C(group)", data=data).fit().get_influence().hat_matrix_diag
N = len(h)

print(f"균형설계 n=(20,20,20)")
print(f"  서로 다른 h = {np.unique(np.round(h, 12))}")
print(f"  sum h = {h.sum():.10f}  (= p = {k})")
print(f"  mean h = {h.mean():.10f}  (= p/N = {k / N})")
print(f"  std h = {h.std(ddof=1):.3e}")

# 불균형 설계에서는 어떻게 되는가
ns = [5, 20, 35]
rg = np.random.default_rng(11)
y2 = np.concatenate([rg.normal(10.0, 1.0, ns[0]),
                     rg.normal(10.8, 1.0, ns[1]),
                     rg.normal(12.0, 1.0, ns[2])])
d2 = pd.DataFrame({"group": np.repeat(["A", "B", "C"], ns), "response": y2})
inf2 = ols("response ~ C(group)", data=d2).fit().get_influence()
h2, r2 = inf2.hat_matrix_diag, inf2.resid_studentized_internal
D2 = inf2.cooks_distance[0]

print(f"\n불균형설계 n=(5,20,35)")
print(f"{'집단':>5}{'n_i':>5}{'h_ii':>10}{'1/n_i':>10}{'n_i*h_ii':>10}")
for i, (g, m) in enumerate(zip("ABC", ns)):
    hi = h2[np.array(d2["group"]) == g][0]
    print(f"{g:>5}{m:>5}{hi:>10.6f}{1 / m:>10.6f}{m * hi:>10.6f}")
print(f"  sum h = {h2.sum():.10f}  (= p = {k})")
print(f"  같은 |r| = 2 라도 Cook 거리는 집단마다 다르다:")
for g, m in zip("ABC", ns):
    hi = h2[np.array(d2["group"]) == g][0]
    print(f"    집단 {g}: D = (2^2/3)*h/(1-h) = {(4 / k) * hi / (1 - hi):.6f}")
print(f"  문턱 4/N = {4 / len(h2):.6f}")

균형설계 n=(20,20,20)
  서로 다른 h = [0.05]
  sum h = 3.0000000000  (= p = 3)
  mean h = 0.0500000000  (= p/N = 0.05)
  std h = 2.099e-17

불균형설계 n=(5,20,35)
   집단  n_i      h_ii     1/n_i  n_i*h_ii
    A    5  0.200000  0.200000  1.000000
    B   20  0.050000  0.050000  1.000000
    C   35  0.028571  0.028571  1.000000
  sum h = 3.0000000000  (= p = 3)
  같은 |r| = 2 라도 Cook 거리는 집단마다 다르다:
    집단 A: D = (2^2/3)*h/(1-h) = 0.333333
    집단 B: D = (2^2/3)*h/(1-h) = 0.070175
    집단 C: D = (2^2/3)*h/(1-h) = 0.039216
  문턱 4/N = 0.066667

import statsmodels.api as sm
from statsmodels.formula.api import ols

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

influence = model.get_influence()
# summary_frame은 진단량을 한 표에 모아 준다.
summary = influence.summary_frame()
print(summary[['hat_diag', 'cooks_d', 'dffits', 'student_resid']].describe())

           hat_diag    cooks_d     dffits  student_resid
count  6.000000e+01  60.000000  60.000000      60.000000
mean   5.000000e-02   0.017544   0.008298       0.036170
std    5.411161e-17   0.065539   0.281456       1.226837
min    5.000000e-02   0.000002  -0.481894      -2.100527
25%    5.000000e-02   0.001276  -0.139287      -0.607140
50%    5.000000e-02   0.005002  -0.015728      -0.068558
75%    5.000000e-02   0.012359   0.097065       0.423097
max    5.000000e-02   0.505775   1.736734       7.570250

import numpy as np
import pandas as pd
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})
inf = ols("response ~ group", data=data).fit().get_influence()
sf = inf.summary_frame()
N = 3 * n
r, t = inf.resid_studentized_internal, inf.resid_studentized_external

print(f"sum r^2 = {(r ** 2).sum():.6f}  (= N = {N})")
print(f"cooks_d 의 평균 = {sf['cooks_d'].mean():.10f},  1/(N-k) = {1 / (N - k):.10f}")

print(f"\nDFFITS = t * sqrt(h/(1-h)) = t / sqrt(n-1) = t / {np.sqrt(n - 1):.6f}")
print(f"  최대 t = {t.max():.6f}  ->  {t.max() / np.sqrt(n - 1):.6f}")
print(f"  summary_frame 의 dffits 최대 = {sf['dffits'].max():.6f}")

thr_d = 2 * np.sqrt(k / N)
print(f"\nDFFITS 문턱 2*sqrt(p/N) = {thr_d:.6f}")
print(f"  |t| 로 옮기면 {thr_d * np.sqrt(n - 1):.6f} = 2*sqrt(1 - p/N)")
print(f"  걸린 관측 = {np.where(np.abs(sf['dffits'].values) > thr_d)[0]}")
print(f"  Cook 4/N 으로 걸린 관측 = {np.where(sf['cooks_d'].values > 4 / N)[0]}")

print(f"\nstudent_resid: 평균 {t.mean():.6f} (0 이 아니다), "
      f"최대 {t.max():.4f}, 최소 {t.min():.4f}")
print(f"  가장 큰 것을 뺀 |t| 의 최대 = {np.sort(np.abs(t))[-2]:.4f}")
print(f"cooks_d: 3사분위 {sf['cooks_d'].quantile(0.75):.6f}, "
      f"최대 {sf['cooks_d'].max():.6f}  ({sf['cooks_d'].max() / sf['cooks_d'].quantile(0.75):.1f} 배)")

sum r^2 = 60.000000  (= N = 60)
cooks_d 의 평균 = 0.0175438596,  1/(N-k) = 0.0175438596

DFFITS = t * sqrt(h/(1-h)) = t / sqrt(n-1) = t / 4.358899
  최대 t = 7.570250  ->  1.736734
  summary_frame 의 dffits 최대 = 1.736734

DFFITS 문턱 2*sqrt(p/N) = 0.447214
  |t| 로 옮기면 1.949359 = 2*sqrt(1 - p/N)
  걸린 관측 = [52 59]
  Cook 4/N 으로 걸린 관측 = [52 59]

student_resid: 평균 0.036170 (0 이 아니다), 최대 7.5702, 최소 -2.1005
  가장 큰 것을 뺀 |t| 의 최대 = 2.1005
cooks_d: 3사분위 0.012359, 최대 0.505775  (40.9 배)

import numpy as np
from scipy import stats

print("자료를 모두 0 으로 두고 첫 관측만 δ 로 바꾼다")
print(f"{'k':>3s} {'n':>4s} {'δ=10':>8s} {'δ=100':>8s} {'δ=10^4':>9s} "
      f"{'극한 F':>9s} {'임계값':>8s} {'유의?':>6s}")
for k, n in [(3, 5), (3, 10), (3, 20), (5, 10), (2, 10)]:
    row = []
    for d in [10, 100, 1e4, 1e12]:
        gs = [np.zeros(n) for _ in range(k)]
        gs[0] = gs[0].copy()
        gs[0][0] = d
        row.append(stats.f_oneway(*gs).statistic)
    crit = stats.f.ppf(0.95, k - 1, k * (n - 1))
    print(f"{k:3d} {n:4d} {row[0]:8.3f} {row[1]:8.3f} {row[2]:9.3f} "
          f"{row[3]:9.3f} {crit:8.3f} {'예' if row[3] > crit else '아니오':>6s}")

자료를 모두 0 으로 두고 첫 관측만 δ 로 바꾼다
  k    n     δ=10    δ=100    δ=10^4      극한 F      임계값    유의?
  3    5    1.000    1.000     1.000     1.000    3.885    아니오
  3   10    1.000    1.000     1.000     1.000    3.354    아니오
  3   20    1.000    1.000     1.000     1.000    3.159    아니오
  5   10    1.000    1.000     1.000     1.000    2.579    아니오
  2   10    1.000    1.000     1.000     1.000    4.414    아니오

import warnings
warnings.filterwarnings("ignore")

rng = np.random.default_rng(10003)
B = 4_000
print("k=3, n=10, σ=1. 한 관측에 δ 를 더한다")
print(f"{'δ':>6s} {'F 의 중앙값':>11s} {'p<0.05 비율':>11s}")
for d in [0, 1, 2, 3, 4, 5, 7, 10, 20, 50, 200]:
    Fs, hit = [], 0
    for _ in range(B):
        gs = [rng.normal(0, 1, 10) for _ in range(3)]
        gs[0] = gs[0].copy()
        gs[0][0] += d
        r = stats.f_oneway(*gs)
        Fs.append(r.statistic)
        hit += r.pvalue < 0.05
    print(f"{d:6d} {np.median(Fs):11.4f} {hit / B:11.4f}")

k=3, n=10, σ=1. 한 관측에 δ 를 더한다
     δ     F 의 중앙값   p<0.05 비율
     0      0.7391      0.0558
     1      0.7081      0.0517
     2      0.7004      0.0510
     3      0.7237      0.0435
     4      0.7492      0.0420
     5      0.7822      0.0382
     7      0.8213      0.0235
    10      0.8993      0.0095
    20      0.9651      0.0000
    50      0.9947      0.0000
   200      0.9992      0.0000

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(10002)
B = 3_000
N = 24
tb = stats.t.ppf(1 - 0.05 / (2 * N), N - 3 - 1)   # 본페로니 임계값
print(f"집단 크기 (8,8,8), 첫 집단에 크기 δ 의 이상점을 1개 또는 2개 넣는다")
print(f"본페로니 임계값 = {tb:.4f}")
print(f"{'δ':>5s} {'이상점 1개 탐지율':>15s} {'2개 중 하나라도':>15s} {'2개 모두':>10s}")
for D in [3.0, 4.0, 5.0, 6.0]:
    a = b = c = 0
    for _ in range(B):
        base = [rng.normal(0, 1, 8) for _ in range(3)]

        gs = [x.copy() for x in base]
        gs[0][0] += D
        df = pd.DataFrame({"y": np.concatenate(gs),
                           "g": np.repeat([0, 1, 2], 8).astype(str)})
        te = (ols("y ~ C(g)", data=df).fit()
              .get_influence().resid_studentized_external)
        a += abs(te[0]) > tb

        gs = [x.copy() for x in base]
        gs[0][0] += D
        gs[0][1] += D
        df = pd.DataFrame({"y": np.concatenate(gs),
                           "g": np.repeat([0, 1, 2], 8).astype(str)})
        te = (ols("y ~ C(g)", data=df).fit()
              .get_influence().resid_studentized_external)
        b += (abs(te[0]) > tb) or (abs(te[1]) > tb)
        c += (abs(te[0]) > tb) and (abs(te[1]) > tb)
    print(f"{D:5.1f} {a / B:15.4f} {b / B:15.4f} {c / B:10.4f}")

집단 크기 (8,8,8), 첫 집단에 크기 δ 의 이상점을 1개 또는 2개 넣는다
본페로니 임계값 = 3.5342
    δ      이상점 1개 탐지율       2개 중 하나라도      2개 모두
  3.0          0.2780          0.1463     0.0000
  4.0          0.5893          0.2540     0.0000
  5.0          0.8500          0.4123     0.0000
  6.0          0.9690          0.5167     0.0000

import numpy as np
from scipy import stats

def trim_welch(gs, tr=0.2):
    """20% 절사평균에 웰치 구조를 씌운 검정."""
    st = []
    for g in gs:
        x = np.sort(g)
        n = len(x)
        k = int(np.floor(tr * n))
        win = np.clip(x, x[k], x[n - k - 1])
        d = win.var(ddof=1) * (n - 1) / ((n - 2 * k - 1) * (n - 2 * k))
        st.append((x[k:n - k].mean(), d, n - 2 * k))
    m = np.array([s[0] for s in st])
    d = np.array([s[1] for s in st])
    h = np.array([s[2] for s in st], float)
    k = len(gs)
    w = 1 / d
    W = w.sum()
    mt = (w * m).sum() / W
    lam = ((1 - w / W)**2 / (h - 1)).sum()
    F = ((w * (m - mt)**2).sum() / (k - 1)
         / (1 + 2 * (k - 2) / (k**2 - 1) * lam))
    return stats.f.sf(F, k - 1, (3 / (k**2 - 1) * lam)**-1)

rng = np.random.default_rng(10004)
B = 4_000
print("k=3, n=20, 참 평균 (0, 0.6, 1.2), 오차 σ=1 에 오염을 섞는다")
print(f"{'오염':>22s} {'F 검정':>8s} {'크러스컬-월리스':>14s} {'20% 절사 웰치':>13s}")
for lab, eps, sc in [("없음", 0.0, 1), ("5% 가 σ=5 인 정규", 0.05, 5),
                     ("10% 가 σ=5 인 정규", 0.10, 5),
                     ("5% 가 σ=10 인 정규", 0.05, 10)]:
    a = b = c = 0
    for _ in range(B):
        gs = []
        for mu in [0, 0.6, 1.2]:
            x = rng.normal(mu, 1, 20)
            m = rng.random(20) < eps
            x[m] = mu + rng.normal(0, sc, m.sum())
            gs.append(x)
        a += stats.f_oneway(*gs).pvalue < 0.05
        b += stats.kruskal(*gs).pvalue < 0.05
        c += trim_welch(gs) < 0.05
    print(f"{lab:>22s} {a / B:8.4f} {b / B:14.4f} {c / B:13.4f}")

k=3, n=20, 참 평균 (0, 0.6, 1.2), 오차 σ=1 에 오염을 섞는다
                    오염     F 검정       크러스컬-월리스     20% 절사 웰치
                    없음   0.9257         0.9062        0.8502
         5% 가 σ=5 인 정규   0.6410         0.8380        0.7977
        10% 가 σ=5 인 정규   0.4990         0.7775        0.7582
        5% 가 σ=10 인 정규   0.4128         0.8200        0.8005

import numpy as np
from scipy import stats

def p_after_removing(gs, howmany):
    """가장 큰 |표준화 잔차| 를 howmany 개 제거한 뒤 F 검정."""
    y = np.concatenate(gs)
    lab = np.repeat([0, 1, 2], [len(g) for g in gs])
    for _ in range(howmany):
        mm = np.array([y[lab == i].mean() for i in range(3)])
        s = np.array([y[lab == i].std(ddof=1) for i in range(3)])
        z = (y - mm[lab]) / s[lab]
        j = np.argmax(np.abs(z))
        y = np.delete(y, j)
        lab = np.delete(lab, j)
    return stats.f_oneway(*[y[lab == i] for i in range(3)]).pvalue

rng = np.random.default_rng(10005)
B = 6_000
print("모든 집단의 참 평균이 같음 (k=3, n=15), 이상점 없음. 명목 0.05")
print(f"{'절차':>34s} {'실제 오류율':>11s}")
for lab, how in [("제거 없음", 0), ("가장 큰 잔차 1개 제거", 1),
                 ("2개 제거", 2), ("3개 제거", 3)]:
    a = sum(p_after_removing([rng.normal(0, 1, 15) for _ in range(3)], how) < 0.05
            for _ in range(B))
    print(f"{lab:>34s} {a / B:11.4f}")

모든 집단의 참 평균이 같음 (k=3, n=15), 이상점 없음. 명목 0.05
                                절차      실제 오류율
                             제거 없음      0.0495
                     가장 큰 잔차 1개 제거      0.0913
                             2개 제거      0.1152
                             3개 제거      0.1538

import warnings
warnings.filterwarnings("ignore")

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

rng = np.random.default_rng(10006)
rows = [pd.DataFrame({"y": np.round(rng.normal(mu, 2.0, n), 2), "g": f"G{g}"})
        for g, (n, mu) in enumerate([(12, 20.0), (12, 22.5), (12, 21.0)])]
df = pd.concat(rows, ignore_index=True)
df.loc[5, "y"] = 31.0                       # G0 에 이상점 하나

def report(d, tag):
    f = ols("y ~ C(g)", data=d).fit()
    t = sm.stats.anova_lm(f, typ=2)
    mse = t.loc["Residual", "sum_sq"] / t.loc["Residual", "df"]
    m = d.groupby("g").y.mean().round(3)
    return (f"{tag:>22s}  F={t.loc['C(g)', 'F']:7.4f}  "
            f"p={t.loc['C(g)', 'PR(>F)']:.4f}  MSE={mse:7.4f}  평균={list(m)}")

print(report(df, "전체 자료"))

inf = ols("y ~ C(g)", data=df).fit().get_influence()
D = inf.cooks_distance[0]
te = inf.resid_studentized_external
order = np.argsort(-D)[:3]
N = len(df)
tb = stats.t.ppf(1 - 0.05 / (2 * N), N - 3 - 1)
print(f"\n쿡 거리 상위 3: 행 {list(order)},  D = {np.round(D[order], 4).tolist()}")
print(f"  해당 스튜던트화 잔차 = {np.round(te[order], 3).tolist()}")
print(f"  본페로니 임계값 = {tb:.4f}  →  "
      f"이상점 판정: {[bool(abs(te[i]) > tb) for i in order]}")
print()
print(report(df.drop(index=order[0]), "최대 D 제거"))
print(report(df.drop(index=list(order[:2])), "상위 2개 제거"))
w = df.copy()
w.loc[order[0], "y"] = df.groupby("g").y.median()[df.loc[order[0], "g"]]
print(report(w, "중앙값으로 대체"))

                 전체 자료  F= 1.3215  p=0.2805  MSE= 7.5736  평균=[20.543, 22.31, 21.025]

쿡 거리 상위 3: 행 [5, 35, 10],  D = [0.4773, 0.102, 0.0812]
  해당 스튜던트화 잔차 = [5.405, -1.907, -1.682]
  본페로니 임계값 = 3.5010  →  이상점 판정: [True, False, False]

               최대 D 제거  F= 5.1898  p=0.0112  MSE= 4.0827  평균=[19.593, 22.31, 21.025]
              상위 2개 제거  F= 6.4731  p=0.0045  MSE= 3.3917  평균=[19.593, 22.31, 21.465]
              중앙값으로 대체  F= 5.5840  p=0.0082  MSE= 3.9590  평균=[19.597, 22.31, 21.025]

관측 #5 (G0, y = 31.0) 이 영향점으로 확인되었다.
  · 쿡 거리 D = 0.477 (다음으로 큰 값은 0.102)
  · 외적 스튜던트화 잔차 t = 5.41, 본페로니 보정 p < 0.01

이 관측의 포함 여부에 따라 결론이 달라진다.
  · 포함: F(2, 33) = 1.32,  p = 0.281
  · 제외: F(2, 32) = 5.19,  p = 0.011

기록을 확인한 결과 [측정 오류로 판명 / 오류를 발견하지 못함].
따라서 [제외한 결과를 / 두 결과를 모두] 보고한다.

import warnings
warnings.filterwarnings("ignore")

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

rng = np.random.default_rng(10006)
rows = [pd.DataFrame({"y": np.round(rng.normal(mu, 2.0, n), 2), "g": f"G{g}"})
        for g, (n, mu) in enumerate([(12, 20.0), (12, 22.5), (12, 21.0)])]
df = pd.concat(rows, ignore_index=True)
df.loc[5, "y"] = 31.0

fit = ols("y ~ C(g)", data=df).fit()
inf = fit.get_influence()
db = inf.dfbetas
names = list(fit.params.index)
order = np.argsort(-inf.cooks_distance[0])[:3]
N = len(df)

print(f"{'행':>4s} " + " ".join(f"{n:>16s}" for n in names))
for i in order:
    print(f"{i:4d} " + " ".join(f"{db[i, j]:16.4f}" for j in range(len(names))))
print(f"\n문턱 2/√N = {2 / np.sqrt(N):.4f}")

   행        Intercept       C(g)[T.G1]       C(g)[T.G2]
   5           1.6297          -1.1524          -1.1524
  35           0.0000          -0.0000          -0.4066
  10          -0.5071           0.3586           0.3586

문턱 2/√N = 0.3333

1. 집단별 상자그림 — 눈으로 먼저
     ↓
2. 쿡 거리를 크기순 정렬 — 뚜렷하게 튀는 점이 있는가
     ↓
3. 스튜던트화(삭제) 잔차 — 본페로니 임계값과 비교
     ↓
4. DFBETAS — 어느 대비가 흔들리는가 (필요하면)
     ↓
5. 민감도 분석 — 빼고/대체하고 다시 적합
     ↓
6. 기록 확인 — 오류인가, 다른 모집단인가, 그냥 극단값인가
     ↓
7. 보고 — 결론이 바뀌면 모든 결과를, 아니면 그 사실을