가정 위반의 처리¶

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

# 이 페이지의 자료를 교체 전 상태로 되돌려 집단 C만 꺼낸다.
rng = np.random.default_rng(42)
n = 20
raw = np.concatenate([
    rng.normal(10.0, 1.0, n),
    rng.normal(10.8, 1.3, n),
    rng.normal(12.0, 1.6, n),
])
yC = raw[2 * n:]

ybar, S = yC.mean(), ((yC - yC.mean()) ** 2).sum()
y_last, d = yC[-1], 20.0 - yC[-1]
print(f"교체 전: ybar = {ybar:.6f}, S = {S:.6f}, s = {np.sqrt(S / (n - 1)):.6f}")
print(f"        y_n = {y_last:.6f},  d = 20 - y_n = {d:.6f}")

# 공식이 예측하는 값
mean_f = ybar + d / n
S_f = S + 2 * d * (y_last - ybar) + (n - 1) / n * d ** 2
print(f"\n  2d(y_n - ybar) = {2 * d * (y_last - ybar):.6f}")
print(f"  (n-1)d^2/n     = {(n - 1) / n * d ** 2:.6f}")
print(f"공식: mean = {mean_f:.6f},  S = {S_f:.6f},  s = {np.sqrt(S_f / (n - 1)):.6f}")

# 실제로 교체해 본 값
yCn = yC.copy()
yCn[-1] = 20.0
print(f"실제: mean = {yCn.mean():.6f},  S = {((yCn - yCn.mean()) ** 2).sum():.6f},  s = {yCn.std(ddof=1):.6f}")

s0, s1 = np.sqrt(S / (n - 1)), yCn.std(ddof=1)
print(f"\ns 가 {s0:.4f} -> {s1:.4f} ({s1 / s0:.3f} 배)")

교체 전: ybar = 12.190922, S = 24.910300, s = 1.145019
        y_n = 13.549245,  d = 20 - y_n = 6.450755

  2d(y_n - ybar) = 17.524419
  (n-1)d^2/n     = 39.531624
공식: mean = 12.513460,  S = 81.966342,  s = 2.077021
실제: mean = 12.513460,  S = 81.966342,  s = 2.077021

s 가 1.1450 -> 2.0770 (1.814 배)

from scipy.stats import kruskal

# 값 대신 순위로 계산하므로 정규성을 요구하지 않는다. 그 대신 정규자료에서는
# 분산분석보다 검정력이 조금 낮고, 견주는 대상이 평균이 아니라 분포의 위치다.
stat, p_value = kruskal(group1, group2, group3)
print(f"Kruskal-Wallis: H = {stat:.4f}, p-value = {p_value:.4f}")

Kruskal-Wallis: H = 26.0698, p-value = 0.0000

import numpy as np
from scipy import stats

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
group1, group2, group3 = response[:n], response[n:2 * n], response[2 * n:]
N = 3 * n

# 전체 순위를 매기고 집단별 평균순위를 본다.
R = stats.rankdata(response)
Rbar = [R[:n].mean(), R[n:2 * n].mean(), R[2 * n:].mean()]
print("동점 없음:", len(np.unique(response)) == N)
print("평균순위:", [f"{r:.3f}" for r in Rbar], " 전체평균 (N+1)/2 =", (N + 1) / 2)

H_formula = 12 / (N * (N + 1)) * sum(n * r ** 2 for r in Rbar) - 3 * (N + 1)
H_scipy, p = stats.kruskal(group1, group2, group3)
print(f"\n공식 H = {H_formula:.6f}")
print(f"scipy H = {H_scipy:.6f},  p = {p:.4e}")
print(f"chi2(2) 꼬리확률 = {stats.chi2.sf(H_formula, 2):.4e}")

# 완전분리일 때의 최대값
Rmax = [10.5, 30.5, 50.5]
print(f"완전분리 H = {12 / (N * (N + 1)) * sum(n * r ** 2 for r in Rmax) - 3 * (N + 1):.4f}")

# 이상점의 크기를 열 배로 키워 본다.
big = response.copy()
big[-1] = 200.0
print(f"\n{'이상점':>8}{'H':>12}{'KW p':>12}{'F':>10}{'ANOVA p':>12}")
for lab, arr in [("20.0", response), ("200.0", big)]:
    gs = (arr[:n], arr[n:2 * n], arr[2 * n:])
    H, ph = stats.kruskal(*gs)
    F, pf = stats.f_oneway(*gs)
    print(f"{lab:>8}{H:>12.6f}{ph:>12.2e}{F:>10.4f}{pf:>12.2e}")

동점 없음: True
평균순위: ['16.900', '29.550', '45.050']  전체평균 (N+1)/2 = 30.5

공식 H = 26.069836
scipy H = 26.069836,  p = 2.1828e-06
chi2(2) 꼬리확률 = 2.1828e-06
완전분리 H = 52.4590

     이상점           H        KW p         F     ANOVA p
    20.0   26.069836    2.18e-06   16.1314    2.81e-06
   200.0   26.069836    2.18e-06    1.3915    2.57e-01

import numpy as np

# 로그 변환은 오른쪽으로 늘어진 자료를 펴 준다. 퍼짐이 평균에 비례해 커지는
# 자료라면 변환 뒤 등분산과 정규성이 함께 좋아지는 일이 흔하다.
# 다만 결과의 해석이 원래 단위가 아니게 된다는 대가가 따른다.
data['log_response'] = np.log(data['response'])
print(data.groupby('group').log_response.agg(['mean', 'std']).round(4))

         mean     std
group                
A      2.2955  0.0900
B      2.3886  0.0918
C      2.5159  0.1460

import numpy as np
import pandas as pd

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})
data["log_response"] = np.log(data["response"])

print(f"{'집단':>4}{'mean':>9}{'sd':>8}{'CV':>9}{'sd(logY)':>10}{'오차':>8}")
for g, sub in data.groupby("group"):
    m, s = sub.response.mean(), sub.response.std(ddof=1)
    cv, sl = s / m, sub.log_response.std(ddof=1)
    print(f"{g:>4}{m:>9.4f}{s:>8.4f}{cv:>9.4f}{sl:>10.4f}{sl / cv - 1:>7.1%}")

print(f"\n{'집단':>4}{'산술평균':>10}{'기하평균':>10}{'GM/AM':>9}{'1-CV^2/2':>10}")
for g, sub in data.groupby("group"):
    am, gm = sub.response.mean(), np.exp(sub.log_response.mean())
    cv = sub.response.std(ddof=1) / am
    print(f"{g:>4}{am:>10.4f}{gm:>10.4f}{gm / am:>9.4f}{1 - cv ** 2 / 2:>10.4f}")

# 분산비가 얼마나 좁혀졌는가
for lab, col in [("원척도", "response"), ("로그척도", "log_response")]:
    sd = data.groupby("group")[col].std(ddof=1)
    print(f"\n{lab}: sd = {sd.round(4).tolist()},  최대/최소 = {sd.max() / sd.min():.3f}")

  집단     mean      sd       CV  sd(logY)      오차
   A   9.9671  0.8702   0.0873    0.0900   3.1%
   B  10.9423  1.0335   0.0945    0.0918  -2.8%
   C  12.5135  2.0770   0.1660    0.1460 -12.1%

  집단      산술평균      기하평균    GM/AM  1-CV^2/2
   A    9.9671    9.9296   0.9962    0.9962
   B   10.9423   10.8978   0.9959    0.9955
   C   12.5135   12.3771   0.9891    0.9862

원척도: sd = [0.8702, 1.0335, 2.077],  최대/최소 = 2.387

로그척도: sd = [0.09, 0.0918, 0.146],  최대/최소 = 1.622

# 제곱근 변환은 로그보다 약하게 편다. 도수 자료처럼 분산이 평균에 비례하는
# 경우에 알맞다.
data['sqrt_response'] = np.sqrt(data['response'])
print(data.groupby('group').sqrt_response.agg(['mean', 'std']).round(4))

         mean     std
group                
A      3.1541  0.1398
B      3.3045  0.1538
C      3.5274  0.2735

import numpy as np
import pandas as pd

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})
data["sqrt_response"] = np.sqrt(data["response"])

print(f"{'집단':>4}{'mean':>9}{'sd':>8}{'sd/(2*sqrt(mean))':>20}{'sd(sqrtY)':>11}{'오차':>8}")
for g, sub in data.groupby("group"):
    m, s = sub.response.mean(), sub.response.std(ddof=1)
    pred, act = s / (2 * np.sqrt(m)), sub.sqrt_response.std(ddof=1)
    print(f"{g:>4}{m:>9.4f}{s:>8.4f}{pred:>20.4f}{act:>11.4f}{act / pred - 1:>7.1%}")

# 세 척도에서 흩어짐의 불균형 정도
data["log_response"] = np.log(data["response"])
print()
for lab, col in [("원척도", "response"), ("제곱근", "sqrt_response"), ("로그", "log_response")]:
    sd = data.groupby("group")[col].std(ddof=1)
    print(f"{lab:>6}: sd 최대/최소 = {sd.max() / sd.min():.3f}")

# lambda = 1 - alpha 규칙: 포아송(alpha=1/2)에서는 제곱근이 맞는다
print(f"\n포아송 자료, 평균을 바꿔 가며 sd(Y^lambda) 를 재 본다 (목표: lambda 열이 평평)")
rg = np.random.default_rng(7)
print(f"{'mu':>6}{'sd(Y)':>9}{'sd(sqrtY)':>11}{'sd(logY)':>10}")
for mu in [5, 20, 80, 320]:
    y = rg.poisson(mu, 200_000).astype(float)
    y[y == 0] = 0.5                                  # 로그를 위해 0만 보정
    print(f"{mu:>6}{y.std(ddof=1):>9.3f}{np.sqrt(y).std(ddof=1):>11.4f}{np.log(y).std(ddof=1):>10.4f}")
print("이론: sd(sqrt Y) -> 1/2 = 0.5000,  sd(log Y) ~ 1/sqrt(mu)")

  집단     mean      sd   sd/(2*sqrt(mean))  sd(sqrtY)      오차
   A   9.9671  0.8702              0.1378     0.1398   1.4%
   B  10.9423  1.0335              0.1562     0.1538  -1.5%
   C  12.5135  2.0770              0.2936     0.2735  -6.8%

   원척도: sd 최대/최소 = 2.387
   제곱근: sd 최대/최소 = 1.956
    로그: sd 최대/최소 = 1.622

포아송 자료, 평균을 바꿔 가며 sd(Y^lambda) 를 재 본다 (목표: lambda 열이 평평)
    mu    sd(Y)  sd(sqrtY)  sd(logY)
     5    2.227     0.5176    0.5340
    20    4.466     0.5047    0.2330
    80    8.957     0.5018    0.1130
   320   17.907     0.5009    0.0561
이론: sd(sqrt Y) -> 1/2 = 0.5000,  sd(log Y) ~ 1/sqrt(mu)

from scipy.stats import boxcox

# boxcox는 양수 자료만 받는다. 0이나 음수가 있으면 상수를 더해야 한다.
transformed_data, best_lambda = boxcox(data['response'])
print(f"Optimal lambda = {best_lambda:.4f}")

Optimal lambda = -1.5414

import numpy as np
from scipy import stats, optimize

rng = np.random.default_rng(42)
n = 20
clean = np.concatenate([
    rng.normal(10.0, 1.0, n), rng.normal(10.8, 1.3, n), rng.normal(12.0, 1.6, n)])
response = clean.copy()
response[-1] = 20.0                       # 이상점
grp = np.repeat([0, 1, 2], n)

def bc(y, lam):
    return np.log(y) if lam == 0 else (y ** lam - 1) / lam

def llf(lam, y, by_group=False):
    """프로파일 로그가능도. by_group=True 면 집단마다 평균을 따로 둔다."""
    z = bc(y, lam)
    mu = np.array([z[grp == g].mean() for g in grp]) if by_group else z.mean()
    r = z - mu
    return -len(y) / 2 * np.log((r ** 2).mean()) + (lam - 1) * np.log(y).sum()

print(f"sum log y = {np.log(response).sum():.4f}")
print(f"\n{'lambda':>9}{'내 식':>12}{'scipy':>12}{'야코비 없이':>14}")
for lam in [-6.0, -3.0, -2.0, -1.5414, -1.0, 0.0, 1.0]:
    z = bc(response, lam)
    print(f"{lam:>9.4f}{llf(lam, response):>12.4f}{stats.boxcox_llf(lam, response):>12.4f}"
          f"{-len(response) / 2 * np.log(z.var()):>14.3f}")

_, lam_scipy = stats.boxcox(response)
mine = optimize.minimize_scalar(lambda L: -llf(L, response), bracket=(-5, 0, 5)).x
lam_grp = optimize.minimize_scalar(
    lambda L: -llf(L, response, by_group=True), bracket=(-5, 0, 5)).x
_, lam_clean = stats.boxcox(clean)
print(f"\nscipy boxcox lambda      = {lam_scipy:.6f}")
print(f"내 식의 argmax           = {mine:.6f}")
print(f"집단평균을 따로 둔 lambda = {lam_grp:.6f}")
print(f"이상점을 빼면 lambda     = {lam_clean:.6f}   (0 = 로그)")

print(f"\n{'척도':>20}{'왜도':>8}{'Shapiro p':>11}{'sd 최대/최소':>12}{'이상점 z':>10}")
for lab, lam in [("원자료 lam=1", 1.0), ("로그 lam=0", 0.0),
                 ("scipy lam=-1.5414", lam_scipy), ("집단별 lam=-1.5262", lam_grp)]:
    z = bc(response, lam)
    sds = [z[grp == g].std(ddof=1) for g in range(3)]
    r = z - np.array([z[grp == g].mean() for g in grp])
    zout = (z[-1] - z[grp == 2].mean()) / sds[2]
    print(f"{lab:>20}{stats.skew(r):>8.3f}{stats.shapiro(r).pvalue:>11.4f}"
          f"{max(sds) / min(sds):>12.3f}{zout:>10.2f}")

sum log y = 143.9986

   lambda         내 식       scipy        야코비 없이
  -6.0000    -41.9710    -41.9710       966.019
  -3.0000    -25.8487    -25.8487       550.146
  -2.0000    -23.7607    -23.7607       408.235
  -1.5414    -23.5122    -23.5122       342.446
  -1.0000    -23.8881    -23.8881       264.109
   0.0000    -26.8420    -26.8420       117.157
   1.0000    -33.4008    -33.4008       -33.401

scipy boxcox lambda      = -1.541353
내 식의 argmax           = -1.541353
집단평균을 따로 둔 lambda = -1.526238
이상점을 빼면 lambda     = -0.024428   (0 = 로그)

                  척도      왜도  Shapiro p    sd 최대/최소     이상점 z
           원자료 lam=1   2.435     0.0000       2.387      3.60
            로그 lam=0   1.154     0.0007       1.622      3.29
   scipy lam=-1.5414  -0.141     0.3817       1.237      2.71
     집단별 lam=-1.5262  -0.132     0.3791       1.235      2.71

from scipy.stats import f_oneway
# Welch 분산분석은 등분산을 가정하지 않는다. 자료를 바꾸지 않고 검정만
# 바꾸는 방법이라 결과를 원래 단위로 읽을 수 있다.
import pingouin as pg

welch_result = pg.welch_anova(dv='response', between='group', data=data)
print(welch_result)

  Source  ddof1      ddof2          F     p_unc       np2
0  group      2  35.386613  14.579033  0.000024  0.361436

import numpy as np
import pandas as pd
from scipy import stats

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
gs = [response[:n], response[n:2 * n], response[2 * n:]]

ni = np.array([len(g) for g in gs], float)
mi = np.array([g.mean() for g in gs])
vi = np.array([g.var(ddof=1) for g in gs])
k = len(gs)

w = ni / vi                                  # 역분산 가중치
W = w.sum()
mt = (w * mi).sum() / W                      # 가중평균
A = (w * (mi - mt) ** 2).sum() / (k - 1)
S = (((1 - w / W) ** 2) / (ni - 1)).sum()
B = 1 + 2 * (k - 2) / (k ** 2 - 1) * S
F_w = A / B
df2 = (k ** 2 - 1) / (3 * S)

print(f"{'집단':>4}{'n':>4}{'mean':>9}{'s^2':>9}{'w_i':>10}{'w_i/W':>8}{'기여':>10}")
for i in range(k):
    print(f"{'ABC'[i]:>4}{ni[i]:>4.0f}{mi[i]:>9.4f}{vi[i]:>9.4f}{w[i]:>10.4f}"
          f"{w[i] / W:>8.4f}{(1 - w[i] / W) ** 2 / (ni[i] - 1):>10.6f}")
print(f"\nW = {W:.6f},  가중평균 = {mt:.6f},  단순평균 = {mi.mean():.6f}")
print(f"A = {A:.6f},  S = {S:.6f},  B = {B:.6f}")
print(f"\nWelch  F = {F_w:.6f},  df1 = {k - 1},  df2 = {df2:.6f},  p = {stats.f.sf(F_w, k - 1, df2):.6e}")
F_c, p_c = stats.f_oneway(*gs)
print(f"고전   F = {F_c:.6f},  df1 = {k - 1},  df2 = {int(ni.sum() - k)},  p = {p_c:.6e}")

# 균형설계에서 df2 의 상한
print(f"\n균형설계 k=3, n=20 에서 df2 의 상한 = k(k+1)(n-1)/(3(k-1)) "
      f"= {k * (k + 1) * (n - 1) / (3 * (k - 1)):.4f}")

  집단   n     mean      s^2       w_i   w_i/W        기여
   A  20   9.9671   0.7572   26.4128  0.5307  0.011593
   B  20  10.9423   1.0682   18.7231  0.3762  0.020482
   C  20  12.5135   4.3140    4.6360  0.0931  0.043283

W = 49.771906,  가중평균 = 10.571115,  단순평균 = 11.140943
A = 14.853695,  S = 0.075358,  B = 1.018840

Welch  F = 14.579033,  df1 = 2,  df2 = 35.386613,  p = 2.408061e-05
고전   F = 16.131362,  df1 = 2,  df2 = 57,  p = 2.807595e-06

균형설계 k=3, n=20 에서 df2 의 상한 = k(k+1)(n-1)/(3(k-1)) = 38.0000

import statsmodels.api as sm

rlm_model = sm.RLM.from_formula('response ~ group', data=data, M=sm.robust.norms.HuberT())
result = rlm_model.fit()
# summary()는 실행 날짜와 시각을 함께 찍으므로 계수 표만 뽑아 본다.
print(result.params.round(4))
print()
print(result.bse.round(4))

Intercept     9.9880
group[T.B]    0.8953
group[T.C]    2.2366
dtype: float64

Intercept     0.2267
group[T.B]    0.3206
group[T.C]    0.3206
dtype: float64

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

rng = np.random.default_rng(42)
n = 20
clean = np.concatenate([
    rng.normal(10.0, 1.0, n), rng.normal(10.8, 1.3, n), rng.normal(12.0, 1.6, n)])
response = clean.copy()
response[-1] = 20.0
data = pd.DataFrame({"group": np.repeat(["A", "B", "C"], n), "response": response})

rlm = sm.RLM.from_formula("response ~ group", data=data,
                          M=sm.robust.norms.HuberT()).fit()
mu = [rlm.params.iloc[0], rlm.params.iloc[0] + rlm.params.iloc[1],
      rlm.params.iloc[0] + rlm.params.iloc[2]]
s, c = rlm.scale, 1.345
print(f"Huber 척도 s = {s:.6f}   (MAD 기준),  c = {c}")

print(f"\n{'집단':>4}{'OLS 평균':>11}{'Huber':>11}{'가중평균':>11}{'최소 w':>9}{'w<1 개수':>10}{'이상점 뺀 평균':>14}")
for i, g in enumerate("ABC"):
    y = response[i * n:(i + 1) * n]
    w = np.minimum(1.0, c / np.abs((y - mu[i]) / s))
    wm = (w * y).sum() / w.sum()
    ref = clean[i * n:(i + 1) * n].mean()
    print(f"{g:>4}{y.mean():>11.6f}{mu[i]:>11.6f}{wm:>11.6f}"
          f"{w.min():>9.4f}{(w < 1).sum():>10d}{ref:>14.6f}")

u_out = (20.0 - mu[2]) / s
print(f"\n이상점: u = (20 - {mu[2]:.4f})/{s:.4f} = {u_out:.4f},  w = c/|u| = {c / abs(u_out):.4f}")

# 표준오차: 샌드위치 식과 statsmodels 의 보정계수
ols_fit = ols("response ~ C(group)", data=data).fit()
u = rlm.resid.values / s
psi, psip = np.clip(u, -c, c), (np.abs(u) <= c).astype(float)
N, p = len(u), 3
sig_h1 = s * np.sqrt((psi ** 2).sum() / (N - p)) / psip.mean()
k = 1 + p / N * psip.var() / psip.mean() ** 2
print(f"\nOLS   잔차 표준편차 = {np.sqrt(ols_fit.mse_resid):.6f},  차의 s.e. = {ols_fit.bse.iloc[1]:.6f}")
print(f"Huber 샌드위치 sigma = {sig_h1:.6f},  보정계수 k = {k:.6f}")
print(f"      k*sigma*sqrt(2/n) = {k * sig_h1 * np.sqrt(2 / n):.6f}"
      f"   statsmodels = {rlm.bse.iloc[1]:.6f}")

Huber 척도 s = 1.130443   (MAD 기준),  c = 1.345

  집단     OLS 평균      Huber       가중평균     최소 w    w<1 개수      이상점 뺀 평균
   A   9.967068   9.987997   9.987997   0.7841         1      9.967068
   B  10.942301  10.883280  10.883280   0.5629         1     10.942301
   C  12.513460  12.224630  12.224630   0.1955         4     12.190922

이상점: u = (20 - 12.2246)/1.1304 = 6.8782,  w = c/|u| = 0.1955

OLS   잔차 표준편차 = 1.430551,  차의 s.e. = 0.452380
Huber 샌드위치 sigma = 1.008230,  보정계수 k = 1.005556
      k*sigma*sqrt(2/n) = 0.320602   statsmodels = 0.320602

import numpy as np
from scipy.stats import f_oneway

# 자료에서 실제로 관측된 F 값
observed_f, _ = f_oneway(group1, group2, group3)

# 귀무가설이 참이라면 집단 이름표는 아무 뜻이 없다. 그래서 이름표를 뒤섞어
# 가며 F 를 다시 구하면, 어떤 분포도 가정하지 않고 귀무분포를 얻을 수 있다.
all_data = np.concatenate([group1, group2, group3])
group_sizes = [len(group1), len(group2), len(group3)]
n_permutations = 10000
perm_f_stats = []

rng = np.random.default_rng(42)
for _ in range(n_permutations):
    shuffled = rng.permutation(all_data)
    g1 = shuffled[:group_sizes[0]]
    g2 = shuffled[group_sizes[0]:group_sizes[0]+group_sizes[1]]
    g3 = shuffled[group_sizes[0]+group_sizes[1]:]
    f_stat, _ = f_oneway(g1, g2, g3)
    perm_f_stats.append(f_stat)

# 뒤섞어 만든 F 중 관측값 이상인 비율이 곧 p-값이다.
p_value = np.mean(np.array(perm_f_stats) >= observed_f)
print(f"Permutation test p-value: {p_value:.4f}")

Permutation test p-value: 0.0000

import numpy as np
from math import comb
from scipy import stats

rng0 = np.random.default_rng(42)
n = 20
clean = np.concatenate([
    rng0.normal(10.0, 1.0, n), rng0.normal(10.8, 1.3, n), rng0.normal(12.0, 1.6, n)])
response = clean.copy()
response[-1] = 20.0

observed_f = stats.f_oneway(response[:n], response[n:2 * n], response[2 * n:])[0]
B = 10_000
print(f"관측 F = {observed_f:.6f}")
print(f"서로 다른 집단배정의 수 60!/(20!)^3 = {comb(60, 20) * comb(40, 20):.3e}")

def perm_F(arr, B, seed):
    rng = np.random.default_rng(seed)
    out = np.empty(B)
    for b in range(B):
        sh = rng.permutation(arr)
        out[b] = stats.f_oneway(sh[:n], sh[n:2 * n], sh[2 * n:])[0]
    return out

fs = perm_F(response, B, 42)
cnt = int((fs >= observed_f).sum())
print(f"\n초과 횟수 = {cnt},  순열 F 의 최대값 = {fs.max():.4f}")
print(f"소박한 p-값   = {cnt / B:.6f}")
print(f"(c+1)/(B+1)  = {(cnt + 1) / (B + 1):.6f}")
print(f"95% 상한 1-0.05^(1/B) = {1 - 0.05 ** (1 / B):.6f}")
print(f"참고: F(2,57) 꼬리확률 = {stats.f.sf(observed_f, 2, 57):.3e}")

# 순열 귀무분포가 F(2,57) 인가
d = stats.f(2, 57)
print(f"\n{'귀무분포':>18}{'평균':>9}{'분산':>9}{'95% 분위':>10}{'최대':>9}")
for lab, arr in [("순열 (이상점 있음)", response), ("순열 (이상점 없음)", clean)]:
    g = perm_F(arr, 20_000, 42)
    print(f"{lab:>18}{g.mean():>9.4f}{g.var(ddof=1):>9.4f}{np.quantile(g, 0.95):>10.4f}{g.max():>9.3f}")
print(f"{'F(2,57)':>18}{d.mean():>9.4f}{d.var():>9.4f}{d.ppf(0.95):>10.4f}{'inf':>9}")

관측 F = 16.131362
서로 다른 집단배정의 수 60!/(20!)^3 = 5.778e+26

초과 횟수 = 0,  순열 F 의 최대값 = 10.8218
소박한 p-값   = 0.000000
(c+1)/(B+1)  = 0.000100
95% 상한 1-0.05^(1/B) = 0.000300
참고: F(2,57) 꼬리확률 = 2.808e-06

              귀무분포       평균       분산    95% 분위       최대
       순열 (이상점 있음)   1.0240   0.9230    2.9202   10.822
       순열 (이상점 없음)   1.0256   1.1375    3.1394   13.571
           F(2,57)   1.0364   1.1551    3.1588      inf

import warnings
warnings.filterwarnings("ignore")

import numpy as np
from scipy import stats

def welch_p(gs):
    n = np.array([len(g) for g in gs], float)
    m = np.array([g.mean() for g in gs])
    v = np.array([g.var(ddof=1) for g in gs])
    k = len(n)
    w = n / v
    W = w.sum()
    mt = (w * m).sum() / W
    lam = ((1 - w / W)**2 / (n - 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)

def perm_p(gs, rng, R=300):
    obs = stats.f_oneway(*gs).statistic
    z = np.concatenate(gs)
    ns = [len(g) for g in gs]
    c = 0
    for _ in range(R):
        p = rng.permutation(z)
        i, parts = 0, []
        for n in ns:
            parts.append(p[i:i + n])
            i += n
        c += stats.f_oneway(*parts).statistic >= obs
    return (c + 1) / (R + 1)

def boot_p(gs, rng, R=300):
    """집단별 중심을 맞춘 뒤 각 집단 안에서 재표집."""
    obs = stats.f_oneway(*gs).statistic
    cen = [g - g.mean() for g in gs]
    c = 0
    for _ in range(R):
        bs = [x[rng.integers(0, len(x), len(x))] for x in cen]
        c += stats.f_oneway(*bs).statistic >= obs
    return (c + 1) / (R + 1)

rng = np.random.default_rng(11002)
B = 1_200
print("정규 자료, σ=(3,1,1,1) — 한 집단만 분산이 큼. 명목 0.05")
print(f"{'설계':>26s} {'표준 F':>8s} {'웰치':>8s} {'순열 F':>8s} {'부트스트랩':>10s}")
for lab, ns, sds, mus in [
    ("귀무 · 역페어링 n=(8,20,20,20)", [8, 20, 20, 20], [3, 1, 1, 1], [0, 0, 0, 0]),
    ("귀무 · 정페어링 n=(20,20,20,8)", [20, 20, 20, 8], [3, 1, 1, 1], [0, 0, 0, 0]),
    ("대립 μ=(0,0,0,1) 역페어링", [8, 20, 20, 20], [3, 1, 1, 1], [0, 0, 0, 1]),
]:
    a = b = c = d = 0
    for _ in range(B):
        gs = [rng.normal(m, s, n) for n, s, m in zip(ns, sds, mus)]
        a += stats.f_oneway(*gs).pvalue < 0.05
        b += welch_p(gs) < 0.05
        c += perm_p(gs, rng) < 0.05
        d += boot_p(gs, rng) < 0.05
    print(f"{lab:>26s} {a / B:8.4f} {b / B:8.4f} {c / B:8.4f} {d / B:10.4f}")

정규 자료, σ=(3,1,1,1) — 한 집단만 분산이 큼. 명목 0.05
                        설계     표준 F       웰치     순열 F      부트스트랩
  귀무 · 역페어링 n=(8,20,20,20)   0.2092   0.0508   0.2133     0.0692
  귀무 · 정페어링 n=(20,20,20,8)   0.0683   0.0592   0.0692     0.0533
       대립 μ=(0,0,0,1) 역페어링   0.6825   0.8075   0.6767     0.2600

1. 집단별 n 과 s 표를 확인한다 (어느 집단에 큰 분산이 붙었는가)
     ↓
2. 웰치 분산분석을 수행한다  (레빈의 p=0.003 은 참고로만)
     ↓
3. 유의하면 게임스-하웰 사후검정
     ↓
4. 효과크기와 신뢰구간을 보고
     ↓
5. 표준 F 의 p=0.008 은 보고하지 않거나 "가정 위반으로 무효"라 밝힌다

import warnings
warnings.filterwarnings("ignore")

import numpy as np
from scipy import stats

def welch_p(gs):
    n = np.array([len(g) for g in gs], float)
    m = np.array([g.mean() for g in gs])
    v = np.array([g.var(ddof=1) for g in gs])
    k = len(n)
    w = n / v
    W = w.sum()
    mt = (w * m).sum() / W
    lam = ((1 - w / W)**2 / (n - 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(11004)
B = 6_000
NS = [12, 15, 10]
print("자료: y = exp(Z),  Z ~ N(μ_g, σ_g²)  — 로그정규, 양수")
print(f"{'상황 (로그 척도의 μ, σ)':>32s} {'표준 F':>8s} {'웰치':>8s} {'KW':>8s} "
      f"{'log 후 F':>9s} {'log 후 웰치':>11s}")
for lab, mus, sds in [("귀무 μ=(0,0,0) σ=(1,1,1)", [0, 0, 0], [1, 1, 1]),
                      ("귀무 μ=(0,0,0) σ=(1.6,1,1)", [0, 0, 0], [1.6, 1, 1]),
                      ("대립 μ=(0,0,1) σ=(1,1,1)", [0, 0, 1], [1, 1, 1])]:
    a = b = c = d = e = 0
    for _ in range(B):
        gs = [np.exp(rng.normal(m, s, n)) for n, m, s in zip(NS, mus, sds)]
        lg = [np.log(g) for g in gs]
        a += stats.f_oneway(*gs).pvalue < 0.05
        b += welch_p(gs) < 0.05
        c += stats.kruskal(*gs).pvalue < 0.05
        d += stats.f_oneway(*lg).pvalue < 0.05
        e += welch_p(lg) < 0.05
    print(f"{lab:>32s} {a / B:8.4f} {b / B:8.4f} {c / B:8.4f} "
          f"{d / B:9.4f} {e / B:11.4f}")

자료: y = exp(Z),  Z ~ N(μ_g, σ_g²)  — 로그정규, 양수
                상황 (로그 척도의 μ, σ)     표준 F       웰치       KW   log 후 F    log 후 웰치
          귀무 μ=(0,0,0) σ=(1,1,1)   0.0377   0.0472   0.0488    0.0490      0.0520
        귀무 μ=(0,0,0) σ=(1.6,1,1)   0.0662   0.0472   0.0542    0.0585      0.0545
          대립 μ=(0,0,1) σ=(1,1,1)   0.5332   0.2728   0.5948    0.6325      0.6097

import warnings
warnings.filterwarnings("ignore")

import numpy as np
from scipy import stats

rng = np.random.default_rng(11005)
B = 8_000
print("자료: y = μ_g + ε,  ε ~ N(0, σ²) — 원 척도에서 이미 정규·등분산")
print("여기에 습관적으로 로그를 씌우면?")
print(f"{'설계':>26s} {'원 척도 F':>10s} {'log 후 F':>9s} {'원 SD비':>8s} {'log SD비':>9s}")
for lab, mus, sd in [("μ=(10,10,10), σ=2", [10, 10, 10], 2.0),
                     ("μ=(3,3,3),  σ=2 (0 에 가까움)", [3, 3, 3], 2.0),
                     ("μ=(3,6,12), σ=2 · 대립", [3, 6, 12], 2.0)]:
    a = b = 0
    r1, r2 = [], []
    for _ in range(B):
        gs = [np.clip(rng.normal(m, sd, 15), 0.05, None) for m in mus]
        lg = [np.log(g) for g in gs]
        a += stats.f_oneway(*gs).pvalue < 0.05
        b += stats.f_oneway(*lg).pvalue < 0.05
        s1 = np.array([g.std(ddof=1) for g in gs])
        s2 = np.array([g.std(ddof=1) for g in lg])
        r1.append(s1.max() / s1.min())
        r2.append(s2.max() / s2.min())
    print(f"{lab:>26s} {a / B:10.4f} {b / B:9.4f} "
          f"{np.median(r1):8.2f} {np.median(r2):9.2f}")

print("\n로그 변환 후 잔차의 치우침 (원 척도에서는 0)")
for lab, mus, sd in [("μ=(10,10,10), σ=2", [10, 10, 10], 2.0),
                     ("μ=(3,3,3), σ=2", [3, 3, 3], 2.0)]:
    sk = []
    for _ in range(2_000):
        gs = [np.clip(rng.normal(m, sd, 15), 0.05, None) for m in mus]
        lg = np.concatenate([np.log(g) - np.log(g).mean() for g in gs])
        sk.append(stats.skew(lg))
    print(f"  {lab:>26s}: 평균 치우침 = {np.mean(sk):+.3f}")

자료: y = μ_g + ε,  ε ~ N(0, σ²) — 원 척도에서 이미 정규·등분산
여기에 습관적으로 로그를 씌우면?
                        설계     원 척도 F   log 후 F    원 SD비   log SD비
         μ=(10,10,10), σ=2     0.0514    0.0516     1.36      1.43
 μ=(3,3,3),  σ=2 (0 에 가까움)     0.0508    0.0430     1.32      1.89
      μ=(3,6,12), σ=2 · 대립     1.0000    1.0000     1.35      7.25

로그 변환 후 잔차의 치우침 (원 척도에서는 0)
           μ=(10,10,10), σ=2: 평균 치우침 = -0.539
              μ=(3,3,3), σ=2: 평균 치우침 = -1.745