이표본 비율 검정¶

import math
from scipy.stats import norm

def test_diff_two_props(k1, n1, k2, n2, delta0=0.0,
                        method="pooled", alt="two-sided", alpha=0.05):
    """귀무가설 H0: p1 - p2 = delta0.

    delta0=0이고 method='pooled'이면 합동 표준오차를 쓴다.
    돌려주는 값은 (z, p, 기각 여부, 이름).
    """
    p1, p2 = k1 / n1, k2 / n2
    d_hat = p1 - p2
    if delta0 == 0.0 and method == "pooled":
        # H0가 "두 비율이 같다"이면 그 공통값의 최선의 추정은 전체를 합친 비율이다.
        # 검정은 H0 아래의 분포를 쓰므로 여기서 합동하는 것이 옳다.
        p_pool = (k1 + k2) / (n1 + n2)
        se = math.sqrt(p_pool * (1 - p_pool) * (1 / n1 + 1 / n2))
        label = "pooled z-test"
    else:
        # delta0이 0이 아니면 "두 비율이 같다"는 가정이 없으므로 합동할 수 없다.
        # 신뢰구간에서 쓰는 것과 같은 표준오차다.
        se = math.sqrt(p1 * (1 - p1) / n1 + p2 * (1 - p2) / n2)
        label = "wald z-test"

    z = (d_hat - delta0) / se
    if alt == "two-sided":
        p = 2 * min(norm.cdf(z), 1 - norm.cdf(z))
    elif alt == "less":
        p = norm.cdf(z)
    else:
        p = 1 - norm.cdf(z)
    return z, p, (p < alpha), label

import numpy as np


def se_pair(k1, n1, k2, n2):
    """합동 SE^2, Wald SE^2, 그리고 닫힌 꼴이 예측하는 차."""
    p1, p2 = k1 / n1, k2 / n2
    N = n1 + n2
    p_pool = (k1 + k2) / N
    se_p2 = p_pool * (1 - p_pool) * (1 / n1 + 1 / n2)
    se_w2 = p1 * (1 - p1) / n1 + p2 * (1 - p2) / n2
    pred = (p1 - p2) * ((n1 - n2) * (1 - p1 - p2)
                        + n1 * n2 * (p1 - p2) / N) / (n1 * n2)
    return se_p2, se_w2, pred


rng = np.random.default_rng(3)
worst = 0.0
for _ in range(200_000):
    a, b = int(rng.integers(2, 301)), int(rng.integers(2, 301))
    u, v = int(rng.integers(0, a + 1)), int(rng.integers(0, b + 1))
    if u + v in (0, a + b):
        continue                      # 합동 비율이 0 또는 1 이면 SE = 0
    sp2, sw2, pred = se_pair(u, a, v, b)
    worst = max(worst, abs(sp2 - sw2 - pred))
print(f"항등식 검사 20만 개: 절대오차의 최대값 = {worst:.3e}")

print("\nn1 = 80, n2 = 60 에서 같은 차이를 위로 옮기면 부호가 뒤집힌다")
print("  k1  k2   p1-p2   p1+p2    SE_p^2 - SE_w^2   큰 쪽")
for k1, k2 in ((30, 18), (44, 29), (48, 32), (50, 33), (52, 34)):
    sp2, sw2, pred = se_pair(k1, 80, k2, 60)
    print(f" {k1:3d} {k2:3d}  {k1/80 - k2/60:6.4f}  {k1/80 + k2/60:6.4f}"
          f"   {sp2 - sw2:+.3e}      {'합동' if sp2 > sw2 else 'Wald'}")
print(f"  부호가 바뀌는 경계: p1+p2 = "
      f"{1 + 80 * 60 / (140 * 20) * 0.075:.6f}  (p1-p2 = 0.075 고정)")

항등식 검사 20만 개: 절대오차의 최대값 = 6.245e-17

n1 = 80, n2 = 60 에서 같은 차이를 위로 옮기면 부호가 뒤집힌다
  k1  k2   p1-p2   p1+p2    SE_p^2 - SE_w^2   큰 쪽
  30  18  0.0750  0.6750   +1.417e-04      합동
  44  29  0.0667  1.0333   +2.249e-05      합동
  48  32  0.0667  1.1333   -5.291e-06      Wald
  50  33  0.0750  1.1750   -1.451e-05      Wald
  52  34  0.0833  1.2167   -2.563e-05      Wald
  부호가 바뀌는 경계: p1+p2 = 1.128571  (p1-p2 = 0.075 고정)

z, p, reject, label = test_diff_two_props(
    k1=30, n1=80, k2=18, n2=60, delta0=0.0, method="pooled"
)
print(label, "z:", z, "p:", p, "reject:", reject)

# 같은 자료를 합동하지 않은 Wald 표준오차로 계산하면 어떻게 다른가?
z2, p2, reject2, label2 = test_diff_two_props(
    k1=30, n1=80, k2=18, n2=60, delta0=0.0, method="wald"
)
print(label2, "z:", z2, "p:", p2, "reject:", reject2)

pooled z-test z: 0.9251909321159419 p: 0.3548665994806586 reject: False
wald z-test z: 0.9353331581027243 p: 0.3496166322568355 reject: False

import numpy as np
from scipy import stats

k1, n1, k2, n2 = 30, 80, 18, 60
p1, p2 = k1 / n1, k2 / n2
N = n1 + n2
p_pool = (k1 + k2) / N
se_p2 = p_pool * (1 - p_pool) * (1 / n1 + 1 / n2)
se_w2 = p1 * (1 - p1) / n1 + p2 * (1 - p2) / n2
pred = (p1 - p2) * ((n1 - n2) * (1 - p1 - p2)
                    + n1 * n2 * (p1 - p2) / N) / (n1 * n2)
print(f"p1 = {p1}, p2 = {p2}, p_pool = {p_pool:.9f}")
print(f"SE_pool = {math.sqrt(se_p2):.9f}   SE_wald = {math.sqrt(se_w2):.9f}")
print(f"SE_p^2 - SE_w^2 = {se_p2 - se_w2:.3e}   닫힌 꼴 = {pred:.3e}")
z_pool = (p1 - p2) / math.sqrt(se_p2)
chi = stats.chi2_contingency([[k1, n1 - k1], [k2, n2 - k2]],
                             correction=False).statistic
print(f"z_pool^2 = {z_pool**2:.12f}   피어슨 카이제곱 = {chi:.12f}")


def exact_level(n1, n2, p, alpha=0.05):
    """p1 = p2 = p 일 때 합동 z 검정의 **정확한** 제1종 오류율.

    (k1, k2) 의 (n1+1)(n2+1) 칸을 모두 훑어 기각 칸의 이항확률을 더한다.
    근사도 모의실험도 아니다.
    """
    zc = stats.norm.ppf(1 - alpha / 2)
    K1, K2 = np.arange(n1 + 1), np.arange(n2 + 1)
    pr = stats.binom.pmf(K1, n1, p)[:, None] * stats.binom.pmf(K2, n2, p)[None, :]
    P1, P2 = K1[:, None] / n1, K2[None, :] / n2
    PP = (K1[:, None] + K2[None, :]) / (n1 + n2)
    with np.errstate(divide="ignore", invalid="ignore"):
        z = (P1 - P2) / np.sqrt(PP * (1 - PP) * (1 / n1 + 1 / n2))
    rej = np.abs(z) > zc
    rej[~np.isfinite(z)] = False        # k1+k2 가 0 또는 n1+n2 인 두 칸
    return float((pr * rej).sum())


print(f"\n정확한 수준 (n1 = {n1}, n2 = {n2}, 명목 0.05)")
print("     p    실제 수준     n2*min(p,1-p)")
for p in (0.01, 0.02, 0.03, 0.05, 0.0833, 0.1, 0.2, 0.343, 0.5):
    print(f"  {p:.4f}   {exact_level(n1, n2, p):.6f}      {n2 * min(p, 1 - p):6.2f}")

ps = np.linspace(0.002, 0.998, 499)
lv = np.array([exact_level(n1, n2, p) for p in ps])
i = int(lv.argmax())
mid = (ps >= 0.03) & (ps <= 0.97)
print(f"\n격자 499 점의 최대 = {lv[i]:.6f} (p = {ps[i]:.3f})")
print(f"0.03 <= p <= 0.97 구간의 범위 = [{lv[mid].min():.6f}, {lv[mid].max():.6f}]")
print(f"0.05 를 넘는 점의 비율 = {float((lv > 0.05).mean()):.3f}")

p1 = 0.375, p2 = 0.3, p_pool = 0.342857143
SE_pool = 0.081064348   SE_wald = 0.080185332
SE_p^2 - SE_w^2 = 1.417e-04   닫힌 꼴 = 1.417e-04
z_pool^2 = 0.855978260870   피어슨 카이제곱 = 0.855978260870

정확한 수준 (n1 = 80, n2 = 60, 명목 0.05)
     p    실제 수준     n2*min(p,1-p)
  0.0100   0.010866        0.60
  0.0200   0.032648        1.20
  0.0300   0.047253        1.80
  0.0500   0.053474        3.00
  0.0833   0.050844        5.00
  0.1000   0.049784        6.00
  0.2000   0.047094       12.00
  0.3430   0.049443       20.58
  0.5000   0.049180       30.00

격자 499 점의 최대 = 0.053604 (p = 0.046)
0.03 <= p <= 0.97 구간의 범위 = [0.045280, 0.053604]
0.05 를 넘는 점의 비율 = 0.365

import numpy as np
from scipy import stats

Z = stats.norm.ppf(0.975)

def wilson(k, n, z=Z):
    p = k / n
    d = 1 + z * z / n
    centre = (p + z * z / (2 * n)) / d
    half = z * np.sqrt(p * (1 - p) / n + z * z / (4 * n * n)) / d
    return centre - half, centre + half

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

def newcombe(k1, n1, k2, n2, z=Z):
    l1, u1 = wilson(k1, n1, z)
    l2, u2 = wilson(k2, n2, z)
    p1, p2 = k1 / n1, k2 / n2
    lo = p1 - p2 - z * np.sqrt(l1 * (1 - l1) / n1 + u2 * (1 - u2) / n2)
    hi = p1 - p2 + z * np.sqrt(u1 * (1 - u1) / n1 + l2 * (1 - l2) / n2)
    return lo, hi

for label, args in [("연습문제 1 (45/200 대 30/200)", (45, 200, 30, 200)),
                    ("연습문제 3 (120/5000 대 95/5000)", (120, 5000, 95, 5000))]:
    lo_w, hi_w = wald(*args)
    lo_n, hi_n = newcombe(*args)
    print(f"{label}")
    print(f"  Wald     ({lo_w:8.5f}, {hi_w:8.5f})")
    print(f"  뉴컴     ({lo_n:8.5f}, {hi_n:8.5f})")

연습문제 1 (45/200 대 30/200)
  Wald     (-0.00115,  0.15115)
  뉴컴     (-0.00172,  0.15099)
연습문제 3 (120/5000 대 95/5000)
  Wald     (-0.00068,  0.01068)
  뉴컴     (-0.00070,  0.01076)

rng = np.random.default_rng(313)
print(f"{'p1':>5s} {'p2':>5s} {'n':>5s} {'Wald':>8s} {'뉴컴':>8s}")
for p1, p2, n in [(0.50, 0.50, 20), (0.10, 0.10, 20), (0.02, 0.02, 50),
                  (0.30, 0.10, 30), (0.05, 0.02, 100)]:
    M = 20_000
    k1, k2 = rng.binomial(n, p1, M), rng.binomial(n, p2, M)
    a = b = 0
    for i in range(M):
        lo, hi = wald(k1[i], n, k2[i], n)
        a += lo <= p1 - p2 <= hi
        lo, hi = newcombe(k1[i], n, k2[i], n)
        b += lo <= p1 - p2 <= hi
    print(f"{p1:5.2f} {p2:5.2f} {n:5d} {a / M:8.4f} {b / M:8.4f}")

   p1    p2     n     Wald       뉴컴
 0.50  0.50    20   0.9210   0.9593
 0.10  0.10    20   0.9591   0.9820
 0.02  0.02    50   0.9850   0.9976
 0.30  0.10    30   0.9352   0.9534
 0.05  0.02   100   0.9345   0.9679

import numpy as np
from scipy import stats

rng = np.random.default_rng(77)
print(f"{'군당 n':>7s} {'p':>5s} {'합동 z':>8s} {'Fisher':>8s} {'χ²(보정)':>10s}")
for n, p in [(10, 0.30), (20, 0.30), (20, 0.10), (50, 0.05), (100, 0.02)]:
    M = 10_000
    k1, k2 = rng.binomial(n, p, M), rng.binomial(n, p, M)
    a = b = c = 0
    for i in range(M):
        x, y = k1[i], k2[i]
        pp = (x + y) / (2 * n)
        if 0 < pp < 1:                       # 둘 다 0이거나 둘 다 n이면 정의 불가
            se = np.sqrt(pp * (1 - pp) * 2 / n)
            a += 2 * stats.norm.sf(abs((x / n - y / n) / se)) < 0.05
        b += stats.fisher_exact([[x, n - x], [y, n - y]]).pvalue < 0.05
        tbl = np.array([[x, n - x], [y, n - y]])
        if tbl.sum(0).min() > 0 and tbl.sum(1).min() > 0:
            c += stats.chi2_contingency(tbl, correction=True)[1] < 0.05
    print(f"{n:7d} {p:5.2f} {a / M:8.4f} {b / M:8.4f} {c / M:10.4f}")

   군당 n     p     합동 z   Fisher     χ²(보정)
     10  0.30   0.0368   0.0119     0.0119
     20  0.30   0.0551   0.0262     0.0229
     20  0.10   0.0392   0.0098     0.0033
     50  0.05   0.0409   0.0064     0.0064
    100  0.02   0.0419   0.0049     0.0049

import numpy as np

def summarise(k1, n1, k2, n2, label):
    p1, p2 = k1 / n1, k2 / n2
    rd = p1 - p2
    rr = p1 / p2
    or_ = (p1 / (1 - p1)) / (p2 / (1 - p2))
    print(f"{label}")
    print(f"  p1={p1:.4f}  p2={p2:.4f}")
    print(f"  위험차 RD = {rd:+.4f}   위험비 RR = {rr:.3f}   "
          f"오즈비 OR = {or_:.3f}   NNT = {1 / abs(rd):.1f}")

summarise(45, 200, 30, 200, "① 임상시험 (연습문제 1)")
summarise(120, 5000, 95, 5000, "② A/B 검정 (연습문제 3)")
summarise(2, 1000, 1, 1000, "③ 드문 부작용")
summarise(900, 1000, 800, 1000, "④ 흔한 결과")

① 임상시험 (연습문제 1)
  p1=0.2250  p2=0.1500
  위험차 RD = +0.0750   위험비 RR = 1.500   오즈비 OR = 1.645   NNT = 13.3
② A/B 검정 (연습문제 3)
  p1=0.0240  p2=0.0190
  위험차 RD = +0.0050   위험비 RR = 1.263   오즈비 OR = 1.270   NNT = 200.0
③ 드문 부작용
  p1=0.0020  p2=0.0010
  위험차 RD = +0.0010   위험비 RR = 2.000   오즈비 OR = 2.002   NNT = 1000.0
④ 흔한 결과
  p1=0.9000  p2=0.8000
  위험차 RD = +0.1000   위험비 RR = 1.125   오즈비 OR = 2.250   NNT = 10.0

import numpy as np
from scipy import stats

def n_two_props(p1, p2, power=0.80, alpha=0.05):
    za, zb = stats.norm.ppf(1 - alpha / 2), stats.norm.ppf(power)
    pbar = (p1 + p2) / 2
    num = (za * np.sqrt(2 * pbar * (1 - pbar))
           + zb * np.sqrt(p1 * (1 - p1) + p2 * (1 - p2)))**2
    return int(np.ceil(num / (p1 - p2)**2))

def n_arcsine(p1, p2, power=0.80, alpha=0.05):
    """코헨의 h 를 쓰는 아크사인 변환 근사."""
    h = 2 * np.arcsin(np.sqrt(p1)) - 2 * np.arcsin(np.sqrt(p2))
    za, zb = stats.norm.ppf(1 - alpha / 2), stats.norm.ppf(power)
    return int(np.ceil(2 * ((za + zb) / h)**2))

print(f"{'p1':>6s} {'p2':>6s} {'차이':>7s} {'위험비':>7s} "
      f"{'군당 n':>8s} {'arcsin':>8s}")
for p1, p2 in [(0.50, 0.40), (0.55, 0.45), (0.30, 0.20),
               (0.10, 0.05), (0.05, 0.025), (0.02, 0.01)]:
    print(f"{p1:6.3f} {p2:6.3f} {p1 - p2:7.3f} {p1 / p2:7.2f} "
          f"{n_two_props(p1, p2):8d} {n_arcsine(p1, p2):8d}")

    p1     p2      차이     위험비     군당 n   arcsin
 0.500  0.400   0.100    1.25      388      388
 0.550  0.450   0.100    1.22      392      392
 0.300  0.200   0.100    1.50      294      292
 0.100  0.050   0.050    2.00      435      424
 0.050  0.025   0.025    2.00      906      882
 0.020  0.010   0.010    2.00     2319     2254

import numpy as np
from scipy import stats

rng = np.random.default_rng(505)

def peek_experiment(peeks, N=5_000, p=0.05, M=5_000):
    """두 변형의 전환율이 실제로 같은데, peeks 번 중간 확인한다."""
    points = [int(N * (i + 1) / peeks) for i in range(peeks)]
    hit = 0
    for _ in range(M):
        a, b = rng.random(N) < p, rng.random(N) < p
        for t in points:
            k1, k2 = a[:t].sum(), b[:t].sum()
            pp = (k1 + k2) / (2 * t)
            if pp <= 0 or pp >= 1:
                continue
            se = np.sqrt(pp * (1 - pp) * 2 / t)
            if 2 * stats.norm.sf(abs((k1 / t - k2 / t) / se)) < 0.05:
                hit += 1
                break                     # 유의해지자마자 중단
    return hit / M

for k in [1, 2, 5, 10, 20]:
    print(f"중간 확인 {k:2d}회 → 제1종 오류율 {peek_experiment(k):.4f}")

중간 확인  1회 → 제1종 오류율 0.0516
중간 확인  2회 → 제1종 오류율 0.0830
중간 확인  5회 → 제1종 오류율 0.1418
중간 확인 10회 → 제1종 오류율 0.1924
중간 확인 20회 → 제1종 오류율 0.2382