Wilcoxon 부호순위검정¶

검산: \(W^+ + W^- = 43 + 12 = 55 = 10(11)/2\). \(\checkmark\)

검정통계량: \(T = \min(43, 12) = 12\).

정규근사:

\(\alpha = 0.05\)에서 \(H_0\)을 기각하지 못한다. 중앙값 반응시간이 250 ms와 다르다는 증거가 충분하지 않다.

import numpy as np
from scipy import stats
d = np.array([268, 241, 273, 230, 255, 262, 289, 237, 295, 275]) - 250
print(stats.wilcoxon(d, method='exact').pvalue)                     # 0.13086
print(stats.wilcoxon(d, method='approx', correction=False).pvalue)  # 0.11413
print(stats.wilcoxon(d, method='approx', correction=True).pvalue)   # 0.12628

0.130859375
0.11412832776223866
0.12627892855469777

(a) \(|d_i|\)를 정렬하면 2, 3, 5, 5, 5, 10, 10, 10, 15, 15이다.

• 2 → 순위 1
• 3 → 순위 2
• 5, 5, 5 → 위치 3, 4, 5 → 중간순위 4
• 10, 10, 10 → 위치 6, 7, 8 → 중간순위 7
• 15, 15 → 위치 9, 10 → 중간순위 9.5

(b)

검산: \(W^+ + W^- = 49 + 6 = 55 = n(n+1)/2 = 10 \times 11/2\).

(c) \(H_0\) 아래에서

보정하지 않은 분산은

이다. 동점 집단은 크기 3(값 5), 크기 3(값 10), 크기 2(값 15)이므로 보정항은

이고 보정된 분산은 \(96.25 - 1.125 = 95.125\), 표준편차는 \(9.753\)이다.

단측 \(p\)값은

import numpy as np
from scipy import stats
d = np.array([10, 5, -5, 15, 2, 5, 10, -3, 15, 10])
print(stats.wilcoxon(d, alternative='greater', method='exact').pvalue)   # 0.013672
print(stats.wilcoxon(d, alternative='greater', method='approx',
                     correction=False).pvalue)                           # 0.013748

출력:

0.013671875
0.01374798396900858

정확 \(p\)값 \(0.01367\)과 거의 일치한다. 동점 보정을 빠뜨리면 \(z = 2.191\), \(p = 0.0142\)가 되어 정확값에서 조금 더 멀어진다.

(d) Wilcoxon 부호순위검정은 \(p \approx 0.0137\)로 \(\alpha = 0.05\)에서 유의한 반면, 부호검정은 \(p \approx 0.055\)로 유의하지 않았다. Wilcoxon이 유의성에 도달한 것은 차이의 크기를 쓰기 때문이다. 큰 양의 차이(10과 15)가 높은 순위를 받아 \(W^+\)에 크게 기여한다. 부호검정은 모든 양의 차이를 똑같이 취급하여 이 정보를 낭비한다.

중앙값이 0인 비대칭 분포를 만들자. \(\text{Exp}(1)\)에서 중앙값 \(\ln 2\)를 빼면 중앙값이 0이면서 오른쪽으로 치우친 분포가 된다.

import numpy as np
from scipy import stats
rng = np.random.default_rng(0)
B = 3000
for n in (10, 20, 50, 100):
    x = rng.exponential(1, (B, n)) - np.log(2)
    pw = np.array([stats.wilcoxon(x[i]).pvalue for i in range(B)])
    ps = np.array([stats.binomtest(int((x[i] > 0).sum()), n).pvalue
                   for i in range(B)])
    print(n, (pw < .05).mean(), (ps < .05).mean())

출력:

10 0.058333333333333334 0.019
20 0.109 0.03966666666666667
50 0.219 0.037
100 0.37566666666666665 0.035

Wilcoxon 검정이 무너진다. 중앙값이 참으로 0인데도 \(n = 100\)에서 표본의 38%를 기각한다. 결정적으로, \(n\)이 커질수록 더 나빠진다. 크기가 명목수준으로 수렴하지 않고 1로 발산한다. 이것은 근사 오차가 아니라 검정이 다른 가설을 검정하고 있다는 신호이다.

이유는 이렇다. Wilcoxon 부호순위검정이 실제로 검정하는 귀무가설은 "중앙값이 0"이 아니라

즉 Hodges-Lehmann 중앙값(모든 쌍 평균 \((X_i + X_j)/2\)의 중앙값, Walsh 평균이라고도 한다)이 0이라는 것이다. 대칭분포에서는 이것이 중앙값과 같지만 비대칭분포에서는 다르다. 중앙값 이동 지수분포에서 Walsh 평균의 중앙값은 약 \(+0.147\)로 0이 아니다.

y = rng.exponential(1, 200000) - np.log(2)
w = (y[:100000] + y[100000:]) / 2
print(np.median(y), np.median(w))     # 0.0001, 0.147

출력:

0.00011479673764763776 0.147140210508641

이 \(+0.147\)이 바로 검정이 잡아내는 것이며, 표본이 커질수록 그 신호가 확실해지므로 기각률이 1로 간다.

부호검정은 이 문제에서 자유롭다. 중앙값을 직접 검정하므로 크기가 \(0.05\) 근처를 유지한다.

실무 지침: 대응차이의 히스토그램이나 Q-Q 그림이 뚜렷한 치우침을 보이면 Wilcoxon 부호순위검정을 쓰지 말라. 검정력이 낮더라도 부호검정이 옳은 답을 준다.

\(H_0\)(대칭성) 아래에서 각 순위 \(i\)의 부호가 독립적으로 \(\pm\) 중 하나를 확률 \(1/2\)씩 갖는다. 지시확률변수

를 두면 \(Z_1, \ldots, Z_n\)은 독립이고 \(Z_i \sim \text{Bernoulli}(1/2)\)이며

이다.

평균. \(E[Z_i] = 1/2\)이므로

분산. \(\text{Var}(Z_i) = 1/4\)이고 \(Z_i\)들이 독립이므로

독립성이 성립하는 이유가 핵심이다. 대칭성 아래에서 \(|D_i|\)와 \(\operatorname{sign}(D_i)\)가 독립이므로, 어떤 절대편차가 어떤 순위를 받았는지와 무관하게 부호가 동전던지기로 결정된다. 이것이 대칭성 가정이 필요한 정확한 지점이며, 연습문제 2에서 이 가정이 깨질 때 검정이 무너지는 이유이기도 하다.

동점이 있으면 순위 \(i\)가 중간순위 \(r_i\)로 대체되어 \(\text{Var}(W^+) = \frac{1}{4}\sum_i r_i^2\)이 된다. 이 합을 정리하면 본문의 보정식이 나온다.

부호순위검정의 신뢰구간은 Walsh 평균에 기반한다. 모든 쌍 \((i \le j)\)에 대해 \((X_i + X_j)/2\)를 계산하면 \(n(n+1)/2 = 55\)개의 값이 나온다. 정렬한 뒤 양 끝에서 \(k\)개씩 잘라내는데, \(k\)는 정확 부호순위 임계값이다.

import numpy as np
from scipy import stats
x = np.array([268, 241, 273, 230, 255, 262, 289, 237, 295, 275])
n = len(x)
walsh = np.array([(x[i] + x[j]) / 2 for i in range(n) for j in range(i, n)])
walsh.sort()
# 양측 95%: W+ <= 8 이면 기각 (2 * P(W+ <= 8) = 0.0488)
k = 8
print(len(walsh), walsh[k], walsh[-(k + 1)])   # 55 248.0 278.5
print(np.median(walsh), x.mean(), np.median(x))  # 263.0 262.5 265.0
print(stats.ttest_1samp(x, 0).confidence_interval())  # (246.90, 278.10)

출력:

55 246.0 278.5
263.0 262.5 265.0
ConfidenceInterval(low=246.89913171804932, high=278.1008682819507)

\(k\)는 \(2P(W^+ \le k) \le 0.05\)를 만족하는 가장 큰 정수이다. \(n = 10\)에서 \(2P(W^+ \le 8) = 0.0488\)이고 \(2P(W^+ \le 9) = 0.0645\)이므로 \(k = 8\)이다.

두 구간이 거의 같다. Wilcoxon 구간의 폭은 \(30.5\), \(t\) 구간의 폭은 \(31.2\)로 오히려 Wilcoxon이 살짝 좁다. 자료가 대칭에 가까워 순위 방법이 잃는 것이 거의 없기 때문이며, \(\text{ARE} = 0.955\)가 실제로 뜻하는 바이다.

두 구간 모두 \(250\)을 포함하므로 앞의 검정 결과(\(p = 0.131\), 기각 못 함)와 일관된다. 검정과 구간이 이렇게 맞물리는 것은 우연이 아니다. 신뢰구간이 검정을 역전시켜 만들어졌으므로 구조적으로 일치한다.

Hodges-Lehmann 추정값 \(263.0\)이 표본평균 \(262.5\)와 표본중앙값 \(265.0\) 사이에 놓인다는 점도 눈여겨보라. Walsh 평균의 중앙값은 평균의 효율성과 중앙값의 로버스트성 사이의 절충으로, 정규분포에서 평균 대비 효율이 \(0.955\)이면서(중앙값은 \(0.637\)) 붕괴점은 \(0.29\)이다(평균은 \(0\)).