μ₁ − μ₂의 신뢰구간¶

(1) 분산이 더해진다. 두 표본이 독립이므로

이고(차를 빼도 분산은 더해진다), 표본값을 넣으면

이다. 표준편차를 더하면 틀린다. 가법성이 있는 것은 분산이지 표준편차가 아니다. 더해 보면

로 참값보다 \(39\%\) 크다. 구간이 그만큼 넓어지므로 틀리되 보수적인 쪽이지만, 틀린 것은 틀린 것이다.

자유도는

이고 \(t_{0.025,\,43.8074} = 2.01562\) 이므로

구간이 \(0\) 을 아슬아슬하게 벗어난다. 점추정값은 차이가 \(10\) 이라 하지만 구간은 \(0.23\) 일 수도 \(19.77\) 일 수도 있다고 말한다. "차이가 있다"까지는 말할 수 있어도 "얼마나 있다"는 거의 말하지 못하는 자료다.

(2) 자유도의 범위. \(a = s_1^2/n_1\), \(b = s_2^2/n_2\) 로 두면

이다.

위 끝 — 코시–슈바르츠. 두 벡터 \(\left(\sqrt{n_1-1},\ \sqrt{n_2-1}\right)\) 과 \(\left(\dfrac{a}{\sqrt{n_1-1}},\ \dfrac{b}{\sqrt{n_2-1}}\right)\) 에 코시–슈바르츠를 쓰면

이므로

이고, 등호는 두 벡터가 평행할 때, 곧

곧 \(s_1/s_2 = 1.2042\) 일 때뿐이다. 웰치의 자유도는 절대로 합동 \(t\) 의 자유도를 넘지 못한다. \(\square\)

아래 끝. \(s_1/s_2 \to 0\) 이면 \(a \to 0\) 이라 \(\nu \to n_2 - 1 = 24\) 이고, \(s_1/s_2 \to \infty\) 면 \(\nu \to n_1 - 1 = 29\) 다. 따라서

이다. 이 자료는 \(s_1/s_2 = 0.75\) 이므로 \(\nu = 43.81\) 로 범위의 가운데쯤에 있다.

뜻. 웰치는 등분산을 가정하지 않는 대가로 자유도를 내놓는다. 그 대가를 수로 보면, 같은 자료의 합동 구간이 \((0.5259,\ 19.4741)\) 로 웰치보다 폭이 \(3.0\%\) 좁다. 등분산이 참이면 이 \(3\%\) 가 공짜로 얻는 정밀도이고, 거짓이면 그만큼 틀린 자리에 놓인 좁음이다. 두 구간 모두 \(0\) 을 배제하므로 여기서는 결론이 갈리지 않는다.

(3) 확인.

import numpy as np
import scipy.stats as stats

n1, n2 = 30, 25
mean1, mean2 = 100, 90
s1, s2 = 15, 20
confidence_level = 0.95

# 두 표본이 독립이므로 분산이 더해진다.
# 표준편차가 아니라 **분산**을 더한 뒤 제곱근이라는 점에 주의.
standard_error = np.sqrt((s1**2 / n1) + (s2**2 / n2))

# Welch-Satterthwaite 자유도. 정수가 아니어도 scipy는 받아 준다.
# 두 표본의 정보량을 하나로 합치는 근사이며, 한쪽 분산이 압도하면
# 그쪽 표본의 자유도 쪽으로 끌려간다.
df = ((s1**2 / n1) + (s2**2 / n2))**2 / (
    ((s1**2 / n1)**2 / (n1 - 1)) + ((s2**2 / n2)**2 / (n2 - 1))
)

t_critical = stats.t.ppf(1 - (1 - confidence_level) / 2, df)
margin_of_error = t_critical * standard_error

confidence_interval = (
    (mean1 - mean2) - margin_of_error,
    (mean1 - mean2) + margin_of_error,
)
print(f"{standard_error = :.5f},  {df = :.4f},  {t_critical = :.5f}")
print(f"{confidence_interval = }")

# (2) 자유도의 위 끝은 n1+n2-2 다. s1/s2 를 훑어 확인한다.
print(f"\n{'s1/s2':>7}{'df':>10}")
for r in (0.01, 0.75, 1.2042, 2.0, 100.0):
    a, b = (r * s2) ** 2 / n1, s2**2 / n2
    print(f"{r:>7.4f}{(a + b) ** 2 / (a**2 / (n1 - 1) + b**2 / (n2 - 1)):>10.4f}")
print(f"n1-1 = {n1 - 1}, n2-1 = {n2 - 1}, n1+n2-2 = {n1 + n2 - 2}")

# 같은 자료의 합동 t 구간
sp2 = ((n1 - 1) * s1**2 + (n2 - 1) * s2**2) / (n1 + n2 - 2)
se_p = np.sqrt(sp2 * (1 / n1 + 1 / n2))
t_p = stats.t.ppf(0.975, n1 + n2 - 2)
print(f"\n합동 t 구간 ({10 - t_p * se_p:.4f}, {10 + t_p * se_p:.4f})"
      f"   폭의 비 {2 * t_p * se_p / (2 * margin_of_error):.4f}")

출력:

standard_error = 4.84768,  df = 43.8074,  t_critical = 2.01562
confidence_interval = (0.22892977648461788, 19.77107022351538)

  s1/s2        df
 0.0100   24.0040
 0.7500   43.8074
 1.2042   53.0000
 2.0000   44.2029
100.0000   29.0070
n1-1 = 29, n2-1 = 24, n1+n2-2 = 53

합동 t 구간 (0.5259, 19.4741)   폭의 비 0.9696

(1)과 맞는다. \(\text{SE} = 4.84768\), \(\nu = 43.8074\), 구간 \((0.2289,\ 19.7711)\) 이다.

(2)도 맞는다. \(s_1/s_2\) 를 훑으면 자유도가 \(s_1/s_2 = 1.2042\) 에서 정확히 \(53.0000 = n_1+n_2-2\) 로 꼭짓점을 찍고, 양쪽으로 가면서 \(24\)(\(s_1/s_2 \to 0\))와 \(29\)(\(s_1/s_2 \to \infty\))로 내려간다. 코시–슈바르츠가 준 위 끝과 등호 조건이 그대로 나타난다.

합동 구간의 폭이 웰치의 \(0.9696\) 배로 \(3.0\%\) 좁고, 두 구간 모두 \(0\) 을 배제한다.

Welch–Satterthwaite 자유도:

\(t_{0.025, 43} \approx 2.017\)이므로:

두 모평균의 참 차이가 0.22와 19.78 사이에 있다고 95% 신뢰한다.

(1) \(100\)회로는 비교가 되지 않는다. 포함률 \(p\) 의 추정값은 되풀이 \(M\) 회에서 표준오차 \(\sqrt{p(1-p)/M}\) 를 갖는다. \(p \approx 0.95\) 이면 \(M = 100\) 에서 \(0.0218\) 이다. 가리려는 차이가 \(0.010\) 이므로 오차가 차이의 두 배다.

차이의 \(95\%\) 구간이 \(0\) 을 벗어나려면, 두 방법을 따로 돌릴 때 차이의 표준오차가

이고 \(2 \times 1.96 \times \text{SE}_{\text{차이}} < 0.010\) 을 풀면

쯤이다. \(100\)회의 \(160\)배가 필요하다.

(2) 같은 표본을 쓰면 짝지은 비교가 된다. 표본 하나에 네 방법을 모두 적용하면, 그 표본이 유난히 극단적이었던 효과가 네 방법에 똑같이 실린다. 그 공통 부분이 차에서 빠지므로 차이의 분산이 크게 줄어든다.

실제로 웰치와 합동 \(t\) 의 포함 지시함수 차 \(\mathbf 1_W - \mathbf 1_P\) 는 두 방법의 판정이 갈릴 때만 \(0\) 이 아닌데, 그 일이 일어나는 비율이

로 합쳐 \(1.1\%\) 뿐이다. 그래서 차의 표준편차가 \(0.1039\) 로 따로 돌릴 때의 \(0.3082\) 보다 세 배 작다. 필요한 되풀이는

으로 약 \(1/10\) 로 줄어든다.

이것이 짝짓기의 일반 원리다. 같은 표본에 여러 방법을 돌리면 표본의 변동이 공통인자가 되어 상쇄된다. 8.4절의 대응표본 설계가 같은 이치이며, 거기서는 짝지은 사람이 공통인자다.

참값은 되풀이 \(200\)만 회로 이렇다. 웰치 \(0.9499\), 합동 \(t\) \(0.9399\), \(z\)-known \(0.9498\), \(z\)-plugin \(0.9317\)(몬테카를로 표준오차 \(0.0002\)). 웰치와 합동 \(t\) 의 차이가 \(0.0100\) 으로 (1)에서 가정한 값과 맞는다.

(3) 그림.

#!/usr/bin/env python3
"""두 평균의 차이에 대한 신뢰구간을 네 방법으로 만들어 포함확률을 비교한다.

Welch 는 두 분산이 다를 수 있다고 보고, 합동(pooled)은 같다고 본다.
분산이 실제로 다를 때 합동 방법의 포함확률이 어떻게 무너지는지가 요점이다.
"""

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import t, norm

rng_seed = 42        # 아래 그림을 재현하려면 고정한다
n_simulations = 100
n1, n2 = 12, 10
mu1, mu2 = 0.0, 0.5
sigma1, sigma2 = 1.0, 1.5
alpha = 0.05
method = "welch"  # 'welch' | 'pooled' | 'z_known' | 'z_plugin'


def main():
    if rng_seed is not None:
        np.random.seed(rng_seed)

    delta_true = mu1 - mu2
    lowers = np.empty(n_simulations)
    uppers = np.empty(n_simulations)
    centers = np.empty(n_simulations)

    for i in range(n_simulations):
        x = np.random.normal(loc=mu1, scale=sigma1, size=n1)
        y = np.random.normal(loc=mu2, scale=sigma2, size=n2)
        xbar, ybar = x.mean(), y.mean()
        s1, s2 = x.std(ddof=1), y.std(ddof=1)
        diff_hat = xbar - ybar
        centers[i] = diff_hat

        if method == "welch":
            se = np.sqrt(s1**2 / n1 + s2**2 / n2)
            num = (s1**2 / n1 + s2**2 / n2) ** 2
            den = (s1**2 / n1) ** 2 / (n1 - 1) + (s2**2 / n2) ** 2 / (n2 - 1)
            df = num / den
            crit = t.ppf(1 - alpha / 2.0, df=df)
        elif method == "pooled":
            df = n1 + n2 - 2
            sp2 = ((n1 - 1) * s1**2 + (n2 - 1) * s2**2) / df
            se = np.sqrt(sp2 * (1.0 / n1 + 1.0 / n2))
            crit = t.ppf(1 - alpha / 2.0, df=df)
        elif method == "z_known":
            se = np.sqrt(sigma1**2 / n1 + sigma2**2 / n2)
            crit = norm.ppf(1 - alpha / 2.0)
        else:  # z_plugin
            se = np.sqrt(s1**2 / n1 + s2**2 / n2)
            crit = norm.ppf(1 - alpha / 2.0)

        lowers[i] = diff_hat - crit * se
        uppers[i] = diff_hat + crit * se

    covered = (lowers <= delta_true) & (delta_true <= uppers)
    n_fail = int((~covered).sum())
    coverage_pct = 100.0 * covered.mean()

    fig, ax = plt.subplots(figsize=(12, 12))
    for i in range(n_simulations):
        color = "k" if covered[i] else "r"
        ax.plot([lowers[i], uppers[i]], [i, i], lw=2, color=color)
        ax.plot(centers[i], i, marker="o", ms=3, color=color)
    ax.axvline(delta_true, linestyle="--", linewidth=1.5)
    ax.set_title(
        f"{n_simulations} Two-Sample Mean CIs ({method}) | n1={n1}, n2={n2}, "
        f"CL={int((1 - alpha) * 100)}% | Fail={n_fail} (Coverage ≈ {coverage_pct:.1f}%)")
    ax.set_yticks([])
    for sp in ["left", "right", "top"]:
        ax.spines[sp].set_visible(False)
    ax.set_xlabel("Δ = μ₁ − μ₂")
    plt.tight_layout()
    plt.show()


if __name__ == "__main__":
    main()

실패가 셋이다. 참값 \(0.9499\) 에서 기대 실패가 다섯인데 셋이 나왔고, \(100\)회의 표준오차 \(0.0218\) 안이다.

폭이 \(1.559\) 에서 \(4.024\) 까지 \(2.58\)배 널뛴다. 일표본 그림보다 심한데, 흔들리는 양이 \(s_1\), \(s_2\) 둘인 데다 임계값까지 흔들리기 때문이다. 새터스웨이트 자유도가 표본마다 달라지므로 \(t_{0.025,\,\nu}\) 가 상수가 아니다. 이 \(100\)개에서 \(\nu\) 가 \(9.49\) 에서 \(19.99\) 까지 오가고, 그에 따라 임계값이 \(2.0861\) 에서 \(2.2446\) 까지 움직인다.

그러나 이 그림은 (1)(2)의 비교에 쓸 수 없다. 그려진 것은 웰치 하나뿐이고, 실패 개수 \(3\) 과 참값 \(0.9499\) 의 거리는 되풀이 \(100\)회의 잡음 안에 있다. 그림은 구간이 어떻게 생겼는지를 보여 주지 포함률이 얼마인지를 말해 주지 않는다.

\(\mathrm{SE} = \sqrt{15^2/36 + 20^2/49} = \sqrt{6.25 + 8.16} = \sqrt{14.41} \approx 3.80\).

\(\mathrm{SE} = \sqrt{64/30 + 100/35} \approx 2.233\). \(n\)이 크므로 \(z\)를 쓴다: 오차한계 = \(1.96 \cdot 2.233 \approx 4.38\).

신뢰구간: \((5 - 4.38, 5 + 4.38) = (0.62, 9.38)\). 0을 포함하지 않는다 — 차이가 있다는 증거이다.

\(\mathrm{SE} = \sqrt{64/30 + 100/35} \approx 2.234\).

Welch 자유도: \(\nu = (2.133 + 2.857)^2/[2.133^2/29 + 2.857^2/34] = 24.9/0.397 \approx 62.7\). \(t_{62}\)를 쓴다.

\(t_{0.975, 62} \approx 2.00\). 신뢰구간: \((78 - 82) \pm 2.00 \cdot 2.234 = -4 \pm 4.47 = (-8.47, 0.47)\).

0을 포함한다 — 5% 수준에서 유의한 차이를 결론지을 수 없다. 대부분 음수여서 방법 B가 더 나을 수 있음을 시사하지만 증거가 결정적이지 않다.

합동 \(t\)-검정: \(\sigma_1 = \sigma_2\)를 가정한다. 두 표본을 합쳐 공통 \(\sigma\)를 추정한다. 자유도 \(= n_1 + n_2 - 2\).

Welch \(t\)-검정: \(\sigma_1 \ne \sigma_2\)를 허용한다. 각각의 분산을 따로 쓴다. Welch-Satterthwaite 공식으로 근사 자유도를 구한다.

합동을 쓸 때: (예컨대 실험설계에 의해) 분산이 같다고 알려져 있을 때. 효율 이득이 크지 않다.

Welch를 쓸 때: 기본. 분산이 같을 필요가 없고, 실제로 분산이 같을 때도 합동에 거의 맞먹는 효율을 낸다.

현대의 권고: 등분산이 구조적으로 보장되지 않는 한 항상 Welch를 쓰라. R의 t.test는 기본이 Welch이다. 분산이 같은데 Welch를 "잘못" 쓰는 대가는 작지만, 다른데 "잘못" 합동하는 대가는 클 수 있다.

아니다. 전후 측정은 같은 피험자를 두 번 잰 것이므로 대응되어 있다. 서로 독립이 아니다. 이표본 신뢰구간을 쓰면 이들을 독립으로 취급하여 피험자 내 상관을 무시하게 된다.

올바른 접근: 피험자마다 차이 \(D_i = \mathrm{post}_i - \mathrm{pre}_i\)를 계산하고, \(\bar D, s_D\)와 자유도 \(n - 1\)로 \(D\)에 대한 일표본 신뢰구간을 구한다.

왜 중요한가: 전후가 (흔히 그렇듯) 양의 상관을 가지면 \(\mathrm{Var}(D) < \mathrm{Var}(\mathrm{pre}) + \mathrm{Var}(\mathrm{post})\)이다. 대응 분석은 표준오차가 작고 신뢰구간이 좁으며 검정력이 크다.

이표본 신뢰구간은 대응에 관한 정보를 잃고 필요 이상으로 넓은 구간을 준다. 분석은 항상 설계에 맞추어야 한다.

연습문제 3에서 합동 표준편차 \(= \sqrt{((29 \cdot 64) + (34 \cdot 100))/63} = \sqrt{(1856 + 3400)/63} = \sqrt{83.4} \approx 9.13\).

Cohen의 \(d = (78 - 82)/9.13 \approx -0.44\).

해석:

• \(|d| = 0.2\): 작은 효과.
• \(|d| = 0.5\): 중간 효과.
• \(|d| = 0.8\): 큰 효과.

\(d = -0.44\)는 "중간에 조금 못 미치는" 효과이다 — 방법 B의 평균이 A의 평균보다 약 0.44 표준편차 위에 있다. 신뢰구간은 유의성에 살짝 못 미쳤지만 효과크기는 무시할 수준이 아니다.

항상 둘 다 보고하라:

• 통계적 유의성 (신뢰구간이 0을 배제, \(p < \alpha\)).
• 실질적 유의성 (효과크기가 의미 있을 만큼 큰가).

효과크기는 무차원이어서 연구와 분야를 넘어 비교할 수 있으므로 메타분석에 유용하다. 유의성 검정은 표본크기에 좌우되지만 효과크기는 그렇지 않다.

공식.

import numpy as np

def satt(n1, n2, s1, s2):
    a, b = s1**2 / n1, s2**2 / n2
    return (a + b)**2 / (a**2 / (n1 - 1) + b**2 / (n2 - 1))

print(f"{'n1':>4s} {'n2':>4s} {'s1':>5s} {'s2':>5s} {'ν':>8s} "
      f"{'하한':>6s} {'상한':>6s}")
for n1, n2, s1, s2 in [(10, 10, 1, 1), (10, 10, 1, 5),
                       (10, 50, 1, 1), (10, 50, 5, 1),
                       (10, 50, 1, 5), (5, 100, 1, 10)]:
    v = satt(n1, n2, s1, s2)
    print(f"{n1:4d} {n2:4d} {s1:5.1f} {s2:5.1f} {v:8.2f} "
          f"{min(n1, n2) - 1:6d} {n1 + n2 - 2:6d}")

  n1   n2    s1    s2        ν     하한     상한
  10   10   1.0   1.0    18.00      9     18
  10   10   1.0   5.0     9.72      9     18
  10   50   1.0   1.0    12.87      9     58
  10   50   5.0   1.0     9.14      9     58
  10   50   1.0   5.0    57.94      9     58
   5  100   1.0  10.0    71.64      4    103

\(\nu\)는 항상 \(\min(n_1,n_2)-1\)과 \(n_1+n_2-2\) 사이에 있다. 이것이 정리로 증명되어 있다.

두 극단.

1. 한쪽 항이 지배하면 \(\nu\)가 그 집단의 자유도에 가까워진다. \((10,50,5,1)\)에서 \(s_1^2/n_1=2.5\)가 \(s_2^2/n_2=0.02\)를 압도해 \(\nu=9.14\approx n_1-1=9\)다. 작은 표본에 큰 분산이면 그 집단의 자유도만 쓰는 셈이다.

2. 다른 쪽이 지배하면 큰 집단의 자유도 쪽으로 간다. \((5,100,1,10)\)에서 \(\nu=71.6\)으로 \(n_2-1=99\)에 가까워진다.

3. 두 항이 비슷하면 합쳐진 자유도에 가까워진다. \((10,10,1,1)\)에서 \(\nu=18=n_1+n_2-2\)로 최댓값이다.

실무적 의미.

• 표본을 어디에 투자할지 알려 준다. 분산이 큰 집단의 표본을 늘려야 \(\nu\)가 오르고 구간이 좁아진다. 분산이 작은 집단을 늘리는 것은 거의 소용이 없다.

• \(\nu\)가 작으면 결과를 조심한다. \(\nu<10\)이면 \(t\) 임계값이 2.2를 넘어 구간이 크게 넓어진다. 이런 상황이면 표본 설계를 재검토할 필요가 있다.

• \(\nu\)는 확률변수다. \(s_1\), \(s_2\)에 의존하므로 표본마다 달라진다. 보고할 때 소수점 한 자리까지 적는 것이 관례다(정수로 반올림하면 정보가 사라진다).

최적 배분과의 연결. 앞서 본 대로 분산이 다르면 \(n_2/n_1=\sigma_2/\sigma_1\)이 최적이다. 그때 \(s_1^2/n_1\approx s_2^2/n_2\)가 되어 \(\nu\)가 최대에 가까워진다. 최적 배분이 자유도도 최대화한다.

원리. 결과와 상관된 공변량 \(X\)를 모형에 넣으면 잔차분산이 \(\sigma^2(1-R^2)\)로 줄어든다.

import numpy as np
from scipy import stats

rng = np.random.default_rng(55)
n, M, tau = 40, 5_000, 1.0                  # 집단당 40명
print(f"{'R²':>6s} {'단순 폭':>9s} {'보정 폭':>9s} {'비':>7s} "
      f"{'유효 n 배율':>12s}")
for R2 in [0.0, 0.1, 0.3, 0.5, 0.7, 0.9]:
    w_raw = w_adj = 0.0
    for _ in range(M // 10):
        g = np.repeat([0, 1], n)
        x = rng.normal(0, 1, 2 * n)
        y = tau * g + np.sqrt(R2) * x + np.sqrt(1 - R2) * rng.normal(0, 1, 2 * n)

        a, b = y[g == 0], y[g == 1]
        sp = np.sqrt((a.var(ddof=1) + b.var(ddof=1)) / 2)
        w_raw += 2 * stats.t.ppf(0.975, 2 * n - 2) * sp * np.sqrt(2 / n)

        X = np.column_stack([np.ones(2 * n), g, x])
        beta, *_ = np.linalg.lstsq(X, y, rcond=None)
        r = y - X @ beta
        df = 2 * n - 3
        se = np.sqrt((r @ r / df) * np.linalg.inv(X.T @ X)[1, 1])
        w_adj += 2 * stats.t.ppf(0.975, df) * se
    w_raw /= M // 10
    w_adj /= M // 10
    print(f"{R2:6.1f} {w_raw:9.4f} {w_adj:9.4f} {w_adj / w_raw:7.4f} "
          f"{(w_raw / w_adj)**2:12.2f}")

    R²      단순 폭      보정 폭       비      유효 n 배율
   0.0    0.8916    0.8978  1.0069         0.99
   0.1    0.8823    0.8411  0.9533         1.10
   0.3    0.8904    0.7456  0.8373         1.43
   0.5    0.8886    0.6293  0.7082         1.99
   0.7    0.8896    0.4906  0.5514         3.29
   0.9    0.8861    0.2848  0.3214         9.68

\(R^2=0.5\)면 폭이 29% 줄고, 표본을 두 배로 늘린 것과 같다. \(R^2=0.9\)면 표본 10배에 가까운 효과다.

\(R^2=0\)이면 손해다. 폭이 0.7% 늘어난다. 자유도를 하나 잃기 때문이다. 무관한 공변량을 넣는 비용은 매우 작지만 이득도 없다.

언제 쓸 가치가 있는가.

무작위 실험에서의 주의점.

1. 공변량은 배정 전에 측정된 것이어야 한다. 사후 변수를 넣으면 효과의 일부를 지운다.
2. 공변량 목록을 사전에 고정한다. 자료를 보고 고르면 선택 편향이 생긴다.
3. 무작위 배정이면 편향은 없다. 공변량을 넣든 안 넣든 추정량은 불편이며, 정밀도만 달라진다. 이 점에서 관측연구와 근본적으로 다르다.

층화 대 보정. 층화(블록화)는 설계 단계에서, 보정은 분석 단계에서 같은 일을 한다.

• 층화는 균형을 보장하고 해석이 쉽다. 다만 변수가 많으면 층이 폭발한다.
• 보정은 연속변수를 그대로 쓸 수 있고 여러 변수를 다루기 쉽다.
• 층화했다면 분석에서도 층을 반영해야 한다. 안 하면 보수적이 된다.

평균 차이의 한계. 두 분포가 같은 평균을 가지면서 전혀 다를 수 있다. 처리가 일부에게만 효과가 있거나 분산을 바꾸는 경우 평균 하나로는 잡히지 않는다.

방법 1 — 분위수 차이. 여러 분위수에서 차이를 재고 각각 구간을 만든다.

import numpy as np

rng = np.random.default_rng(19)
x = rng.normal(10, 2, 60)                   # 대조군
y = rng.normal(10, 4, 60)                   # 처리군: 평균 같고 분산만 큼
B = 20_000
bx = x[rng.integers(0, len(x), (B, len(x)))]
by = y[rng.integers(0, len(y), (B, len(y)))]

print(f"{'분위수':>7s} {'차이':>8s} {'95% 구간':>20s}")
for q in [10, 25, 50, 75, 90]:
    d = np.percentile(y, q) - np.percentile(x, q)
    bd = np.percentile(by, q, axis=1) - np.percentile(bx, q, axis=1)
    lo, hi = np.percentile(bd, [2.5, 97.5])
    print(f"{q:6d}% {d:8.3f}   ({lo:7.3f}, {hi:7.3f})")

d0 = y.mean() - x.mean()
bd0 = by.mean(1) - bx.mean(1)
print(f"\n평균 차이 {d0:.3f}   "
      f"({np.percentile(bd0, 2.5):.3f}, {np.percentile(bd0, 97.5):.3f})")

    분위수       차이               95% 구간
    10%   -3.179   ( -5.003,  -1.031)
    25%   -1.515   ( -2.904,  -0.058)
    50%    0.978   ( -1.710,   2.846)
    75%    2.723   (  1.396,   3.443)
    90%    3.411   (  1.747,   5.185)

평균 차이 0.431   (-0.805, 1.676)

평균 차이는 0에 가깝고 유의하지 않은데, 분위수를 보면 분포가 확연히 다르다. 아래 꼬리는 낮아지고 위 꼬리는 높아졌다 — 정확히 분산이 커진 모습이다. 분위수 분석이 평균이 놓친 것을 잡아낸다.

방법 2 — 이동-척도 그림(Q-Q plot). 두 표본의 분위수를 서로 대응시켜 그린다. 직선이면 이동만 있고, 기울기가 1이 아니면 척도가 다르며, 휘어 있으면 모양이 다르다.

방법 3 — 콜모고로프-스미르노프 통계량. 두 경험분포함수의 최대 차이에 대한 구간. 분포가 같은지에 대한 전체적 판단을 준다. 다만 어디가 다른지는 말해 주지 않는다.

방법 4 — 확률적 우세. \(P(Y>X)\)를 추정한다. 만-휘트니 통계량을 \(n_1n_2\)로 나눈 것이 불편추정량이다.

from scipy import stats
u = stats.mannwhitneyu(y, x, alternative="two-sided")
print(f"P(Y > X) 추정 {u.statistic / (len(x) * len(y)):.4f}")

P(Y > X) 추정 0.5494

0.55로 0.5에 가깝다. 확률적 우세로도 두 분포의 차이가 잡히지 않는다. 이 지표는 위치 차이에는 민감하지만 척도 차이에는 무디다.

방법 5 — 분위수회귀. 공변량이 있는 경우로 확장한다. 각 분위수에서 처리효과를 추정하고 구간을 준다.

권고. 그림을 먼저 그린다. 두 집단의 밀도나 상자그림을 나란히 놓으면 어떤 요약이 필요한지 보인다. 평균 차이 하나만 보고하는 관행은 분포가 이동만 하는 경우에나 충분하다.

구간 자체는 인과적 주장을 하지 않는다. \(\bar X_1-\bar X_2\)의 구간은 두 모집단 평균의 차이를 담을 뿐이다. 그것을 "처리의 효과"로 읽으려면 추가 조건이 필요하다.

필요한 가정 — 넷.

1. 교환가능성(무교란). 측정한 공변량 \(L\)을 고정하면 처리 배정이 잠재결과와 독립이다.

$$ Y^{(a)}\perp A\mid L $$

검증할 수 없다. 무작위 배정에서는 설계로 보장되지만, 관측연구에서는 영역 지식에 기댄 가정이다.

1. 양(positivity). 모든 \(L\) 값에서 두 처리가 모두 일어날 확률이 양수다. 특정 집단이 언제나 한 처리만 받으면 그 부분에서는 비교가 불가능하다.

2. 일관성. 관측된 처리에 대응하는 잠재결과가 관측값과 같다. "처리"가 잘 정의되어야 한다 — "운동"처럼 모호한 노출은 여러 버전이 섞여 있어 위험하다.

3. 간섭 없음(SUTVA). 한 개체의 처리가 다른 개체의 결과에 영향을 주지 않는다. 감염병, 교육, 네트워크 효과에서 자주 깨진다.

무엇을 해야 하는가.

흔한 실수 — 넷.

1. 중간변수를 보정. 처리의 결과인 변수를 넣으면 효과의 일부가 지워진다. 앞서 본 과잉짝짓기와 같은 문제다.

2. 충돌자를 보정. 처리와 결과의 공통 결과를 넣으면 없던 연관이 생긴다. "모든 변수를 넣으면 안전하다"는 직관이 틀린 이유다.

3. "보정했으니 인과"라고 말하기. 보정은 관측된 교란만 처리한다. 측정되지 않은 교란은 그대로다.

4. 구간의 폭을 불확실성 전부로 착각. 구간은 표집 불확실성만 담는다. 인과 가정의 불확실성은 그 안에 없다. 실제 불확실성은 훨씬 크다.

정직한 보고의 예.

보정 후 두 군의 평균 차이는 2.3단위였다(95% 신뢰구간 0.8~3.8). 이를 인과효과로 해석하려면 나이·성별·기저질환으로 조정한 뒤 잔여 교란이 없다는 가정이 필요하다. 로젠바움 민감도 분석에서, 측정되지 않은 요인이 노출 오즈를 1.4배 이상 바꾸면 결론이 유지되지 않는다. 따라서 이 결과는 인과관계를 시사하되 확정하지 못한다.

한 문장 요약. 신뢰구간은 "얼마나 다른가"에 답하고, "왜 다른가"에는 답하지 않는다. 후자는 설계와 가정의 몫이다.