Fisher의 정확검정 (2×2)¶

import numpy as np
from scipy import stats

# 2x2 분할표
#               성공   실패
# 처치군          1      5
# 대조군          8      2
# 기대도수가 5에 못 미치는 칸이 있어 카이제곱 근사를 쓸 수 없는 상황이다.
observed = np.array([[1, 5],
                     [8, 2]])

print("Observed contingency table:")
print(observed)

# Fisher의 정확검정. 근사가 전혀 없고 초기하분포로 확률을 직접 더한다.
# 주변 합계를 고정한 채 가능한 표들을 모두 나열할 수 있기 때문에 가능한 일이다.
odds_ratio, p_value = stats.fisher_exact(observed)

print(f"Odds ratio : {odds_ratio:.4f}")
print(f"p-value    : {p_value:.4f}")

if p_value < 0.05:
    print("Reject H0: significant association (alpha = 0.05).")
else:
    print("Fail to reject H0: no significant association (alpha = 0.05).")

# (1) 일곱 표의 확률을 유리수로
from fractions import Fraction
from math import comb

n, r1, c1 = 16, 6, 9
denom = comb(n, c1)
prob = {a: Fraction(comb(r1, a) * comb(n - r1, c1 - a), denom) for a in range(r1 + 1)}
print(f"\n분모 C(16,9) = {denom}")
print("  a  분자   P(a)")
for a, pr in prob.items():
    print(f"  {a}  {comb(r1, a) * comb(n - r1, c1 - a):5d}   {float(pr):.5f}")
print(f"  분자의 합 {sum(comb(r1, a) * comb(n - r1, c1 - a) for a in prob)}"
      f"   확률의 합 {sum(prob.values())}")

# (2) 양측 p-값. 기준은 확률 자체.
extreme = [a for a in prob if prob[a] <= prob[1]]
p_exact = sum(prob[a] for a in extreme)
print(f"\n극단적인 표 {extreme}   p = {p_exact} = {float(p_exact):.6f}")

# (3) 카이제곱과 견주기
exp = np.outer(observed.sum(1), observed.sum(0)) / n
print(f"\n기대도수\n{exp}   최소 {exp.min():.3f}" + ("  < 5" if exp.min() < 5 else ""))
for name, corr in [("비보정", False), ("예이츠", True)]:
    x, pv = stats.chi2_contingency(observed, correction=corr)[:2]
    print(f"{name}  chi2 = {x:.5f}  p = {pv:.6f}   정확값의 {pv / float(p_exact):.2f} 배"
          f"   {'기각' if pv < 0.05 else '기각 못 함'}")

# 명목 5% 검정의 실제 크기
size = sum(prob[a] for a in prob
           if sum(prob[k] for k in prob if prob[k] <= prob[a]) <= Fraction(1, 20))
print(f"\n명목 5% 검정의 실제 크기 {size} = {float(size):.6f}")

Observed contingency table:
[[1 5]
 [8 2]]
Odds ratio : 0.0500
p-value    : 0.0350
Reject H0: significant association (alpha = 0.05).

분모 C(16,9) = 11440
  a  분자   P(a)
  0     10   0.00087
  1    270   0.02360
  2   1800   0.15734
  3   4200   0.36713
  4   3780   0.33042
  5   1260   0.11014
  6    120   0.01049
  분자의 합 11440   확률의 합 1

극단적인 표 [0, 1, 6]   p = 5/143 = 0.034965

기대도수
[[3.375 2.625]
 [5.625 4.375]]   최소 2.625  < 5
비보정  chi2 = 6.11217  p = 0.013425   정확값의 0.38 배   기각
예이츠  chi2 = 3.80952  p = 0.050962   정확값의 1.46 배   기각 못 함

명목 5% 검정의 실제 크기 5/143 = 0.034965

odds_ratio, p_value = stats.fisher_exact([[3, 1], [1, 3]])
print(f"OR = {odds_ratio:.4f}, p = {p_value:.4f}")

OR = 9.0000, p = 0.4857

import numpy as np
from scipy import stats

rng = np.random.default_rng(1234)
M = 20_000

print(f"{'n1,n2':>8s} {'p':>5s} {'Fisher':>8s} {'중간 p':>8s} "
      f"{'χ²(보정없음)':>13s} {'χ²(Yates)':>11s}")
for (n1, n2), p in [((10, 10), 0.30), ((15, 15), 0.30), ((20, 20), 0.10),
                    ((30, 30), 0.50), ((50, 50), 0.05)]:
    k1, k2 = rng.binomial(n1, p, M), rng.binomial(n2, p, M)
    f = mid = cn = cy = 0
    for i in range(M):
        t = np.array([[k1[i], n1 - k1[i]], [k2[i], n2 - k2[i]]])
        f += stats.fisher_exact(t).pvalue < 0.05

        # 중간 p 값: 관측된 표의 확률을 절반만 더한다
        N, K, n = n1 + n2, n1, k1[i] + k2[i]
        xs = np.arange(0, min(K, n) + 1)
        pmf = stats.hypergeom.pmf(xs, N, K, n)
        obs = stats.hypergeom.pmf(k1[i], N, K, n)
        mid += pmf[pmf <= obs * (1 + 1e-9)].sum() - obs / 2 < 0.05

        if t.sum(0).min() > 0 and t.sum(1).min() > 0:
            cn += stats.chi2_contingency(t, correction=False)[1] < 0.05
            cy += stats.chi2_contingency(t, correction=True)[1] < 0.05
    print(f"{f'{n1},{n2}':>8s} {p:5.2f} {f / M:8.4f} {mid / M:8.4f} "
          f"{cn / M:13.4f} {cy / M:11.4f}")

   n1,n2     p   Fisher     중간 p      χ²(보정없음)   χ²(Yates)
   10,10  0.30   0.0117   0.0199        0.0371      0.0117
   15,15  0.30   0.0160   0.0255        0.0530      0.0142
   20,20  0.10   0.0131   0.0131        0.0389      0.0050
   30,30  0.50   0.0262   0.0273        0.0508      0.0262
   50,50  0.05   0.0077   0.0211        0.0440      0.0077

import numpy as np
from scipy import stats

rng = np.random.default_rng(999)
M, n1, n2 = 3_000, 15, 15
for label, p1, p2 in [("H0 (0.3 대 0.3)", 0.3, 0.3),
                      ("H1 (0.7 대 0.3)", 0.7, 0.3),
                      ("H1 (0.6 대 0.2)", 0.6, 0.2)]:
    k1, k2 = rng.binomial(n1, p1, M), rng.binomial(n2, p2, M)
    f = b = c = 0
    for i in range(M):
        t = np.array([[k1[i], n1 - k1[i]], [k2[i], n2 - k2[i]]])
        f += stats.fisher_exact(t).pvalue < 0.05
        b += stats.barnard_exact(t).pvalue < 0.05
        if t.sum(0).min() > 0 and t.sum(1).min() > 0:
            c += stats.chi2_contingency(t, correction=False)[1] < 0.05
    print(f"{label:>16s}:  Fisher {f / M:.4f}   Barnard {b / M:.4f}   "
          f"χ²(보정없음) {c / M:.4f}")

  H0 (0.3 대 0.3):  Fisher 0.0157   Barnard 0.0400   χ²(보정없음) 0.0517
  H1 (0.7 대 0.3):  Fisher 0.4183   Barnard 0.5853   χ²(보정없음) 0.5920
  H1 (0.6 대 0.2):  Fisher 0.4503   Barnard 0.6140   χ²(보정없음) 0.6407

import numpy as np
from scipy import stats

def two_sided_variants(table):
    a = table[0, 0]
    R1, R2 = table[0].sum(), table[1].sum()
    C1, N = table[:, 0].sum(), table.sum()
    xs = np.arange(max(0, C1 - R2), min(R1, C1) + 1)
    pmf = stats.hypergeom.pmf(xs, N, R1, C1)
    obs = stats.hypergeom.pmf(a, N, R1, C1)

    p_prob = pmf[pmf <= obs * (1 + 1e-9)].sum()              # 확률법
    p_double = min(1.0, 2 * min(stats.hypergeom.cdf(a, N, R1, C1),
                                stats.hypergeom.sf(a - 1, N, R1, C1)))
    mean = R1 * C1 / N
    p_central = pmf[np.abs(xs - mean) >= abs(a - mean) - 1e-9].sum()
    return p_prob, p_double, p_central, p_prob - obs / 2

for t in [np.array([[1, 5], [8, 2]]), np.array([[3, 1], [1, 3]]),
          np.array([[7, 2], [3, 8]]), np.array([[2, 10], [8, 4]])]:
    a, b, c, m = two_sided_variants(t)
    print(f"{str(t.tolist()):>22s}  확률법 {a:.4f} "
          f"(scipy {stats.fisher_exact(t).pvalue:.4f})  "
          f"두 배법 {b:.4f}  중심법 {c:.4f}  중간p {m:.4f}")

      [[1, 5], [8, 2]]  확률법 0.0350 (scipy 0.0350)  두 배법 0.0490  중심법 0.0350  중간p 0.0232
      [[3, 1], [1, 3]]  확률법 0.4857 (scipy 0.4857)  두 배법 0.4857  중심법 0.4857  중간p 0.3714
      [[7, 2], [3, 8]]  확률법 0.0698 (scipy 0.0698)  두 배법 0.0698  중심법 0.0698  중간p 0.0537
     [[2, 10], [8, 4]]  확률법 0.0361 (scipy 0.0361)  두 배법 0.0361  중심법 0.0361  중간p 0.0277

import numpy as np
from scipy import stats
from scipy.optimize import brentq
from scipy.special import comb

def nchg(a, N, K, n, psi):
    """비중심 초기하분포의 지지집합과 확률질량."""
    xs = np.arange(max(0, n - (N - K)), min(K, n) + 1)
    w = comb(K, xs) * comb(N - K, n - xs) * psi**xs
    return xs, w / w.sum()

def cond_mle(a, N, K, n):
    lo, hi = max(0, n - (N - K)), min(K, n)
    if a == lo:
        return 0.0
    if a == hi:
        return np.inf
    def gap(log_psi):
        xs, pm = nchg(a, N, K, n, np.exp(log_psi))
        return (xs * pm).sum() - a
    return np.exp(brentq(gap, -20, 20))

def exact_ci(a, N, K, n, alpha=0.05):
    lo_b, hi_b = max(0, n - (N - K)), min(K, n)
    def upper_tail(log_psi):          # P(X ≥ a) = α/2  → 하한
        xs, pm = nchg(a, N, K, n, np.exp(log_psi))
        return pm[xs >= a].sum() - alpha / 2
    def lower_tail(log_psi):          # P(X ≤ a) = α/2  → 상한
        xs, pm = nchg(a, N, K, n, np.exp(log_psi))
        return pm[xs <= a].sum() - alpha / 2
    lo = 0.0 if a == lo_b else np.exp(brentq(upper_tail, -30, 30))
    hi = np.inf if a == hi_b else np.exp(brentq(lower_tail, -30, 30))
    return lo, hi

for t in [np.array([[1, 5], [8, 2]]), np.array([[3, 1], [1, 3]]),
          np.array([[7, 2], [3, 8]])]:
    a, N = t[0, 0], t.sum()
    K, n = t[0].sum(), t[:, 0].sum()
    sample_or = (t[0, 0] * t[1, 1]) / (t[0, 1] * t[1, 0])
    lo, hi = exact_ci(a, N, K, n)
    print(f"{str(t.tolist()):>22s}")
    print(f"   표본 OR {sample_or:8.4f}   조건부 MLE {cond_mle(a, N, K, n):8.4f}"
          f"   정확 95% CI ({lo:.4f}, {hi:.4f})")

      [[1, 5], [8, 2]]
   표본 OR   0.0500   조건부 MLE   0.0647   정확 95% CI (0.0009, 0.9912)
      [[3, 1], [1, 3]]
   표본 OR   9.0000   조건부 MLE   6.4083   정확 95% CI (0.2117, 626.2435)
      [[7, 2], [3, 8]]
   표본 OR   9.3333   조건부 MLE   8.1538   정확 95% CI (0.8821, 127.0090)

2×2 표를 분석한다
    │
    ├─ 자료가 대응(짝지어진)인가
    │      └─ 예 ──→ McNemar 검정 (이 장의 다른 절)
    │
    └─ 독립인 두 군
         │
         ├─ 설계에서 무엇이 고정되었는가
         │      ├─ 두 군의 크기만 ──→ 바너드 또는 카이제곱
         │      ├─ 전체 n 만 ────────→ 카이제곱
         │      └─ 두 주변합 모두 ───→ 피셔
         │
         ├─ 기대도수를 확인한다
         │      ├─ 모두 5 이상 ──→ 카이제곱(보정 없음)
         │      └─ 5 미만 있음 ──→ 정확검정
         │
         └─ 언제나: 오즈비·위험비·위험차 + 신뢰구간을 함께