X̄₁ - X̄₂의 표본분포 (이분산과 불균형)¶

(1) 해석적으로. \(c = 1/n_1 + 1/n_2\), \(d_i = n_i-1\)로 두자. 분자는 이분산이어도 정확히 정규이므로

이다. 분모의 제곱 \(V_p = S_p^2 c\)는 서로 독립인 두 카이제곱의 척도가 다른 합이다.

이 합은 카이제곱이 아니다. 그래서 Satterthwaite의 방법을 합동분산 쪽에 그대로 적용한다. 평균과 분산을 계산하면

이고, \(V_p \approx \theta\chi^2_m/m\)에 두 적률을 맞추면 \(\theta = E[V_p]\)와

를 얻는다. \(c\)가 약분되어 사라지는 것에 주목하라. \(m\)은 분산비와 자유도만으로 정해진다. 등분산이면 \(m = d_1+d_2\)로 합동 자유도와 같아진다.

이제 \(R = E[V_p]/\tau^2\)이 본문의 배율이므로 \(V_p \approx \tau^2 R\,\chi^2_m/m\)이고, 분자와 분모가 독립이므로

이다. 따라서 명목 5% 검정의 실제 오류율은

가 된다. 고장 난 손잡이가 둘이라는 것이 이 식의 요점이다. \(R \ne 1\)이면 통계량의 척도가 틀리고, \(m < d_1+d_2\)이면 합동 자유도가 주장하는 것보다 분모가 더 흔들려 꼬리가 무거워진다. 본문이 쓴 거친 어림 \(P(|Z| > t_{0.975}\sqrt R)\)은 뒤의 손잡이를 빠뜨린 것이며, 그래서 체계적으로 작게 나온다.

격자에 넣어 두 어림을 나란히 둔다.

거친 어림은 등분산 세 줄에서 이미 \(0.036\)과 \(0.044\)를 주어 \(0.05\)를 맞히지 못한다. 아무것도 고장 나지 않은 자리에서 틀리는 어림이니 고장 난 자리의 숫자도 믿을 수 없다. 다듬은 어림은 그 세 줄에서 \(0.0500\)을 정확히 주고, 균형 설계의 \(R = 1\) 줄에서도 \(0.0554\)와 \(0.0617\)로 \(m\)이 줄어든 몫만큼 부푼다.

(2) 모의실험.

import numpy as np
from scipy import stats

rng = np.random.default_rng(1)
B = 50_000          # 반복 횟수. 오류율 추정의 표준오차는 약 0.001이다.
alpha = 0.05


def error_rates(n1, n2, s1, s2):
    """H0(mu1 = mu2)가 참인 자료를 B번 만들어 두 검정의 기각 비율을 센다."""
    x1 = rng.normal(0, s1, size=(B, n1))
    x2 = rng.normal(0, s2, size=(B, n2))
    d = x1.mean(axis=1) - x2.mean(axis=1)
    v1, v2 = x1.var(axis=1, ddof=1), x2.var(axis=1, ddof=1)

    # 합동 t: 자유도가 n1+n2-2로 고정된다.
    sp2 = ((n1 - 1) * v1 + (n2 - 1) * v2) / (n1 + n2 - 2)
    t_pool = d / np.sqrt(sp2 * (1 / n1 + 1 / n2))
    rej_pool = np.mean(np.abs(t_pool) > stats.t(n1 + n2 - 2).ppf(1 - alpha / 2))

    # Welch t: 표본마다 자유도가 달라진다.
    a, b = v1 / n1, v2 / n2
    t_welch = d / np.sqrt(a + b)
    nu = (a + b) ** 2 / (a ** 2 / (n1 - 1) + b ** 2 / (n2 - 1))
    rej_welch = np.mean(np.abs(t_welch) > stats.t(nu).ppf(1 - alpha / 2))
    return rej_pool, rej_welch, nu.mean()


print("sigma1/sigma2  (n1, n2)   pooled t   Welch t   mean nu")
for r in (1, 2, 4):
    for n1, n2 in ((10, 10), (10, 40), (40, 10)):
        rp, rw, nu = error_rates(n1, n2, r, 1.0)
        print(f"{r:>13}  ({n1:>2}, {n2:>2})     {rp:.3f}      {rw:.3f}     {nu:5.1f}")

출력:

sigma1/sigma2  (n1, n2)   pooled t   Welch t   mean nu
            1  (10, 10)     0.051      0.049      16.5
            1  (10, 40)     0.051      0.051      15.4
            1  (40, 10)     0.049      0.050      15.4
            2  (10, 10)     0.054      0.050      13.6
            2  (10, 40)     0.186      0.052      10.5
            2  (40, 10)     0.005      0.048      31.1
            4  (10, 10)     0.060      0.049      10.4
            4  (10, 40)     0.291      0.050       9.4
            4  (40, 10)     0.001      0.049      46.1

첫 세 줄이 등분산 기준선이다. 두 검정이 모두 0.05를 지킨다. 아래로 내려가면서 합동 \(t\)의 값이 좌우로 벌어지는데, Welch 열은 아홉 칸 모두 0.048에서 0.052 사이에 머문다.

아홉 칸이 모두 맞는다. 반복 5만 회에서 \(p = 0.05\) 둘레의 몬테카를로 오차가 \(0.0010\), \(p = 0.29\) 둘레에서 \(0.0020\)이다.

아홉 줄의 어긋남이 모두 몬테카를로 오차의 두 배 안이다. 표의 모든 수가 두 손잡이 \(R\)과 \(m\)에서 나왔다.

두 손잡이가 각각 무엇을 하는지 줄별로 읽힌다. \((10,10)\) 세 줄은 \(R = 1\)이라 척도는 멀쩡한데 \(m\)이 \(18.0 \to 13.2 \to 10.1\)로 줄어 오류율이 \(0.050 \to 0.055 \to 0.062\)로 민다. 분산비 하나만으로 생기는 어긋남은 이만큼이고, 실무에서 감당할 만하다. \((10,40)\) 줄은 \(R\)이 \(0.46\)과 \(0.29\)로 무너지면서 \(0.183\)과 \(0.294\)가 된다. 같은 분산비인데 표본크기의 배치만 바꿔 어긋남이 다섯 배가 되었다. \((40,10)\) 줄은 \(R > 1\)이라 반대로 \(0.005\)와 \(0.0007\)로 주저앉는다.

Welch 열은 왜 아홉 칸 모두 0.05인가. 두 손잡이를 다 잠갔기 때문이다. 분모를 \(S_1^2/n_1+S_2^2/n_2\)로 바꾸면 \(E[S_i^2] = \sigma_i^2\)이므로 겨누는 값이 참 분산과 정확히 같아 \(R \equiv 1\)이 된다. 그리고 자유도를 고정하지 않고 같은 적률맞춤으로 \(\nu\)를 자료마다 다시 계산하므로, 위에서 합동 자유도가 틀렸던 \(m\)의 몫까지 제자리를 찾는다. 본문의 표에서 \(\nu\)가 \(\sigma_1/\sigma_2 = 4\), \((10,40)\)에서 \(9.4\)까지 내려가는 것이 그 값이다. 관측값 50개를 썼지만 정보는 10개짜리 표본만큼이라고 정직하게 세는 것이고, 합동 \(t\)는 같은 자리에서 48개라고 주장한다.

(1) 그림이 될 모양은 표에 이미 적혀 있다. 보기 1의 아홉 칸을 분산비로 묶으면 이렇게 된다.

• 왼쪽 패널(\(\sigma_1/\sigma_2 = 1\)). 여섯 막대가 모두 \(0.049 \sim 0.051\)이므로 점선 위에 나란히 선다. 높이차가 눈에 보이지 않아야 한다.
• 가운데와 오른쪽 패널. 파란 막대(Welch)는 여전히 여섯 개 모두 \(0.048 \sim 0.052\)라 점선에 붙어 있고, 빨간 막대(합동)만 가운데 칸에서 솟고(\(0.186 \to 0.291\)) 오른쪽 칸에서 주저앉는다(\(0.005 \to 0.001\)).
• 왼쪽 칸의 빨간 막대는 조금씩만 오른다. \(0.051 \to 0.054 \to 0.060\)이다. 균형 설계에서는 \(R = 1\)이라 합동 \(t\)가 버티기 때문이고, 이분산 하나만으로는 큰일이 나지 않는다는 보기 1의 결론이 여기에 그림으로 남는다.

한 문장으로 줄이면, 빨간 막대는 패널을 지날수록 좌우로 찢어지고 파란 막대는 아홉 칸 내내 같은 높이여야 한다.

(2) 그려서 확인한다.

import matplotlib.pyplot as plt
import numpy as np

rng = np.random.default_rng(1)   # 보기 1과 같은 난수열을 다시 쓴다

sizes = [(10, 10), (10, 40), (40, 10)]
ratios = [1, 2, 4]
labels = [f"$n_1={a},\\ n_2={b}$" for a, b in sizes]

fig, axes = plt.subplots(1, 3, figsize=(12, 3.8), sharey=True)
xs = np.arange(len(sizes))
for ax, r in zip(axes, ratios):
    pool, welch = [], []
    for n1, n2 in sizes:
        rp, rw, nu = error_rates(n1, n2, r, 1.0)     # 보기 1의 함수
        pool.append(rp)
        welch.append(rw)
    ax.bar(xs - 0.19, pool, 0.36, color="#c44e52", edgecolor="white", label="pooled $t$")
    ax.bar(xs + 0.19, welch, 0.36, color="#4c72b0", edgecolor="white", label="Welch $t$")
    for x, v in zip(xs - 0.19, pool):
        ax.text(x, v + 0.006, f"{v:.3f}", ha="center", fontsize=7)
    for x, v in zip(xs + 0.19, welch):
        ax.text(x, v + 0.006, f"{v:.3f}", ha="center", fontsize=7)
    ax.axhline(0.05, color="black", ls="--", lw=1, label="nominal 0.05")
    ax.set_xticks(xs)
    ax.set_xticklabels(labels, fontsize=8)
    ax.set_title(rf"$\sigma_1/\sigma_2 = {r}$")
    ax.set_ylim(0, 0.34)

axes[0].set_ylabel("actual type I error rate")
axes[2].legend(fontsize=8)
plt.tight_layout()
plt.show()

왼쪽 패널(등분산)에서는 여섯 개의 막대가 모두 점선 위에 나란히 있다. 오른쪽으로 갈수록 빨간 막대만 튀어 오르거나 주저앉고 파란 막대는 점선에 붙어 있다. 파란 막대가 아홉 칸 내내 같은 높이라는 것이 이 그림의 요점이다.

예측한 모양 그대로다. 막대 위에 찍힌 수가 보기 1의 표와 한 자리도 다르지 않다. 같은 씨앗으로 같은 격자를 다시 돈 것이니 그래야 맞다.

이 그림이 가리는 것이 하나 있다. 세로축이 \(0\)에서 \(0.34\)까지의 선형 눈금이므로 \(0.005\)와 \(0.001\)이 바닥에 깔려 구별되지 않는다. 오른쪽 패널의 \((40,10)\) 칸에서 빨간 막대는 아예 보이지 않는다. 그런데 거기서 일어난 일은 "별일 없음"이 아니라 명목 5% 검정의 실제 오류율이 0.1%로 내려앉은 것이고, 그만큼 검정력을 통째로 잃었다는 뜻이다. 보수적인 고장은 그림에서 지워지고 관대한 고장만 솟아 보인다.

막대그림이 보이지 않는 것이 하나 더 있다. 몬테카를로 오차다. 5만 번에서 \(0.05\) 둘레의 오차가 \(\pm 0.001\)인데 막대 두께에 묻혀 사라진다. 파란 막대들이 "정확히 같은 높이"로 보이는 것은 사실 \(\pm 0.002\) 안에서 같다는 뜻이며, 그 구별은 보기 1의 표로 돌아가야 읽힌다. 그림은 패턴을 보이고 표는 크기를 보인다.

(1) 극한에 \(N\)이 남지 않는다. \(n_i\)가 크면 \(n_i - 1 \approx n_i\)로 두어도 극한이 바뀌지 않으므로

이고 \(n_i = c_iN\)을 넣으면 앞 분수에서 \(N\)이 약분되고 뒤 괄호에서 \(1/N\)이 떨어져 나온다. 참 분산 \(\frac{1}{N}\!\left(\frac{\sigma_1^2}{c_1}+\frac{\sigma_2^2}{c_2}\right)\)에도 같은 \(1/N\)이 있으므로 비를 잡으면

이 남는다. \(N\)이 사라졌다. \(R\)을 정하는 것은 두 비 \(c_1:c_2\)와 \(\sigma_1^2:\sigma_2^2\)뿐이다.

\(\sigma_1^2:\sigma_2^2 = 16:1\)을 넣어 두 배치를 계산한다.

표본이 크면 자유도도 함께 커져 \(t\)가 \(Z\)로 가므로, 극한 오류율은

이다. 어느 쪽도 \(0.05\)로 가지 않는다.

(2) 유한한 \(N\)에서의 예측. 보기 1에서 쓴 \(T_{\text{pool}} \approx t_m/\sqrt R\)을 그대로 적용하면 네 줄의 예측값이 나온다(\(m\)은 합동분산을 맞춘 카이제곱의 자유도다).

\(R\)이 극한 \(0.3077\)과 \(3.2500\)을 향해 옆걸음으로 다가가고, 오류율은 \(0.2942 \to 0.2786\)과 \(0.0007 \to 0.0004\)로 움직인다. 줄어드는 것이 아니라 극한값에 가 닿는 중이다.

이제 모의실험.

import numpy as np
from scipy import stats

rng = np.random.default_rng(1)
B = 50_000
alpha = 0.05

# error_rates 는 보기 1과 같다(생략).

print("sigma1/sigma2 = 4,  n1 : n2 = 1 : 4 를 유지하며 표본을 키운다")
print(" (n1,  n2)    pooled t   Welch t   mean nu")
for n1, n2 in ((10, 40), (20, 80), (40, 160), (100, 400)):
    rp, rw, nu = error_rates(n1, n2, 4.0, 1.0)
    print(f"({n1:>3}, {n2:>3})     {rp:.3f}      {rw:.3f}    {nu:6.1f}")

print()
print("sigma1/sigma2 = 4,  n1 : n2 = 4 : 1 을 유지하며 표본을 키운다")
print(" (n1,  n2)    pooled t   Welch t   mean nu")
for n1, n2 in ((40, 10), (80, 20), (160, 40), (400, 100)):
    rp, rw, nu = error_rates(n1, n2, 4.0, 1.0)
    print(f"({n1:>3}, {n2:>3})     {rp:.3f}      {rw:.3f}    {nu:6.1f}")

출력:

sigma1/sigma2 = 4,  n1 : n2 = 1 : 4 를 유지하며 표본을 키운다
 (n1,  n2)    pooled t   Welch t   mean nu
( 10,  40)     0.292      0.050       9.4
( 20,  80)     0.286      0.050      19.7
( 40, 160)     0.282      0.050      40.3
(100, 400)     0.275      0.050     102.2

sigma1/sigma2 = 4,  n1 : n2 = 4 : 1 을 유지하며 표본을 키운다
 (n1,  n2)    pooled t   Welch t   mean nu
( 40,  10)     0.001      0.051      46.1
( 80,  20)     0.001      0.049      96.0
(160,  40)     0.001      0.052     196.0
(400, 100)     0.000      0.050     496.0

표본을 10배로 늘려도 합동 \(t\)의 오류율이 0.292에서 0.275로 움직일 뿐이다. 줄어드는 것이 아니라 극한값으로 다가가는 중이다. \(n_1 = c_1N\), \(n_2 = c_2N\)으로 두고 \(N \to \infty\)를 취하면

인데 이 값에 \(N\)이 없다. \(c_1:c_2 = 1:4\), \(\sigma_1^2:\sigma_2^2 = 16:1\)이면 \(R \to 0.308\)이고, 대응하는 오류율의 극한은 \(P(|Z| > 1.96\sqrt{0.308}) = 0.277\)이다. 표에서 실제로 그 값으로 수렴하고 있다.

반대 방향도 마찬가지다. \(4:1\) 배치의 오류율은 0.001 근처에 머물다 극한 0.0004로 간다. 검정이 사실상 아무것도 기각하지 않는다.

예측과 모의가 네 줄씩 맞는다. 소수 넷째 자리까지 다시 찍으면 위 블록이 \(0.2924,\ 0.2856,\ 0.2819,\ 0.2748\)을, 아래 블록이 \(0.0010,\ 0.0007,\ 0.0005,\ 0.0004\)를 준다.

반복 5만 회에서 \(p \approx 0.28\)의 몬테카를로 오차가 \(\sqrt{0.28 \times 0.72/50000} = 0.0020\)이므로 왼쪽 네 줄은 모두 \(2\) 오차 안이다. 오른쪽은 극한 \(0.0004\)에 자리까지 맞춰 내려간다.

읽어야 할 것은 수렴하는 방향이다. 표본을 10배로 늘려 \(0.292\)가 \(0.275\)가 되었다. 이것을 "좋아지고 있다"로 읽으면 안 된다. \(0.05\)가 아니라 \(0.277\)로 가는 중이고, 100배를 더 써도 \(0.277\)에서 멈춘다. 극한식에 \(N\)이 들어 있지 않기 때문이다. 반대 배치는 \(0.0004\)로 내려가 검정이 사실상 아무것도 기각하지 않는 자리로 간다.

이것이 다음 쪽과 갈리는 지점이다. 비정규성이라면 중심극한정리가 표본크기로 구해 준다. 설계의 불균형은 표본크기가 손댈 수 없다. 고치려면 설계를 바꾸거나(\(n_1 = n_2\)로 맞추거나) 방법을 바꿔야(Welch를 쓰거나) 한다. Welch 열이 여덟 줄 내내 \(0.049\)에서 \(0.052\) 사이에 머무는 것이 뒤쪽 답이다.

(1) \(\nu\)가 갇히는 구간. \(n_1 = n_2 = n\)이고 \(\sigma_1 = \sigma_2\)이면 \(a = S_1^2/n\)과 \(b = S_2^2/n\)이 독립이고 같은 분포를 따른다. \(d = n-1\)로 두고 \(u = a/(a+b)\)라 하면

이다. \(a\)와 \(b\)가 각각 척도가 같은 \(\chi^2_d\)에 비례하므로 \(u \sim \text{Beta}(d/2,\, d/2)\)이고, 척도는 \(u\)에서 약분되어 사라진다. 분모 \(u^2+(1-u)^2 = 1 - 2u(1-u)\)는 \(u \in [0,1]\)에서 \(1/2\)(가운데)과 \(1\)(끝점) 사이를 움직이므로

이다. 위쪽 끝은 두 표본분산이 정확히 같을 때, 아래쪽 끝은 한쪽이 전부를 지배할 때다. 연습문제 4의 일반 한계 \(\min(n_i-1) \le \nu \le n_1+n_2-2\)가 균형 설계에서 \([n-1,\ 2n-2]\)로 좁혀진 꼴이다.

\(E[\nu] = d\,E\!\left[1/(1-2u(1-u))\right]\)은 초기하함수로 적을 수 있지만, 베타밀도를 한 번 적분해 수로 내는 편이 빠르다.

import numpy as np
from scipy import integrate, stats

# 등분산·균형(n1 = n2 = n)에서 웰치 자유도는 u ~ Beta(d/2, d/2) 하나로 줄어든다.
#   nu = (n-1) / (u^2 + (1-u)^2),   u = a/(a+b)
# 분모가 1/2 과 1 사이이므로 nu 는 n-1 과 2(n-1) 사이에 갇힌다.
print(f"{'n':>4}{'하한 n-1':>10}{'E[nu] 이론':>12}{'상한 2(n-1)':>12}{'t_.975(E nu)':>14}{'t_.975(2n-2)':>14}")
for n in (5, 10, 20, 50):
    d = n - 1
    f = lambda u: stats.beta(d / 2, d / 2).pdf(u) / (u ** 2 + (1 - u) ** 2)
    val, _ = integrate.quad(f, 0, 1)
    Enu = d * val
    print(f"{n:>4}{d:>10}{Enu:>12.2f}{2 * d:>12}"
          f"{stats.t(Enu).ppf(0.975):>14.4f}{stats.t(2 * n - 2).ppf(0.975):>14.4f}")

출력:

   n    하한 n-1    E[nu] 이론   상한 2(n-1)  t_.975(E nu)  t_.975(2n-2)
   5         4        6.85           8        2.3752        2.3060
  10         9       16.54          18        2.1143        2.1009
  20        19       36.32          38        2.0275        2.0244
  50        49       96.14          98        1.9849        1.9845

\(E[\nu]\)가 상한 \(2(n-1)\)에 바짝 붙어 있다. \(n = 10\)에서 \(16.54\) 대 \(18\), \(n = 50\)에서 \(96.14\) 대 \(98\)이다. 등분산이면 \(u\)가 \(1/2\) 둘레에 몰려 \(\nu\)가 위쪽 끝으로 쏠리기 때문이다. 그 결과 임계값의 차이도 작다. \(n = 5\)에서 \(2.3752\) 대 \(2.3060\)으로 \(3.0\%\), \(n = 10\)에서 \(0.6\%\), \(n = 50\)에서 \(0.02\%\)다. Welch가 치르는 비용은 이 임계값 차이가 전부이므로 표본과 함께 빠르게 사라질 것이라 예측된다.

(2) 합동 \(t\)의 검정력은 정확히 계산된다. 등분산이 실제로 맞으므로 \(H_1\) 아래에서 합동 \(t\)는 비중심 \(t\)를 따르고, 비중심모수는 \(\mathrm{ncp} = 1/\sqrt{2/n}\)이다. Welch 쪽은 \(\nu\)가 확률변수라 정확한 식이 없으므로, \(E[\nu]\)를 꽂아 넣은 어림으로 둔다.

이제 모의실험.

import numpy as np
from scipy import stats

rng = np.random.default_rng(1)
B = 50_000
alpha = 0.05

print("등분산 정규모집단에서 Welch를 쓰는 값: 크기와 검정력")
print("  n     size(pool) size(Welch)   power(pool) power(Welch)  mean nu")
for n in (5, 10, 20, 50):
    x1 = rng.normal(0, 1, size=(B, n))
    x2 = rng.normal(0, 1, size=(B, n))
    v1, v2 = x1.var(axis=1, ddof=1), x2.var(axis=1, ddof=1)
    sp2, a, b = (v1 + v2) / 2, v1 / n, v2 / n
    nu = (a + b) ** 2 / (a ** 2 / (n - 1) + b ** 2 / (n - 1))
    c_pool = stats.t(2 * n - 2).ppf(1 - alpha / 2)
    c_welch = stats.t(nu).ppf(1 - alpha / 2)

    # H0가 참인 경우(크기)와 참 차이가 1인 경우(검정력)를 같은 표본으로 비교한다.
    d0 = x1.mean(axis=1) - x2.mean(axis=1)
    out = []
    for d in (d0, d0 + 1.0):
        out += [np.mean(np.abs(d / np.sqrt(sp2 * 2 / n)) > c_pool),
                np.mean(np.abs(d / np.sqrt(a + b)) > c_welch)]
    print(f"{n:>4}      {out[0]:.3f}       {out[1]:.3f}         "
          f"{out[2]:.3f}       {out[3]:.3f}    {nu.mean():5.1f}")

출력:

등분산 정규모집단에서 Welch를 쓰는 값: 크기와 검정력
  n     size(pool) size(Welch)   power(pool) power(Welch)  mean nu
   5      0.048       0.042         0.287       0.265      6.9
  10      0.049       0.048         0.563       0.557     16.5
  20      0.051       0.051         0.868       0.867     36.3
  50      0.049       0.049         0.998       0.998     96.1

등분산이 실제로 맞는 상황에서 Welch가 잃는 것을 잰 표다. 집단당 5개일 때 검정력이 0.287에서 0.265로 2.2%포인트 떨어지고, 10개면 0.6%포인트, 20개 이상이면 차이가 사라진다.

잃는 것은 이만큼이고 얻는 것은 0.291 대 0.050이다. 저울이 한쪽으로 완전히 기운다.

합동 열은 정확히 맞고 Welch 열은 한 자리에서 어긋난다.

합동 열은 네 줄 모두 소수 셋째 자리까지 맞는다. 등분산·정규에서 그 분포가 정확히 비중심 \(t\)이므로 당연하다.

Welch 열은 \(n = 5\)에서 \(0.2762\)를 예측했으나 모의는 \(0.265\)를 주어 \(0.011\) 어긋난다. 나머지 세 줄은 \(0.001\) 안에서 맞는다. 까닭은 같은 표에 이미 적혀 있다. \(n = 5\)에서 Welch의 크기가 \(0.042\)로 명목 \(0.05\)에 못 미친다. 보수적인 검정은 검정력도 함께 잃는다. 어림은 "크기가 \(0.05\)인 \(t_{6.85}\) 검정"을 가정했는데 실제 검정은 그렇지 않은 것이다. 더 근본적으로 \(\nu\)는 확률변수이고 분모 \(S_1^2/n+S_2^2/n\)과 같은 자료에서 나와 서로 얽혀 있으므로, 평균 \(E[\nu]\)를 꽂아 넣는 어림이 자유도가 한 자릿수일 때 깨진다. \(n = 10\)부터는 그 얽힘이 보이지 않을 만큼 작아진다.

저울의 두 접시를 같은 표에서 읽는다. 등분산이 실제로 맞을 때 Welch가 잃는 것은 집단당 5개에서 검정력 \(2.2\)%포인트, 10개에서 \(0.6\)%포인트, 20개 이상에서 \(0\)이다. 등분산이 틀렸을 때 합동 \(t\)가 잃는 것은 보기 1의 \(0.291\) 대 \(0.050\), 곧 오류율 자체다. 앞은 정도의 문제이고 뒤는 타당성의 문제다.

\(a = s_1^2/n_1 = 1.6\), \(b = s_2^2/n_2 = 0.025\)로 두면 \(a + b = 1.625\)이므로 분자는

이고 분모는

이다. 따라서

이다. 임계값은 다음과 같다.

Welch 자유도가 9.28로 보수적 자유도 9에 거의 붙어 있다. 관측값 50개를 썼는데 정보는 10개짜리 표본만큼이라는 뜻이다. 두 번째 집단이 40개나 되지만 분산이 16분의 1이라 차의 분산에 \(0.025/1.625 = 1.5\%\)밖에 기여하지 않기 때문이다.

이것이 합동 \(t\)가 오류율 0.291을 내는 이유이기도 하다. 합동 \(t\)는 이 상황에서 자유도가 48이라고 주장한다. 없는 정보를 있다고 세는 것이다.

참 분산은

이고, 합동분산이 겨누는 값은

이다. 비는 \(R = 0.4766/1.625 = 0.293\)이고 \(\sqrt R = 0.542\)이다.

분모가 참값의 0.542배이므로 \(t\) 통계량은 평균적으로 \(1/0.542 = 1.85\)배로 부푼다. 임계값 2.011을 넘기가 그만큼 쉬워진다. 실제로 기각률을 근사하면

이고 모의실험 값 0.291과 가깝다. 근사가 조금 작게 나오는 것은 \(S_p^2\) 자체의 흔들림을 무시했기 때문이다.

관측값 50개 중 40개가 분산이 작은 집단에 있다는 사실이 합동분산을 작게 만든다. 그런데 차의 분산은 표본이 적은 집단이 지배한다. 이 어긋남이 전부다.

\(n_1 = n_2 = n\)을 넣으면

이고 이것이 참 분산이다. \(\square\) 자유도 가중과 \(1/n_i\) 가중이 둘 다 대칭이 되어 서로 상쇄된다.

기댓값이 맞는다고 분포가 같은 것은 아니다. 합동 \(t\)가 정확히 \(t_{2n-2}\)가 되려면 분모가 참 분산에 비례하는 카이제곱이어야 하는데, 이분산이면

처럼 척도가 다른 두 카이제곱의 합이 되어 카이제곱이 아니다. 평균은 맞지만 분산이 카이제곱보다 크다. 분모가 더 흔들리면 통계량의 꼬리가 무거워지고 기각이 조금 늘어난다.

분산비가 커질수록 이 효과가 커진다. 표에서 \(\sigma_1/\sigma_2\)가 1, 2, 4일 때 \((10,10)\)의 오류율이 0.051, 0.054, 0.060으로 서서히 오르는 것이 그것이다. 그러나 이 정도는 실무에서 감당할 만한 어긋남이며, 부등한 표본크기가 만드는 0.291과는 종류가 다르다.

\(a = S_1^2/n_1\), \(b = S_2^2/n_2\), \(d_i = n_i - 1\)로 두면

이다.

위쪽 한계. 코시-슈바르츠 부등식을 쓴다.

양변을 분모로 나누면 \(\nu \le d_1+d_2 = n_1+n_2-2\)이다. 등호는 \(a/d_1 = b/d_2\)일 때, 즉 두 집단이 분산과 자유도 면에서 균형을 이룰 때다.

아래쪽 한계. 일반성을 잃지 않고 \(d_1 \le d_2\)라 하자. 그러면

이 성립하므로(\(b^2/d_2 \le b^2/d_1 \le (b^2+2ab)/d_1\)) \(\nu \ge d_1 = \min(n_1-1, n_2-1)\)이다. \(\square\) 등호는 한쪽이 전부를 지배할 때, 즉 \(b \to 0\)일 때 극한으로 다가간다.

연습문제 1이 그 극한에 가까운 예다. \(b/a = 0.0156\)이라 \(\nu = 9.28\)로 하한 9에 거의 붙었다.

정규모집단에서 \((n_i-1)S_i^2/\sigma_i^2 \sim \chi^2_{n_i-1}\)이므로

이다. \(a_i = \sigma_i^2/n_i\)로 두면 두 표본이 독립이므로

이다.

이제 \(V \approx \theta\chi^2_\nu/\nu\)로 두자. 이 근사분포는

를 갖는다. 평균을 맞추면 \(\theta = a_1+a_2\)이고, 분산을 맞추면

이므로

이다. \(\square\) 실제로는 \(\sigma_i^2\)을 모르므로 \(S_i^2\)을 대입하며, 그래서 \(\nu\)가 확률변수가 된다.

검산. \(\sigma_1 = \sigma_2\)이고 \(n_1 = n_2 = n\)이면 \(a_1 = a_2 = a\)이므로

로 합동 자유도와 같아진다. 등분산·균형 설계에서 Welch가 합동 \(t\)로 되돌아온다는 뜻이다.

\(n_i\)가 크면 \(n_i - 1 \approx n_i\), \(n_1+n_2-2 \approx n_1+n_2\)로 두어도 극한은 바뀌지 않는다. 그러면

이고 \(n_i = c_i N\)을 넣으면 \(N\)이 분자와 분모에서 한 번씩 약분되어

이 된다. 참 분산도 \(\frac{1}{N}\left(\frac{\sigma_1^2}{c_1}+\frac{\sigma_2^2}{c_2}\right)\)이므로 비를 잡으면 \(1/N\)이 약분되어

이 남는다. 표본크기 \(N\)이 사라졌다. \(\square\) \(R\)을 정하는 것은 두 비 \(c_1:c_2\)와 \(\sigma_1:\sigma_2\)뿐이다.

수치. \(c_1=1\), \(c_2=4\), \(\sigma_1^2=16\), \(\sigma_2^2=1\)을 넣으면

이다. 표본이 크면 \(t\)가 \(Z\)로 가므로 오류율의 극한은

이다. 보기 3의 표에서 0.292, 0.286, 0.282, 0.275로 이 값 근처로 다가가고 있다. 0.05로 가지 않는다.

반대 배치 \(c_1:c_2 = 4:1\)이면 \(R \to 3.25\)이고 오류율의 극한은 \(P(|Z| > 1.96\times1.803) = 0.0004\)다. 표의 0.001, 0.001, 0.001, 0.000과 맞는다.

이 성질이 이 쪽의 핵심이다. 비정규성이라면 중심극한정리가 표본크기로 구해 주지만, 설계의 불균형은 표본크기로 고칠 수 없다. 고치려면 설계를 바꾸거나(\(n_1 = n_2\)로 맞추거나) 방법을 바꿔야(Welch를 쓰거나) 한다. 다음 쪽에서 이 대조를 정면으로 다룬다.

import numpy as np
from scipy import stats

rng = np.random.default_rng(0)
REP, d = 40_000, 0.8
print(f"{'n1,n2':>10}{'합동 t':>12}{'웰치':>10}{'손실':>10}")
for n1, n2 in [(10,10), (20,20), (10,30), (5,45)]:
    a = rng.normal(0, 1, (REP, n1))
    b = rng.normal(d, 1, (REP, n2))
    p1 = stats.ttest_ind(a, b, axis=1, equal_var=True).pvalue
    p2 = stats.ttest_ind(a, b, axis=1, equal_var=False).pvalue
    r1, r2 = (p1 < 0.05).mean(), (p2 < 0.05).mean()
    print(f"{f'{n1},{n2}':>10}{r1:>12.4f}{r2:>10.4f}{r1-r2:>+10.4f}")

출력:

     n1,n2        합동 t        웰치        손실
     10,10      0.3968    0.3917   +0.0051
     20,20      0.6964    0.6953   +0.0011
     10,30      0.5690    0.5374   +0.0315
      5,45      0.3838    0.3001   +0.0837

표본크기가 같으면 손실이 사실상 없다. \((10,10)\)에서 \(0.005\), \((20,20)\)에서 \(0.001\)이다. 연습문제 3에서 본 대로 \(n_1=n_2\)이면 두 방법의 통계량이 같고 자유도만 조금 다르기 때문이다.

불균형하면 손실이 커진다. \((5,45)\)에서 \(0.084\)로, 검정력의 \(22\%\)를 잃는다. 웰치–새터스웨이트 자유도가 작은 쪽으로 끌려가(\(\nu \approx 4.9\)) 임계값이 커지기 때문이다(연습문제 4).

그래도 웰치를 기본으로 쓰는 이유.

잃는 것과 지키는 것의 크기가 다르다. 등분산일 때 웰치가 잃는 것은 검정력 \(0.005 \sim 0.08\)이고, 이분산일 때 합동이 잃는 것은 유의수준 자체다(다음 연습문제에서 \(0.21\)까지 치솟는다). 전자는 정도의 문제이고 후자는 타당성의 문제다.

불균형이 심하면 설계를 고치는 편이 낫다. \((5,45)\)처럼 극단적인 배분은 검정력 면에서도 나쁘다(연습문제 10의 최적 배분). 웰치의 손실을 걱정하기 전에 배분을 균등에 가깝게 만드는 것이 먼저다.

소수 자유도는 근사의 산물이다. 연습문제 5에서 본 대로 \(V = S_1^2/n_1 + S_2^2/n_2\)의 분포를 척도화 카이제곱 \(\theta\chi^2_\nu/\nu\)로 적률맞춤했는데, 평균과 분산 두 개를 맞추다 보면 \(\nu\)가 실수로 나온다. 실제로 \(V\)가 카이제곱을 따르는 것이 아니라 가장 닮은 카이제곱을 고른 것이다.

\(t\) 분포는 실수 자유도에서도 잘 정의된다. 밀도함수가 감마함수로 쓰이므로 \(\nu = 12.7\) 같은 값을 그대로 넣을 수 있다. scipy.stats.t.ppf(0.975, 12.7)이 문제없이 계산된다.

import numpy as np
from scipy import stats

n1, n2, s1, s2 = 10, 40, 4.0, 1.0
v1, v2 = s1**2/n1, s2**2/n2
nu = (v1 + v2)**2 / (v1**2/(n1-1) + v2**2/(n2-1))
print(f"  웰치 자유도 nu = {nu:.4f}")
for lab, d in [("그대로", nu), ("내림", np.floor(nu)),
               ("보수적 min(n-1)", min(n1-1, n2-1)), ("합동", n1+n2-2)]:
    print(f"    {lab:>16}: df={d:>7.3f}  t_0.975={stats.t.ppf(0.975, d):.4f}")

출력:

  웰치 자유도 nu = 9.2829
                 그대로: df=  9.283  t_0.975=2.2517
                  내림: df=  9.000  t_0.975=2.2622
        보수적 min(n-1): df=  9.000  t_0.975=2.2622
                  합동: df= 48.000  t_0.975=2.0106

내림과 그대로 쓰기의 차이가 미미하다. 임계값이 \(2.2517\) 대 \(2.2622\)로 \(0.5\%\) 차이다. 자유도가 클수록 차이는 더 줄어든다.

내림은 표를 보던 시절의 유물이다. 자유도 표에 정수만 있었으므로 보수적인 쪽(작은 자유도)으로 내려 읽었다. 컴퓨터로 계산하는 지금은 그대로 쓰는 것이 표준이며 scipy, R, SAS 모두 그렇게 한다.

보수적 자유도 \(\min(n_1-1, n_2-1)\)는 다른 물건이다. 이것은 근사가 아니라 하한(연습문제 4)이라 언제나 안전하지만 검정력을 더 잃는다. 이 예에서는 \(\nu = 9.28\)이 이미 하한 \(9\)에 거의 붙어 있어 둘이 같은 \(2.2622\)를 주지만, 분산비가 덜 극단적이면 둘이 벌어진다. 계산기 없이 손으로 빠르게 판단해야 할 때의 어림이며, 정식 보고에는 쓰지 않는다.

합동 자유도 \(48\)을 쓰면 임계값이 \(2.0106\)으로 크게 작아진다. 이 자료는 \(\sigma_1 = 4\sigma_2\)로 분산이 크게 다르므로 이 임계값은 부당하게 관대하다. 연습문제 2에서 본 오류율 폭증이 바로 이 차이에서 온다.

import numpy as np
from scipy import stats

rng = np.random.default_rng(0)
REP = 40_000
print(f"{'n1,n2':>9}{'s1/s2':>7}{'항상 합동':>12}{'항상 웰치':>12}{'사전검정 후':>14}")
for n1, n2, r in [(10,30,1.0), (10,30,3.0), (30,10,3.0)]:
    a = rng.normal(0, r, (REP, n1))
    b = rng.normal(0, 1, (REP, n2))
    pp = stats.ttest_ind(a, b, axis=1, equal_var=True).pvalue
    pw = stats.ttest_ind(a, b, axis=1, equal_var=False).pvalue
    F = a.var(1, ddof=1) / b.var(1, ddof=1)
    lo, hi = stats.f.ppf(0.025, n1-1, n2-1), stats.f.ppf(0.975, n1-1, n2-1)
    two = np.where((F >= lo) & (F <= hi), pp, pw)     # 사전검정 후 선택
    print(f"{f'{n1},{n2}':>9}{r:>7.1f}{(pp<0.05).mean():>12.4f}"
          f"{(pw<0.05).mean():>12.4f}{(two<0.05).mean():>14.4f}")

출력:

    n1,n2  s1/s2       항상 합동       항상 웰치        사전검정 후
    10,30    1.0      0.0488      0.0497        0.0498
    10,30    3.0      0.2114      0.0520        0.0549
    30,10    3.0      0.0042      0.0503        0.0494

항상 합동이 최악이다. 작은 표본 쪽 분산이 크면(\(10,30\)에 \(\sigma\)비 \(3\)) 오류율이 \(0.21\)로 명목의 네 배다. 반대로 큰 표본 쪽 분산이 크면(\(30,10\)) \(0.0042\)로 지나치게 보수적이다. 어느 쪽이든 명목 수준을 지키지 못한다.

항상 웰치가 세 상황 모두에서 \(0.05\) 근처다(\(0.0497 \sim 0.0520\)).

사전검정 절차는 웰치보다 나을 것이 없다. \(0.0549\)로 약간 부풀었고, 등분산일 때도 웰치와 사실상 같은 \(0.0498\)이다. 번거로움만 늘고 얻는 것이 없다.

왜 사전검정이 도움이 안 되는가.

• 사전검정의 검정력이 낮다. 5.9절에서 본 대로 \(F\) 검정은 표본이 작을 때 분산비 \(3\)도 절반쯤밖에 못 잡는다. 못 잡으면 합동 \(t\)로 넘어가 오류율이 오른다.
• 조건부 분포가 달라진다. "\(F\) 검정을 통과했다"는 조건 아래에서 \(S_p^2\)의 분포는 원래 분포가 아니다. 선택 자체가 자료에 의존하므로 두 번째 단계의 이론이 깨진다. 9장의 선택 후 추론 문제와 같은 구조다.
• \(F\) 검정 자체가 정규성에 취약하다. 비정규 자료에서는 사전검정이 엉뚱한 판정을 내린다(5.9절).

결론은 단순하다. 그냥 웰치를 쓴다. 등분산일 때 잃는 검정력은 미미하고(연습문제 7), 이분산일 때 지키는 것은 검정의 타당성 자체다. R의 t.test가 웰치를 기본값으로 두는 이유이며, 15장에서 이 논점을 등분산 검정 전반으로 확장한다.

통상의 \(F\)는 모든 집단이 같은 \(\sigma^2\)을 갖는다고 가정한다. 집단 내 제곱합을 합쳐 하나의 \(S_p^2\)을 만들고 그것으로 나눈다. 집단마다 분산이 다르면 이 합치기가 부당해진다.

웰치 분산분석은 가중치를 분산의 역수로 준다.

로 정밀한 집단에 큰 가중치를 주어 전체 평균을 잡고, 통계량은

형태가 된다. 분모의 \(\Lambda\)는 표본크기와 가중치의 불균형을 반영하는 보정항이고, 자유도도 새터스웨이트 방식으로 근사한다.

두 표본에서 웰치 \(t\)로 환원된다. \(k=2\)면 위 식이 \(t_W^2\)이 되며, 연습문제 9에서 본 \(F = t^2\) 관계의 웰치판이다.

무엇이 달라지는가.

실무 권장도 두 표본에서와 같다. 집단별 분산이 비슷하다는 확신이 없으면 웰치 분산분석을 쓴다. R의 oneway.test(var.equal=FALSE)가 이것이며, 11장에서 자세히 다룬다.

다만 사후비교가 까다로워진다. 통상의 투키 방법은 등분산을 전제하므로, 웰치 틀에서는 게임스–하웰 같은 별도 절차를 써야 한다. 가정을 느슨하게 하면 뒤따르는 도구도 전부 바뀐다는 점이 이분산 문제의 성가신 면이다.