대응 순열검정¶

\(\alpha = 0.05\)에서 \(H_0\)을 기각하지 못한다. 5쌍만으로는 중간 정도의 효과를 탐지할 검정력이 부족하다.

import numpy as np, itertools
from scipy import stats
D = np.array([6, -3, 7, 3, 4])

a = np.abs(D); T = D.sum()
allT = np.array([np.dot(s, a) for s in itertools.product([-1, 1], repeat=5)])
print((np.abs(allT) >= abs(T)).mean())        # 0.1875
print(stats.wilcoxon(D, method='exact').pvalue)  # 0.1875
print(stats.ttest_1samp(D, 0).pvalue)            # 0.1239
print(stats.binomtest(4, 5).pvalue)              # 0.3750

출력:

0.1875
0.1875
0.12385684836470104
0.375

순열검정과 Wilcoxon이 소수점 이하 전부 일치한다. 우연이 아니라 \(n = 5\)에서 순열분포가 아주 성기기 때문이다.

\(p\)값은 \(\{|T(\mathbf{s})| \ge |T_{\text{obs}}|\}\)인 부호 배정의 개수로 결정된다. Wilcoxon은 같은 논리를 \(|D|\) 대신 순위 \((1, \ldots, 5)\)에 적용한다. 이 자료의 절대차이 \((6, 3, 7, 3, 4)\)가 순위로는 \((4, 1.5, 5, 1.5, 3)\)인데, 두 벡터의 순서 구조가 같으므로 극단적인 배정 6개가 정확히 같은 6개가 된다.

두 검정이 언제나 일치하는 것은 아니다. 절대차이 사이의 비율이 순위 사이의 비율과 크게 다르면 갈라진다. 예를 들어 \(|D| = (1, 2, 3, 4, 100)\)이면 순열검정은 \(100\)의 압도적 크기를 반영하지만 Wilcoxon에게는 그냥 순위 5일 뿐이다.

두 경우를 보자. \(D_A = (1, 2, 3, 4, 100)\)(모두 양수)과 \(D_B = (1, 2, 3, 4, -100)\)(가장 큰 것만 음수).

import numpy as np, itertools
from scipy import stats

def perm_p(D):
    a = np.abs(D); T = D.sum()
    allT = np.array([np.dot(s, a) for s in itertools.product([-1, 1], repeat=len(D))])
    return (np.abs(allT) >= abs(T)).mean()

for D in (np.array([1, 2, 3, 4, 100]), np.array([1, 2, 3, 4, -100])):
    print(D, perm_p(D),
          stats.wilcoxon(D, method='exact').pvalue,
          stats.ttest_1samp(D, 0).pvalue)

출력:

[  1   2   3   4 100] 0.0625 0.0625 0.32245920880686274
[   1    2    3    4 -100] 1.0 0.625 0.4296209319991539

첫 번째 자료에서는 세 검정이 갈린다. 순열검정과 Wilcoxon은 \(n = 5\)의 최소값 \(2/32 = 0.0625\)를 준다. 모든 차이가 양수이므로 가능한 한 극단적이다. 반면 \(t\) 검정은 \(0.32\)로 아무것도 못 찾는다. 분자 \(\bar{D} = 22\)가 크지만 \(100\)이 표본표준편차를 \(43.6\)까지 부풀려 \(t = 1.13\)밖에 안 되기 때문이다.

여기서 순위 기반 검정과 순열검정의 로버스트성이 드러난다. 두 검정 모두 "다섯 개가 모두 같은 방향"이라는 사실만 보고 판정하는 반면, \(t\) 검정은 이상치 하나에 무너진다.

두 번째 자료에서 순열검정이 정확히 \(p = 1\)을 준다. \(T = -90\)인데, \(|D| = (1,2,3,4,100)\)의 어떤 부호 배정에서도 \(100\)이 나머지 합 \(10\)을 압도하므로 모든 \(|T(\mathbf{s})|\)가 \(90\) 이상이다(\(90\)에서 \(110\) 사이). 관측값이 가장 덜 극단적인 배정 중 하나인 셈이다.

Wilcoxon은 \(W^+ = 1+2+3+4 = 10\), \(W^- = 5\)로 순위상 어느 정도 불균형이 있어 \(p = 0.625\)이다. 순위를 쓰면 \(100\)이 그저 "가장 큰 것"일 뿐이므로 다른 넷이 반대 방향이라는 사실이 약한 신호로 남는다.

핵심 차이: 순열검정은 크기의 비율에 민감하고 Wilcoxon은 그렇지 않다. \(100\)을 \(1000\)으로 바꾸어도 Wilcoxon \(p\)값은 \(0.625\) 그대로지만, 순열검정에서는 \(100\)이 지배적일수록 \(p\)가 1에 붙는다. 이것이 항상 좋은 성질은 아니다. 지배적인 이상치 하나가 순열검정을 완전히 무력화할 수 있다.

\(H_0\)이 참인 자료(\(\mathcal{N}(0,1)\), \(n = 15\))를 생성하여 두 방식의 실제 크기를 비교한다.

import numpy as np
rng = np.random.default_rng(0)
n, M = 15, 40000

for B in (20, 100, 200, 1000):
    raw = adj = 0
    for _ in range(M):
        D = rng.normal(0, 1, n)
        a = np.abs(D); T = abs(D.sum())
        s = rng.choice([-1, 1], size=(B, n))
        c = np.sum(np.abs(s @ a) >= T)
        raw += (c / B) <= 0.05
        adj += ((c + 1) / (B + 1)) <= 0.05
    print(B, raw / M, adj / M)

출력:

20 0.0942 0.04645
100 0.0588 0.049675
200 0.05315 0.04795
1000 0.05105 0.0499

\(B = 20\)에서 보정하지 않은 검정의 실제 크기가 \(0.094\)로 명목값의 거의 두 배이다. \(B\)가 커지면 차이가 줄어들지만 \(B = 1000\)에서도 여전히 \(0.0510\) 대 \(0.0499\)로 보정한 쪽이 정확하다.

이유는 이렇다. \(c/B \le 0.05\)가 되려면 \(B = 20\)에서 \(c \le 1\)이면 된다. 관측값보다 극단적인 순열이 하나까지 허용되므로 기준이 느슨하다. 반면 \((c+1)/21 \le 0.05\)는 \(c = 0\)을 요구한다.

개념적으로도 보정이 옳다. 귀무가설 아래에서 관측된 자료는 다른 모든 순열과 교환 가능하다. 따라서 \(B\)개의 무작위 순열과 관측값을 합쳐 \(B+1\)개를 만들고 관측값이 그중 상위 \(\alpha\)에 드는지를 물어야 한다. 분자와 분모의 \(+1\)이 바로 "관측값 자신"을 세는 것이다.

이 때문에 실무에서 \(B\)를 \(999\)나 \(9999\)로 잡는 관행이 생겼다. \(B + 1\)이 깔끔한 수가 되어 \(\hat{p}\)의 분모가 \(1000\)이나 \(10000\)이 된다.

import numpy as np
from scipy import stats
rng = np.random.default_rng(9)

def perm_p_mc(D, B=4999, rng=rng):
    a = np.abs(D); T = abs(D.sum())
    signs = rng.choice([-1, 1], size=(B, len(D)))
    Ts = np.abs(signs @ a)
    return (np.sum(Ts >= T) + 1) / (B + 1)

for n in (5, 10, 20, 50):
    pp, pt = [], []
    for _ in range(300):
        D = rng.normal(0.5, 1, n)
        pp.append(perm_p_mc(D)); pt.append(stats.ttest_1samp(D, 0).pvalue)
    pp, pt = np.array(pp), np.array(pt)
    print(n, round(np.corrcoef(pp, pt)[0, 1], 4),
          round(np.max(np.abs(pp - pt)), 4))

출력:

5 0.9752 0.4711
10 0.999 0.0784
20 0.9996 0.0391
50 0.9997 0.0108

상관이 \(n = 10\)에서 이미 \(0.999\)이다. 최대 차이는 \(n = 5\)에서 \(0.47\)로 매우 크지만(\(2^5 = 32\)가지 순열밖에 없어 \(p\)값이 \(1/32\) 단위로 계단을 이룬다) \(n = 20\)부터는 \(0.04\) 아래로 떨어진다.

이유는 두 검정이 결국 같은 통계량 \(\sum D_i\)를 쓰기 때문이다. 차이는 귀무분포를 어떻게 얻느냐뿐이다. \(t\) 검정은 정규이론으로, 순열검정은 열거로 얻는다. 자료가 정규이고 \(n\)이 커지면 순열분포가 \(t\) 분포로 수렴한다.

실무적 함의는 양방향이다.

• 자료가 정규에 가까우면 순열검정을 써도 잃는 것이 거의 없다. 즉 순열검정을 기본으로 삼아도 비용이 없다.
• 반대로 자료가 정규이면 굳이 순열검정의 계산 비용을 치를 필요도 없다.

순열검정의 진짜 가치는 자료가 정규가 아닐 때 드러난다. 연습문제 2의 \((1, 2, 3, 4, 100)\)에서 \(t\) 검정은 \(0.347\), 순열검정은 \(0.0625\)였다. 여기서는 두 검정이 전혀 수렴하지 않는다.