비율에 대한 Z-검정¶

\(\hat p = 0.59\). \(z = (0.59 - 0.5)/\sqrt{0.25/100} = 0.09/0.05 = 1.8\).

단측 p-값: \(P(Z > 1.8) \approx 0.036\). \(0.036 < 0.05\)이므로 \(H_0\)을 기각한다.

동전이 앞면을 50% 넘게 보인다는 증거이며, 치우친 동전과 부합한다.

\(H_0: p = 0.80\) 대 \(H_1: p \ne 0.80\).

\(z = (0.74 - 0.80)/\sqrt{0.80 \cdot 0.20/100} = -0.06/0.04 = -1.5\). \(|z| < 1.96\). 기각하지 못한다.

80%를 반박할 증거가 부족하다.

\(\hat p = 28/500 = 0.056\). \(z = (0.056 - 0.04)/\sqrt{0.04 \cdot 0.96/500} = 0.016/0.00877 \approx 1.83\).

양측 p-값 \(= 2 \cdot P(Z > 1.83) \approx 0.067\). \(\alpha = 0.05\)에서 기각하지 못한다.

조건 확인: \(np_0 = 20 \ge 10\), \(n(1-p_0) = 480 \ge 10\). 타당하다.

\(z = 1.83\), 단측 p-값 \(= P(Z > 1.83) \approx 0.034\). \(\alpha = 0.05\)에서 기각한다.

판정이 달라진다: 양측은 아니라고 하고 단측은 그렇다고 한다. 증가만 의심할 사전 근거가 있었다면 단측이 적절하다.

주의: 자료를 본 뒤에 단측을 고르면 제1종 오류가 부풀려진다. 방향은 미리 지정하라.

조건: \(np_0 \ge 10\)이고 \(n(1 - p_0) \ge 10\). 분포의 모양만 보는 느슨한 기준은 문턱값이 5지만, 검정의 오류율을 명목 수준에 맞추려는 자리이므로 보수적 기준 쪽을 진술한다.

연습문제 1: \(np_0 = 100 \cdot 0.5 = 50 \ge 10\), \(n(1-p_0) = 50 \ge 10\). 타당하다.

조건이 깨지면(\(n\)이 작거나 \(p_0\)이 극단적이면) 정확한 이항검정을 쓴다. scipy의 binomtest를 쓸 수 있다.

공식: \(n = (z_{\alpha/2}\sqrt{p_0(1-p_0)} + z_\beta\sqrt{p_1(1-p_1)})^2/(p_1 - p_0)^2\).

\(= (1.96 \sqrt{0.24} + 0.84 \sqrt{0.2436})^2/(0.02)^2\)

\(= (1.96 \cdot 0.490 + 0.84 \cdot 0.494)^2/0.0004\)

\(= (0.960 + 0.415)^2/0.0004 = 1.89/0.0004 \approx 4730\).

2퍼센트포인트 차이를 검정력 80%로 탐지하려면 관측값이 약 4700개 필요하다. 작은 차이에는 큰 표본이 든다.

수정된 통계량.

import numpy as np
from scipy import stats

def level(n, p0, cc=False, alpha=0.05):
    k = np.arange(n + 1)
    ph = k / n
    num = np.abs(ph - p0) - (1 / (2 * n) if cc else 0)
    num = np.maximum(num, 0)
    z = num / np.sqrt(p0 * (1 - p0) / n)
    rej = z > stats.norm.ppf(1 - alpha / 2)
    return stats.binom.pmf(k, n, p0)[rej].sum()

print(f"{'n':>5s} {'p0':>6s} {'점수':>9s} {'점수+수정':>11s} {'정확':>9s}")
for n, p0 in [(20, 0.5), (30, 0.5), (50, 0.3), (100, 0.1),
              (200, 0.05)]:
    k = np.arange(n + 1)
    ex = np.array([stats.binomtest(int(j), n, p0).pvalue
                   for j in k]) <= 0.05
    exact = stats.binom.pmf(k, n, p0)[ex].sum()
    print(f"{n:5d} {p0:6.2f} {level(n, p0):9.4f} "
          f"{level(n, p0, True):11.4f} {exact:9.4f}")

    n     p0        점수       점수+수정        정확
   20   0.50    0.0414      0.0414    0.0414
   30   0.50    0.0428      0.0428    0.0428
   50   0.30    0.0433      0.0433    0.0433
  100   0.10    0.0636      0.0284    0.0443
  200   0.05    0.0328      0.0328    0.0328

대부분의 경우 수정이 아무 차이를 만들지 않는다. 다섯 줄 중 네 줄에서 수정 전후가 같다. 이산성 때문에 \(\frac1{2n}\)을 빼도 기각역이 바뀌지 않기 때문이다.

차이가 나는 경우에는 과잉보정이다. \(n=100\), \(p_0=0.1\)에서 수정 후 0.0284로, 정확검정(0.0443)보다도 훨씬 보수적이 된다. 정확검정을 흉내 내려다 지나쳐 버린다.

왜 그런가. 연속성 수정은 이항의 계단함수를 정규곡선의 중간점에서 맞추는 장치인데, 그 보정량이 꼬리 확률의 실제 오차와 정확히 맞지 않는다. 방향은 맞지만 크기가 과하다.

판단 — 목적에 달렸다.

현대적 관점. 연속성 수정은 정확검정을 손으로 계산할 수 없던 시절의 근사다. 지금은 정확검정이 즉시 계산되므로, 수정의 유일한 용도는 "왜 정확검정이 보수적인가"를 설명하는 것이다.

두 표본의 경우는 다르다. \(2\times2\) 표에서 예이츠의 연속성 수정은 오랜 논쟁의 대상이며, 지나치게 보수적이라는 비판이 강하다. 표본이 작지 않으면 쓰지 않는 것이 권장된다.

교과서마다 다른 조건.

어느 것이 맞는가 — 실제 수준으로 확인한다.

import numpy as np
from scipy import stats

def level_score(n, p0, alpha=0.05):
    k = np.arange(n + 1)
    z = np.abs(k / n - p0) / np.sqrt(p0 * (1 - p0) / n)
    rej = z > stats.norm.ppf(1 - alpha / 2)
    return stats.binom.pmf(k, n, p0)[rej].sum()

print(f"{'p0':>6s} {'n':>6s} {'np0':>6s} {'실제 수준':>10s} {'평가':>8s}")
for p0 in [0.5, 0.2, 0.1, 0.05]:
    for n in [20, 50, 100, 200, 500]:
        lvl = level_score(n, p0)
        tag = "양호" if 0.025 <= lvl <= 0.07 else "주의"
        print(f"{p0:6.2f} {n:6d} {n * p0:6.1f} {lvl:10.4f} {tag:>8s}")
    print()

    p0      n    np0      실제 수준       평가
  0.50     20   10.0     0.0414       양호
  0.50     50   25.0     0.0649       양호
  0.50    100   50.0     0.0569       양호
  0.50    200  100.0     0.0560       양호
  0.50    500  250.0     0.0544       양호

  0.20     20    4.0     0.0437       양호
  0.20     50   10.0     0.0493       양호
  0.20    100   20.0     0.0595       양호
  0.20    200   40.0     0.0415       양호
  0.20    500  100.0     0.0502       양호

  0.10     20    2.0     0.0432       양호
  0.10     50    5.0     0.0297       양호
  0.10    100   10.0     0.0636       양호
  0.10    200   20.0     0.0439       양호
  0.10    500   50.0     0.0437       양호

  0.05     20    1.0     0.0755       주의
  0.05     50    2.5     0.0378       양호
  0.05    100    5.0     0.0341       양호
  0.05    200   10.0     0.0328       양호
  0.05    500   25.0     0.0501       양호

관찰 셋.

1. \(np_0\) 하나로 예측되지 않는다. \(np_0=2\)인데도(\(p_0=0.1\), \(n=20\)) 0.043으로 양호하고, \(np_0=10\)에서(\(p_0=0.1\), \(n=100\)) 0.064로 더 어긋난다.

2. 톱니 때문이다. 앞서 여러 번 본 이산성의 진동이며, \(n\)과 \(p_0\)의 조합마다 다르게 나온다.

3. 정말 위험한 것은 \(np_0<2\)쯤이다. 표에서 유일하게 "주의"인 행이 \(p_0=0.05\), \(n=20\)(\(np_0=1\))로 0.076이다. 검정의 종류에 따라 요구 수준이 다르다는 것이 여기서 드러난다. 표준오차에 \(\hat p\)를 넣는 왈드는 보수적 기준 \(np_0\ge10\)을 채워야 오류율이 명목값 근처에 오지만, \(H_0\) 아래의 분산 \(p_0(1-p_0)/n\)을 쓰는 점수 검정은 그보다 훨씬 덜 요구한다. 위 표대로라면 느슨한 기준 \(np_0\ge5\)조차 점수 검정에는 과한 셈이다. 점수 검정이 왈드보다 훨씬 강건하기 때문이다.

실무 권고.

• 규칙을 외우기보다 실제 수준을 계산한다. 위 코드는 몇 줄이면 되고 정확하다.
• \(\min(np_0,n(1-p_0))\)이 5 아래면 정확검정을 쓴다. 계산 비용이 없으니 의심스러우면 그냥 정확검정으로 가도 좋다. 위 표에서 점수 검정이 정말 흔들리는 것은 \(np_0\)이 2 아래로 내려갈 때였다.
• 교과서 규칙은 대략적 지침으로만 받아들인다.

왜 조건이 필요한가 — 왜도 때문이다. 이항의 왜도가

이고, 정규근사의 오차가 \(\gamma_1\)에 비례한다. \(p_0=0.5\)면 왜도가 0이라 근사가 좋고, \(p_0\)가 극단이면 나쁘다. \(n\gamma_1^2\)가 작아야 한다는 것이 네 번째 조건의 근거다.

\(p\) 관리도. 각 시점 \(i\)에서 \(n_i\)개를 검사해 불량률 \(\hat p_i\)를 그린다.

import numpy as np
from scipy import stats

pbar, n = 0.04, 200
se = np.sqrt(pbar * (1 - pbar) / n)
lcl, ucl = max(pbar - 3 * se, 0), pbar + 3 * se
print(f"중심선 {pbar:.4f},  관리한계 ({lcl:.4f}, {ucl:.4f})")
print(f"불량 개수로는 ({lcl * n:.1f}, {ucl * n:.1f})개")

# 실제 오경보율(이산성 반영)
k = np.arange(n + 1)
alarm = (k / n < lcl) | (k / n > ucl)
a1 = stats.binom.pmf(k, n, pbar)[alarm].sum()
print(f"\n정상 공정의 오경보율 {a1:.5f}  "
      f"→ 평균 실행 길이 {1 / a1:.0f}점")

# 공정이 나빠졌을 때 탐지 확률
print(f"\n{'참 p':>7s} {'점당 탐지확률':>12s} {'평균 탐지까지':>12s}")
for p1 in [0.04, 0.05, 0.06, 0.08, 0.10]:
    d = stats.binom.pmf(k, n, p1)[alarm].sum()
    print(f"{p1:7.3f} {d:12.4f} {1 / d:12.1f}점")

중심선 0.0400,  관리한계 (0.0000, 0.0816)
불량 개수로는 (0.0, 16.3)개

정상 공정의 오경보율 0.00303  → 평균 실행 길이 330점

    참 p      점당 탐지확률      평균 탐지까지
  0.040       0.0030        330.1점
  0.050       0.0238         42.0점
  0.060       0.0946         10.6점
  0.080       0.4340          2.3점
  0.100       0.7925          1.3점

하한이 음수가 되어 0으로 잘렸다. \(\bar p=0.04\), \(n=200\)이면 \(3\) 표준오차가 0.0416이라 아래쪽 한계가 존재하지 않는다. 불량률이 낮아지는 것은 이 차트로 탐지할 수 없다.

오경보율이 0.00303이다. 정규 \(3\sigma\)의 0.0027보다 크며, 한쪽 꼬리만 작동하는데도 그렇다. 이산성과 왜도 때문이다.

평균 실행 길이(ARL)가 330점. 정규 양측 \(3\sigma\)의 370점과 비슷하다.

불량률이 4%에서 6%로 오르면 평균 11점 만에 탐지한다. 8%면 2.3점이다.

설계 시 고려 넷.

1. 한계선을 정확 이항으로 정한다. 정규근사 대신 \(\text{Bin}(n,\bar p)\)의 분위수를 쓰면 원하는 ARL을 정확히 맞출 수 있다.

2. \(np\)가 작으면 \(p\) 관리도가 부적절하다. 위 예에서 \(n\bar p=8\)인데도 하한이 사라졌다. \(np<5\)면 불량 개수를 세는 \(c\) 관리도(포아송)나 불량 사이의 간격을 보는 방법이 낫다.

3. \(\bar p\)를 추정한다는 사실. 초기 단계에서 \(\bar p\)가 부정확하면 한계선이 틀린다. 위상 I(한계선 설정)과 위상 II(감시)를 구분한다.

4. 작은 이동에는 CUSUM이나 EWMA. 슈하트 관리도는 큰 이동에 강하고 작은 이동에 약하다. 1\(\sigma\) 수준의 이동을 잡으려면 누적합 기반 차트가 훨씬 빠르다.

검정과의 연결. 관리도의 각 점은 \(H_0:p=\bar p\)의 검정이다. 앞서 본 대로 반복 검정이므로 \(3\sigma\)라는 보수적 한계를 쓰며, 그것이 ARL 설계의 다른 표현이다.

보고할 것 일곱.

1. 분자와 분모. "200명 중 18명(9.0%)". 비율만 적으면 안 된다.
2. 비교 기준 \(p_0\)와 그 근거.
3. 검정 방법. 정확/점수/왈드, 단측/양측.
4. 검정통계량과 \(p\)-값.
5. 신뢰구간. 방법을 명시한다(윌슨/클로퍼-피어슨).
6. 절대 차이와 상대 차이.
7. 표본설계. 단순임의추출인가, 군집·가중이 있는가.

좋은 보고의 예.

검사한 제품 200개 중 18개(9.0%)가 불량이었다. 목표 불량률 5%와 비교하는 양측 정확 이항검정에서 \(p=0.032\)였다. 불량률의 95% 신뢰구간은 5.4%~13.8%(윌슨 방법)로 목표값 5%를 간신히 벗어난다. 절대 차이는 4.0%포인트, 상대적으로는 1.8배다. 표본은 하루 생산분에서 단순임의추출했다.

흔한 잘못 여섯.

유효숫자의 기준. 표준오차가 \(\sqrt{0.09\times0.91/200}=0.020\)이므로, 소수 첫째 자리(%)까지가 적당하다. "9.0%"는 좋고 "9.00%"는 과잉이다.

희귀사건의 특별한 경우. \(k=0\)이면

300명 중 이상반응이 0건이었다. 이상반응률의 95% 단측 상한은 1.0%다(클로퍼-피어슨, "3의 법칙" \(3/n=1.0\%\)와 일치).

처럼 상한을 보고한다. "이상반응이 없었다"만 적으면 "발생률이 0"이라는 오해를 준다.

여러 비율을 보고할 때. 표로 정리하고 각 행에 \(k/n\), 백분율, 구간을 함께 싣는다. 비교가 목적이면 차이와 그 구간을 별도 열로 둔다.

한 문장. 분자, 분모, 구간 — 이 셋이 있으면 독자가 무엇이든 다시 계산할 수 있다.