p̂의 표본분포 (희귀사건)¶

(1) 이론값. \(n\)은 고정이고 \(p\)만 바뀐다. 왜도는 \(\gamma_1 = (1-2p)/\sqrt{npq}\)이므로

표본크기는 셋 다 \(100\)인데 왜도가 \(0\)에서 \(0.69\)까지 간다. \(p\)가 작으면 \(\gamma_1 \approx 1/\sqrt{np}\)이고(\(p = 0.02\)에서 \(0.6857\) 대 \(0.7071\)로 \(3\%\) 차이), 그래서 \(n\)은 홀로 등장하지 않고 늘 \(np\)로만 나타난다.

마지막 열이 다른 종류의 문제다. \(p = 0.02\)에서 \(0.98^{100} = 0.1326\)이다. 여덟 번에 한 번꼴로 성공이 한 번도 나오지 않는다. 정규분포는 연속이라 어떤 한 점에도 확률을 주지 못하므로 이 \(13\%\)를 담을 자리가 없다. 대신 왼쪽으로 새어 나간다.

존재할 수 없는 영역에 \(7.66\%\)를 흘리는 것이고, 이는 점질량 \(0.1326\)의 \(58\%\)에 해당한다. 포아송은 이산이므로 같은 자리에 \(e^{-2} = 0.1353\)을 제대로 얹는다.

(2) 모의실험.

import matplotlib.pyplot as plt
import numpy as np
from scipy import stats

n = 100

print("   p    np   왜도    1/sqrt(np)   P(p^=0) 정확 / 포아송   정규가 음수에 흘리는 확률")
fig, axes = plt.subplots(1, 3, figsize=(12, 3.8))
for ax, p in zip(axes, (0.5, 0.1, 0.02)):
    se = np.sqrt(p * (1 - p) / n)
    k = np.arange(n + 1)
    pmf = stats.binom(n, p).pmf(k)
    x = k / n

    print(f"{p:>5} {n * p:>5.0f} {(1 - 2 * p) / np.sqrt(n * p * (1 - p)):>6.4f} "
          f"{1 / np.sqrt(n * p):>11.4f} {pmf[0]:>13.4f} / {stats.poisson(n * p).pmf(0):.4f}"
          f"{stats.norm.cdf((0 - p) / se):>22.4f}")

    lo, hi = p - 4.5 * se, p + 4.5 * se
    if lo < 0.02:            # 0 의 점질량이 보이도록 왼쪽을 0 아래까지 연다
        lo = -0.004
    m = (x >= lo) & (x <= hi)
    ax.vlines(x[m], 0, pmf[m], color="tab:blue", lw=1.2, alpha=0.75,
              label=r"exact PMF of $\hat p$")
    ax.plot(x[m], pmf[m], "o", color="tab:blue", ms=4)

    g = np.linspace(lo, hi, 400)
    ax.plot(g, stats.norm(p, se).pdf(g) / n, "--r", lw=1.8,
            label=r"$N(p,\,pq/n)$ / $n$")

    if p == 0.02:
        # 이 패널에서만 포아송 근사를 함께 올린다. 이산분포라 0 을 다룰 수 있다.
        ax.plot(x[m], stats.poisson(n * p).pmf(k[m]), "x", color="tab:green",
                ms=6, label=r"Poisson($np$)")
        ax.annotate(f"$P(\\hat p = 0) = {(1-p)**n:.4f}$",
                    xy=(0, pmf[0]), xytext=(0.012, pmf[0] * 0.97),
                    fontsize=9, arrowprops=dict(arrowstyle="->", lw=1))

    ax.axvline(0, color="gray", lw=0.8, ls=":")
    ax.set_title(f"p = {p},  np = {n*p:g},  "
                 f"skew = {(1-2*p)/np.sqrt(n*p*(1-p)):.3f}", fontsize=10)
    ax.set_xlim(lo, hi)
    ax.set_ylim(bottom=0)
    ax.set_xlabel(r"$\hat p$")
    ax.legend(fontsize=8)

axes[0].set_ylabel("Probability")
plt.tight_layout()
plt.show()

출력:

   p    np   왜도    1/sqrt(np)   P(p^=0) 정확 / 포아송   정규가 음수에 흘리는 확률
  0.5    50 0.0000      0.1414        0.0000 / 0.0000                0.0000
  0.1    10 0.2667      0.3162        0.0000 / 0.0000                0.0004
 0.02     2 0.6857      0.7071        0.1326 / 0.1353                0.0766

\(p = 0.5\) 패널은 교과서 그림이다. 왜도가 정확히 \(0\)이고 점들이 곡선 위에 얹혀 좌우가 대칭이다.

\(p = 0.1\) 패널에서 왜도가 보인다. \(\gamma_1 = 0.2667\)이고, 그림에서는 왼쪽 꼬리의 점들이 곡선보다 낮고 오른쪽 꼬리의 점들이 곡선보다 높은 모습으로 나타난다. 오른쪽으로 끌린 것이다. 그래도 아직 "조금 어긋난 정규"라 부를 만하다.

\(p = 0.02\) 패널은 종류가 다르다. 왜도가 \(0.6857\)로 커진 것은 사실이지만, 그것만으로 설명이 끝나지 않는다. 가능한 값이 실질적으로 \(\hat p = 0, 0.01, \ldots, 0.07\)의 여덟 개뿐이고, 그 가운데 왼쪽 끝 하나에 \(0.1326\)이 쌓여 있다. 붉은 곡선은 그 점을 지나 점선(0의 위치) 왼쪽으로 계속 뻗어 존재하지 않는 영역에 \(7.66\%\)를 둔다.

그러므로 답은 "왜도가 커진 탓"이 아니다. 근사의 지지집합이 틀렸다. 왜도는 모양에 관한 흠이라 에지워스 보정이나 연속성 수정으로 깎을 수 있지만, 연속분포로 점질량을 표현하는 일은 어떤 보정으로도 되지 않는다. 가로축 \(\hat p < 0\)에 밀도가 놓여 있다는 사실 자체가 고칠 수 없는 결함이다.

초록색 \(\times\)가 대안을 보여 준다. 포아송은 이산분포이고 지지집합이 \(\{0, 1, 2, \ldots\}/n\)으로 맞으므로, \(\hat p = 0\)에 \(0.1353\)을 얹어 정확값 \(0.1326\)을 \(0.003\) 차이로 맞힌다. 점질량 위에 거의 그대로 겹치는 것이 그림에서 보인다. \(p\)가 작을 때 바꿀 것은 보정이 아니라 근사하는 분포다.

(1) 두 어림값.

정규 쪽. 앞 쪽에서 얻은 분해를 그대로 쓴다.

\(n = 100\)에서 \(\sqrt{npq}\)가 \(5,\ 3,\ 1.4\)이므로

가 된다. 두 항이 모두 커지는데 격자 항이 더 빨리 커진다. \(p\)를 줄이면 \(\sqrt{npq}\)가 작아져 격자가 거칠어지기 때문이고, \(p = 0.02\)에서는 격자 항이 왜도 항의 세 배다.

포아송 쪽. 르캉 부등식이 총변동거리를 억제한다.

누적분포의 최대 차는 총변동거리 이하이므로 이 값이 그대로 상한이 된다. \(n = 100\)에서 \(np^2\)는 \(25\), \(1\), \(0.04\)다. 앞의 둘은 \(1\) 이상이라 아무것도 말해 주지 않는 수이고, \(p = 0.02\)에서만 \(0.04\)라는 쓸모 있는 보장을 준다.

(2) 수치로.

import numpy as np
from scipy import stats

n = 100
k = np.arange(n + 1)

print("   p    np      SE    skew   P(p^=0)  Normal  Poisson  sup|F-Phi|  sup|F-Pois|  에지워스 예측  르캉 상한")
for p in (0.5, 0.1, 0.02):
    q = 1 - p
    se = np.sqrt(p * q / n)
    skew = (1 - 2 * p) / np.sqrt(n * p * q)

    binom = stats.binom(n, p)
    pois = stats.poisson(n * p)

    # p^ = 0 이 될 확률을 세 방법으로 구한다.
    exact0 = binom.pmf(0)                     # (1-p)^n
    normal0 = stats.norm.cdf((0 - p) / se)    # 정규근사가 주는 P(p^ <= 0)
    pois0 = pois.pmf(0)                       # e^{-np}

    # 격자점에서 CDF 의 최대 어긋남
    F = binom.cdf(k)
    d_norm = np.max(np.abs(F - stats.norm.cdf((k / n - p) / se)))
    d_pois = np.max(np.abs(F - pois.cdf(k)))

    # 어림값 두 개. 정규 쪽은 격자항 + 왜도항, 포아송 쪽은 르캉 부등식의 상한이다.
    edgeworth = stats.norm.pdf(0) / 2 / np.sqrt(n * p * q) + stats.norm.pdf(0) / 6 * skew
    lecam = n * p ** 2

    print(f"{p:>5} {n*p:>5.0f} {se:>7.4f} {skew:>7.4f} {exact0:>9.4f} "
          f"{normal0:>7.4f} {pois0:>8.4f} {d_norm:>11.4f} {d_pois:>12.4f} "
          f"{edgeworth:>13.4f} {lecam:>10.4f}")

출력:

   p    np      SE    skew   P(p^=0)  Normal  Poisson  sup|F-Phi|  sup|F-Pois|  에지워스 예측  르캉 상한
  0.5    50  0.0500  0.0000    0.0000  0.0000   0.0000      0.0398       0.0854        0.0399    25.0000
  0.1    10  0.0300  0.2667    0.0000  0.0004   0.0000      0.0832       0.0142        0.0842     1.0000
 0.02     2  0.0140  0.6857    0.1326  0.0766   0.1353      0.1767       0.0027        0.1881     0.0400

에지워스 어림이 맞는다. \(p = 0.5\)에서 \(0.0399\) 대 실제 \(0.0398\), \(p = 0.1\)에서 \(0.0842\) 대 \(0.0832\)로 \(1\%\) 차이다. \(p = 0.02\)에서는 \(0.1881\) 대 \(0.1767\)로 \(6\%\) 크게 어림하는데, 그 자리에서는 \(np = 2\)여서 전개의 다음 항까지 무시할 수 없기 때문이다. \(np = 2\)에서 에지워스가 흔들린다는 것 자체가 "정규근사를 쓰지 말라"는 신호다.

르캉 상한은 옳지만 느슨하다. \(p = 0.02\)에서 실제 오차가 \(0.0027\)인데 상한은 \(0.0400\)으로 \(15\)배 크다. 부등식이 깨진 적은 없고(세 줄 모두 실제 \(\le\) 상한), 다만 \(p = 0.5\)와 \(0.1\)에서는 상한이 \(25\)와 \(1\)이라 쓸모가 없다. 보장의 쓸모는 그것이 \(1\)보다 작아지는 자리에서 생긴다.

자리를 바꾸는 지점. 정규근사의 오차는 \(0.0398 \to 0.0832 \to 0.1767\)로 커지고 포아송은 \(0.0854 \to 0.0142 \to 0.0027\)로 작아진다. \(p = 0.5\)에서는 정규가 두 배 낫고, \(p = 0.1\)에서는 벌써 포아송이 \(5.9\)배 낫다. \(p = 0.02\)에 이르면 \(0.1767/0.0027 = 65\)배다. 역전은 \(p = 0.1\)보다 위에서 이미 일어나 있다.

역전점이 정확히 어디인지는 \(p\)를 훑어 재면 된다. \(n = 100\)에서

이므로 역전은 \(p = 0.30\)과 \(0.35\) 사이에서 일어난다. \(p = 0.1\)은 역전점이 아니라 이미 포아송이 한참 앞선 자리다. 르캉 상한으로는 이 경계를 찾을 수 없다. 그 자리에서 \(np^2\)가 \(9\)에서 \(12\) 사이여서 아무 정보가 없기 때문이고, 부등식은 "안전하다"를 말할 수 있어도 "어느 쪽이 더 나은가"를 말해 주지 못한다.

마지막으로 \(P(\hat p = 0)\) 열을 다시 보라. \(p = 0.02\)에서 정확값 \(0.1326\)을 포아송은 \(0.1353\)으로 맞히고 정규는 음수 영역의 \(0.0766\)으로 "대신한다". \(\sup\lvert F - \Phi \rvert = 0.1767\)이라는 수의 대부분이 바로 이 한 점에서 생긴다.

\(\hat p = 0\)은 100번 모두 실패한다는 뜻이므로

이다. \(e^{-np} = e^{-2} = 0.1353\)과 가깝다(\(\log 0.98 = -0.0202\)이므로 \(0.98^{100} = e^{-2.02}\)다).

정규근사는 이 확률을 다루지 못한다. 연속분포이므로 한 점의 확률이 0이다. 굳이 대응시키자면 근사분포가 \(\hat p \le 0\)에 주는 확률

이 그 자리를 대신하지만, 정확값 0.1326의 6할에도 못 미치고 그마저도 음수 영역에 놓인다.

현실적 함의가 있다. 100명을 검사해 아무도 양성이 아니었다는 결과는 \(p = 0.02\)인 모집단에서 13%의 확률로 벌어지는 일이다. "아무도 없었으니 \(p = 0\)"이라고 결론지을 수 없다.

두 조건 중 더 까다로운 쪽이 결정한다.

따라서 \(n \ge 250\)이다. \(p\)가 작을 때 걸리는 조건은 언제나 \(np \ge 5\)다.

\(p = 0.5\)에서는 두 조건이 모두 \(n \ge 10\)이다. 25배 차이다. 일반적으로 \(n \ge 5/\min(p, 1-p)\)이며 \(p\)가 0이나 1에 가까울수록 발산한다.

다만 \(n = 250\)은 느슨한 기준을 겨우 채운 값일 뿐 좋은 근사라는 보장은 아니다. \(np = 5\)에서 왜도는 여전히 0.44다. 신뢰구간을 세울 작정이라면 보수적 기준 \(np \ge 10\)을 써야 하고, 그러면 \(p = 0.02\)에서 \(n \ge 500\)으로 요구가 두 배가 된다.

\(\gamma_1 = (1-2p)/\sqrt{npq}\)에 \(n = 100\)을 넣는다.

\(\gamma_1 < 0.1\)을 \(n\)에 대해 풀면

이므로 \(n \ge 4703\)이 필요하다.

앞 페이지의 \(p = 0.3\)에서는 \(n = 77\)이면 충분했다. 같은 대칭성을 얻는 데 표본이 61배 더 든다. \(np\)로 바꿔 읽으면 두 경우 모두 \(np \approx 23 \sim 94\)쯤에서 왜도가 0.1 아래로 내려가며, 결국 하나의 기준 \(np\)로 통일된다는 것을 알 수 있다.

총변동거리는 정의상 1을 넘지 않으므로 상한 25나 1은 아무 정보가 없다. 부등식이 쓸모 있는 것은 \(np^2 < 1\), 즉 \(p < 1/\sqrt n\)일 때다. \(n = 100\)이면 \(p < 0.1\)이다.

\(p = 0.02\)에서 상한 0.04는 실제 0.0027의 15배쯤으로 느슨하지만, 자료를 보지 않고 미리 보장되는 값이라는 점이 중요하다. 정규근사에는 이런 형태의 보장이 없다. 베리-에센 상한은 \(C\gamma_1\)의 꼴이어서 \(p\)가 작으면 함께 커진다.

\(np^2 = \lambda p\)로 쓰면 부등식의 의미가 잘 보인다. 포아송 근사의 오차는 \(\lambda = np\)가 아니라 개별 성공확률 \(p\)가 작은지에 달려 있다. \(\lambda\)는 아무리 커도 좋고, \(p\)만 작으면 된다.

\(\hat p \ge 0.05\)는 \(X \ge 5\)와 같다.

포아송은 정확값과 0.002 차이지만 정규근사는 3배 이상 작다.

결론이 뒤집힌다. \(H_0: p = 0.02\)를 단측 유의수준 0.05로 검정한다고 하자.

• 정확한 계산: \(p\text{-값} = 0.0508 > 0.05\)이므로 기각하지 못한다.
• 정규근사: \(p\text{-값} = 0.0161 < 0.05\)이므로 기각한다.

같은 자료에서 정반대의 결론이 나온다. 정규근사가 오른쪽 꼬리를 과소평가하기 때문이다. 왜도가 양수인 분포의 오른쪽 꼬리는 정규보다 두껍다.

방향이 위험한 쪽이다. 희귀사건 자료에서 정규근사를 쓰면 \(p\)값이 실제보다 작게 나와 없는 유의성을 만들어 낸다. 이것이 5.5절의 규칙을 지켜야 하는 실무적 이유다.

(나) 정확한 방법부터. 상한 \(p_U\)는 "\(p = p_U\)였다면 지금처럼 0을 볼 확률이 겨우 5%"가 되는 값으로 정한다.

풀면

이다. 이것이 \(x = 0\)에서의 클로퍼-피어슨 상한이다.

(가) 삼의 법칙. 위 식의 로그를 취한다.

\(p_U\)가 작으면 \(\log(1-p_U) \approx -p_U\)이므로

이다. \(\square\) \(n = 100\)이면 \(3/100 = 0.03\)이다.

두 값의 비교. 정확값 0.0295와 어림값 0.0300의 차이는 0.0005다. 근사에 쓰인 것은 \(\log(1-p) \approx -p\)뿐이고 그 오차가 \(O(p^2)\)이므로, \(p_U\)가 작을수록, 즉 \(n\)이 클수록 더 잘 맞는다. 이름의 3은 \(-\log 0.05 = 2.996\)에서 왔다.

참고로 같은 상황에서 다음 페이지에 나올 윌슨 구간의 상한은

이고, 왈드 구간은 \(\hat p \pm 1.96\sqrt{0 \times 1/100} = [0, 0]\)으로 붕괴한다. 성공이 0이라는 자료에서 왈드 구간은 "\(p\)는 정확히 0"이라고 주장한다.

읽는 법. 성공이 하나도 없다는 자료가 말해 주는 것은 "\(p = 0\)"이 아니라 "\(p\)는 대략 \(3/n\)보다 작다"다. 약물 부작용, 제품 결함, 희귀질환 검사처럼 0을 관측하는 일이 흔한 분야에서 이 어림법이 널리 쓰인다.

포아송 구간은 카이제곱으로 적힌다. \(K \sim \text{Poisson}(n\lambda)\)일 때

이다(\(k=0\)이면 하한은 \(0\)).

import numpy as np
from scipy import stats

n = 500
print(f"{'k':>4}{'이항 CP 구간':>28}{'포아송 구간':>30}")
for k in (0, 1, 3, 10):
    lo_b = 0 if k == 0 else stats.beta.ppf(0.025, k, n - k + 1)
    hi_b = stats.beta.ppf(0.975, k + 1, n - k)
    lo_p = 0 if k == 0 else stats.chi2.ppf(0.025, 2*k) / 2 / n
    hi_p = stats.chi2.ppf(0.975, 2*k + 2) / 2 / n
    print(f"{k:>4}   [{lo_b:.6f}, {hi_b:.6f}]      [{lo_p:.6f}, {hi_p:.6f}]")

출력:

   k                    이항 CP 구간                        포아송 구간
   0   [0.000000, 0.007351]      [0.000000, 0.007378]
   1   [0.000051, 0.011092]      [0.000051, 0.011143]
   3   [0.001239, 0.017434]      [0.001237, 0.017535]
  10   [0.009631, 0.036472]      [0.009591, 0.036781]

소수점 셋째 자리까지 같다. \(k=10\)에서 상한이 \(0.03647\) 대 \(0.03678\)로 \(1\%\) 차이 나는 것이 가장 큰 불일치다.

이유는 이 페이지의 주제 그대로다. \(p\)가 작고 \(n\)이 크면 이항이 포아송에 수렴하며, 르캉 부등식(연습문제 4)이 그 거리를 \(np^2\)으로 묶는다. 여기서는 \(p \lesssim 0.02\)이므로 \(np^2 \lesssim 0.2\)이고, 실제 차이는 그보다 훨씬 작다.

포아송 쪽이 편리한 자리가 있다.

• \(n\)을 모를 때. "이 도시에서 작년에 희귀질환 \(7\)건"처럼 노출 인구·기간만 알고 개별 시행 수를 세지 않는 자료에서는 이항을 쓸 수 없다. 포아송은 발생률 \(\lambda\)(단위 시간·인구당)를 직접 다룬다.
• 노출량이 관측마다 다를 때. 사람–년(person-year)을 더해 총 노출로 나누는 방식이 자연스럽다.
• 합치기가 쉽다. 포아송은 가법적이라 여러 지역의 건수를 더하면 그대로 포아송이다.

그래서 역학에서는 발생률 구간에 포아송을 쓴다. 이항 정확구간과 사실상 같은 답을 주면서 자료 구조에 더 잘 맞기 때문이다. \(k=0\)일 때 상한이 \(\chi^2_{0.975}(2)/(2n) = 3.689/n\)인데, 이것이 연습문제 6의 삼의 법칙 \(3/n\)의 정확판이다(양측 \(95\%\) 기준).

\(\operatorname{sd}(\hat p)/p = \sqrt{(1-p)/(np)}\)를 \(0.20\) 이하로 만든다.

import numpy as np

print(f"{'p':>9}{'필요한 n':>14}{'기대 사건 수 np':>18}")
for p in (0.1, 0.01, 0.001, 0.0001):
    n = int(np.ceil((1 - p) / (p * 0.20**2)))
    print(f"{p:>9}{n:>14,}{n*p:>18.1f}")

출력:

        p         필요한 n        기대 사건 수 np
      0.1           225              22.5
     0.01         2,475              24.8
    0.001        24,975              25.0
   0.0001       249,975              25.0

기대 사건 수가 모두 \(25\)로 수렴한다. \(p\)가 \(1{,}000\)배 작아지면 \(n\)도 \(1{,}000\)배 커지는데, 그 곱인 사건 수는 그대로다.

이것이 희귀사건 설계의 핵심 규칙이다.

필요한 것은 표본 크기가 아니라 사건 수다. 상대오차 \(20\%\)면 \(25\)건, \(10\%\)면 \(100\)건, \(5\%\)면 \(400\)건이다. \(p\)는 식에서 사라진다.

실무에서 바로 쓰인다.

• 임상시험에서 "사건 수 기반 설계"를 한다. 환자 수가 아니라 목표 사건 수에 도달할 때까지 추적하며, 21장 생존분석의 표준 설계다.
• A/B 테스트에서 전환율이 낮으면 트래픽이 아니라 전환 건수를 기준으로 실험 기간을 정한다.
• 결함 검사에서 결함률이 \(10^{-5}\)면 상대오차 \(20\%\)에 \(250\)만 개를 검사해야 한다. 현실적으로 불가능하므로 가속 시험이나 물리 모형으로 우회한다.

연습문제 6의 삼의 법칙과 이어진다. 사건이 \(0\)건이면 상대오차가 정의조차 되지 않는다. 사건 수가 적을 때 우리가 할 수 있는 말은 "\(p\)가 얼마인가"가 아니라 "\(p\)는 이보다 작다"뿐이다.

from scipy import stats

n = 100
print(f"  n={n}, k=0 일 때 p 의 95% 상한")
print(f"    삼의 법칙  3/n         = {3/n:.5f}")
print(f"    CP  단측 95%           = {stats.beta.ppf(0.95, 1, n):.5f}")
print(f"    CP  양측 95% 의 상한   = {stats.beta.ppf(0.975, 1, n):.5f}")
print(f"    제프리스 단측 95%      = {stats.beta.ppf(0.95, 0.5, n + 0.5):.5f}")

출력:

  n=100, k=0 일 때 p 의 95% 상한
    삼의 법칙  3/n         = 0.03000
    CP  단측 95%           = 0.02951
    CP  양측 95% 의 상한   = 0.03622
    제프리스 단측 95%      = 0.01898

삼의 법칙과 단측 CP가 거의 같다(\(0.0300\) 대 \(0.0295\)). 연습문제 6에서 본 대로 \((1-p)^n = 0.05\)를 풀면 \(p = 1-0.05^{1/n} \approx -\ln(0.05)/n = 2.996/n\)이기 때문이다.

셋이 다른 이유는 각각 다른 질문에 답하기 때문이다.

양측 CP가 더 넓은 것은 당연하다. 하한 쪽에 \(2.5\%\)를 떼어 두었는데 \(k=0\)이면 하한이 어차피 \(0\)이라, 그 \(2.5\%\)가 낭비된다. 상한만 궁금하다면 단측을 쓰는 것이 맞다.

제프리스가 가장 좁은 것이 흥미롭다. 사전분포가 \(p\)를 \(0\) 쪽으로 약간 끌어당기기 때문이다. \(\text{Beta}(0.5, 0.5)\)는 양 끝에서 밀도가 높은 U자 모양이라, \(k=0\)이라는 자료와 합쳐지면 사후분포가 \(0\) 근처에 더 몰린다.

어느 것을 쓸 것인가.

• 안전성 보고(규제)에는 보수적인 쪽을 쓴다. 상한을 낮게 잡으면 위험을 과소평가하게 되므로 CP가 선호된다.
• 여러 연구를 합치거나 사후 확률을 말하고 싶으면 베이즈가 자연스럽다.
• 어림셈에는 삼의 법칙이면 충분하다. \(n=100\)이면 "\(3\%\) 이하"다.

공통된 교훈은 \(k=0\)이 \(p=0\)을 뜻하지 않는다는 것이다. 세 방법 모두 \(0.019 \sim 0.036\)의 상한을 준다. \(100\)명에게서 부작용이 없었다는 것이 "\(100\)명 중 \(2\)명꼴의 부작용률"과 모순되지 않는다. \(2\%\) 부작용률이라면 \(100\)명 중 아무도 겪지 않을 확률이 \(0.98^{100} = 13\%\)나 된다.

점추정은 둘 다 \(0\)이다. \(\hat p_1 = \hat p_2 = 0\)이므로 위험차는 \(0\), 위험비는 \(0/0\)으로 정의되지 않는다. 그러나 이것이 "위험이 같다"는 뜻은 아니다.

import numpy as np
from scipy import stats

n1 = n2 = 50
# 각 집단의 단측 95% 상한 (삼의 법칙 / 정확)
u = stats.beta.ppf(0.95, 1, n1)
print(f"  각 집단의 p 상한(단측 95%) = {u:.4f}  (삼의 법칙 {3/n1:.4f})")

# 위험차의 구간: 두 p 가 각각 [0, u] 안에 있을 때 차의 범위
print(f"  위험차가 놓일 수 있는 범위 ≈ [{-u:.4f}, {u:.4f}]")

# 피셔 정확검정
table = [[0, n1], [0, n2]]
print(f"  피셔 정확검정 p값 = {stats.fisher_exact(table)[1]:.4f}")

# 0 이 하나라도 있으면 흔히 0.5 를 더한다(하네스-헐 보정)
a, b, c, d = 0.5, n1 + 0.5, 0.5, n2 + 0.5
lor = np.log((a*d)/(b*c))
se = np.sqrt(1/a + 1/b + 1/c + 1/d)
print(f"  보정 오즈비 = {np.exp(lor):.4f}, "
      f"95% CI = [{np.exp(lor-1.96*se):.4f}, {np.exp(lor+1.96*se):.4f}]")

출력:

  각 집단의 p 상한(단측 95%) = 0.0582  (삼의 법칙 0.0600)
  위험차가 놓일 수 있는 범위 ≈ [-0.0582, 0.0582]
  피셔 정확검정 p값 = 1.0000
  보정 오즈비 = 1.0000, 95% CI = [0.0195, 51.3856]

피셔 검정의 \(p\)값이 정확히 \(1\)이다. 자료가 두 집단을 구별할 어떤 증거도 주지 않는다는 뜻이다. 그러나 이것은 "위험이 같다"의 증거가 아니라 "아무 말도 할 수 없다"의 표현이다.

보정 오즈비의 구간이 \([0.02,\ 51.4]\)다. 한쪽 집단의 위험이 다른 쪽의 \(51\)배일 가능성과 \(1/51\)일 가능성이 모두 배제되지 않는다. \(2{,}600\)배에 걸친 구간이며 실질적으로 아무것도 좁히지 못했다.

왜 이런가 — 사건이 없으면 정보가 없다. 연습문제 8에서 본 대로 정밀도는 사건 수가 결정하는데 여기서는 \(0\)건이다. \(n\)을 아무리 키워도 사건이 나오지 않으면 비교의 정밀도는 오르지 않는다. 다만 두 집단의 위험 상한은 함께 내려간다.

실무에서 흔한 오독. "두 군 모두 부작용 \(0\)건, 유의차 없음(\(p=1.0\))"이라는 보고를 "안전성이 동등함이 확인되었다"로 읽으면 안 된다. 올바른 서술은

두 군 모두 사건이 관측되지 않았다. 각 군의 발생률 상한은 \(5.8\%\)(단측 \(95\%\))이며, 이 표본크기로는 군 간 차이를 평가할 수 없다.

동등성을 주장하려면 동등성 검정을 설계해야 한다. 허용 차이를 미리 정하고(예: 위험차 \(2\)%포인트 이내), 그것을 배제할 수 있을 만큼의 표본을 확보해야 한다. 9장에서 다룰 비열등성·동등성 설계의 문제이며, "유의하지 않음"과 "차이 없음"이 다르다는 원칙의 가장 극단적인 사례다. \(\square\)