McNemar 검정 (대응된 이진 자료)¶

import numpy as np
from scipy import stats

def mcnemar_test(table):
    """
    Perform McNemar's test on a 2x2 table of paired counts.

    Parameters
    ----------
    table : array-like, shape (2, 2)
        Contingency table where off-diagonal cells (b, c)
        represent discordant pairs:
            [[a, b],
             [c, d]]

    Returns
    -------
    statistic : float   McNemar chi-square statistic (continuity-corrected)
    p_value   : float   Two-sided p-value from chi-square(1)
    """
    table = np.asarray(table)
    # 대각선의 a와 d는 아예 쓰이지 않는다.
    # 두 시점에서 상태가 **바뀐** 쌍만이 변화의 증거이기 때문이다.
    # 그래서 표본이 500쌍이어도 실제 정보량은 b + c에 달려 있다.
    b = table[0, 1]
    c = table[1, 0]
    # 연속성 보정: |b - c|에서 1을 뺀다.
    # 이산인 이항분포를 연속인 카이제곱으로 근사하는 데서 오는 편향을 줄인다.
    chi2 = (abs(b - c) - 1) ** 2 / (b + c)
    p_value = stats.chi2(1).sf(chi2)
    return chi2, p_value

# 치료 전후의 질병 상태. 같은 사람을 두 번 관찰한 대응자료다.
#                  치료후+  치료후-
# 치료전+            101      121
# 치료전-             59       33
table = np.array([[101, 121],
                  [ 59,  33]])

chi2, p = mcnemar_test(table)

print(f"McNemar chi2 = {chi2:.4f}")
print(f"p-value      = {p:.4e}")

if p < 0.05:
    print("Reject H0: significant change after treatment (alpha = 0.05).")
else:
    print("Fail to reject H0: no significant change (alpha = 0.05).")

# (1) 비보정판과 z
(a, b), (c, d) = table
m, n = b + c, table.sum()
raw = (b - c) ** 2 / m
print(f"\n비보정 chi2 = {raw:.5f}   z = {(b - c) / np.sqrt(m):.5f}"
      f"   p = {stats.chi2(1).sf(raw):.4e}")

# (2) 정확 이항검정
exact = stats.binomtest(int(b), int(m), 0.5).pvalue
print(f"정확 이항 p = {exact:.4e}")
print(f"  비보정/정확 = {stats.chi2(1).sf(raw) / exact:.2f} 배")
print(f"  보정  /정확 = {p / exact:.2f} 배")

# (3) 효과크기와 신뢰구간
p_before, p_after = (a + b) / n, (a + c) / n
d_hat = (b - c) / n
var_paired = (m - (b - c) ** 2 / n) / n ** 2
se_paired = np.sqrt(var_paired)
se_indep = np.sqrt(p_before * (1 - p_before) / n + p_after * (1 - p_after) / n)
cov = (a / n - p_before * p_after) / n
phi = (a * d - b * c) / np.sqrt((a + b) * (c + d) * (a + c) * (b + d))
print(f"\n전 {p_before:.5f}  후 {p_after:.5f}  차 {d_hat:.5f}")
print(f"대응 SE {se_paired:.6f}   독립이라 보면 {se_indep:.6f}")
print(f"cov {cov:.8f}   Var_독립 - 2cov = {se_indep**2 - 2 * cov:.8f}"
      f"   Var_대응 = {var_paired:.8f}")
print(f"phi = {phi:.4f}")
print(f"95% 구간 [{d_hat - 1.96 * se_paired:.4f}, {d_hat + 1.96 * se_paired:.4f}]")

McNemar chi2 = 20.6722
p-value      = 5.4501e-06
Reject H0: significant change after treatment (alpha = 0.05).

비보정 chi2 = 21.35556   z = 4.62121   p = 3.8151e-06
정확 이항 p = 4.4344e-06
  비보정/정확 = 0.86 배
  보정  /정확 = 1.23 배

전 0.70701  후 0.50955  차 0.19745
대응 SE 0.041249   독립이라 보면 0.038152
cov -0.00012294   Var_독립 - 2cov = 0.00170147   Var_대응 = 0.00170147
phi = -0.1697
95% 구간 [0.1166, 0.2783]

from scipy import stats

b, c = 121, 59            # 위 보기의 불일치 쌍
n_discordant = b + c
p_exact = stats.binomtest(b, n_discordant, 0.5).pvalue
print(f"exact p = {p_exact:.4e}")

exact p = 4.4344e-06

import numpy as np
from scipy import stats

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

def four_tests(b, c, alpha=0.05):
    """불일치 도수 b, c 에 대해 네 가지 검정의 기각 여부."""
    n_disc = b + c
    if n_disc == 0:
        return (False, False, False, False)
    chi_corr = (abs(b - c) - 1)**2 / n_disc if abs(b - c) >= 1 else 0.0
    chi_raw = (b - c)**2 / n_disc
    k = min(b, c)
    p_exact = min(1.0, 2 * stats.binom.cdf(k, n_disc, 0.5))
    p_mid = p_exact - stats.binom.pmf(k, n_disc, 0.5)
    return (stats.chi2.sf(chi_corr, 1) < alpha,
            stats.chi2.sf(chi_raw, 1) < alpha,
            p_exact < alpha, p_mid < alpha)

for title, split in [("① 수준 (두 방향이 같음)", (0.5, 0.5)),
                     ("② 검정력 (한 방향이 70%)", (0.7, 0.3))]:
    print(title)
    print(f"{'n':>5s} {'불일치율':>9s} {'보정 χ²':>9s} {'무보정 χ²':>10s} "
          f"{'정확':>8s} {'중간 p':>8s}")
    for n, p_disc in [(50, 0.20), (100, 0.20), (50, 0.40), (200, 0.10)]:
        cnt = np.zeros(4)
        for _ in range(M):
            m = rng.multinomial(n, [1 - p_disc,
                                    p_disc * split[0], p_disc * split[1]])
            cnt += np.array(four_tests(m[1], m[2]), float)
        print(f"{n:5d} {p_disc:9.2f} " + " ".join(f"{v / M:9.4f}" for v in cnt))
    print()

① 수준 (두 방향이 같음)
    n      불일치율     보정 χ²     무보정 χ²       정확     중간 p
   50      0.20    0.0217    0.0417    0.0220    0.0405
  100      0.20    0.0275    0.0530    0.0305    0.0469
   50      0.40    0.0267    0.0513    0.0296    0.0453
  200      0.10    0.0284    0.0544    0.0315    0.0487

② 검정력 (한 방향이 70%)
    n      불일치율     보정 χ²     무보정 χ²       정확     중간 p
   50      0.20    0.1461    0.2184    0.1475    0.2152
  100      0.20    0.3407    0.4531    0.3567    0.4295
   50      0.40    0.3405    0.4502    0.3549    0.4258
  200      0.10    0.3398    0.4484    0.3535    0.4253

import numpy as np
from scipy import stats

rng = np.random.default_rng(975)
M, n = 20_000, 200

print(f"{'일치 비율':>10s} {'McNemar':>9s} {'독립 χ²(잘못)':>14s}")
for pa, pb, pc, pd in [(0.40, 0.10, 0.10, 0.40),
                       (0.25, 0.25, 0.25, 0.25),
                       (0.45, 0.05, 0.05, 0.45)]:
    mc = ind = 0
    for _ in range(M):
        a, b, c, d = rng.multinomial(n, [pa, pb, pc, pd])
        if b + c > 0:
            mc += 2 * min(1.0, stats.binom.cdf(min(b, c),
                                               b + c, 0.5)) < 0.05
        # 잘못된 분석: 전/후를 독립인 두 표본으로 본다
        wrong = np.array([[a + b, c + d], [a + c, b + d]], float)
        if wrong.sum(0).min() > 0 and wrong.sum(1).min() > 0:
            ind += stats.chi2_contingency(wrong,
                                          correction=False)[1] < 0.05
    print(f"{pa + pd:10.2f} {mc / M:9.4f} {ind / M:14.4f}")

     일치 비율   McNemar      독립 χ²(잘못)
      0.80    0.0357         0.0022
      0.50    0.0423         0.0534
      0.90    0.0294         0.0001

print(f"{'일치 비율':>10s} {'McNemar':>9s} {'독립 χ²(잘못)':>14s}")
for pa, pb, pc, pd in [(0.40, 0.14, 0.06, 0.40),
                       (0.45, 0.08, 0.02, 0.45)]:
    mc = ind = 0
    for _ in range(M):
        a, b, c, d = rng.multinomial(n, [pa, pb, pc, pd])
        if b + c > 0:
            mc += 2 * min(1.0, stats.binom.cdf(min(b, c),
                                               b + c, 0.5)) < 0.05
        wrong = np.array([[a + b, c + d], [a + c, b + d]], float)
        if wrong.sum(0).min() > 0 and wrong.sum(1).min() > 0:
            ind += stats.chi2_contingency(wrong,
                                          correction=False)[1] < 0.05
    print(f"{pa + pd:10.2f} {mc / M:9.4f} {ind / M:14.4f}")

     일치 비율   McNemar      독립 χ²(잘못)
      0.80    0.6770         0.2837
      0.90    0.7362         0.0444

import numpy as np
from scipy import stats

def n_mcnemar(p_disc, pi, power=0.80, alpha=0.05):
    """필요한 (전체 쌍 수, 불일치 쌍 수)."""
    za, zb = stats.norm.ppf(1 - alpha / 2), stats.norm.ppf(power)
    n_disc = ((za + zb * 2 * np.sqrt(pi * (1 - pi))) / (2 * pi - 1))**2
    return int(np.ceil(n_disc / p_disc)), int(np.ceil(n_disc))

print(f"{'불일치율':>9s} {'한 방향 비율 π':>14s} {'필요 쌍':>9s} "
      f"{'필요 불일치 쌍':>14s}")
for p_disc in [0.10, 0.20, 0.30]:
    for pi in [0.65, 0.70, 0.80]:
        n_total, n_disc = n_mcnemar(p_disc, pi)
        print(f"{p_disc:9.2f} {pi:14.2f} {n_total:9d} {n_disc:14d}")

     불일치율      한 방향 비율 π      필요 쌍       필요 불일치 쌍
     0.10           0.65       849             85
     0.10           0.70       467             47
     0.10           0.80       193             20
     0.20           0.65       425             85
     0.20           0.70       234             47
     0.20           0.80        97             20
     0.30           0.65       283             85
     0.30           0.70       156             47
     0.30           0.80        65             20

b, c = 30, 10
n_disc = b + c
lo, hi = stats.beta.ppf([0.025, 0.975],
                        [b, b + 1], [n_disc - b + 1, n_disc - b])
print(f"b={b}, c={c}:  OR = b/c = {b / c:.4f}")
print(f"  π 의 정확 95% CI ({lo:.4f}, {hi:.4f})")
print(f"  OR 의 95% CI ({lo / (1 - lo):.4f}, {hi / (1 - hi):.4f})")

b=30, c=10:  OR = b/c = 3.0000
  π 의 정확 95% CI (0.5880, 0.8731)
  OR 의 95% CI (1.4274, 6.8793)

import numpy as np
from scipy import stats

def bowker(table):
    """보커의 대칭성 검정 (k×k 대응 표)."""
    table = np.asarray(table, float)
    k = table.shape[0]
    stat = 0.0
    for i in range(k):
        for j in range(i + 1, k):
            s = table[i, j] + table[j, i]
            if s > 0:
                stat += (table[i, j] - table[j, i])**2 / s
    df = k * (k - 1) // 2
    return stat, df, stats.chi2.sf(stat, df)

# 치료 전(행) × 치료 후(열): 경·중·중증
tbl = np.array([[40, 12,  3],
                [ 6, 30,  8],
                [ 2,  5, 24]], float)
stat, df, p = bowker(tbl)
print(f"보커 대칭성:  χ² = {stat:.4f},  df = {df},  p = {p:.4f}")
print(f"  전 주변분포 {np.round(tbl.sum(1) / tbl.sum(), 4).tolist()}")
print(f"  후 주변분포 {np.round(tbl.sum(0) / tbl.sum(), 4).tolist()}")

# 주변분포가 같지만 대칭이 아닌 표
tbl2 = np.array([[40, 15,  5],
                 [ 5, 40, 15],
                 [15,  5, 40]], float)
stat2, df2, p2 = bowker(tbl2)
print(f"\n순환 구조 표")
print(f"  보커 대칭성:  χ² = {stat2:.4f},  df = {df2},  p = {p2:.4f}")
print(f"  전 주변분포 {np.round(tbl2.sum(1) / tbl2.sum(), 4).tolist()}")
print(f"  후 주변분포 {np.round(tbl2.sum(0) / tbl2.sum(), 4).tolist()}")
print("  → 주변분포는 같은데 대칭성은 기각된다")

보커 대칭성:  χ² = 2.8923,  df = 3,  p = 0.4085
  전 주변분포 [0.4231, 0.3385, 0.2385]
  후 주변분포 [0.3692, 0.3615, 0.2692]

순환 구조 표
  보커 대칭성:  χ² = 15.0000,  df = 3,  p = 0.0018
  전 주변분포 [0.3333, 0.3333, 0.3333]
  후 주변분포 [0.3333, 0.3333, 0.3333]
  → 주변분포는 같은데 대칭성은 기각된다

이진 결과를 비교한다
    │
    ├─ 같은 대상을 두 번 측정했는가
    │      ├─ 예 ──→ McNemar (대응)
    │      └─ 아니오 ──→ 보통의 χ² / 피셔 (독립)
    │
    ├─ 조건이 셋 이상인가
    │      └─ 예 ──→ 코크런 Q (다음 절)
    │
    └─ 범주가 셋 이상인가
           └─ 예 ──→ 보커 대칭성 / 주변 동질성 (연습문제 9)