대응자료에 대한 순열검정¶

(1) 해석적으로. 독립인 부호 \(\varepsilon_i \in \{-1, +1\}\)을 각각 확률 \(1/2\)로 붙이므로 \(E[\varepsilon_i] = 0\), \(\operatorname{Var}(\varepsilon_i) = 1\)이고

이다. \(\sum_i d_i^2 = 449\)이므로

지지집합. 삼각부등식에서 \(\lvert \bar d^{*}\rvert \le \frac1n\sum_i \lvert d_i\rvert = \overline{\lvert d\rvert}\)이고, 등호는 \(\varepsilon_i d_i\)의 부호가 모두 같을 때 성립한다. \(\sum_i \lvert d_i \rvert = 63\)이므로 지지집합은 \([-5.25,\ 5.25]\)이고 양 끝은 실제로 도달한다. 순열분포는 유계다. \(t\) 검정이 쓰는 \(t_{11}\)에는 그런 경계가 없다.

(2) 해석적으로. \(\varepsilon_i d_i = \pm\lvert d_i\rvert\)이고 부호가 자유로우므로, 부호 배정을 "어느 \(\lvert d_i\rvert\)를 음수로 둘 것인가"로 바꾸어 적을 수 있다. 음수로 둔 자리의 집합을 \(A\), 그 합을 \(s = \sum_{i \in A}\lvert d_i\rvert\)라 하면

이다. 관측값이 \(\bar d_{\text{obs}} = -4.25\), 곧 합으로는 \(-51\)이므로

이다. \(s \ge 57\)은 여집합의 합이 \(6\) 이하라는 뜻이므로 \(A \mapsto A^c\)가 두 쪽을 일대일로 맞바꾼다. 한쪽만 세면 된다.

\(\lvert d\rvert = (8, 3, 2, 12, 5, 1, 4, 7, 6, 2, 9, 4)\)에서 합이 \(6\) 이하인 부분집합을 센다. \(6\)을 넘는 \(7, 8, 9, 12\)는 아예 들어갈 수 없고, 값 \(2\)와 값 \(4\)는 두 개씩 있다는 데 주의한다.

합이 \(25\)다. 양쪽을 더하면

수치적으로. \(n = 12\)이므로 \(2^{12} = 4{,}096\)가지 부호 배정을 모두 열거한다.

import numpy as np, itertools
from scipy import stats

d = np.array([-8, -3, 2, -12, -5, -1, 4, -7, -6, -2, -9, -4])
print(d.mean(), round(d.std(ddof=1), 4))     # -4.25  4.595

S = np.array(list(itertools.product([1, -1], repeat=12)))
m = (S * d).mean(axis=1)
count = (np.abs(m) >= abs(d.mean()) - 1e-12).sum()
print(count, len(m), count / len(m))         # 50  4096  0.012207

print(stats.ttest_1samp(d, 0))               # t = -3.204, p = 0.00839
print(stats.wilcoxon(d))                     # W = 8.0,   p = 0.012207

출력:

-4.25 4.595
50 4096 0.01220703125
TtestResult(statistic=-3.2040407749287168, pvalue=0.008394085815475436, df=11)
WilcoxonResult(statistic=8.0, pvalue=0.01220703125)

(1)과 (2)의 수를 모두 확인한다. 위 블록의 변수를 그대로 이어 쓴다.

import collections

n = len(d)
print(f"sum|d| = {np.abs(d).sum():.0f},  mean|d| = {np.abs(d).mean():.2f},"
      f"  sum d^2 = {(d ** 2).sum():.0f}")
print(f"귀무 SD 닫힌 꼴 = {np.sqrt((d ** 2).sum()) / n:.6f},  열거 = {m.std():.6f}")
print(f"지지집합의 끝   = {np.abs(m).max():.2f}")
print(f"t 검정의 척도   = {d.std(ddof=1) / np.sqrt(n):.6f}")

# |합| >= 51 은 뒤집는 쪽 |d_i| 들의 합 s 가 6 이하이거나 57 이상인 것과 같다.
cnt = collections.Counter()
for k in range(n + 1):
    for c in itertools.combinations(range(n), k):
        s = int(sum(abs(d[i]) for i in c))
        if s <= 6:
            cnt[s] += 1
print(f"합이 s 인 부분집합의 수: {sorted(cnt.items())}"
      f"  ->  s <= 6 은 {sum(cnt.values())} 가지")

# 두 통계량이 고르는 배정이 같은 집합인가
r = stats.rankdata(np.abs(d))
Wp = ((S * d > 0) * r).sum(axis=1)
W_obs = ((d > 0) * r).sum()
sel_mean = np.abs(m) >= abs(d.mean()) - 1e-12
sel_w = (Wp <= W_obs + 1e-12) | (Wp >= r.sum() - W_obs - 1e-12)
print(f"W+ obs = {W_obs:.1f},  순위합 전체 = {r.sum():.0f},"
      f"  Wilcoxon 이 고른 배정 = {sel_w.sum()}")
print(f"두 집합이 같은가: {np.array_equal(sel_mean, sel_w)}")

출력:

sum|d| = 63,  mean|d| = 5.25,  sum d^2 = 449
귀무 SD 닫힌 꼴 = 1.765802,  열거 = 1.765802
지지집합의 끝   = 5.25
t 검정의 척도   = 1.326450
합이 s 인 부분집합의 수: [(0, 1), (1, 1), (2, 2), (3, 3), (4, 4), (5, 6), (6, 8)]  ->  s <= 6 은 25 가지
W+ obs = 8.0,  순위합 전체 = 78,  Wilcoxon 이 고른 배정 = 50
두 집합이 같은가: True

손으로 센 것이 그대로 맞는다. 합 \(0\)부터 \(6\)까지의 개수가 \(1, 1, 2, 3, 4, 6, 8\)로 표와 한 칸도 다르지 않고, 합이 \(25\)다. 따라서 \(4{,}096\)가지 중 \(50\)가지가 걸리고 \(p_{\text{정확}} = 50/4096 = 0.012207\)이다. 닫힌 꼴 표준편차 \(1.765802\)와 지지집합의 끝 \(5.25\)도 열거와 일치한다.

Wilcoxon 부호순위검정이 정확히 같은 답을 주는 까닭도 드러난다. \(W^+ \le 8\) 또는 \(W^+ \ge 70\)인 배정이 \(50\)개인데, 코드가 확인하듯 그것이 평균 통계량이 고른 \(50\)개와 같은 집합이다. 두 검정이 같은 \(2^{12}\)가지 위에서 정의되고 이 자료에서는 극단 배정의 순서까지 겹치기 때문이다. 통계량이 달라도 기각역이 일치할 수 있다.

대응 \(t\) 검정만 다른 답을 준다. \(t = -4.25/1.326450 = -3.204\)에서 \(p = 0.0084\)로 순열검정의 \(0.0122\)보다 작다. 이 자료의 왜도는 \(0.197\)로 거의 대칭이므로 원인은 치우침이 아니다. 두 가지가 겹친다.

• 척도가 다르다. 순열은 \(\sqrt{\sum d_i^2}/n = 1.766\)으로 재고 \(t\) 검정은 \(s_d/\sqrt n = 1.326\)으로 잰다. 순열검정 보기 3에서 본 대로 비가 \(\sqrt{(n-1+t^2)/n} = \sqrt{(11 + 10.266)/12} = 1.331\)이다.
• 꼬리가 다르다. 순열분포는 \(\pm 5.25\)에서 끊기는데 \(t_{11}\)은 끝없이 뻗는다.

더 넓고 유계인 귀무분포에서 같은 관측값을 재니 \(p\)값이 커진다. 순열검정의 \(0.0122\)는 이 자료에서 정확한 값이고, \(t\) 검정의 \(0.0084\)는 정규성을 빌려 온 근사다.

import numpy as np, itertools
from scipy import stats

before = np.array([148, 142, 136, 155, 129, 161, 138, 144])
after  = np.array([140, 138, 132, 147, 131, 152, 135, 139])
d = after - before
print(d)                      # [-8 -4 -4 -8  2 -9 -3 -5]
print(d.mean(), round(d.std(ddof=1), 4))   # -4.875  3.5632

S = np.array(list(itertools.product([1, -1], repeat=8)))
m = (S * d).mean(axis=1)
cnt = (np.abs(m) >= abs(d.mean()) - 1e-12).sum()
print(cnt, len(m), cnt / len(m))           # 4  256  0.015625

print(stats.ttest_rel(after, before))
print(stats.wilcoxon(d, method='exact'))
print(stats.binomtest((d > 0).sum(), 8, 0.5).pvalue)

출력:

[-8 -4 -4 -8  2 -9 -3 -5]
-4.875 3.5632
4 256 0.015625
TtestResult(statistic=-3.8697136256193168, pvalue=0.006133954352036696, df=7)
WilcoxonResult(statistic=1.0, pvalue=0.015625)
0.0703125

(a) 차이는 \((-8, -4, -4, -8, 2, -9, -3, -5)\)이고 \(\bar{d} = -4.875\)이다. \(256\)가지 부호 배정 중 \(|\bar{d}^*| \ge 4.875\)를 만족하는 것은 \(4\)가지이므로

(b) 네 검정의 비교

순열검정과 Wilcoxon이 또 일치한다. 본문 보기와 같은 현상이다. 두 통계량이 같은 \(4\)개의 배정을 극단으로 지목한다.

\(t\) 검정이 \(2.5\)배 작은 \(p\)값을 준다(\(0.0061\) 대 \(0.0156\)). 본문 보기와 방향이 같다. \(n = 8\)에서 \(t_7\)의 꼬리 근사가 순열분포보다 관대하다.

부호검정은 \(p = 0.070\)으로 기각하지 못한다. 부호검정은 차이의 크기를 완전히 버리고 방향만 세기 때문이다. 여덟 명 중 일곱 명이 내려갔다는 사실만으로는 \(\alpha = 0.05\)에 못 미친다. 반면 순열검정과 Wilcoxon은 "내려간 사람들의 하락폭이 유일하게 올라간 사람의 상승폭(\(+2\))보다 훨씬 크다"는 정보를 쓴다.

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

def signflip_p(d, stat, B=499):
    obs = abs(stat(d))
    S = rng.choice([-1, 1], size=(B, len(d)))
    v = np.abs(np.apply_along_axis(stat, 1, S * d))
    return ((v >= obs - 1e-12).sum() + 1) / (B + 1)

def wil(a):
    r = stats.rankdata(np.abs(a))
    return (r * (a > 0)).sum() - (r * (a < 0)).sum()
sign = lambda a: (a > 0).sum() - (a < 0).sum()
trim = lambda a: stats.trim_mean(a, 0.2)

제1종 오류율 (이동 = 0)

네 순열검정 모두 크기를 지킨다. 부호 통계량이 \(0.032\)로 보수적인데, 이는 부호검정의 귀무분포가 매우 이산적이기 때문이다(\(n = 15\)에서 이항분포의 눈금).

검정력 (이동 = 0.8)

정규 차이에서는 평균이 최선이고 \(t\) 검정과 사실상 같다(\(0.811\) 대 \(0.823\)). Wilcoxon은 \(0.795\)로 \(2\)%p만 뒤진다. 이는 이론이 예측하는 바이다. 정규분포에서 Wilcoxon 부호순위검정의 \(t\) 검정 대비 점근상대효율은 \(3/\pi \approx 0.955\)이다.

부호 통계량은 정규자료에서 \(0.595\)로 크게 뒤진다. 크기 정보를 버리는 대가이다. 점근상대효율 \(2/\pi \approx 0.637\)과 부합한다.

오염 차이에서 순위가 뒤집힌다. 절사평균이 \(0.644\)로 최선이고, Wilcoxon \(0.604\), 평균 \(0.540\), \(t\) 검정 \(0.493\) 순이다. 평균 기반 검정이 \(0.811 \to 0.540\)으로 \(33\)%의 검정력을 잃는 동안 Wilcoxon은 \(0.795 \to 0.604\)로 \(24\)%만 잃는다.

주목할 점은 부호 통계량이 오염자료에서 \(0.507\)로 거의 손실이 없다는 것이다(\(0.595 \to 0.507\)). 크기를 아예 보지 않으므로 이상값의 영향을 받지 않는다. 그럼에도 절사평균에 크게 뒤지는 것은 버리는 정보가 너무 많기 때문이다.

결론: 자료의 꼬리를 모를 때는 절사평균이나 Wilcoxon이 합리적인 기본값이다. 두 상황 모두에서 최선에 가깝다.

import numpy as np
rng = np.random.default_rng(21)

def signflip_p(d, B=199, plus1=True):
    obs = abs(d.mean())
    S = rng.choice([-1, 1], size=(B, len(d)))
    v = np.abs((S * d).mean(axis=1))
    if plus1:
        return ((v >= obs - 1e-12).sum() + 1) / (B + 1)
    return (v >= obs - 1e-12).mean()

a = b = 0
M = 4000
for _ in range(M):
    d = rng.normal(0, 1, 12)
    a += signflip_p(d, plus1=True) < 0.05
    b += signflip_p(d, plus1=False) < 0.05
print(a/M, b/M)

출력:

0.04575 0.05075

보정 없는 정의가 명목수준을 넘는다(\(0.0512 > 0.05\)). 보정을 넣으면 \(0.0458\)로 \(\alpha\) 이하가 된다.

차이가 \(0.005\)로 작아 보이지만 방향이 중요하다. 보정 없는 검정은 원리적으로 크기를 보장하지 못한다. 보정된 검정은 어떤 \(B\)에서도 \(P(\hat{p} \le \alpha) \le \alpha\)를 만족한다.

왜 그런가. \(H_0\) 아래에서 관측된 배열은 \(B\)개의 무작위 배열과 교환 가능하다. 즉 관측 통계량의 순위는 \(B+1\)개 값 중에서 균등분포를 따른다. 따라서

가 정확히 성립한다. 여기서 \(R\)은 관측값 이상인 순열의 개수이다. 관측값을 세지 않으면 이 논증이 무너진다.

\(B\)가 클수록 차이가 줄어든다. \(B = 199\)에서 \(0.005\)인 차이가 \(B = 9999\)에서는 \(0.0001\) 수준이다. 그래도 보정을 넣는 데 드는 비용이 \(0\)이므로 항상 넣는 것이 맞다.

import numpy as np
from scipy import stats
rng = np.random.default_rng(21)
# signflip_p 는 연습문제 2 의 것과 같고, 통계량으로 평균을 쓴다

def size(gen, n, M=3000):
    a = b = 0
    for _ in range(M):
        d = gen(n)
        a += signflip_p(d, np.mean, B=499) < 0.05
        b += stats.ttest_1samp(d, 0).pvalue < 0.05
    return round(a/M, 3), round(b/M, 3)

exp_c = lambda k: rng.exponential(1, k) - 1.0
ln_c  = lambda k: rng.lognormal(0, 1, k) - np.exp(0.5)

부호 뒤집기 검정의 크기가 무너진다. 지수 차이에서 \(0.092\), 로그정규 차이에서 \(0.146\)으로 명목값의 \(2\)--\(3\)배이다.

"정확검정"이 정확성을 잃은 것이 아니다. 다른 가설을 검정하고 있는 것이다. 부호 뒤집기 검정의 실제 귀무가설은

이며, \(E[d_i] = 0\)이지만 치우친 분포에서 이 가설은 거짓이다. 검정이 참인 대칭성 위반을 정확히 탐지하고 있으며, 우리가 그것을 "평균이 0이 아니다"로 잘못 읽는 것이다.

\(t\) 검정도 똑같이 나쁘다. \(0.095\), \(0.150\)으로 사실상 구별되지 않는다. 이 상황에서 순열검정을 쓴다고 해서 더 나빠지지는 않는다.

\(n\)을 늘리면 개선된다. 지수 차이에서 \(n = 15\)의 \(0.092\)가 \(n = 40\)에서 \(0.071\)로 줄어든다. 중심극한정리가 작동하기 때문이다. 다만 수렴이 느리다. 왜도가 \(\gamma_1\)일 때 \(t\) 통계량의 왜도는 대략 \(\gamma_1/\sqrt{n}\)로 줄어들므로, \(\text{Exp}(1)\)의 \(\gamma_1 = 2\)를 \(0.2\) 아래로 낮추려면 \(n \ge 100\)이 필요하다.

어떻게 해야 하는가.