대응 부호검정¶

1단계. 위에서 차이를 계산했다. 동점 하나(\(D_6 = 0\))를 제외한다.

2단계. 남은 것: \(n = 9\), \(n_+ = 7\), \(n_- = 2\).

3단계. \(H_0\) 아래에서 \(n_+ \sim \text{Binomial}(9, 0.5)\)이다.

4단계. 오른쪽 꼬리 \(p\)값(\(H_a\): 프로그램이 심박수를 낮추므로 \(D_i > 0\)이 기대된다):

\(\alpha = 0.05\)에서 \(H_0\)을 기각하지 못한다. 동점이 아닌 9쌍 중 7쌍이 개선을 보였는데도 표본이 너무 작아 부호검정이 유의성에 이르지 못했다. 소표본에서 부호검정의 검정력이 낮음을 보여 주는 예이다.

import numpy as np, itertools
from scipy import stats
d = np.array([4, 5, -2, 2, 7, 0, 8, 2, 3, -2])
nz = d[d != 0]          # 동점 제외, n = 9

print(stats.binomtest(int((nz > 0).sum()), len(nz),
                      alternative='greater').pvalue)          # 0.0898
print(stats.wilcoxon(nz, alternative='greater',
                     method='exact').pvalue)                  # 0.0195
print(stats.ttest_1samp(nz, 0, alternative='greater').pvalue) # 0.0165

a = np.abs(nz); T = nz.sum()
allT = np.array([np.dot(s, a) for s in itertools.product([-1, 1], repeat=9)])
print((allT >= T).mean())                                     # 0.0254

출력:

0.08984375
0.01953125
0.01652701624852805
0.025390625

부호검정만 결론이 다르다. 나머지 셋은 모두 유의하다고 판정한다.

이유는 이 자료의 구조에 있다. 음의 차이 두 개는 모두 \(-2\)로 절댓값이 가장 작다. 양의 차이는 \(4, 5, 2, 7, 8, 2, 3\)으로 \(8\)까지 올라간다. 크기를 쓰는 세 검정은 "반대 방향 증거가 있긴 하나 모두 미미하다"를 읽어 내지만, 부호검정에게는 \(-2\)가 \(-200\)과 같은 한 개의 음수일 뿐이다.

순열검정(\(0.0254\))이 Wilcoxon(\(0.0195\))보다 \(p\)값이 살짝 크다는 점도 흥미롭다. 순열검정은 원래 값을 쓰므로 이 자료에서 순위보다 정보를 더 쓰는데도 그렇다. \(n = 9\)에서 순열분포가 \(2^9 = 512\)개의 값밖에 갖지 못해 이산성이 크기 때문이며, 이 정도 차이는 잡음 범위이다.

from scipy import stats
# n+ = 7 은 고정, 동점만 늘려 전체 쌍 수를 키운다
for extra in range(0, 5):
    n = 9                # 0이 아닌 쌍의 수는 변하지 않는다
    print(10 + extra, n, stats.binomtest(7, n, alternative='greater').pvalue)

출력:

10 9 0.08984375
11 9 0.08984375
12 9 0.08984375
13 9 0.08984375
14 9 0.08984375

출력은 모두 \(0.0898\)이다. 동점을 아무리 많이 붙여도 부호검정의 \(p\)값은 전혀 변하지 않는다.

이것이 부호검정의 동점 처리가 가진 근본적 성질이자 위험이다. 참가자 100명 중 91명이 변화가 없고 7명이 개선, 2명이 악화했다고 하자. 부호검정은 여전히 \(p = 0.0898\)을 준다. 하지만 "100명 중 7명만 개선"이라는 사실은 실질적으로 프로그램의 효과가 미미함을 강하게 시사한다.

부호검정이 검정하는 귀무가설이 "\(D \ne 0\)인 쌍들 중에서" \(P(D > 0) = P(D < 0)\)이기 때문이다. 변화가 없는 쌍은 이 조건부 확률에 대해 아무 정보도 주지 않으므로 옳게 버려지는 것이다.

실무 권고: 동점 비율이 높으면 부호검정 \(p\)값과 함께 동점 비율을 반드시 함께 보고해야 한다. 검정 결과만 보고하면 효과의 실질적 크기를 완전히 오도한다. 효과 크기가 관심사라면 \(\hat{p} = n_+/(n_+ + n_-)\)가 아니라 \(n_+/N\)을 보는 것이 낫다.

부호검정을 역전시키면 중앙값의 신뢰구간이 나온다. 차이를 정렬한 순서통계량 \(D_{(1)} \le \cdots \le D_{(n)}\)에 대해 구간은

이며, \(k\)는 \(2P(S \le k) \le \alpha\)를 만족하는 가장 큰 정수이다(\(S \sim \text{Binomial}(n, 0.5)\)).

import numpy as np
from scipy import stats
d = np.array([4, 5, -2, 2, 7, 0, 8, 2, 3, -2])
nz = np.sort(d[d != 0]); n = len(nz)     # n = 9
for k in range(0, 4):
    print(k, 2 * stats.binom.cdf(k, n, 0.5))
# k=0: 0.0039, k=1: 0.0391, k=2: 0.1797  →  k = 1
print(nz)                     # [-2 -2  2  2  3  4  5  7  8]
print(nz[1], nz[n - 2])       # -2, 7

출력:

0 0.00390625
1 0.0390625
2 0.1796875
3 0.5078125
[-2 -2  2  2  3  4  5  7  8]
-2 7

\(k = 1\)이므로 95% 신뢰구간은 \(\bigl(D_{(2)}, D_{(8)}\bigr) = (-2, \; 7)\)이다.

구간이 \(0\)을 포함하므로 검정 결과(\(p = 0.0898\), 기각 못 함)와 일관된다. 구간과 검정이 이렇게 맞물리는 것은 구간이 검정을 역전시켜 만들어졌기 때문이다.

구간이 매우 넓다는 점에 주목하라. 폭이 \(9\) bpm으로, 자료 전체 범위(\(-2\)에서 \(8\))와 거의 같다. \(n = 9\)에서 부호검정 기반 구간이 쓸 수 있는 것은 \(k = 1\)뿐이므로 양 끝에서 하나씩만 잘라 낼 수 있다. 실제 포함확률은 \(1 - 0.0391 = 0.961\)로 명목값보다 높다. 소표본에서 부호검정 구간이 지나치게 보수적인 이유이다.

import numpy as np
from scipy import stats
before = np.array([72, 80, 68, 76, 85, 74, 90, 78, 82, 70])
after  = np.array([68, 75, 70, 74, 78, 74, 82, 76, 79, 72])
print(stats.mannwhitneyu(before, after, alternative='greater'))
print(np.corrcoef(before, after)[0, 1])

출력:

MannwhitneyuResult(statistic=61.0, pvalue=0.2128050157999556)
0.9232351791538461

대응 구조를 무시하면 \(p\)값이 \(0.213\)으로 크게 나빠진다. 대응 Wilcoxon의 \(0.0195\)와 비교하면 11배이다. 이유는 두 측정값의 상관이 \(0.923\)으로 매우 높기 때문이다.

직관은 이렇다. 참가자들의 기저 심박수가 \(68\)에서 \(90\)까지 넓게 퍼져 있다. 두 집단을 독립으로 보면 이 개인차가 모두 잡음으로 들어간다. 대응 분석은 각 참가자 안에서 차이를 계산하여 개인차를 통째로 소거하고, 남은 것은 처치 효과뿐이다.

이것이 대응설계의 존재 이유이며, 대응자료를 독립표본으로 분석하는 것이 흔하면서도 값비싼 실수인 이유이다. 반대 방향의 실수(독립자료를 억지로 짝지어 대응 분석) 역시 위험하다. 이 경우 자유도를 절반으로 잃으면서 아무것도 얻지 못한다.