p̂₁ - p̂₂의 표본분포 (작은 기대도수)¶

(1) 붕괴 조건은 간단하다. 왈드의 표준오차 \(\widehat{\text{SE}} = \sqrt{\hat p_1\hat q_1/n_1 + \hat p_2\hat q_2/n_2}\)는 두 항이 모두 \(0\)일 때만 \(0\)이 된다. \(\hat p_i\hat q_i = 0\)은 \(X_i \in \{0, n_i\}\)를 뜻하므로, 붕괴하는 칸은 네 모퉁이 \((0,0)\), \((0,n)\), \((n,0)\), \((n,n)\)뿐이다. \(p_1 = p_2 = p\)이고 \(n_1 = n_2 = n\)이면

이고, \(p\)가 작으면 첫 항만 남는다. \(n = 30\)에서

\(p = 0.02\)에서 열 번에 세 번꼴이다. 드문 사고가 아니다.

\((0,0)\)에서 아그레스티-카포의 반폭은 손으로 나온다. \(X_1 = X_2 = 0\)이면

이고 두 항이 같으므로

이다. \(\tilde p_i\)는 \(0\)이 될 수 없으므로 붕괴가 원천적으로 불가능하다. 이것이 "각 표본에 성공 하나와 실패 하나를 더한다"가 하는 일의 전부다.

(2) 전수 계산.

import numpy as np
from scipy import stats

n1 = n2 = 30
z975 = stats.norm.ppf(0.975)

k1 = np.arange(n1 + 1)
k2 = np.arange(n2 + 1)
X1, X2 = np.meshgrid(k1, k2, indexing="ij")
ph1, ph2 = X1 / n1, X2 / n2
d = ph1 - ph2

def weight(pa, pb):
    """(X1, X2) = (i, j) 가 나올 확률을 담은 31 x 31 행렬."""
    return np.outer(stats.binom.pmf(k1, n1, pa), stats.binom.pmf(k2, n2, pb))

# --- 왈드 ---
se_w = np.sqrt(ph1 * (1 - ph1) / n1 + ph2 * (1 - ph2) / n2)
wald_lo, wald_hi = d - z975 * se_w, d + z975 * se_w

# --- 아그레스티-카포: 각 칸에 1 씩 더한 뒤 왈드 공식을 그대로 ---
t1, t2 = (X1 + 1) / (n1 + 2), (X2 + 1) / (n2 + 2)
se_ac = np.sqrt(t1 * (1 - t1) / (n1 + 2) + t2 * (1 - t2) / (n2 + 2))
ac_lo, ac_hi = (t1 - t2) - z975 * se_ac, (t1 - t2) + z975 * se_ac

# --- 윌슨 기반(뉴컴 혼합 점수 구간) ---
def wilson(x, n, z=z975):
    m = n + z ** 2
    center = (x + z ** 2 / 2) / m
    half = z / m * np.sqrt(x * (n - x) / n + z ** 2 / 4)
    return center - half, center + half

l1, u1 = wilson(X1, n1)
l2, u2 = wilson(X2, n2)
nc_lo = d - np.sqrt((ph1 - l1) ** 2 + (u2 - ph2) ** 2)
nc_hi = d + np.sqrt((u1 - ph1) ** 2 + (ph2 - l2) ** 2)

print(f"P(X = 0) = 0.95^30 = {0.95 ** 30:.4f}")
print(f"두 표본 모두 0    = {0.95 ** 60:.4f}")
print()
print("p1 = p2 = p 일 때 왈드 구간의 폭이 0 이 될 확률")
for p in (0.02, 0.05, 0.10, 0.20):
    print(f"  p = {p:.2f}:  {weight(p, p)[se_w == 0].sum():.4f}")
print()
print("관측값별 구간 (n1 = n2 = 30)")
for a, b in [(0, 0), (1, 0), (2, 1)]:
    print(f"  (x1, x2) = ({a}, {b}):"
          f"  왈드 [{wald_lo[a,b]:+.4f}, {wald_hi[a,b]:+.4f}]"
          f"  AC [{ac_lo[a,b]:+.4f}, {ac_hi[a,b]:+.4f}]"
          f"  윌슨 [{nc_lo[a,b]:+.4f}, {nc_hi[a,b]:+.4f}]")

출력:

P(X = 0) = 0.95^30 = 0.2146
두 표본 모두 0    = 0.0461

p1 = p2 = p 일 때 왈드 구간의 폭이 0 이 될 확률
  p = 0.02:  0.2976
  p = 0.05:  0.0461
  p = 0.10:  0.0018
  p = 0.20:  0.0000

관측값별 구간 (n1 = n2 = 30)
  (x1, x2) = (0, 0):  왈드 [+0.0000, +0.0000]  AC [-0.0853, +0.0853]  윌슨 [-0.1135, +0.1135]
  (x1, x2) = (1, 0):  왈드 [-0.0309, +0.0976]  AC [-0.0720, +0.1345]  윌슨 [-0.0834, +0.1667]
  (x1, x2) = (2, 1):  왈드 [-0.0766, +0.1433]  AC [-0.1000, +0.1625]  윌슨 [-0.1085, +0.1824]

\(p = 0.02\)이면 붕괴 확률이 30%에 이른다. \(p = 0.05\)에서도 4.6%다. 드문 사고가 아니라 흔한 결과다.

\((0,0)\) 칸을 보면 세 방법의 성격이 한눈에 드러난다. 왈드는 폭 0, 아그레스티-카포는 \(\pm 0.085\), 윌슨은 \(\pm 0.114\)다. 뒤의 두 값은 "각 군 30명으로 알 수 있는 정도"와 잘 맞는다. 한 표본에서 \(0/30\)이면 참 비율의 95% 상한이 대략 \(3/30 = 0.10\)이라는 3의 법칙과도 일관된다.

(1)의 두 값이 그대로 나왔다. 붕괴 확률이 \(p = 0.02\)에서 \(0.2976\), \(p = 0.05\)에서 \(0.0461\)이고, \((0,0)\)에서 아그레스티-카포 구간이 \([-0.0853,\ +0.0853]\)이다. 손으로 낸 \(0.0853\)과 같다.

세 구간이 같은 자료에 전혀 다른 답을 준다.

왈드의 \(0\)이 뜻하는 바를 말로 옮기면 터무니없다. 각 군 \(30\)명에서 사건을 하나도 보지 못하고 "두 비율의 차가 정확히 \(0\)이라고 \(95\%\) 확신한다"고 말하는 셈이다. 자료가 가장 적게 말해 주는 순간에 구간이 가장 확신하는 구조이며, \(\hat p(1-\hat p)\)를 표준오차의 재료로 쓴 데서 오는 필연적 결과다.

뒤의 두 값은 3의 법칙과 맞는다. 한 표본에서 \(0/n\)을 관측했을 때 참 비율의 \(95\%\) 상한이 대략 \(3/n\)이라는 어림이 있다. \(n = 30\)이면 \(3/30 = 0.10\)이다. 실제로 윌슨 구간을 한 표본에 적용하면 \(X = 0\), \(n = 30\)에서 \([0,\ 0.1135]\)가 나오고, 어림 \(0.10\)과 같은 자릿수다. 그러므로 각 비율이 \(0\)에서 \(0.11\) 사이 어디든 될 수 있고, 따라서 차도 \(\pm 0.11\) 안 어디든 될 수 있다. 윌슨 기반 구간의 \(\pm 0.1135\)가 바로 그 말이다. 왈드만 이 상식에서 벗어나 있다.

\((1,0)\)과 \((2,1)\) 줄을 보면 어긋남이 작은 관측값 전반에 걸쳐 있음도 보인다. \((1,0)\)에서 왈드의 폭이 \(0.1285\)인데 아그레스티-카포는 \(0.2065\), 윌슨은 \(0.2501\)이다. 붕괴는 가장 눈에 띄는 증상일 뿐 병 자체는 아니다.

(1) 같은 칸이 참값에 따라 공과 과를 바꾼다. 보기 1에서 \(X_1 = X_2 = 0\)이면 왈드 구간이 한 점 \([0,0]\)이었다. 이 칸의 운명은 담아야 할 참값 \(\delta = p_1-p_2\)가 무엇이냐로 갈린다.

• \(\delta = 0\)이면 \([0,0]\)이 \(0\)을 담으므로 성공이다.
• \(\delta \ne 0\)이면 폭이 \(0\)인 구간이 \(\delta\)를 담을 길이 없으므로 반드시 실패다.

\((0.05,\ 0.05)\)는 앞의 경우다. 붕괴 확률 \(0.95^{60} = 0.0461\)이 통째로 성공으로 세어지므로 포함률이 명목 위로 밀려 올라간다. 왈드가 멀쩡해 보이는 것은 고장이 상을 받고 있기 때문이다.

\((0.05,\ 0.02)\)는 뒤의 경우다. 붕괴 확률이

이고 담아야 할 값이 \(0.03\)이므로 이 \(0.1171\)이 고스란히 실패다. 다른 칸들이 완벽하더라도 포함률은 \(1 - 0.1171 = 0.8829\)를 넘을 수 없고, 나머지 칸의 평범한 손실까지 더하면 \(0.85\) 언저리로 내려갈 것이라 예측된다.

아그레스티-카포는 이 칸에서 \([-0.0853,\ +0.0853]\)을 주므로 \(0.03\)을 담는다. 붕괴가 없으니 손실도 없고, 오히려 작은 기대도수에서 구간이 넉넉해 명목을 꽤 넘을 것이다.

(2) 전수 계산.

print("명목 95% 구간의 정확 포함률")
print(f"{'':>14}{'왈드':>9}{'AC':>9}{'윌슨':>9}")
for pa, pb in [(0.05, 0.05), (0.05, 0.02), (0.10, 0.05), (0.15, 0.05),
               (0.20, 0.20), (0.30, 0.10), (0.50, 0.20)]:
    W = weight(pa, pb)
    t = pa - pb          # 담아야 할 참값
    cw = W[(wald_lo <= t) & (t <= wald_hi)].sum()
    ca = W[(ac_lo <= t) & (t <= ac_hi)].sum()
    cn = W[(nc_lo <= t) & (t <= nc_hi)].sum()
    print(f"  ({pa:.2f}, {pb:.2f}) {cw:9.4f}{ca:9.4f}{cn:9.4f}")

출력:

명목 95% 구간의 정확 포함률
                     왈드       AC       윌슨
  (0.05, 0.05)    0.9716   0.9929   0.9929
  (0.05, 0.02)    0.8482   0.9976   0.9975
  (0.10, 0.05)    0.9298   0.9833   0.9807
  (0.15, 0.05)    0.9260   0.9707   0.9748
  (0.20, 0.20)    0.9396   0.9547   0.9547
  (0.30, 0.10)    0.9374   0.9593   0.9552
  (0.50, 0.20)    0.9412   0.9515   0.9484

가장 극적인 줄은 \((p_1, p_2) = (0.05,\ 0.02)\)다. 왈드의 포함률이 \(0.8482\)로 명목 \(0.95\)에서 6.5퍼센트포인트나 내려간다. 같은 자리에서 아그레스티-카포는 \(0.9976\)이다.

왜 이렇게 벌어지는가. 이 설정에서 두 표본이 모두 0일 확률은

이다. 이때 왈드 구간은 \([0,0]\)인데 담아야 할 참값은 \(0.03\)이므로 반드시 실패한다. 붕괴 확률 \(0.1171\)이 그대로 포함률 손실로 돌아온다.

반면 첫 줄 \((0.05,\ 0.05)\)에서는 왈드가 \(0.9716\)으로 오히려 명목을 넘는다. \([0,0]\)이 참값 \(0\)을 담기 때문이다. 같은 붕괴가 참값이 0이면 성공으로, 0이 아니면 실패로 세어진다. 왈드의 결함을 \(p_1 = p_2\)만 보고 판단하면 안 되는 이유다.

아그레스티-카포와 윌슨은 \(p_1 = p_2\)인 두 줄에서 값이 완전히 같다. 우연이 아니다. 목표값이 \(0\)일 때 두 방법의 채택 영역이 \(961\)개 칸 전부에서 일치한다(왈드는 \(353\)칸, 나머지 둘은 같은 \(369\)칸을 채택한다). \(p_1 \ne p_2\)인 줄에서는 갈라진다.

예측이 둘 다 맞는다. \((0.05,\ 0.05)\)에서 왈드가 \(0.9716\)으로 명목을 넘고, \((0.05,\ 0.02)\)에서 \(0.8482\)로 \(6.5\)퍼센트포인트 내려간다. 뒤 칸의 예측 상한 \(0.8829\) 아래이며, 붕괴 \(0.1171\)에 다른 칸의 손실 \(0.0347\)이 더해져 실패 확률이 \(0.1518\)이 된 것이다(보기 5에서 칸별로 센다).

\(0.1171\) 하나가 \(0.1518\)의 \(77\%\)다. 왈드의 과소포함은 여러 칸에 고루 흩어진 것이 아니라 한 칸에서 온다. 그래서 그 한 칸만 고치면 거의 다 고쳐진다.

\((0.10,\ 0.05)\)에서는 붕괴 확률이 \(0.0091\)에 그치는데도 왈드가 \(0.9298\)이다. 붕괴를 모두 빼도 \(0.991\)까지는 갈 수 있어야 하는데 그보다 \(6\)퍼센트포인트 낮다. 붕괴가 왈드의 유일한 병은 아니다. \(\hat p_i\)가 \(0\)은 아니지만 작을 때 \(\widehat{\text{SE}}\)가 참값보다 작게 나오는 경향이 남아 있고, 그 몫이 기대도수가 커질 때까지 따라다닌다. 마지막 줄 \((0.20,0.20)\)의 \(0.9396\)이 그 잔여분이다.

아그레스티-카포와 윌슨이 \(p_1 = p_2\)인 두 줄에서 값이 완전히 같다. 우연이 아니다. 목표값이 \(0\)일 때 두 방법의 채택 영역이 \(961\)칸 전부에서 일치한다(왈드는 \(353\)칸, 나머지 둘은 같은 \(369\)칸을 채택한다). 두 방법은 서로 다른 공식에서 출발했지만 \(\delta = 0\)을 담느냐 마느냐라는 물음에는 똑같이 답한다. \(p_1 \ne p_2\)인 줄에서는 \(0.9833\) 대 \(0.9807\)처럼 갈라진다.

기대도수가 커지면 셋이 모여든다. 마지막 줄 \((0.50,\ 0.20)\)에서 \(0.9412 / 0.9515 / 0.9484\)로 셋 다 \(0.95\) 둘레다. 이 쪽에서 재는 것은 기대도수가 작을 때의 차이이고, 그 차이가 가장 크게 벌어지는 자리가 \((0.05,\ 0.02)\)의 \(0.8482\) 대 \(0.9976\)이다.

(1) 기각역이 격자 위에 있다. \(961\)칸 각각에서 합동 \(z\) 통계량을 계산해 \(|z| > 1.96\)인 칸을 모으면 \(600\)칸이 기각역이다. 수로는 많아 보이지만, 문제는 그 칸들이 어디 있는가다. \(x_2\)를 고정하고 기각되는 가장 작은 \(x_1\)을 찾으면 이렇다.

두 표본이 모두 거의 비어 있으면 아무것도 기각할 수 없다. \(x_2 = 0\)일 때 \(x_1\)이 \(4\)는 되어야 하고, \(x_2 = 1\)이면 \(6\)이 필요하다. 그런데 \(p = 0.02\)에서는 \(X_1 + X_2 \sim \text{Bin}(60, 0.02)\)의 평균이 \(1.2\)라 그만한 격차 자체가 나오지 않는다.

크기의 상한은 바로 나온다. 기각하려면 적어도 \(X_1 + X_2 \ge 4\)여야 하므로

이고, 실제로는 \(X_1+X_2 = 4\) 가운데 \((4,0)\)과 \((0,4)\) 두 칸만 기각되므로 그보다 훨씬 작다. 명목 \(0.05\)는 애초에 도달할 수 없는 값이다.

반대쪽도 미리 보인다. \(p\)가 커지면 도달 가능한 칸이 많아지고 격자가 촘촘해져 \(0.05\)에 가까워지는데, 격자가 성긴 중간 영역에서는 넘어설 수도 있다. 기각역의 경계가 \(0.05\)에 딱 맞지 않고 칸 단위로만 움직이기 때문이다.

(2) 전수 계산.

ppool = (X1 + X2) / (n1 + n2)
se0 = np.sqrt(ppool * (1 - ppool) * (1 / n1 + 1 / n2))
zstat = np.divide(d, se0, out=np.zeros_like(d), where=se0 > 0)
rej_z = np.abs(zstat) > z975

# 961 개 칸마다 피셔 정확검정의 p-값을 미리 계산해 둔다.
fisher = np.empty_like(d)
for a in k1:
    for b in k2:
        fisher[a, b] = stats.fisher_exact([[a, n1 - a], [b, n2 - b]])[1]
rej_f = fisher <= 0.05

print("명목 유의수준 5% 의 실제 1종오류율")
for p in (0.02, 0.05, 0.10, 0.20, 0.30, 0.50):
    W = weight(p, p)
    print(f"  p1 = p2 = {p:.2f}:  z 검정 {W[rej_z].sum():.4f}   피셔 {W[rej_f].sum():.4f}")

출력:

명목 유의수준 5% 의 실제 1종오류율
  p1 = p2 = 0.02:  z 검정 0.0032   피셔 0.0000
  p1 = p2 = 0.05:  z 검정 0.0284   피셔 0.0015
  p1 = p2 = 0.10:  z 검정 0.0544   피셔 0.0094
  p1 = p2 = 0.20:  z 검정 0.0490   피셔 0.0222
  p1 = p2 = 0.30:  z 검정 0.0487   피셔 0.0261
  p1 = p2 = 0.50:  z 검정 0.0519   피셔 0.0274

\(z\) 검정의 실제 크기가 \(p\)에 따라 \(0.0032\)에서 \(0.0554\)까지 움직인다(\(p\)를 \(0.02\)에서 \(0.50\)까지 훑었을 때의 최솟값과 최댓값이며, 최댓값은 \(p = 0.116\) 부근에서 나온다). \(p\)가 아주 작으면 극단적으로 보수적이고, \(p = 0.1\) 근처에서는 명목 \(0.05\)를 넘는다. 한쪽에서는 검정력을 버리고 다른 쪽에서는 거짓 유의성을 낸다.

피셔 정확검정은 언제나 명목보다 작다. 최댓값이 \(0.0274\)로 명목의 절반을 겨우 넘는다. "정확검정"이라는 이름은 \(p\)-값 계산이 정확하다는 뜻이지 크기가 \(0.05\)에 맞는다는 뜻이 아니다. 주변합을 모두 고정하는 조건화와 이산성 때문에 달성 가능한 \(p\)-값이 띄엄띄엄하고, 그래서 실제 크기가 명목에 결코 도달하지 못한다. 안전하지만 검정력이 낮다.

\(z\) 검정의 크기가 \(p\)에 따라 널뛴다. \(p = 0.02\)에서 \(0.0032\), \(p = 0.10\)에서 \(0.0544\)다. \(p\)를 \(0.02\)에서 \(0.50\)까지 \(0.0005\) 간격으로 훑으면 최솟값이 \(0.0032\)(\(p = 0.02\)), 최댓값이 \(0.0554\)(\(p = 0.116\))다. 한쪽에서는 검정력을 버리고 다른 쪽에서는 거짓 유의성을 낸다.

\(p = 0.02\)의 \(0.0032\)는 (1)의 상한 \(0.0322\)의 10분의 1이다. 기각역에서 확률이 가장 큰 두 칸 \((4,0)\)과 \((0,4)\)가 각각 \(0.0014\)이고 둘을 더하면 \(0.0028\)이니, 크기의 \(88\%\)가 이 두 칸에서 나온다. 나머지는 사실상 도달 불가능한 칸들이다.

피셔 정확검정은 언제나 명목보다 작다. 최댓값이 \(0.0274\)로 명목의 절반을 겨우 넘는다. "정확"이라는 이름은 \(p\)-값 계산에 근사가 없다는 뜻이지 크기가 \(0.05\)에 맞는다는 뜻이 아니다. 까닭은 둘이다.

• 주변합을 모두 고정한다. \(H_0\) 아래에서 \(X_1+X_2\)를 조건으로 잡아 초기하분포를 쓰므로, 검정이 보는 표본공간이 전체 \(961\)칸이 아니라 대각선 하나로 줄어든다. 그 안에서 달성 가능한 \(p\)-값이 더 성겨진다.
• 이산성. 달성 가능한 \(p\)-값이 띄엄띄엄하므로 "\(p \le 0.05\)"가 걸러 내는 확률은 \(0.05\)에 못 미치는 어떤 값에서 멈춘다. 넘어서려면 다음 칸을 통째로 들여야 하는데 그러면 \(0.05\)를 훌쩍 넘는다.

그래서 셋 가운데 고르는 일이 맞바꿈이 된다. \(z\) 검정은 평균적으로 \(0.05\) 둘레이지만 어떤 \(p\)에서는 넘는다. 피셔는 절대 넘지 않지만 검정력을 크게 잃는다. 중간을 노리는 것이 중간-\(p\)값이나 조건화를 풀고 보는 비조건부 정확검정이며, 어느 쪽을 쓸지는 거짓 유의성과 놓친 발견 가운데 무엇이 더 비싼가로 정한다.

(1) 포함률은 칸의 확률을 더한 값이고, 무엇이 움직이느냐가 다르다. 포함률은

로 쓰인다. 모수를 움직이면 가중값 \(P(X_1=i)P(X_2=j)\)와 채택 집합 \(\mathcal A\)가 함께 움직이는데, 앞의 것은 \(p\)에 대해 매끄러운 함수이고 뒤의 것은 유한집합이라 계단으로만 변한다.

• 왼쪽 그림(\(p_1 = p_2 = p\)). 담아야 할 값이 \(\delta = 0\)으로 고정되어 있다. \(\mathcal A(0)\)이 바뀌지 않으므로 움직이는 것은 가중값뿐이고, 가중값은 \(p\)의 매끄러운 함수다. 곡선이 매끄러워야 한다.
• 오른쪽 그림(\(\delta\)를 움직임). \(\delta\)가 어떤 칸의 구간 끝 \(L_{ij}\)나 \(U_{ij}\)를 지나는 순간 그 칸이 \(\mathcal A\)에서 통째로 빠지거나 들어온다. 그 칸의 확률만큼 곡선이 수직으로 떨어지거나 솟는다. 톱니가 되어야 한다.

톱니의 이빨 크기도 짐작할 수 있다. 가장 큰 이빨은 확률이 가장 큰 칸이 빠질 때 생긴다. \(p_2 = 0.05\) 둘레에서 가장 무거운 칸들의 확률이 \(0.1\) 안팎이므로(보기 5에서 \((0,0)\)이 \(0.117\)), 이빨 하나가 수 퍼센트포인트에 이를 수 있다.

왼쪽 그림에서 왈드의 모양도 예측된다. 보기 1에서 붕괴 확률이 \([(1-p)^{30}+p^{30}]^2\)이었고, 붕괴 칸 \([0,0]\)은 \(\delta = 0\)을 담는다. \(p\)가 작으면 이 공짜 성공이 커서 포함률을 명목 위로 밀어 올리고, \(p\)가 커지면 그 몫이 사라져 왈드 본래의 과소포함이 드러난다. 작은 \(p\)에서 명목을 넘고 큰 \(p\)에서 명목 아래가 될 것이다.

(2) 그려서 확인한다.

포함률을 한두 점이 아니라 연속적으로 훑어 본다. 왼쪽은 \(p_1 = p_2 = p\)로 두고 \(p\)를 움직인 것이고, 오른쪽은 \(p_2 = 0.05\)로 고정한 채 담아야 할 참값 \(\delta = p_1 - p_2\)를 움직인 것이다.

import matplotlib.pyplot as plt

# 왼쪽: 참값을 0 에 고정하고 공통 비율 p 를 움직인다.
g = np.arange(0.02, 0.5001, 0.001)
cw = [weight(p, p)[(wald_lo <= 0) & (0 <= wald_hi)].sum() for p in g]
ca = [weight(p, p)[(ac_lo <= 0) & (0 <= ac_hi)].sum() for p in g]

# 오른쪽: p2 를 0.05 에 고정하고 담아야 할 참값 delta 를 움직인다.
D = np.arange(0.0, 0.4501, 0.0005)
cw2 = [weight(0.05 + t, 0.05)[(wald_lo <= t) & (t <= wald_hi)].sum() for t in D]
ca2 = [weight(0.05 + t, 0.05)[(ac_lo <= t) & (t <= ac_hi)].sum() for t in D]

fig, axes = plt.subplots(1, 2, figsize=(12, 4.2))
panels = [(g, cw, ca, r"$p$   ($p_1=p_2=p$,  target $p_1-p_2=0$)",
           r"$n_1=n_2=30$:  coverage vs. the common $p$", 1.3),
          (D, cw2, ca2, r"$\delta = p_1 - p_2$   ($p_2 = 0.05$ fixed)",
           "coverage vs. the target difference:  sawtooth", 1.0)]
for ax, (xs, a, b, xlab, title, lw) in zip(axes, panels):
    ax.plot(xs, a, color="C3", lw=lw, label="Wald")
    ax.plot(xs, b, color="C0", lw=lw, label="Agresti-Caffo")
    ax.axhline(0.95, color="k", ls="--", lw=1, label="nominal 0.95")
    ax.set_xlabel(xlab)
    ax.set_title(title)
    ax.set_ylim(0.88, 1.005)
    ax.legend(fontsize=8, loc="lower right")
axes[0].set_ylabel("exact coverage")

plt.tight_layout()
plt.show()

왼쪽 그림. \(p_1 = p_2\)인 축을 따라가면 곡선이 매끄럽다. 담아야 할 참값이 \(0\)에 고정되어 있어 칸의 포함 여부가 바뀌지 않고, 게다가 두 이항의 차는 확률이 더 많은 값에 나뉘어 퍼지므로 한 표본일 때만큼 들쭉날쭉하지 않다(연습문제 6). 왈드는 \(p < 0.083\)에서 붕괴 덕에 명목을 넘고(\(p = 0.02\)에서 \(0.9968\)), 그 뒤로는 계속 명목 아래에 머문다(최저 \(0.9396\), \(p = 0.21\) 부근). 아그레스티-카포는 작은 \(p\)에서 보수적이고 큰 \(p\)에서 \(0.95\)에 붙는다(최저 \(0.9481\)).

오른쪽 그림. 참값 \(\delta\)를 움직이면 톱니가 드러난다. 구간의 끝이 놓일 수 있는 자리가 유한하므로, \(\delta\)가 어떤 칸의 구간 경계를 지나는 순간 그 칸의 확률이 통째로 빠지고 곡선이 수직으로 떨어진다. 왈드는 \(0.8946\)에서 \(0.9716\) 사이를 오가며 평균 \(0.9351\), 아그레스티-카포는 \(0.9428\)에서 \(0.9952\) 사이를 오가며 평균 \(0.9606\)이다.

포함률이 참값의 위치에 따라 이렇게 뛴다는 것은 "명목 95%"가 하나의 수가 아니라 참값에 따라 달라지는 함수임을 뜻한다. 이산자료에서는 피할 수 없는 성질이며, 좋은 방법은 진동의 중심을 \(0.95\)에 두고 아래쪽으로 깊이 파이지 않게 하는 방법이다.

예측이 둘 다 맞았다.

왼쪽 그림은 매끄럽다. 담아야 할 값이 \(0\)에 고정되어 채택 집합이 바뀌지 않기 때문이다. 왈드는 \(p = 0.02\)에서 \(0.9968\)로 출발해 \(p < 0.083\) 구간에서 명목을 넘고, 그 뒤로는 계속 명목 아래에 머문다(최저 \(0.9396\), \(p = 0.208\) 부근). 붕괴가 공짜 성공으로 세어지는 구간과 그렇지 않은 구간이 \(p \approx 0.083\)에서 갈린다. 아그레스티-카포는 작은 \(p\)에서 보수적이고 큰 \(p\)에서 \(0.95\)에 붙으며 최저가 \(0.9481\)이다.

오른쪽 그림은 톱니다. 두 방법의 진동 폭과 평균은 이렇다.

읽을 것은 세 가지다. 첫째, 두 방법 모두 진동한다. 이산자료에서 포함률은 하나의 수가 아니라 참값에 따라 달라지는 함수이며, 이것은 방법의 흠이 아니라 자료의 성질이다. 둘째, 진동의 중심이 다르다. 왈드의 평균이 \(0.9351\)로 명목보다 \(1.5\)퍼센트포인트 낮고 아그레스티-카포는 \(0.9606\)으로 조금 높다. 셋째, 아래로 파인 깊이가 다르다. 왈드는 \(0.8946\)까지 떨어지는 반면 아그레스티-카포는 \(0.9428\)에서 멈춘다.

좋은 구간의 기준이 이 그림에서 정해진다. 진동을 없앨 수는 없으니, 진동의 중심을 \(0.95\)에 두고 아래로 깊이 파이지 않게 하는 것이 목표다. 아그레스티-카포는 평균을 \(0.96\)으로 조금 올려 깊은 골을 메우는 쪽을 택한 것이고, 그 대가로 구간이 평균적으로 약간 넓다. 명목을 지키는 일과 구간을 좁게 만드는 일이 맞바꿈이며, 기대도수가 작을 때는 앞쪽을 택하는 것이 표준이다.

(1) 세 격자는 서로 다른 것을 보인다. 왼쪽은 어디에 확률이 있는가, 가운데와 오른쪽은 어디가 성공이고 어디가 실패인가다. 포함률은 이 둘의 곱을 더한 것이므로, 실패 영역이 아무리 넓어도 거기에 확률이 없으면 손실이 없고 반대로 실패 칸 하나에 확률이 몰려 있으면 그 한 칸이 전부를 결정한다.

\((0.05, 0.02)\)에서 확률이 어디 있는지는 손으로 셀 수 있다.

이므로 \((0,0)\) 한 칸이 \(0.2146 \times 0.5455 = 0.1171\)이다. 두 번째 집단은 절반 넘는 확률이 \(X_2 = 0\)에 쏠려 있고, 그래서 \(X_2 = 0\)인 세로줄 하나가 그림 전체의 절반이 넘는다. 앞 쪽의 매끄러운 종 모양과는 전혀 다른 그림이 될 것이다.

그리고 \((0,0)\) 칸에서 왈드 구간은 \([0,0]\)인데 담아야 할 값은 \(0.03\)이므로 반드시 실패한다. 보기 2가 이미 그 대가를 \(0.95 - 0.8482 = 0.1018\)로 보였고, 붕괴 확률 \(0.1171\)이 그 손실의 뼈대다. 실패 확률 전체에서 이 한 칸이 얼마를 차지하는지 세어 둔다.

# (0.05, 0.02) 에서 담아야 할 참값은 delta = 0.03 이다.
W = weight(0.05, 0.02)
t = 0.03
ok_w = (wald_lo <= t) & (t <= wald_hi)
ok_a = (ac_lo <= t) & (t <= ac_hi)
print(f"확률이 큰 세 칸:  (0,0) {W[0,0]:.4f}   (1,0) {W[1,0]:.4f}   (2,0) {W[2,0]:.4f}")
print(f"왈드 실패 확률 합 = {W[~ok_w].sum():.4f}"
      f"   그중 (0,0) 한 칸의 몫 = {W[0,0] / W[~ok_w].sum():.1%}")
print(f"AC   실패 확률 합 = {W[~ok_a].sum():.4f}"
      f"   왈드의 {W[~ok_w].sum() / W[~ok_a].sum():.0f}분의 1")

출력:

확률이 큰 세 칸:  (0,0) 0.1171   (1,0) 0.1849   (2,0) 0.1411
왈드 실패 확률 합 = 0.1518   그중 (0,0) 한 칸의 몫 = 77.1%
AC   실패 확률 합 = 0.0024   왈드의 62분의 1

왈드의 실패 확률 \(0.1518\) 가운데 \(77\%\)가 \((0,0)\) 한 칸에서 나온다. 나머지 \(960\)칸을 다 합쳐도 \(0.0347\)뿐이다. 그림에서 보아야 할 것이 이 한 칸으로 좁혀졌다.

(2) 그림을 읽는다.

\((p_1, p_2) = (0.05,\ 0.02)\)에서 각 칸의 확률과, 그 칸에서 구간이 참값 \(0.03\)을 담는지를 나란히 본다.

왼쪽 격자를 보면 확률이 왼쪽 아래 구석에 몰려 있다. \((0,0)\) 한 칸이 \(0.117\), \((1,0)\)이 \(0.185\), \((2,0)\)이 \(0.141\)이다. 앞 페이지의 매끄러운 종 모양과는 전혀 다른 그림이다.

가운데가 왈드다. 파란 테두리를 두른 \((0,0)\) 칸이 붕괴 지점이며, 여기서만 확률 \(0.117\)이 통째로 실패로 넘어간다. 실패 확률의 합은 \(0.1518\)이고 그중 77%가 이 한 칸에서 나온다.

오른쪽이 아그레스티-카포다. \((0,0)\) 칸이 성공으로 바뀌었고, 실패 영역 전체가 확률이 거의 없는 바깥쪽으로 밀려났다. 실패 확률의 합이 \(0.0024\)로 62분의 1이다.

세 격자가 (1)의 계산을 그대로 보인다. 왼쪽에서 확률이 왼쪽 아래 구석에 몰려 있고, 가운데에서 그 구석의 \((0,0)\) 칸만 실패로 칠해져 있으며, 오른쪽에서는 그 칸이 성공으로 바뀌고 실패 영역이 확률이 거의 없는 바깥쪽으로 밀려났다. 왈드와 아그레스티-카포의 차이가 격자 한 칸의 색으로 나타난 셈이다.

이 그림이 가리는 것 하나 — 칸의 색은 확률의 크기를 말해 주지 않는다. 가운데 격자에서 실패로 칠해진 칸은 여럿이지만 그중 \((0,0)\) 하나가 실패 확률의 \(77\%\)다. 색만 세면 "실패 칸이 꽤 많다"로 읽히고, 확률을 세면 "실패는 사실상 한 칸에서 일어난다"가 된다. 왼쪽 격자를 함께 보지 않으면 가운데 격자는 오해를 부른다.

가리는 것 둘 — 이 그림은 \((p_1,p_2) = (0.05, 0.02)\) 한 점의 사진이다. 참값이 \(0\)인 \((0.05,0.05)\)로 옮기면 똑같은 \((0,0)\) 칸이 성공으로 뒤집히고 왈드의 포함률이 \(0.9716\)으로 명목을 넘는다. 같은 붕괴가 같은 자리에서 일어나는데 참값의 위치에 따라 공과 과가 뒤바뀐다. 한 장의 그림으로는 그 뒤바뀜을 볼 수 없고, 그래서 보기 4가 참값 \(\delta\)를 연속으로 훑어 톱니를 그린 것이다.

고치는 값은 한 줄이다. \(\hat p_i\)를 \(\tilde p_i = (X_i+1)/(n_i+2)\)로 바꾸면 \((0,0)\) 칸의 구간이 \([0,0]\)에서 \(\pm 0.0853\)으로 벌어지고, 실패 확률이 \(0.1518\)에서 \(0.0024\)로 62분의 1이 된다. 자료를 더 모은 것도 아니고 공식을 바꾼 것도 아니며, 각 표본에 성공 하나와 실패 하나를 보탠 것이 전부다.

한 표본에서 사건이 하나도 없을 확률은

이다. 두 표본이 독립이므로

곧 약 4.6%다. 이때 \(\hat p_1 = \hat p_2 = 0\)이므로

이고 구간은 \(0 \pm 1.96 \times 0 = [0,\ 0]\)이다.

\(p\)가 더 작으면 훨씬 자주 일어난다. \(p = 0.02\)이면 \(0.98^{60} = 0.2976\)으로 세 번에 한 번 가까이 붕괴한다.

각 표본에 성공 하나와 실패 하나를 더한다.

차의 중심은 \(\tilde p_1 - \tilde p_2 = 0\)이고 표준오차는

이므로

이다. 보기 1의 출력과 일치한다.

해석. "참 차가 대략 \(\pm 8.5\)퍼센트포인트 안에 있다"는 말이며, 각 군 30명으로 알 수 있는 정도와 잘 맞는다. 폭 0인 왈드 구간과 비교하면 차이가 분명하다.

차이를 만드는 것은 붕괴 칸이 성공으로 세어지느냐 실패로 세어지느냐다.

\((0.05, 0.05)\)에서 담아야 할 참값은 \(0\)이다. \(X_1 = X_2 = 0\)이면 구간은 \([0,0]\)이고 이 구간은 \(0\)을 담는다. 확률 \(0.0461\)이 성공으로 들어간다.

\((0.05, 0.02)\)에서 담아야 할 참값은 \(0.03\)이다. 같은 칸의 구간 \([0,0]\)은 \(0.03\)을 담지 못한다. 이 경우 붕괴 확률은

로 더 크고, 그것이 통째로 실패로 들어간다. 보기 5의 그림에서 이 한 칸이 전체 실패 확률 \(0.1518\)의 77%를 차지하는 것을 볼 수 있다.

교훈. \(p_1 = p_2\)에서만 포함률을 재면 왈드가 멀쩡해 보인다. 신뢰구간은 모든 \((p_1, p_2)\)에서 명목 수준을 지켜야 하는 물건이므로, 귀무가설 위에서만 평가하는 것은 의미가 없다. 신뢰구간의 존재 이유가 \(p_1 - p_2\)가 0이 아닐 가능성을 다루는 데 있기 때문이다.

첫째, 경계에서 밀어낸다. \(\tilde p_i = (X_i+1)/(n_i+2)\)는 \(X_i\)가 무엇이든 \(1/(n_i+2)\)와 \((n_i+1)/(n_i+2)\) 사이에 있다. 결코 0이나 1이 되지 않으므로 \(\tilde p_i\tilde q_i > 0\)이고 표준오차가 0이 될 수 없다. 붕괴가 구조적으로 불가능하다.

둘째, 왈드의 편향을 상쇄한다. \(\hat p\hat q\)는 평균적으로 \(pq\)보다 작다. 실제로

이고, \(\hat p\)가 경계에 가까우면 어긋남이 특히 심하다(\(\hat p = 0\)이면 \(0\)). 그래서 왈드 구간이 좁아진다. \(\tilde p\)는 \(\hat p\)를 \(0.5\) 쪽으로 당기므로(\(\tilde p - \hat p = \frac{n-2X}{n(n+2)}\)) \(\tilde p\tilde q \ge \hat p\hat q\)가 되어 구간이 적절히 넓어진다.

셋째, \(1\)인 이유. 한 표본의 \(p\)에 대한 95% 윌슨 구간의 중심은

이며 \(z = 1.96\)이면 \(z^2/2 = 1.92\), \(z^2 = 3.84\)다. 곧 "성공 약 2개와 실패 약 2개를 더한" 자리다. 이것이 한 표본의 더하기 4 규칙이다.

두 표본의 차에서는 사정이 다르다. 두 표본 각각에 2씩 더하면 보정이 겹쳐 지나치게 보수적이 된다. 아그레스티와 카포는 각 표본에 1씩만 더했을 때 포함률이 \(0.95\)에 가장 가깝게 유지된다는 것을 보였다. 합치면 성공 2개와 실패 2개가 더해지는 셈이니, 한 표본의 더하기 4 규칙과 총량이 맞는다.

대가. 구간이 넓어진다. \(p_1 = p_2 = 0.05\)에서 평균 폭이 왈드 \(0.2055\), 아그레스티-카포 \(0.2573\), 윌슨 \(0.2876\)이다. 그러나 \(p_1 = p_2 = 0.20\)에서는 각각 \(0.3960\), \(0.3982\), \(0.3991\)로 거의 같아진다. 작은 \(p\)에서만 값을 치르고, 그 자리가 정확히 왈드가 실패하는 자리다.

이름의 뜻. "정확"은 \(p\)-값을 근사 없이 계산한다는 뜻이다. 주변합을 모두 고정하면 \(H_0\) 아래에서 \(X_1\)이 초기하분포를 따르고, 그 확률을 그대로 더해 \(p\)-값을 얻는다. 정규근사가 개입하지 않는다.

왜 크기가 맞지 않는가. 검정통계량이 이산이므로 달성 가능한 \(p\)-값이 띄엄띄엄하다. \(p\)-값이 \(0.05\) 바로 아래에 오는 자료 구성은 몇 개 없고, 대부분은 \(0.05\)를 훌쩍 넘거나 훨씬 밑돈다. "\(p \le 0.05\)"라는 규칙이 실제로 뽑아내는 확률은 \(0.05\)가 아니라 그보다 작은 어떤 값이다.

여기에 조건화가 더해진다. 주변합을 둘 다 고정하는 것은 실제 설계와 맞지 않는 경우가 많다. 보통은 \(n_1\)과 \(n_2\)만 고정되고 성공의 총합은 고정되지 않는다. 필요 이상으로 조건을 걸면 쓸 수 있는 자료 구성이 줄어들어 검정이 더 보수적이 된다. 보기 3에서 최댓값이 \(0.0274\)이고 \(p = 0.05\)에서는 \(0.0015\)로 명목의 3%에 불과한 이유다.

실무적 함의. 피셔 검정은 1종오류를 절대 초과하지 않지만 그 대가로 검정력을 잃는다. 검정력이 중요하면 주변합을 하나만 고정하는 무조건부 정확검정(바너드 검정)을 쓰거나, \(p\)-값의 덩어리를 절반만 세는 중간 \(p\)-값을 쓴다. 어느 쪽도 이산성을 없애지는 못하고 완화할 뿐이다.

톱니의 기원. 포함률은 참값 \(\delta\)의 함수로

이다. 여기서 \((L_{ij}, U_{ij})\)는 자료 \((i,j)\)에서 만든 구간이며, \(961\)개의 칸마다 하나씩, 모두 유한개다.

\(\delta\)를 연속으로 움직이면 각 칸의 조건 \(L_{ij} \le \delta \le U_{ij}\)는 \(\delta\)가 \(L_{ij}\)나 \(U_{ij}\)를 지나는 그 순간에만 바뀌고, 그때 그 칸의 확률이 합에서 통째로 빠지거나 들어온다. 곧 \(C(\delta)\)는 조각마다 매끄럽고 유한개의 점에서 뛰는 함수이며, 뛰는 폭이 그 칸의 확률이다. 작은 \(p\)에서 진폭이 큰 이유도 여기 있다. 확률이 몇 칸에 몰려 있어 한 칸이 빠질 때의 낙차가 크다. \(\square\)

왼쪽 그림이 매끄러운 이유는 두 가지다. 첫째, 왼쪽에서는 참값이 \(\delta = 0\)으로 고정되어 있다. 뛰게 만드는 것은 \(\delta\)의 이동이므로, \(\delta\)가 멈춰 있으면 칸의 포함 여부도 바뀌지 않는다. 움직이는 것은 확률 \(P(X_1=i)P(X_2=j)\)뿐이고 이것은 \(p\)의 다항식이므로 \(C(p)\)가 연속이다. 둘째, 두 이항의 차는 확률이 더 고르게 퍼진다. 격자 간격 자체는 한 표본과 같은 \(1/30\)이지만, 한 표본의 \(\hat p\)가 \(31\)개 값에 확률을 싣는 데 비해 \(\hat p_1 - \hat p_2\)는 \(61\)개 값에 싣고 각 값에 여러 칸이 대응해 덩어리가 잘게 쪼개진다. 한 표본의 포함률 곡선에서 보이는 거친 톱니가 두 표본에서 상당 부분 완화되는 까닭이다.

진동의 원인은 자료의 이산성이지만, 드러나는 방향은 무엇을 축으로 놓느냐에 달려 있다. 참값을 축으로 놓으면 톱니가 보이고 공통 비율을 축으로 놓으면 숨는다. 한 방향만 보고 "이 방법은 괜찮다"고 결론 내리면 안 된다.

각 비율의 윌슨 구간 \([l_i,\ u_i]\)를 구한 뒤

로 둔다. 하한에는 \(l_1\)과 \(u_2\)를, 상한에는 \(u_1\)과 \(l_2\)를 쓴다는 점이 핵심이다. 차가 가장 작아지는 조합과 가장 커지는 조합을 각각 쓰는 것이다.

import numpy as np
from scipy import stats

z = 1.959964

def wilson(k, n):
    ph = k/n; c = 1 + z**2/n
    ctr = (ph + z**2/(2*n))/c
    half = z*np.sqrt(ph*(1-ph)/n + z**2/(4*n**2))/c
    return ctr-half, ctr+half

def wald(k1,n1,k2,n2):
    p1,p2 = k1/n1, k2/n2
    s = z*np.sqrt(p1*(1-p1)/n1 + p2*(1-p2)/n2)
    return p1-p2-s, p1-p2+s

def ac(k1,n1,k2,n2):
    p1,p2 = (k1+1)/(n1+2), (k2+1)/(n2+2)
    s = z*np.sqrt(p1*(1-p1)/(n1+2) + p2*(1-p2)/(n2+2))
    return p1-p2-s, p1-p2+s

def newcombe(k1,n1,k2,n2):
    l1,u1 = wilson(k1,n1); l2,u2 = wilson(k2,n2)
    d = k1/n1 - k2/n2
    return (d - z*np.sqrt(l1*(1-l1)/n1 + u2*(1-u2)/n2),
            d + z*np.sqrt(u1*(1-u1)/n1 + l2*(1-l2)/n2))

def cover(n1,n2,p1,p2,f):
    tot = 0.0; d0 = p1-p2
    for k1 in range(n1+1):
        w1 = stats.binom.pmf(k1,n1,p1)
        for k2 in range(n2+1):
            lo,hi = f(k1,n1,k2,n2)
            if lo <= d0 <= hi:
                tot += w1*stats.binom.pmf(k2,n2,p2)
    return tot

print(f"{'p1,p2':>12}{'Wald':>9}{'AC':>9}{'Newcombe':>11}")
for p1,p2 in [(0.05,0.05),(0.05,0.02),(0.2,0.1),(0.5,0.5)]:
    print(f"{f'{p1},{p2}':>12}" +
          "".join(f"{cover(30,30,p1,p2,f):>9.4f}" for f in (wald,ac,newcombe)))

출력:

       p1,p2     Wald       AC   Newcombe
   0.05,0.05   0.9716   0.9929   0.9929
   0.05,0.02   0.8482   0.9976   0.9975
     0.2,0.1   0.9425   0.9622   0.9631
     0.5,0.5   0.9478   0.9481   0.9481

왈드가 무너지는 자리에서 둘 다 살아난다. \((0.05,\ 0.02)\)에서 왈드 \(0.8482\)인데 AC와 뉴콤은 \(0.9975\) 안팎이다.

AC와 뉴콤이 거의 같다. 네 상황 모두에서 소수 셋째 자리까지 일치한다. 이 표본크기에서는 둘 중 아무거나 써도 좋다는 뜻이다.

그러면 왜 뉴콤을 따로 아는가. 표본크기가 크게 불균형하거나 \(p\)가 극단으로 갈 때 뉴콤이 조금 더 안정적이고, 무엇보다 윌슨 구간이라는 하나의 원리에서 일관되게 유도된다. AC의 "\(+1\)씩 더한다"는 어림에 가까운 반면, 뉴콤은 5.5절 포함률 문서 연습문제 9의 점수 검정 논리를 두 표본으로 확장한 것이다.

둘 다 보수적이라는 점은 유의하라. \(p\)가 작을 때 포함률이 \(0.99\)를 넘는다. 명목 \(0.95\)보다 안전하지만 구간이 그만큼 넓다. 이산성 때문에 정확히 \(0.95\)를 맞추는 것은 불가능하며(5.5절 포함률 문서 연습문제 10), 어느 방향으로 틀릴지 고르는 문제다.

차이는 무엇을 고정하는가에 있다.

피셔는 열합까지 고정해 조건부 분포를 쓰므로 계산이 간단하지만, 실제로는 열합이 고정된 설계가 아닌 경우가 대부분이다. 두 군에 각각 \(n_1\), \(n_2\)명을 배정했을 뿐 전체 성공 수를 정해 둔 것이 아니다. 불필요하게 조건을 걸어 정보를 버리는 셈이라 보수적이 된다.

바너드는 조건을 걸지 않는 대신 성가신 모수 \(p\)를 남긴다. 이를 처리하려고 \(p\)에 대해 최댓값을 취한다(최대화 \(p\)값).

import numpy as np
from scipy import stats

n1 = n2 = 12
tab = [[9, 3], [4, 8]]      # 관측: 1군 9/12 성공, 2군 4/12 성공
print(f"  피셔  양측 p = {stats.fisher_exact(tab)[1]:.4f}")
try:
    print(f"  바너드 양측 p = {stats.barnard_exact(tab).pvalue:.4f}")
    print(f"  보슈루프 p    = {stats.boschloo_exact(tab).pvalue:.4f}")
except AttributeError:
    print("  (scipy 1.7+ 필요)")

출력:

  피셔  양측 p = 0.0995
  바너드 양측 p = 0.0485
  보슈루프 p    = 0.0494

같은 자료에서 피셔는 \(0.0995\), 바너드는 \(0.0485\)다. \(\alpha=0.05\) 기준으로 결론이 갈린다(피셔는 기각 못 하고 바너드는 기각한다). 피셔의 보수성이 실질적인 검정력 손실을 낳는다는 것이 이 한 예로 드러난다.

그러면 왜 피셔가 표준처럼 쓰이는가.

• 역사와 관성. 손으로 계산할 수 있었고 교과서에 먼저 실렸다.
• 계산 비용. 바너드는 \(p\)에 대한 최대화를 해야 해서 표가 커지면 느리다.
• 조건부 논증에 대한 철학적 선호. 열합이 성가신 모수에 대한 보조통계량이므로 조건을 거는 것이 옳다는 견해도 있다.

오늘날의 권고는 바뀌고 있다. \(2\times2\) 표에서 두 군의 크기만 고정된 설계라면 바너드(또는 보슈루프)가 더 강력하면서도 크기를 지킨다는 것이 다수 의견이다. scipy.stats.barnard_exact, boschloo_exact로 바로 쓸 수 있다.

정의. 통상의 \(p\)값은 "관측값만큼 또는 그보다 극단적인" 확률인데, 중간 \(p\)값은 관측값 자체의 확률을 절반만 센다.

import numpy as np
from scipy import stats

def fisher_mid(k1, n1, k2, n2):
    """단측 중간 p값 (초기하)."""
    m, nn, N = k1 + k2, n1, n1 + n2
    rv = stats.hypergeom(N, m, nn)
    return rv.sf(k1) + 0.5 * rv.pmf(k1)

rng = np.random.default_rng(0)
n1 = n2 = 15
REP = 40_000
for p in (0.3, 0.5):
    a = rng.binomial(n1, p, REP); b = rng.binomial(n2, p, REP)
    pf = np.array([stats.fisher_exact([[x, n1-x], [y, n2-y]],
                                      alternative='greater')[1]
                   for x, y in zip(a[:4000], b[:4000])])
    pm = np.array([fisher_mid(x, n1, y, n2)
                   for x, y in zip(a[:4000], b[:4000])])
    print(f"  p={p}: 피셔 단측 크기 {np.mean(pf<0.05):.4f}   "
          f"중간 p값 {np.mean(pm<0.05):.4f}")

출력:

  p=0.3: 피셔 단측 크기 0.0213   중간 p값 0.0428
  p=0.5: 피셔 단측 크기 0.0217   중간 p값 0.0490

피셔가 명목 \(0.05\)의 절반도 안 된다(\(0.021 \sim 0.022\)). 중간 \(p\)값을 쓰면 \(0.043 \sim 0.049\)로 명목에 훨씬 가까워진다.

왜 듣는가. 이산분포에서는 "\(p\)값 \(\le 0.05\)"를 만족하는 표가 듬성듬성해서, 실제 기각 확률이 \(0.05\)에 한참 못 미친다. 관측값의 확률을 절반만 세면 경계에 있는 표들이 기각역에 포함되어 틈이 메워진다.

대가는 보장을 잃는다는 것이다.

여기서는 두 \(p\) 모두에서 \(0.05\)를 넘지 않았지만 보장은 없다. 다른 \(p\)나 표본크기에서는 명목을 넘을 수 있으며, 그 순간 엄밀한 의미의 "수준 \(\alpha\) 검정"이 아니게 된다.

그래서 쓰임이 갈린다. 규제·안전성 판정처럼 보장이 필요하면 피셔를, 역학이나 유전체 연구처럼 여러 검정을 합쳐 평균적 성질이 중요하면 중간 \(p\)값을 쓴다. 특히 다중검정 보정과 결합할 때 피셔의 보수성이 누적되면 검정력이 거의 사라지므로, 중간 \(p\)값이 실질적인 대안이 된다.

먼저 설계를 확인한다. 어떤 주변합이 설계에 의해 고정되었는지가 방법을 정한다.

1. 두 군의 크기 n1, n2 가 설계로 정해졌는가?
     예  -> 행합 고정. 바너드/보슈루프 또는 근사 검정
     아니오(전체 N 만 정함) -> 다항 모형
2. 전체 성공 수도 설계로 고정되었는가?
     예  -> 피셔가 자연스럽다(드물다)
     아니오 -> 피셔는 보수적임을 감수하는 선택

그다음 표본크기와 기대도수를 본다.

쓰지 말아야 할 것도 분명하다.

• 왈드 구간. \(p\)가 극단이면 포함률이 \(0.85\)까지 떨어진다(연습문제 3). 대안이 있는데 쓸 이유가 없다.
• 연속성 수정한 \(Z\) 검정. 앞 쪽 연습문제 7에서 보았듯 과보정으로 크기가 절반이 된다.
• 기대도수가 작은데 쓰는 카이제곱 근사. 예이츠 보정으로 가리려 하지 말고 정확검정으로 옮긴다.

보고할 때 함께 적을 것.

- 네 칸의 도수 (비율만 적지 않는다)
- 쓴 방법의 이름 (피셔/바너드/카이제곱, 보정 여부)
- 효과크기와 구간 (차 또는 오즈비, 어느 쪽인지 명시)
- 표본이 작으면 그 사실과 한계

마지막으로 검정보다 구간이 낫다. \(p\)값은 "차이가 있는가"에만 답하지만 구간은 "얼마나 다른가"와 "얼마나 모르는가"를 함께 말한다. 5.5절 희귀사건 문서 연습문제 10에서 두 군 모두 \(0\)건일 때 \(p\)값이 \(1.000\)이었지만 구간은 \(2{,}600\)배에 걸쳐 있었다. \(p\)값만 보았다면 "차이 없음"으로 오독했을 자리였다. \(\square\)