비율 차이에 대한 이표본 Z-검정¶

import numpy as np
from scipy import stats

def two_proportion_z_test(x1, n1, x2, n2, alternative="two-sided"):
    """두 비율의 차이에 대한 이표본 z 검정.

    Parameters
    ----------
    x1, x2 : int
        Number of successes in each sample.
    n1, n2 : int
        Sample sizes.
    alternative : str
        'two-sided', 'less', or 'greater'.

    Returns
    -------
    z_stat : float
        Test statistic.
    p_value : float
        P-value.
    """
    p1_hat = x1 / n1
    p2_hat = x2 / n2
    # H0가 "두 비율이 같다"이므로 그 공통값을 전체를 합쳐 추정한다.
    # 신뢰구간을 만들 때는 합동하지 않는다. 목적이 다르기 때문이다.
    p_hat = (x1 + x2) / (n1 + n2)

    se = np.sqrt(p_hat * (1 - p_hat) * (1 / n1 + 1 / n2))
    z_stat = (p1_hat - p2_hat) / se

    if alternative == "two-sided":
        p_value = 2 * (1 - stats.norm.cdf(abs(z_stat)))
    elif alternative == "less":
        p_value = stats.norm.cdf(z_stat)
    elif alternative == "greater":
        p_value = 1 - stats.norm.cdf(z_stat)
    else:
        raise ValueError("alternative must be 'two-sided', 'less', or 'greater'")

    return z_stat, p_value


# 보기 1 의 A/B 검정
z, p = two_proportion_z_test(45, 500, 58, 480, alternative="two-sided")
print(f"z = {z:.3f}, p-value = {p:.3f}")

z = -1.573, p-value = 0.116

x1, n1, x2, n2 = 45, 500, 58, 480
p1, p2 = x1 / n1, x2 / n2
pb = (x1 + x2) / (n1 + n2)
d = p1 - p2
se_p = np.sqrt(pb * (1 - pb) * (1 / n1 + 1 / n2))
se_u = np.sqrt(p1 * (1 - p1) / n1 + p2 * (1 - p2) / n2)
print(f"p1 = {p1:.7f}   p2 = {p2:.7f}   차이 = {d:.7f}")
print(f"합동 p = {pb:.7f}   합동 SE = {se_p:.7f}")
print(f"손계산의 차이 -0.031 →  z = {-0.031 / se_p:.6f}")
print(f"정확한 차이          →  z = {d / se_p:.6f}"
      f"   p = {2 * stats.norm.sf(abs(d / se_p)):.6f}")

print("\n자료를 k 배로 늘릴 때")
print("   k      z          p            sqrt(k)*z_1")
z1 = d / se_p
for k in (1, 10, 100):
    z, _ = two_proportion_z_test(k * x1, k * n1, k * x2, k * n2)
    p = 2 * stats.norm.sf(abs(z))
    print(f"{k:4d}  {z:9.6f}  {p:.3e}   {np.sqrt(k) * z1:9.6f}")

print(f"\n비합동 SE = {se_u:.7f}   |z| = {abs(d / se_u):.6f}")
print(f"합동   SE = {se_p:.7f}   |z| = {abs(d / se_p):.6f}")
print(f"→ 이 불균형 자료에서는 {'합동' if se_p > se_u else '비합동'} 쪽 SE 가 크다")

# 균형 설계에서의 항등식: SE_p^2 - SE_u^2 = d^2/(2n)
rng = np.random.default_rng(0)
worst = 0.0
for _ in range(20_000):
    n = int(rng.integers(20, 2000))
    a, b = int(rng.integers(1, n)), int(rng.integers(1, n))
    q1, q2, qb = a / n, b / n, (a + b) / (2 * n)
    sp2 = qb * (1 - qb) * (2 / n)
    su2 = q1 * (1 - q1) / n + q2 * (1 - q2) / n
    worst = max(worst, abs(sp2 - su2 - (q1 - q2) ** 2 / (2 * n)))
print(f"\n균형 항등식 SE_p^2 - SE_u^2 = d^2/(2n) 의 최대 오차 = {worst:.2e}"
      f"  (무작위 표 20000 개)")

# 총 980 명을 어떻게 나누면 부호가 뒤집히는가
print("\n  n1   n2    합동 SE   비합동 SE   큰 쪽")
for m in (400, 450, 490, 495, 500, 560):
    k = 980 - m
    PB = (m * p1 + k * p2) / 980
    sp = np.sqrt(PB * (1 - PB) * (1 / m + 1 / k))
    su = np.sqrt(p1 * (1 - p1) / m + p2 * (1 - p2) / k)
    print(f"{m:4d} {k:4d}   {sp:.7f}  {su:.7f}   {'합동' if sp > su else '비합동'}")

# z_p^2 은 2x2 피어슨 카이제곱과 같고 N 을 넘지 못한다
chi = stats.chi2_contingency([[x1, n1 - x1], [x2, n2 - x2]], correction=False)
print(f"\nz_p^2 = {z1**2:.9f}   피어슨 카이제곱 = {chi.statistic:.9f}")
worst = 0.0
for _ in range(20_000):
    a, b = int(rng.integers(5, 500)), int(rng.integers(5, 500))
    u, v = int(rng.integers(0, a + 1)), int(rng.integers(0, b + 1))
    if u + v in (0, a + b):
        continue
    P1, P2, PB = u / a, v / b, (u + v) / (a + b)
    z = (P1 - P2) / np.sqrt(PB * (1 - PB) * (1 / a + 1 / b))
    worst = max(worst, z**2 / (a + b))
print(f"z_p^2 / N 의 최대값 = {worst:.9f}  (1 을 넘지 못한다)")

p1 = 0.0900000   p2 = 0.1208333   차이 = -0.0308333
합동 p = 0.1051020   합동 SE = 0.0195975
손계산의 차이 -0.031 →  z = -1.581838
정확한 차이          →  z = -1.573333   p = 0.115642

자료를 k 배로 늘릴 때
   k      z          p            sqrt(k)*z_1
   1  -1.573333  1.156e-01   -1.573333
  10  -4.975317  6.514e-07   -4.975317
 100  -15.733333  8.938e-56   -15.733333

비합동 SE = 0.0196244   |z| = 1.571171
합동   SE = 0.0195975   |z| = 1.573333
→ 이 불균형 자료에서는 비합동 쪽 SE 가 크다

균형 항등식 SE_p^2 - SE_u^2 = d^2/(2n) 의 최대 오차 = 6.94e-18  (무작위 표 20000 개)

  n1   n2    합동 SE   비합동 SE   큰 쪽
 400  580   0.0201930  0.0196954   합동
 450  530   0.0197882  0.0195560   합동
 490  490   0.0196192  0.0195945   합동
 495  485   0.0196073  0.0196084   비합동
 500  480   0.0195975  0.0196244   비합동
 560  420   0.0196385  0.0199796   비합동

z_p^2 = 2.475377593   피어슨 카이제곱 = 2.475377593
z_p^2 / N 의 최대값 = 1.000000000  (1 을 넘지 못한다)

import numpy as np
from scipy import stats

z95 = stats.norm.ppf(0.975)
n1, k1 = 30, 3         # 10%
n2, k2 = 100, 1        # 1%
p1, p2 = k1 / n1, k2 / n2
pp = (k1 + k2) / (n1 + n2)
se_pooled = np.sqrt(pp * (1 - pp) * (1 / n1 + 1 / n2))
se_unpool = np.sqrt(p1 * (1 - p1) / n1 + p2 * (1 - p2) / n2)

print(f"p̂1 = {p1:.4f} ({k1}/{n1}),  p̂2 = {p2:.4f} ({k2}/{n2}),  "
      f"합동 p̂ = {pp:.4f}\n")
print(f"합동   SE={se_pooled:.5f}  z={(p1 - p2) / se_pooled:.4f}  "
      f"p={2 * stats.norm.sf(abs((p1 - p2) / se_pooled)):.4f}  ← 기각")
print(f"비합동 SE={se_unpool:.5f}  z={(p1 - p2) / se_unpool:.4f}  "
      f"p={2 * stats.norm.sf(abs((p1 - p2) / se_unpool)):.4f}")
print(f"Wald 95% CI ({p1 - p2 - z95 * se_unpool:.5f}, "
      f"{p1 - p2 + z95 * se_unpool:.5f})  ← 0 을 포함")

def wilson(k, n, z=z95):
    p = k / n
    d = 1 + z * z / n
    c = (p + z * z / (2 * n)) / d
    h = z * np.sqrt(p * (1 - p) / n + z * z / (4 * n * n)) / d
    return c - h, c + h

l1, u1 = wilson(k1, n1)
l2, u2 = wilson(k2, n2)
lo = p1 - p2 - z95 * np.sqrt(l1 * (1 - l1) / n1 + u2 * (1 - u2) / n2)
hi = p1 - p2 + z95 * np.sqrt(u1 * (1 - u1) / n1 + l2 * (1 - l2) / n2)
print(f"뉴컴 95% CI ({lo:.5f}, {hi:.5f})  ← 0 을 포함하지 않음")

print(f"\n피셔 정확검정 p = "
      f"{stats.fisher_exact([[k1, n1 - k1], [k2, n2 - k2]]).pvalue:.4f}")
print(f"기대도수의 최솟값 = "
      f"{min(n1 * pp, n1 * (1 - pp), n2 * pp, n2 * (1 - pp)):.2f}")

p̂1 = 0.1000 (3/30),  p̂2 = 0.0100 (1/100),  합동 p̂ = 0.0308

합동   SE=0.03595  z=2.5036  p=0.0123  ← 기각
비합동 SE=0.05567  z=1.6167  p=0.1059
Wald 95% CI (-0.01911, 0.19911)  ← 0 을 포함
뉴컴 95% CI (0.01090, 0.24643)  ← 0 을 포함하지 않음

피셔 정확검정 p = 0.0382
기대도수의 최솟값 = 0.92

import numpy as np
from scipy import stats

z95 = stats.norm.ppf(0.975)

def compare(k1, n1, k2, n2, label):
    p1, p2 = k1 / n1, k2 / n2
    pp = (k1 + k2) / (n1 + n2)
    z_pooled = (p1 - p2) / np.sqrt(pp * (1 - pp) * (1 / n1 + 1 / n2))
    se_u = np.sqrt(p1 * (1 - p1) / n1 + p2 * (1 - p2) / n2)

    rr = p1 / p2
    se_lrr = np.sqrt(1 / k1 - 1 / n1 + 1 / k2 - 1 / n2)
    or_ = (k1 / (n1 - k1)) / (k2 / (n2 - k2))
    se_lor = np.sqrt(1 / k1 + 1 / (n1 - k1) + 1 / k2 + 1 / (n2 - k2))

    print(f"{label}:  p̂1={p1:.4f}  p̂2={p2:.4f}")
    print(f"  합동 z = {z_pooled:.4f},  "
          f"p = {2 * stats.norm.sf(abs(z_pooled)):.4f}")
    print(f"  위험차 {p1 - p2:+.5f}  95% CI "
          f"({p1 - p2 - z95 * se_u:+.5f}, {p1 - p2 + z95 * se_u:+.5f})")
    print(f"  RR     {rr:.4f}     95% CI "
          f"({np.exp(np.log(rr) - z95 * se_lrr):.4f}, "
          f"{np.exp(np.log(rr) + z95 * se_lrr):.4f})")
    print(f"  OR     {or_:.4f}     95% CI "
          f"({np.exp(np.log(or_) - z95 * se_lor):.4f}, "
          f"{np.exp(np.log(or_) + z95 * se_lor):.4f})")

compare(60, 200, 54, 180, "연습문제 1 (금연)")
print()
compare(240, 2000, 270, 2000, "연습문제 2 (A/B 전환)")

연습문제 1 (금연):  p̂1=0.3000  p̂2=0.3000
  합동 z = 0.0000,  p = 1.0000
  위험차 +0.00000  95% CI (-0.09228, +0.09228)
  RR     1.0000     95% CI (0.7352, 1.3601)
  OR     1.0000     95% CI (0.6444, 1.5518)

연습문제 2 (A/B 전환):  p̂1=0.1200  p̂2=0.1350
  합동 z = -1.4222,  p = 0.1550
  위험차 -0.01500  95% CI (-0.03567, +0.00567)
  RR     0.8889     95% CI (0.7556, 1.0457)
  OR     0.8737     95% CI (0.7254, 1.0525)

import numpy as np
from scipy import stats

z95 = stats.norm.ppf(0.975)
k_v, n_v = 8, 21_720        # 백신군: 21,720명 중 8명 발병
k_p, n_p = 162, 21_728      # 위약군: 21,728명 중 162명 발병
p_v, p_p = k_v / n_v, k_p / n_p
print(f"백신군 {k_v}/{n_v:,} = {p_v * 100:.4f}%")
print(f"위약군 {k_p}/{n_p:,} = {p_p * 100:.4f}%")

# ── 방법 ①: log RR 의 델타법 구간
rr = p_v / p_p
se = np.sqrt(1 / k_v - 1 / n_v + 1 / k_p - 1 / n_p)
lo, hi = np.exp(np.log(rr) - z95 * se), np.exp(np.log(rr) + z95 * se)
print(f"\n① 로그 RR 델타법")
print(f"   RR = {rr:.5f}  95% CI ({lo:.5f}, {hi:.5f})")
print(f"   VE = {(1 - rr) * 100:.2f}%  95% CI "
      f"({(1 - hi) * 100:.2f}%, {(1 - lo) * 100:.2f}%)")

# ── 방법 ②: 전체 사례 수를 고정한 조건부 이항
n_case = k_v + k_p
theta = k_v / n_case                   # 사례 중 백신군의 비율
ratio = n_v / n_p                      # 배정 비율

def wilson(k, n, z=z95):
    p = k / n
    d = 1 + z * z / n
    c = (p + z * z / (2 * n)) / d
    h = z * np.sqrt(p * (1 - p) / n + z * z / (4 * n * n)) / d
    return c - h, c + h

def ve(th):
    return 1 - (th / (1 - th)) / ratio

t_lo, t_hi = wilson(k_v, n_case)
print(f"\n② 조건부 이항 (사례 {n_case}건 중 백신군 {k_v}건)")
print(f"   θ̂ = {theta:.5f}  윌슨 95% CI ({t_lo:.5f}, {t_hi:.5f})")
print(f"   VE = {ve(theta) * 100:.2f}%  95% CI "
      f"({ve(t_hi) * 100:.2f}%, {ve(t_lo) * 100:.2f}%)")

rd = p_p - p_v
print(f"\n위험차 {rd * 1e5:.2f}건 / 10만명")
print(f"1명의 발병을 막는 데 필요한 접종 수 = {1 / rd:.0f}명")
pool = (k_v + k_p) / (n_v + n_p)
zs = (p_v - p_p) / np.sqrt(pool * (1 - pool) * (1 / n_v + 1 / n_p))
print(f"합동 z = {zs:.4f},  p = {2 * stats.norm.sf(abs(zs)):.3g}")

백신군 8/21,720 = 0.0368%
위약군 162/21,728 = 0.7456%

① 로그 RR 델타법
   RR = 0.04940  95% CI (0.02430, 0.10044)
   VE = 95.06%  95% CI (89.96%, 97.57%)

② 조건부 이항 (사례 170건 중 백신군 8건)
   θ̂ = 0.04706  윌슨 95% CI (0.02404, 0.09010)
   VE = 95.06%  95% CI (90.09%, 97.54%)

위험차 708.75건 / 10만명
1명의 발병을 막는 데 필요한 접종 수 = 141명
합동 z = -11.8320,  p = 2.67e-32