p₁ − p₂의 신뢰구간¶

import numpy as np
import scipy.stats as stats

n1, n2 = 200, 250
x1, x2 = 120, 130
confidence_level = 0.95

p1 = x1 / n1
p2 = x2 / n2

# 두 표본이 독립이므로 분산이 더해진다.
# 여기서는 두 비율을 **따로** 추정해 넣는다. 검정에서 쓰는 합동비율은
# "두 비율이 같다"는 귀무가설 아래의 계산이라 신뢰구간에는 맞지 않는다.
standard_error = np.sqrt((p1 * (1 - p1) / n1) + (p2 * (1 - p2) / n2))
z_critical = stats.norm.ppf(1 - (1 - confidence_level) / 2)
margin_of_error = z_critical * standard_error

confidence_interval = ((p1 - p2) - margin_of_error, (p1 - p2) + margin_of_error)
print(f"{confidence_interval = }")

confidence_interval = (-0.011897024279429055, 0.171897024279429)

import numpy as np
from scipy import stats

n_1 = 120  # north
n_2 = 140  # south
p_1_hat = 54 / n_1
p_2_hat = 77 / n_2

confidence_level = 0.90
alpha = 1 - confidence_level
# 앞에서는 ppf(1 - alpha/2)를 썼고 여기서는 -ppf(alpha/2)를 썼다.
# 표준정규가 0을 중심으로 대칭이라 두 값이 같다.
z_star = -stats.norm().ppf(alpha / 2)
margin_of_error = z_star * np.sqrt(
    p_1_hat * (1 - p_1_hat) / n_1 + p_2_hat * (1 - p_2_hat) / n_2
)
print(f"90% CI: {p_1_hat - p_2_hat:.4f} ± {margin_of_error:.4f}")

90% CI: -0.1000 ± 0.1018

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

Wald 는 식이 가장 간단하지만 비율이 0 이나 1 에 가까우면 포함확률이
명목수준에 한참 못 미친다. Newcombe 는 그 약점을 고친 방법이다.
"""

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

rng_seed = 42        # 아래 그림을 재현하려면 고정한다
n_simulations = 100
n1, n2 = 50, 40
p1_true, p2_true = 0.60, 0.50
alpha = 0.05
method = "newcombe"  # 'newcombe' | 'wald' | 'cp'


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

    delta_true = p1_true - p2_true
    z = norm.ppf(1 - alpha / 2.0)
    lowers = np.empty(n_simulations)
    uppers = np.empty(n_simulations)
    centers = np.empty(n_simulations)

    for i in range(n_simulations):
        k1 = np.random.binomial(n1, p1_true)
        k2 = np.random.binomial(n2, p2_true)
        p1hat, p2hat = k1 / n1, k2 / n2
        centers[i] = p1hat - p2hat

        if method == "wald":
            se = np.sqrt(p1hat * (1 - p1hat) / n1 + p2hat * (1 - p2hat) / n2)
            lo, hi = centers[i] - z * se, centers[i] + z * se
        elif method == "newcombe":
            # 집단마다 Wilson 구간을 구한다.
            denom1 = 1 + z**2 / n1
            center1 = (p1hat + z**2 / (2 * n1)) / denom1
            half1 = z * np.sqrt(p1hat * (1 - p1hat) / n1 + z**2 / (4 * n1**2)) / denom1
            L1, U1 = center1 - half1, center1 + half1
            denom2 = 1 + z**2 / n2
            center2 = (p2hat + z**2 / (2 * n2)) / denom2
            half2 = z * np.sqrt(p2hat * (1 - p2hat) / n2 + z**2 / (4 * n2**2)) / denom2
            L2, U2 = center2 - half2, center2 + half2
            # 두 구간을 **제곱합**으로 합친다. L1 - U2 처럼 양끝을 그냥 빼면
            # 두 집단이 동시에 최악으로 어긋나는 경우를 가정하는 셈이라
            # 구간이 지나치게 넓어진다(아래 설명 참조).
            lo = (p1hat - p2hat) - np.sqrt((p1hat - L1)**2 + (U2 - p2hat)**2)
            hi = (p1hat - p2hat) + np.sqrt((U1 - p1hat)**2 + (p2hat - L2)**2)
        elif method == "cp":
            L1 = 0.0 if k1 == 0 else beta.ppf(alpha / 2.0, k1, n1 - k1 + 1)
            U1 = 1.0 if k1 == n1 else beta.ppf(1 - alpha / 2.0, k1 + 1, n1 - k1)
            L2 = 0.0 if k2 == 0 else beta.ppf(alpha / 2.0, k2, n2 - k2 + 1)
            U2 = 1.0 if k2 == n2 else beta.ppf(1 - alpha / 2.0, k2 + 1, n2 - k2)
            lo, hi = L1 - U2, U1 - L2

        lowers[i] = max(-1.0, lo)
        uppers[i] = min(1.0, hi)

    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, color="r")
    ax.set_title(
        f"{n_simulations} Δ=p1−p2 CIs ({method.title()}) | 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("Δ = p1 − p2")
    plt.tight_layout()
    plt.show()


if __name__ == "__main__":
    main()

import numpy as np
from scipy import stats

z = stats.norm.ppf(0.975)

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

def newcombe(k1, n1, k2, n2):
    p1, p2 = k1 / n1, k2 / n2
    l1, u1 = wilson(k1, n1)
    l2, u2 = wilson(k2, n2)
    D = p1 - p2
    return (D - np.sqrt((p1 - l1)**2 + (u2 - p2)**2),
            D + np.sqrt((u1 - p1)**2 + (p2 - l2)**2))

cases = [(150, 200, 120, 180), (5, 50, 0, 50), (0, 20, 0, 20),
         (20, 20, 15, 20), (3, 100, 1, 100)]
for k1, n1, k2, n2 in cases:
    p1, p2 = k1 / n1, k2 / n2
    D = p1 - p2
    se = np.sqrt(p1 * (1 - p1) / n1 + p2 * (1 - p2) / n2)
    lo, hi = newcombe(k1, n1, k2, n2)
    print(f"{k1:3d}/{n1:3d} vs {k2:3d}/{n2:3d}: "
          f"왈드 ({D - z * se:7.4f}, {D + z * se:7.4f})  "
          f"뉴콤 ({lo:7.4f}, {hi:7.4f})")

150/200 vs 120/180: 왈드 (-0.0080,  0.1747)  뉴콤 (-0.0079,  0.1737)
  5/ 50 vs   0/ 50: 왈드 ( 0.0168,  0.1832)  뉴콤 ( 0.0090,  0.2136)
  0/ 20 vs   0/ 20: 왈드 ( 0.0000,  0.0000)  뉴콤 (-0.1611,  0.1611)
 20/ 20 vs  15/ 20: 왈드 ( 0.0602,  0.4398)  뉴콤 ( 0.0378,  0.4687)
  3/100 vs   1/100: 왈드 (-0.0187,  0.0587)  뉴콤 (-0.0287,  0.0751)

import numpy as np
from scipy import stats

rng = np.random.default_rng(31)
M = 20_000
z = stats.norm.ppf(0.975)
print(f"{'p1':>5s} {'p2':>5s} {'n':>5s} {'왈드':>8s} {'뉴콤':>8s} "
      f"{'AC':>8s}")
for p1, p2, n in [(0.50, 0.40, 50), (0.10, 0.05, 50), (0.05, 0.02, 100),
                  (0.02, 0.01, 50), (0.90, 0.80, 30)]:
    k1, k2 = rng.binomial(n, p1, M), rng.binomial(n, p2, M)
    ph1, ph2 = k1 / n, k2 / n
    D, true = ph1 - ph2, p1 - p2

    se = np.sqrt(ph1 * (1 - ph1) / n + ph2 * (1 - ph2) / n)
    cw = np.mean((D - z * se <= true) & (true <= D + z * se))

    d = 1 + z**2 / n
    def wil(k):
        ph = k / n
        c = (ph + z**2 / (2 * n)) / d
        h = z / d * np.sqrt(ph * (1 - ph) / n + z**2 / (4 * n**2))
        return c - h, c + h
    l1, u1 = wil(k1)
    l2, u2 = wil(k2)
    lo = D - np.sqrt((ph1 - l1)**2 + (u2 - ph2)**2)
    hi = D + np.sqrt((u1 - ph1)**2 + (ph2 - l2)**2)
    cn = np.mean((lo <= true) & (true <= hi))

    nt = n + 2
    pt1, pt2 = (k1 + 1) / nt, (k2 + 1) / nt
    sea = np.sqrt(pt1 * (1 - pt1) / nt + pt2 * (1 - pt2) / nt)
    ca = np.mean((pt1 - pt2 - z * sea <= true) & (true <= pt1 - pt2 + z * sea))

    print(f"{p1:5.2f} {p2:5.2f} {n:5d} {cw:8.4f} {cn:8.4f} {ca:8.4f}")

   p1    p2     n       왈드       뉴콤       AC
 0.50  0.40    50   0.9469   0.9476   0.9476
 0.10  0.05    50   0.9478   0.9702   0.9692
 0.05  0.02   100   0.9396   0.9694   0.9737
 0.02  0.01    50   0.7785   0.9991   0.9998
 0.90  0.80    30   0.9425   0.9648   0.9629

import numpy as np
from scipy import stats

a, n1 = 30, 200      # 처리군: 사건 30건
b, n2 = 50, 200      # 대조군: 사건 50건
z = stats.norm.ppf(0.975)
p1, p2 = a / n1, b / n2

# 위험차
se_d = np.sqrt(p1 * (1 - p1) / n1 + p2 * (1 - p2) / n2)
print(f"위험차 {p1 - p2:7.4f}  ({p1 - p2 - z * se_d:7.4f}, "
      f"{p1 - p2 + z * se_d:7.4f})")

# 위험비 (로그 척도)
rr = p1 / p2
se_lr = np.sqrt(1 / a - 1 / n1 + 1 / b - 1 / n2)
print(f"위험비 {rr:7.4f}  ({rr * np.exp(-z * se_lr):7.4f}, "
      f"{rr * np.exp(z * se_lr):7.4f})")

# 오즈비 (로그 척도, 울프의 방법)
c, d = n1 - a, n2 - b
orr = (a * d) / (b * c)
se_lo = np.sqrt(1 / a + 1 / b + 1 / c + 1 / d)
print(f"오즈비 {orr:7.4f}  ({orr * np.exp(-z * se_lo):7.4f}, "
      f"{orr * np.exp(z * se_lo):7.4f})")

lo_d, hi_d = p1 - p2 - z * se_d, p1 - p2 + z * se_d
print(f"\n치료필요수 NNT = 1/|위험차| = {1 / abs(p1 - p2):.1f}")
print(f"  95% 구간 {1 / abs(lo_d):.1f} ~ {1 / abs(hi_d):.1f}")

위험차 -0.1000  (-0.1778, -0.0222)
위험비  0.6000  ( 0.3990,  0.9023)
오즈비  0.5294  ( 0.3201,  0.8755)

치료필요수 NNT = 1/|위험차| = 10.0
  95% 구간 5.6 ~ 45.0

import numpy as np
from scipy import stats

# 층별 (a, n1, b, n2) : 처리군 사건/총, 대조군 사건/총
strata = [(20, 100, 30, 100),      # 저위험군
          (40, 120, 55, 110),      # 중위험군
          (35,  60, 42,  55)]      # 고위험군
z = stats.norm.ppf(0.975)

print(f"{'층':>3s} {'p1':>7s} {'p2':>7s} {'위험차':>8s} {'위험비':>7s}")
for i, (a, n1, b, n2) in enumerate(strata, 1):
    print(f"{i:3d} {a / n1:7.4f} {b / n2:7.4f} "
          f"{a / n1 - b / n2:8.4f} {(a / n1) / (b / n2):7.4f}")

# 역분산 가중 (위험차)
ws, ds = [], []
for a, n1, b, n2 in strata:
    p1, p2 = a / n1, b / n2
    v = p1 * (1 - p1) / n1 + p2 * (1 - p2) / n2
    ws.append(1 / v)
    ds.append(p1 - p2)
ws, ds = np.array(ws), np.array(ds)
d_iv = (ws * ds).sum() / ws.sum()
se_iv = np.sqrt(1 / ws.sum())
print(f"\n역분산 가중 위험차 {d_iv:.4f}  "
      f"({d_iv - z * se_iv:.4f}, {d_iv + z * se_iv:.4f})")

# 이질성 검정 (코크런 Q)
Q = (ws * (ds - d_iv)**2).sum()
print(f"코크런 Q = {Q:.3f}  df = {len(strata) - 1}  "
      f"p = {stats.chi2.sf(Q, len(strata) - 1):.4f}")
print(f"I² = {max(0, (Q - (len(strata) - 1)) / Q) * 100:.1f}%")

# 뭉뚱그린 분석 (층 무시)
A = sum(a for a, _, _, _ in strata)
N1 = sum(n1 for _, n1, _, _ in strata)
B = sum(b for _, _, b, _ in strata)
N2 = sum(n2 for _, _, _, n2 in strata)
p1, p2 = A / N1, B / N2
se = np.sqrt(p1 * (1 - p1) / N1 + p2 * (1 - p2) / N2)
print(f"\n층 무시 위험차 {p1 - p2:.4f}  "
      f"({p1 - p2 - z * se:.4f}, {p1 - p2 + z * se:.4f})")

  층      p1      p2      위험차     위험비
  1  0.2000  0.3000  -0.1000  0.6667
  2  0.3333  0.5000  -0.1667  0.6667
  3  0.5833  0.7636  -0.1803  0.7639

역분산 가중 위험차 -0.1418  (-0.2187, -0.0648)
코크런 Q = 0.824  df = 2  p = 0.6622
I² = 0.0%

층 무시 위험차 -0.1400  (-0.2218, -0.0581)

strata2 = [(20, 100, 6, 20), (40, 120, 55, 110), (9, 20, 42, 55)]
ws2, ds2 = [], []
for a, n1, b, n2 in strata2:
    p1, p2 = a / n1, b / n2
    ws2.append(1 / (p1 * (1 - p1) / n1 + p2 * (1 - p2) / n2))
    ds2.append(p1 - p2)
ws2, ds2 = np.array(ws2), np.array(ds2)
d2 = (ws2 * ds2).sum() / ws2.sum()
se2 = np.sqrt(1 / ws2.sum())
A = sum(x[0] for x in strata2); N1 = sum(x[1] for x in strata2)
B = sum(x[2] for x in strata2); N2 = sum(x[3] for x in strata2)
q1, q2 = A / N1, B / N2
sep = np.sqrt(q1 * (1 - q1) / N1 + q2 * (1 - q2) / N2)
print(f"역분산 {d2:.4f}  ({d2 - z * se2:.4f}, {d2 + z * se2:.4f})")
print(f"층 무시 {q1 - q2:.4f}  ({q1 - q2 - z * sep:.4f}, {q1 - q2 + z * sep:.4f})")

역분산 -0.1766  (-0.2760, -0.0773)
층 무시 -0.2693  (-0.3609, -0.1776)

import numpy as np
from scipy import stats

z = stats.norm.ppf(0.975)

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

for k1, n1, k2, n2 in [(9, 10, 1, 10), (10, 12, 0, 8)]:
    p1, p2 = k1 / n1, k2 / n2
    D = p1 - p2
    se = np.sqrt(p1 * (1 - p1) / n1 + p2 * (1 - p2) / n2)
    nt1, nt2 = n1 + 2, n2 + 2
    pt1, pt2 = (k1 + 1) / nt1, (k2 + 1) / nt2
    sea = np.sqrt(pt1 * (1 - pt1) / nt1 + pt2 * (1 - pt2) / nt2)
    l1, u1 = wilson(k1, n1)
    l2, u2 = wilson(k2, n2)
    lo = D - np.sqrt((p1 - l1)**2 + (u2 - p2)**2)
    hi = D + np.sqrt((u1 - p1)**2 + (p2 - l2)**2)
    print(f"{k1}/{n1} vs {k2}/{n2}")
    print(f"  왈드 ({D - z * se:7.4f}, {D + z * se:7.4f})")
    print(f"  AC   ({pt1 - pt2 - z * sea:7.4f}, {pt1 - pt2 + z * sea:7.4f})")
    print(f"  뉴콤 ({lo:7.4f}, {hi:7.4f})")

9/10 vs 1/10
  왈드 ( 0.5370,  1.0630)
  AC   ( 0.3685,  0.9649)
  뉴콤 ( 0.3699,  0.9161)
10/12 vs 0/8
  왈드 ( 0.6225,  1.0442)
  AC   ( 0.4015,  0.9699)
  뉴콤 ( 0.4039,  0.9530)