비모수 검정 모음¶

(1) 순위합검정을 쓰면 안 되고, 나머지 둘의 순서는 ARE 가 보장하지 않는다.

자료는 같은 학생 15명을 처치 전후로 잰 대응자료다. ranksums 는 두 열을 서로 독립인 표본으로 보므로 짝 구조를 버린다. 그 \(p = 0.1409\) 는 "보수적인 결과" 가 아니라 잘못 적용한 검정의 값이고, 다른 둘과 견줄 수 있는 수가 아니다.

남은 둘의 순서는 조금 더 조심해야 한다. ARE 가 뜻하는 것은

정규모집단에서, 귀무가설로 수렴하는 국소 대립가설 열에 대해, \(n \to \infty\) 에서 같은 검정력을 얻는 데 필요한 표본크기의 비

다. ARE \(= 0.955\) 는 "윌콕슨이 \(t\) 검정의 검정력에 맞추려면 관측값이 \(1/0.955 = 1.047\) 배 필요하다" 는 뜻이고, \(0.637\) 은 부호검정에 \(1/0.637 = 1.571\) 배가 필요하다는 뜻이다. 세 조건이 모두 붙어 있다.

그러므로 ARE 는 한 자료에서 어느 검정의 p-값이 작을지 보장하지 않는다. 이 자료에서 부호순위검정이 이긴 것은 음의 차이 둘이 하필 가장 작은 쪽(\(-2\) 와 \(-3\))이어서 \(W^- = 11\) 로 작게 나왔기 때문이다. 같은 부호 구성에서 크기만 뒤바꾸면 부호검정의 p-값은 \(0.0209\) 그대로인데 부호순위검정은 \(0.1095\) 로 올라가 순서가 뒤집힌다(일표본 비모수 검정 보기 5 가 그 자료를 만든다).

(2) 분포를 바꾸면 ARE 가 1 을 넘는다. 네 대칭분포에서 두 공식을 계산하면

정규분포는 \(t\) 검정의 고향이므로 둘 다 1 아래다. 그런데 꼬리를 무겁게 하면

• 로지스틱에서 윌콕슨 ARE 가 \(\pi^2/9 = 1.0966\) 으로 1 을 넘는다. 정규보다 꼬리가 조금 두꺼운 것만으로도 순위검정이 \(t\) 검정을 앞선다.
• 이중지수(라플라스)에서는 윌콕슨 \(1.5\), 부호검정 \(2.0\) 으로 둘 다 1 을 넘고, 순서까지 뒤집힌다. 부호검정이 윌콕슨보다 나아진다 — 이 분포에서 \(f(0)\) 이 뾰족해 중앙값 근처의 정보가 가장 값지기 때문이다.
• 반대로 균등분포에서는 꼬리가 아예 없어 부호검정이 \(1/3\) 로 최악이 되고 윌콕슨은 정확히 \(1\) 이다.

이 네 줄이 비모수 검정을 쓰는 참된 이유를 말해 준다. "정규성이 성립할 때 5% 손해" 는 최악의 경우이고, 꼬리가 조금만 무거워지면 손해가 이익으로 바뀐다. 반면 \(t\) 검정이 이기는 범위는 정규분포 주변으로 좁다.

수치적으로.

import numpy as np
from scipy import stats

paired_data = np.array([
    [93, 76], [70, 72], [81, 75], [65, 68], [79, 65],
    [54, 54], [94, 88], [91, 81], [77, 65], [65, 57],
    [95, 86], [89, 87], [78, 78], [80, 77], [76, 76]
])

post, pre = paired_data[:, 0], paired_data[:, 1]

# 부호검정(정규근사)
diffs = post - pre
nonzero = diffs[diffs != 0]
n_plus = (nonzero > 0).sum()
n = len(nonzero)
z_sign = (2 * n_plus - n) / np.sqrt(n)
p_sign = 2 * stats.norm.cdf(-abs(z_sign))
print(f"Sign test:          Z = {z_sign:.4f}, p = {p_sign:.4f}")
# Sign test:          Z = 2.3094, p = 0.0209

# Wilcoxon 부호순위검정
stat_sr, p_sr = stats.wilcoxon(post, pre, alternative="two-sided",
                                method="approx", zero_method="pratt")
print(f"Signed-rank test:   W = {stat_sr}, p = {p_sr:.4f}")
# Signed-rank test:   W = 11.0, p = 0.0086

# Wilcoxon 순위합검정 — 이 대응자료에는 맞지 않는다. 견주어 보려고 싣는다
stat_rs, p_rs = stats.ranksums(post, pre)
print(f"Rank-sum test:      Z = {stat_rs:.4f}, p = {p_rs:.4f}")
# Rank-sum test:      Z = 1.4725, p = 0.1409

출력:

Sign test:          Z = 2.3094, p = 0.0209
Signed-rank test:   W = 11.0, p = 0.0086
Rank-sum test:      Z = 1.4725, p = 0.1409

대응자료에 적절한 두 검정 중에서는 부호순위검정이 정보를 더 많이 쓰므로 \(p\)값이 작다.

from scipy import integrate

cases = {
    "정규 N(0,1)": stats.norm(),
    "로지스틱": stats.logistic(),
    "이중지수": stats.laplace(),
    "균등 U(-.5,.5)": stats.uniform(loc=-0.5, scale=1),
}
print("분포              sigma^2      f(0)   int f^2   부호 ARE  윌콕슨 ARE")
for name, dist in cases.items():
    var, f0 = dist.var(), dist.pdf(0)
    I, _ = integrate.quad(lambda x: dist.pdf(x) ** 2, -60, 60, limit=400)
    print(f"{name:<14} {var:>9.6f} {f0:>9.6f} {I:>9.6f}"
          f" {4 * var * f0 ** 2:>10.6f} {12 * var * I ** 2:>11.6f}")

print(f"\n이론값: 2/pi = {2 / np.pi:.6f},  3/pi = {3 / np.pi:.6f},"
      f"  pi^2/12 = {np.pi ** 2 / 12:.6f},  pi^2/9 = {np.pi ** 2 / 9:.6f}")

출력:

분포              sigma^2      f(0)   int f^2   부호 ARE  윌콕슨 ARE
정규 N(0,1)       1.000000  0.398942  0.282095   0.636620    0.954930
로지스틱            3.289868  0.250000  0.166667   0.822467    1.096623
이중지수            2.000000  0.500000  0.250000   2.000000    1.500000
균등 U(-.5,.5)    0.083333  1.000000  1.000000   0.333333    1.000000

이론값: 2/pi = 0.636620,  3/pi = 0.954930,  pi^2/12 = 0.822467,  pi^2/9 = 1.096623

rng = np.random.default_rng(0)
B, m = 2000, 20
delta = 0.5
print(f"n = {m},  delta = {delta},  반복 {B} 회,  alpha = 0.05")
print("분포        t 검정   윌콕슨   부호검정")
for name, draw in (("정규", lambda: rng.normal(0, 1, m)),
                   ("이중지수", lambda: rng.laplace(0, 1, m))):
    p_t, p_w, p_s = [], [], []
    for _ in range(B):
        d = draw() + delta
        p_t.append(stats.ttest_1samp(d, 0).pvalue)
        p_w.append(stats.wilcoxon(d).pvalue)
        p_s.append(stats.binomtest(int((d > 0).sum()), m).pvalue)
    rate = lambda p: np.mean(np.array(p) < 0.05)
    print(f"{name:<9} {rate(p_t):>7.3f} {rate(p_w):>8.3f} {rate(p_s):>10.3f}")
print(f"몬테카를로 오차 <= {np.sqrt(0.25 / B):.3f}")

출력:

n = 20,  delta = 0.5,  반복 2000 회,  alpha = 0.05
분포        t 검정   윌콕슨   부호검정
정규          0.566    0.552      0.392
이중지수        0.348    0.406      0.385
몬테카를로 오차 <= 0.011

네 분포의 ARE 가 이론값 \(2/\pi\), \(3/\pi\), \(\pi^2/12\), \(\pi^2/9\), 그리고 \(2.0\), \(1.5\), \(1/3\), \(1.0\) 과 소수 여섯째 자리까지 모두 맞는다. 수치적분으로 구한 \(\int f^2\) 가 로지스틱에서 \(1/6\), 이중지수에서 \(1/4\), 균등에서 \(1\) 로 떨어지는 것도 확인된다.

둘째 모의실험이 유한 \(n\) 에서 그 역전을 보여 준다. 정규에서는 $t > $ 윌콕슨 \(>\) 부호(\(0.566 / 0.552 / 0.392\))인데 이중지수에서는 \(t\) 가 꼴찌로 내려간다(\(0.348\)). 다만 역전이 절반만 일어났다 — ARE 가 예측하는 "부호 \(>\) 윌콕슨" 까지는 가지 않고 \(0.385\) 대 \(0.406\) 이다(몬테카를로 오차 \(\pm 0.011\)). 이상한 일이 아니다. ARE 는 \(\delta \to 0\), \(n \to \infty\) 의 극한값이고 여기서는 \(n = 20\), \(\delta = 0.5\) 로 둘 다 극한과 멀기 때문이다. 유한 \(n\) 에서는 부호검정의 이산성 손실(기각역을 촘촘히 채울 수 없는 것)이 ARE 의 이득을 상당 부분 먹는다.

부호검정을 써야 한다. Wilcoxon 부호순위검정은 \(H_0\) 아래에서 차이의 분포가 중앙값을 중심으로 대칭임을 요구한다. 이 가정이 깨지면 \(W^+\)의 귀무분포가 틀리고 제1종 오류가 부풀려질 수 있다.

맞교환은 검정력이다. 부호검정은 차이의 부호만 쓰고 크기 정보를 버리므로, 대칭성이 실제로 성립한다면 검정력이 낮다. 그러나 이 경우에는 연구자가 대칭성을 확인할 수 없으므로 제1종 오류를 통제하는 안전성이 검정력보다 중요하다.

표본이 더 크다면 차이의 히스토그램으로 대칭성을 평가한 뒤 부호순위검정으로 바꿀 수 있다. 다만 \(n = 12\)로는 히스토그램으로 대칭성을 판정하기 어렵다. 맥락 지식(예: "이 측정값은 바닥효과가 있어 차이가 오른쪽으로 치우친다")이 있다면 그것이 더 믿을 만한 근거이다. \(\square\)

ARE는 두 검정의 효능(efficacy) 제곱의 비이다. 검정통계량 \(T_n\)의 효능은

으로 정의된다.

\(t\) 검정. \(T_n = \bar{D}\)이고 \(\operatorname{E}_\delta[\bar D] = \delta\)이므로 미분값이 \(1\)이다. \(\operatorname{Var}_0(\bar D) = \sigma^2/n\)이므로

부호검정. \(T_n = \hat{p} = n_+/n\)이라 하자. \(D_i \sim \mathcal{N}(\delta, \sigma^2)\)이면

이고, \(\delta = 0\)에서 미분하면

이다(\(f\)는 \(D_i\)의 밀도함수). \(\operatorname{Var}_0(\hat p) = 1/(4n)\)이므로

ARE.

일반적으로 밀도 \(f\)와 표준편차 \(\sigma\)를 갖는 대칭분포에서 \(\text{ARE}(\text{부호}, t) = 4\sigma^2 f(0)^2\)이다. 이 공식이 유용한 것은 다른 분포에서도 바로 계산할 수 있기 때문이다.

라플라스에서 부호검정이 \(t\) 검정보다 두 배 효율적임에 주목하라. 꼬리가 두꺼울수록 \(f(0)\)이 커지고 부호검정이 유리해진다.

from scipy import stats

a = [5, 8, 12, 15]
b = [7, 11, 14, 18, 20]
c = [3, 6, 9]

print(stats.kruskal(a, b, c))
# KruskalResult(statistic=4.0295, pvalue=0.13335)

r = stats.median_test(a, b, c)
print(r.statistic, r.pvalue, r.median)
print(r.table)
# 4.8  0.09072  10.0
# [[2 4 0]
#  [2 1 3]]

출력:

KruskalResult(statistic=4.029487179487184, pvalue=0.13335459246779172)
4.8 0.0907179532894125 10.0
[[2 4 0]
 [2 1 3]]

놀랍게도 Mood 중앙값검정의 \(p\)값이 더 작다. "Kruskal--Wallis가 언제나 더 강력하다"는 통념이 이 자료에서는 성립하지 않는다.

이유는 이 자료의 구조에 있다. 전체 중앙값은 \(10\)이고, 집단 C는 값이 \(3, 6, 9\)로 전부 중앙값 아래이다. Mood 검정에게 이는 매우 강한 신호이다(\(0\) 대 \(3\)). 반면 Kruskal--Wallis는 집단 C의 순위가 \(1, 2, 4\)로 낮긴 하지만 집단 A의 \(3, 5, 8, 10\)과 상당히 겹친다는 것도 함께 본다.

더 근본적으로, 표본이 \(4, 5, 3\)으로 매우 작아 두 검정 모두 \(\chi^2\) 근사가 믿을 만하지 않다. 이 크기에서 \(p\)값의 순서를 두고 검정력을 논하는 것은 무리이다.

교훈: "검정 A가 검정 B보다 강력하다"는 진술은 언제나 분포와 대립가설의 종류에 대한 평균을 두고 하는 말이다. 개별 자료에서는 어느 쪽이든 이길 수 있다. 자료를 본 뒤 \(p\)값이 작은 검정을 고르면 실제 유의수준이 명목값을 크게 넘는다. \(\square\)

Kruskal--Wallis 검정의 귀무가설은 \(H_0{:}\; F_1 = F_2 = \cdots = F_k\), 즉 모든 집단이 같은 분포에서 왔다는 것이다. 이는 분포의 동일성에 대한 검정이지 중앙값이 같은지에 대한 검정이 아니다. \(H_0\)을 기각한다는 것은 집단들이 위치, 척도, 모양 또는 그 조합에서 다르다는 뜻일 수 있다.

분포들이 위치만 다르다는 위치이동 가정(\(F_j(x) = F(x - \mu_j)\), \(F\)는 공통 모양)을 추가하면 분포가 다를 수 있는 유일한 통로가 평균/중앙값이 되므로, 이 경우 Kruskal--Wallis가 중앙값(또는 평균) 동일성 검정이 된다.

위치이동 가정이 없으면 유의한 Kruskal--Wallis 결과가 중심이 아니라 산포의 차이를 반영할 수 있고, 이를 중앙값 검정으로 해석하는 연구자를 오도한다.

Kruskal-Wallis 연습문제 3에서 이를 정량적으로 확인했다. 중앙값이 모두 정확히 0이지만 왜도가 다른 세 집단에서 \(n = 100\)일 때 기각률이 \(0.614\)까지 올라간다. \(\square\)

Dunn 검정은 원래의 Kruskal--Wallis 합친 순위에서 얻은 평균순위를 이용해 \(\binom{k}{2}\)개의 모든 집단 쌍을 비교한다. 평균순위가 \(\bar{R}_i\), \(\bar{R}_j\)인 집단 \(i\)와 \(j\)에 대해

이며 \(N\)은 전체 표본크기이다. 각 \(Z_{ij}\)를 표준정규분포와 비교한다.

다중비교 보정(Bonferroni, Holm, Benjamini--Hochberg 등)이 필요한 이유는 \(\binom{k}{2}\)개의 검정을 수행하면 집단별 오류율이 부풀려지기 때문이다. \(m\)개의 쌍을 각각 유의수준 \(\alpha\)로 검정하면 검정들이 독립일 때 적어도 하나가 거짓 기각될 확률이 \(1 - (1 - \alpha)^m\)까지 올라간다.

\(k = 8\)이면 모든 집단이 동일해도 76%의 확률로 "유의한" 쌍을 하나 이상 찾게 된다.

Bonferroni 보정은 각 쌍별 검정에 \(\alpha/m\)을 쓰므로 전체 집단별 오류율을 \(\alpha\) 이하로 보장한다. Holm 절차는 FWER을 똑같이 통제하면서 균일하게 더 강력하므로 대체로 더 낫다.