비율의 대응 구간 (McNemar)¶

점추정값은:

표준오차는:

95% 신뢰구간은:

구간이 0을 포함하므로 95% 수준에서 주변비율이 다르다고 결론지을 충분한 증거가 없다.

불일치 짝의 개수는 \(b + c = 25 + 15 = 40\)으로 최소 기준 10을 훨씬 넘는다. 따라서 McNemar 신뢰구간의 정규근사가 적절하다.

일치 짝(\(a\)와 \(d\))은 두 조건에서 결과가 같으므로 변화의 증거를 주지 않는다. 두 번 모두 양성이거나 두 번 모두 음성인 피험자는 어느 조건이 더 높은 양성 비율을 내는지 가려내는 데 도움이 되지 않는다.

조건 사이에서 결과가 바뀐 불일치 짝만이 차이에 관한 정보를 담는다. \(b > c\)이면 양성에서 음성으로 바뀐 피험자가 그 반대보다 많다는 뜻이므로 \(p_1 > p_2\)를 시사한다. McNemar 검정과 신뢰구간은 오직 비 \(b/(b+c)\)와 그것이 0.5와 다른지에 집중하여 이를 형식화한다.

불일치 짝의 개수가 \(b + c = 3 + 2 = 5\)로, 정규근사에 권장되는 최소 기준 10에 못 미친다. 불일치 짝이 이렇게 적으면 Wald 신뢰구간을 믿을 수 없다 — 개수가 작을 때 이항분포에 대한 정규근사가 나쁘기 때문이다.

권장되는 대안은 확률 \(\pi = P(\text{불일치 짝이 } b \text{ 유형}) = b/(b+c)\)에 대한 정확한 이항 신뢰구간이다. \(H_0\): \(p_1 = p_2\) 아래에서 \(b\) 유형 불일치 짝의 개수는 Binomial\((b+c, 0.5)\)을 따른다. \(\pi\)에 대한 정확한 신뢰구간을 구한 뒤 \(p_1 - p_2\)의 신뢰구간으로 되돌릴 수 있다.

조건부 논증. 불일치 짝의 총수 \(m=b+c\)를 주어진 것으로 보면, \(H_0\)(주변동질성) 아래에서

다. 일치 짝 \(a\), \(d\)는 이 조건부 분포에 전혀 들어오지 않는다. 따라서 맥니마 문제는 "\(m\)번 던진 동전이 공정한가"로 환원된다.

import numpy as np
from scipy import stats

b, c, n = 25, 15, 150
m = b + c

print(f"정확 이항 양측 p = {2 * stats.binom.cdf(min(b, c), m, 0.5):.4f}")
chi = (b - c)**2 / m
print(f"맥니마 χ² = {chi:.3f}  p = {stats.chi2.sf(chi, 1):.4f}")
chi_cc = (abs(b - c) - 1)**2 / m
print(f"연속성 수정 χ² = {chi_cc:.3f}  p = {stats.chi2.sf(chi_cc, 1):.4f}")

# π = b/(b+c) 의 클로퍼-피어슨 구간 → 차이의 구간으로 변환
lo = stats.beta.ppf(0.025, b, m - b + 1)
hi = stats.beta.ppf(0.975, b + 1, m - b)
print(f"π = b/(b+c) 구간 ({lo:.4f}, {hi:.4f})")
print(f"p1 - p2 구간 ({(2 * lo - 1) * m / n:.4f}, {(2 * hi - 1) * m / n:.4f})")

정확 이항 양측 p = 0.1539
맥니마 χ² = 2.500  p = 0.1138
연속성 수정 χ² = 2.025  p = 0.1547
π = b/(b+c) 구간 (0.4580, 0.7727)
p1 - p2 구간 (-0.0224, 0.1455)

변환의 근거. \(\pi=P(B\text{쪽}\mid\text{불일치})\)라 두면

이므로 \(\pi\)의 구간에 이 단조변환을 적용하면 된다. 다만 \(m/n\)을 고정으로 취급하는 근사가 들어 있다.

세 \(p\)-값의 관계.

연속성 수정이 정확값을 잘 흉내 낸다. 이것이 수정을 넣는 이유다. 다만 앞서 비율 구간에서 보았듯, 수정이 항상 좋은 것은 아니다.

언제 정확 방법이 필수인가. \(m=b+c<20\)이면 카이제곱 근사가 믿을 수 없다. 이 예는 \(m=40\)이라 근사가 쓸 만하지만, 정확 방법의 계산 비용이 0이므로 그냥 정확 방법을 쓰는 것이 낫다.

착안. 두 주변비율 각각에 윌슨 구간을 만들고, 둘을 상관 \(\phi\)로 결합한다. 왈드의 경계 문제를 피하면서 대응 구조를 반영한다.

여기서 \((l_j,u_j)\)는 \(\hat p_j\)의 윌슨 구간, \(\hat\phi\)는 \(2\times2\) 표의 파이 계수다.

import numpy as np
from scipy import stats

a, b, c, d = 60, 25, 15, 50
n = a + b + c + d
z = stats.norm.ppf(0.975)

def wilson(k, nn):
    ph = k / nn
    dd = 1 + z**2 / nn
    cc = (ph + z**2 / (2 * nn)) / dd
    h = z / dd * np.sqrt(ph * (1 - ph) / nn + z**2 / (4 * nn**2))
    return cc - h, cc + h

p1, p2 = (a + b) / n, (a + c) / n
l1, u1 = wilson(a + b, n)
l2, u2 = wilson(a + c, n)
phi = (a * d - b * c) / np.sqrt((a + b) * (c + d) * (a + c) * (b + d))
D = p1 - p2

lo_n = D - np.sqrt((p1 - l1)**2 - 2 * phi * (p1 - l1) * (u2 - p2) + (u2 - p2)**2)
hi_n = D + np.sqrt((u1 - p1)**2 - 2 * phi * (u1 - p1) * (p2 - l2) + (p2 - l2)**2)

se = np.sqrt(b + c - (b - c)**2 / n) / n
print(f"p1 = {p1:.4f}  p2 = {p2:.4f}  차이 {D:.4f}   φ = {phi:.4f}")
print(f"왈드   ({D - z * se:.4f}, {D + z * se:.4f})   폭 {2 * z * se:.4f}")
print(f"뉴콤   ({lo_n:.4f}, {hi_n:.4f})   폭 {hi_n - lo_n:.4f}")

p1 = 0.5667  p2 = 0.5000  차이 0.0667   φ = 0.4709
왈드   (-0.0153, 0.1486)   폭 0.1639
뉴콤   (-0.0151, 0.1468)   폭 0.1619

이 자료에서는 거의 같다. \(n=150\)이고 비율이 0.5 근처라 근사가 잘 듣는 영역이기 때문이다.

뉴콤이 빛나는 곳. \(b\)나 \(c\)가 아주 작거나 0일 때다.

• \(b=c=0\)이면 왈드의 표준오차가 0이 되어 구간이 점 \(\{0\}\)으로 퇴화한다. 뉴콤은 유한한 폭을 준다.
• \(b=5\), \(c=0\)이면 왈드는 여전히 작동하지만 하한이 경계를 넘을 수 있다.
• \(p_1\)이나 \(p_2\)가 0이나 1 근처면 왈드 구간이 \([-1,1]\)을 벗어난다. 뉴콤은 윌슨을 바탕으로 하므로 구조적으로 그럴 수 없다.

\(\phi\)의 역할. 상관이 클수록 두 윌슨 구간의 변동이 연동되어 차이의 불확실성이 줄어든다. 여기서 \(\phi=0.47\)이다. \(\phi=0\)이면 독립표본 구간과 같아지고, \(\phi\to1\)이면 구간이 좁아진다.

권고. 모의실험 연구에서 뉴콤의 방법이 포함확률과 폭의 균형에서 가장 낫다고 보고된다. 왈드는 \(b+c\)가 크고 비율이 중간일 때만 안전하다.

정의. 불일치 짝 중 "\(+\to-\)"의 오즈다.

\(\pi=b/(b+c)\)와는 \(\text{OR}=\pi/(1-\pi)\)의 관계에 있다.

import numpy as np
from scipy import stats

b, c = 25, 15
m = b + c
z = stats.norm.ppf(0.975)
OR = b / c
print(f"조건부 오즈비 {OR:.4f}")

# 정확 구간: π의 CP 구간을 오즈로 변환
lo = stats.beta.ppf(0.025, b, m - b + 1)
hi = stats.beta.ppf(0.975, b + 1, m - b)
print(f"정확(CP 기반)  ({lo / (1 - lo):.4f}, {hi / (1 - hi):.4f})")

# 로그 척도 왈드
se = np.sqrt(1 / b + 1 / c)
print(f"로그 왈드      ({OR * np.exp(-z * se):.4f}, {OR * np.exp(z * se):.4f})")

조건부 오즈비 1.6667
정확(CP 기반)  (0.8451, 3.4002)
로그 왈드      (0.8787, 3.1612)

두 구간 모두 1을 담는다. 주변비율 차이의 구간이 0을 담았던 것과 일관된다.

차이 대 오즈비 — 무엇이 다른가.

핵심 차이 — 일치 짝의 역할. 오즈비는 \(a\)와 \(d\)를 완전히 무시한다. 같은 \(b=25\), \(c=15\)라도

• \(n=150\)(\(a+d=110\))이면 차이는 \(10/150=0.067\),
• \(n=1000\)(\(a+d=960\))이면 차이는 \(10/1000=0.010\)

으로 차이는 크게 달라지지만 오즈비는 1.667로 같다.

어느 것을 보고할 것인가.

• 공중보건이나 정책에서는 절대 차이가 중요하다. "전체 인구의 몇 %가 영향을 받는가"가 관심사다.
• 기전이나 개별 위험에서는 오즈비가 유용하다. 기저율이 다른 집단 간 비교에 안정적이다.
• 둘 다 보고하는 것이 가장 좋다. 하나만 보면 오해가 생긴다.

조건부 논리의 다른 예. 이 구조는 조건부 로지스틱회귀로 일반화된다. 짝마다 층을 만들고 개체효과를 조건부로 소거하면, 정확히 이 오즈비가 최대가능도추정값으로 나온다.

두 질문.

• 주변동질성: \(p_1=p_2\)인가? 즉 전체 양성률이 두 시점에서 같은가. 맥니마가 답한다.
• 일치도: 개별 환자의 두 결과가 얼마나 잘 맞는가. 카파가 답한다.

import numpy as np

a, b, c, d = 60, 25, 15, 50
n = a + b + c + d

po = (a + d) / n                                    # 관측 일치율
pe = ((a + b) * (a + c) + (c + d) * (b + d)) / n**2  # 우연 일치율
kappa = (po - pe) / (1 - pe)
print(f"관측 일치율 {po:.4f}   우연 일치율 {pe:.4f}   카파 {kappa:.4f}")
print(f"주변비율 p1 = {(a + b) / n:.4f}   p2 = {(a + c) / n:.4f}   "
      f"차이 {(b - c) / n:.4f}")

관측 일치율 0.7333   우연 일치율 0.5000   카파 0.4667
주변비율 p1 = 0.5667   p2 = 0.5000   차이 0.0667

두 값이 전혀 다른 것을 말한다.

• 주변비율 차이 0.067은 작고 유의하지 않다. 전체 양성률은 거의 그대로다.
• 카파 0.47은 "보통" 수준의 일치다. 개별 환자 40명(27%)의 결과가 바뀌었다.

극단적인 예로 명확히 하자.

셋 다 주변동질성이 완벽히 성립하지만(\(b=c\)), 일치도는 완전 일치에서 완전 불일치까지 걸쳐 있다. 맥니마 검정은 이 셋을 구별하지 못한다.

반대 방향. \(a=90\), \(b=10\), \(c=0\), \(d=0\)이면 일치율이 0.9로 높지만 주변비율이 \(1.0\) 대 \(0.9\)로 다르고 맥니마가 유의하다.

실무적 함의.

1. 검사법 비교에서는 둘 다 봐야 한다. "새 검사가 기존 검사와 같은 비율로 양성을 낸다"와 "같은 사람을 같게 판정한다"는 다른 요구다.

2. 진단 일치도 연구에서 카파만 보고하고 주변 차이를 무시하면, 한 평가자가 체계적으로 후하게 주는 경향을 놓친다.

3. 카파의 역설. 유병률이 극단적이면(\(p\)가 0이나 1 근처) 일치율이 높아도 카파가 낮게 나온다. 우연 일치율 \(p_e\)가 커지기 때문이다. 카파를 해석할 때 주변분포를 함께 제시해야 한다.

문제. \(k\times k\) 표에서 \(H_0:\ p_{i\cdot}=p_{\cdot i}\)(\(i=1,\dots,k\))를 검정한다. 제약이 \(k-1\)개이므로 자유도가 \(k-1\)이다.

스튜어트-맥스웰 검정. 주변도수 차이 벡터 \(\mathbf d=(d_1,\dots,d_{k-1})^\top\), \(d_i=n_{i\cdot}-n_{\cdot i}\)에 대해

import numpy as np
from scipy import stats

# 세 등급(개선/유지/악화)으로 두 번 평가한 120명
N = np.array([[30,  12,  6],
              [ 8,  25, 10],
              [ 4,   9, 16]], dtype=float)
k = N.shape[0]
n = N.sum()
row, col = N.sum(1), N.sum(0)
print("행 주변", row, " 열 주변", col)

d = (row - col)[:k - 1]
V = np.zeros((k - 1, k - 1))
for i in range(k - 1):
    V[i, i] = row[i] + col[i] - 2 * N[i, i]
    for j in range(k - 1):
        if i != j:
            V[i, j] = -(N[i, j] + N[j, i])
chi = d @ np.linalg.solve(V, d)
print(f"스튜어트-맥스웰 χ² = {chi:.4f}  df = {k - 1}  "
      f"p = {stats.chi2.sf(chi, k - 1):.4f}")

# 쌍별 맥니마 (참고용, 다중비교 보정 필요)
for i in range(k):
    for j in range(i + 1, k):
        bb, cc = N[i, j], N[j, i]
        if bb + cc > 0:
            pv = 2 * stats.binom.cdf(min(bb, cc), int(bb + cc), 0.5)
            print(f"  ({i+1},{j+1}) b={bb:.0f} c={cc:.0f}  정확 p = {min(pv, 1):.4f}")

행 주변 [48. 43. 29.]  열 주변 [42. 46. 32.]
스튜어트-맥스웰 χ² = 1.2390  df = 2  p = 0.5382
  (1,2) b=12 c=8  정확 p = 0.5034
  (1,3) b=6 c=4  정확 p = 0.7539
  (2,3) b=10 c=9  정확 p = 1.0000

주변동질성을 기각하지 못한다(\(p=0.538\)). 쌍별 맥니마도 모두 유의하지 않아 결론이 일관된다.

주의할 점.

1. \(V\)가 특이할 수 있다. 어떤 범주의 불일치가 전혀 없으면 역행렬이 없다. 그 범주를 합치거나 일반화역행렬을 쓴다.

2. 쌍별 맥니마를 그냥 여러 번 하면 안 된다. \(k=3\)이면 비교가 3개이므로 다중비교 보정이 필요하고, 게다가 서로 독립이 아니다. 전체 검정을 먼저 하고, 유의하면 사후비교로 넘어간다.

3. 순서형 범주라면 순서를 활용하는 검정이 훨씬 강력하다. 주변동질성의 대립가설을 "한쪽으로 이동"으로 좁히면 자유도가 1로 줄어 검정력이 오른다.

4. 보와커 검정과 구별. 보와커는 대칭성(\(p_{ij}=p_{ji}\) for all \(i,j\))을 검정하며, 자유도가 \(k(k-1)/2\)다. 대칭성은 주변동질성보다 강한 조건이다. \(2\times2\)에서는 둘이 같아진다.

구간은? \(k>2\)에서 "차이"가 벡터이므로 신뢰영역이 된다. 실무에서는 각 범주별 주변 차이에 개별 구간을 만들고 본페로니로 보정하는 것이 보통이다.

핵심 — 검정력을 결정하는 것은 불일치 짝의 수다. 앞서 보았듯 맥니마 검정은 \(B\mid m\sim\text{Bin}(m,\pi)\)에 대한 검정이므로, \(m=b+c\)가 실질적인 표본크기다.

두 단계.

1단계. \(H_0:\pi=1/2\)를 대립 \(\pi_1\)에서 검정하는 데 필요한 불일치 짝 수

2단계. 불일치 비율 \(p_{\text{disc}}=(b+c)/n\)으로 나누어 총 쌍 수를 얻는다.

import numpy as np
from scipy import stats

za, zb = stats.norm.ppf(0.975), stats.norm.ppf(0.80)
print(f"{'불일치율':>8s} {'π1':>6s} {'불일치 짝':>10s} {'총 쌍':>8s}")
for p_disc, pi1 in [(0.20, 0.70), (0.30, 0.65), (0.40, 0.625), (0.10, 0.70)]:
    m = (za + zb)**2 / (4 * (pi1 - 0.5)**2)
    print(f"{p_disc:8.2f} {pi1:6.3f} {np.ceil(m):10.0f} "
          f"{np.ceil(m / p_disc):8.0f}")

    불일치율     π1      불일치 짝      총 쌍
    0.20  0.700         50      246
    0.30  0.650         88      291
    0.40  0.625        126      314
    0.10  0.700         50      491

미리 알아야 할 두 가지.

1. 불일치 비율 \(p_{\text{disc}}\). 두 측정이 얼마나 일치하는지. 예비연구나 문헌에서 얻어야 한다. 이것을 모르면 표본크기를 정할 수 없다.

2. 불일치의 방향 비 \(\pi_1\). 탐지하고자 하는 효과의 크기. "불일치 중 70%가 한 방향"이 대립가설이다.

직관적 함의.

• 일치도가 높을수록 더 많은 쌍이 필요하다. 불일치율 0.10이면 491쌍, 0.40이면 314쌍이다. 역설적이지만 정보는 불일치 짝에서만 나오기 때문이다.
• 다만 일치도가 높으면 대체로 \(\pi_1\)을 크게 잡을 수 있어, 실제로는 상쇄되는 경우가 많다.

주변 차이로 환산. \(p_1-p_2=p_{\text{disc}}(2\pi_1-1)\)이다. 위 첫 줄은 \(0.20\times0.40=0.08\), 즉 8%포인트 차이를 탐지하는 설계다.

독립표본과 비교. 같은 8%포인트 차이를 \(p\approx0.5\)인 두 독립집단에서 탐지하려면 집단당 약 613명, 총 1,226명이 필요하다. 대응 설계가 246쌍(492 관측)으로 해결한다. 절반 이하다.

놓치기 쉬운 것. 계획 단계에서 \(p_{\text{disc}}\)를 과대추정하면 표본이 모자라게 된다. 중간분석에서 관측된 불일치 수를 확인하고 필요하면 표본을 늘리는 설계(내부 예비추정)를 고려한다. 다만 이 경우 \(\alpha\) 보정이 필요하다.