사후 Dunn 검정¶

쌍별 비교:

\(\alpha = 0.05\)에서의 해석:

• B 대 C: \(p^* = 0.0018 < 0.05\). 비료 B와 C가 유의하게 다르다. 비료 B가 더 큰 식물을 낸다(평균순위가 높다).
• A 대 B: \(p^* = 0.183 > 0.05\). A와 B 사이에 유의한 차이가 없다.
• A 대 C: \(p^* = 0.359 > 0.05\). A와 C 사이에 유의한 차이가 없다.

Kruskal-Wallis의 기각은 주로 비료 B와 비료 C의 큰 차이가 이끌어 낸 것이다.

보정하지 않은 \(p\)값을 작은 것부터 정렬한다. \(m = 3\)이다.

Bonferroni: 모두 \(\times 3\).

Holm: \(r\)번째 작은 값에 \((m - r + 1)\)을 곱하고, 앞선 값보다 작아지지 않도록 누적 최댓값을 취한다.

Benjamini-Hochberg: \(r\)번째 값에 \(m/r\)을 곱하고, 뒤에서부터 누적 최솟값을 취한다.

import numpy as np
from scipy import stats
from statsmodels.stats.multitest import multipletests

p = np.array([0.000604, 0.061, 0.1198])   # B-C, A-B, A-C
for m in ("bonferroni", "holm", "fdr_bh"):
    print(m, multipletests(p, method=m)[1].round(5))

출력:

bonferroni [0.00181 0.183   0.3594 ]
holm [0.00181 0.122   0.122  ]
fdr_bh [0.00181 0.0915  0.1198 ]

결론은 세 방법 모두 같다. B 대 C만 유의하다.

그러나 \(p\)값의 크기는 뚜렷이 다르다. A 대 C에서 Bonferroni는 \(0.359\), Holm은 \(0.122\), BH는 \(0.120\)이다. Bonferroni가 가장 보수적이고 BH가 가장 관대하다.

Holm이 Bonferroni를 지배하는 이유. Holm은 가장 작은 \(p\)값에만 \(m\)을 곱하고, 그다음부터는 \(m-1\), \(m-2\)로 줄여 간다. 이미 기각된 가설이 있으면 남은 비교에 대한 "예산"이 줄어들 필요가 없다는 논리이다. 어떤 상황에서도 Bonferroni보다 작거나 같은 보정값을 주므로 무조건 더 강력하다.

BH는 다른 것을 통제한다. FWER(거짓 발견을 하나라도 낼 확률)이 아니라 FDR(발견 중 거짓의 기대 비율)을 통제한다. 비교가 수십 개 이상이고 탐색적 분석이라면 BH가 적절하지만, 확증적 연구에서는 FWER 통제가 표준이다.

import numpy as np, itertools
from scipy import stats
from statsmodels.stats.multitest import multipletests

A = [12, 15, 14, 10, 13]; B = [20, 18, 22, 17, 19]; C = [8, 11, 9, 7, 10]
groups = {"A": A, "B": B, "C": C}
pairs = list(itertools.combinations(groups, 2))
pv = [stats.mannwhitneyu(groups[i], groups[j], method='exact').pvalue
      for i, j in pairs]
print(pairs)
print(np.round(pv, 5))
print(np.round(multipletests(pv, method='bonferroni')[1], 4))

출력:

[('A', 'B'), ('A', 'C'), ('B', 'C')]
[0.00794 0.03175 0.00794]
[0.0238 0.0952 0.0238]

결론이 갈린다. Dunn 검정은 B 대 C만 유의하다고 하지만, 쌍별 Mann-Whitney는 A 대 B와 B 대 C를 유의하다고 한다.

차이의 원인은 두 가지이다.

1. 순위를 다시 매기는가. Mann-Whitney는 각 쌍에 대해 \(10\)개 관측값만 놓고 순위를 다시 매긴다. A와 B는 완전히 분리되어 있으므로(\(A\)의 최댓값 \(15 < B\)의 최솟값 \(17\)) \(U = 0\)이 되어 가능한 최소 \(p\)값 \(2/\binom{10}{5} = 2/252 = 0.0079\)가 나온다. B와 C도 마찬가지로 완전 분리이다. Dunn 검정은 전역 순위를 쓰므로 A의 순위(\(4.5\)--\(10\))와 B의 순위(\(11\)--\(15\))가 겹치지 않는다는 사실이 평균순위 차이 \(5.3\)으로만 반영된다.

2. 정확검정 대 정규근사. Mann-Whitney는 \(n = 5, 5\)에서 정확 \(p\)값을 쓰지만 Dunn 검정은 정규근사를 쓴다. \(N = 15\)는 정규근사에 작다.

어느 쪽이 옳은가? 정답은 없지만 주의할 점이 있다.

• 쌍별 Mann-Whitney는 각 쌍에서 다른 순위 체계를 쓰므로 결과가 서로 모순될 수 있다(비이행성).
• Dunn 검정은 전역 순위를 쓰므로 일관되지만, 관계없는 세 번째 집단이 두 집단의 비교에 영향을 준다. 여기서 집단 C의 존재가 A와 B의 순위 간격을 넓혀 놓았다.

실무 권고: 전체검정과 사후검정의 논리적 일관성을 중시하면 Dunn을, 각 쌍의 정확한 비교를 중시하면 쌍별 Mann-Whitney를 쓴다. 둘 다 돌려 보고 유의한 쪽을 고르는 것은 명백한 오용이다.

귀무가설이 완전히 참인(\(k\)개 집단이 모두 같은 분포) 상황에서, 어느 쌍이라도 유의하다고 선언할 확률을 센다.

import numpy as np, itertools
from scipy import stats
from statsmodels.stats.multitest import multipletests
rng = np.random.default_rng(0)
B, n = 1500, 15

def dunn_p(groups):
    N = sum(len(g) for g in groups)
    r = stats.rankdata(np.concatenate(groups))
    idx, mr = 0, []
    for g in groups:
        mr.append(r[idx:idx + len(g)].mean()); idx += len(g)
    out = []
    for i, j in itertools.combinations(range(len(groups)), 2):
        sd = np.sqrt(N * (N + 1) / 12 * (1/len(groups[i]) + 1/len(groups[j])))
        out.append(2 * stats.norm.sf(abs(mr[i] - mr[j]) / sd))
    return np.array(out)

for k in (3, 5, 8):
    gated = unrestricted = 0
    for _ in range(B):
        gs = [rng.normal(0, 1, n) for _ in range(k)]
        pw = multipletests(dunn_p(gs), method='holm')[1]
        any_sig = (pw < 0.05).any()
        unrestricted += any_sig
        gated += any_sig and (stats.kruskal(*gs).pvalue < 0.05)
    print(k, gated / B, unrestricted / B)

출력:

3 0.03333333333333333 0.034666666666666665
5 0.034 0.04
8 0.02266666666666667 0.03333333333333333

두 절차 모두 \(\alpha = 0.05\) 아래로 FWER을 유지한다. Holm 보정만으로도 이미 통제되기 때문이다.

그렇다면 전체검정 게이트는 왜 필요한가? 세 가지 이유가 있다.

1. 보정을 하지 않는 경우를 막는다. 실무에서 사후 비교의 다중성 보정을 잊는 일이 흔하다. 게이트가 최소한의 방어선이 된다.
2. 보정된 검정을 더 보수적으로 만든다. 위 표에서 게이트 있는 쪽이 언제나 더 낮고, \(k = 8\)에서는 \(0.023\) 대 \(0.033\)으로 격차가 벌어진다.
3. 논리적 일관성. 전체검정이 기각하지 못했는데 어떤 쌍이 유의하다고 보고하면 독자가 혼란스럽다.

다만 게이트에도 대가가 있다. 전체검정이 놓치는 상황 --- 예를 들어 \(k = 8\) 중 두 집단만 크게 다른 경우 --- 에서 게이트가 참인 발견을 막을 수 있다. Kruskal-Wallis는 모든 집단의 차이를 평균하므로 국소적인 큰 차이를 희석하기 때문이다.

\(H_0\) 아래에서 \(1, \ldots, N\)의 순위가 집단들에 무작위로 배정된다. 이는 크기 \(N\)인 유한모집단 \(\{1, \ldots, N\}\)에서 비복원추출하는 것과 같다.

유한모집단의 모수. 순위 전체의 평균과 분산은

이다.

평균순위의 분산. 크기 \(n_i\)인 비복원 표본평균의 분산은 유한모집단 보정을 포함하여

이다.

공분산. 두 집단이 같은 순위 집합을 나누어 가지므로 음의 상관이 있다. 비복원추출에서

이다.

차이의 분산.

import numpy as np
rng = np.random.default_rng(0)
N, ni, nj = 15, 5, 5
ranks = np.arange(1, N + 1)
d = []
for _ in range(200000):
    perm = rng.permutation(ranks)
    d.append(perm[:ni].mean() - perm[ni:ni + nj].mean())
print(np.var(d), N * (N + 1) / 12 * (1 / ni + 1 / nj))   # 8.049, 8.0

출력:

8.049307944999 8.0

모의실험 분산 \(8.049\)가 공식값 \(8.0\)과 일치한다.

유도에서 공분산 항이 결정적이다. 이를 빠뜨리고 두 집단이 독립이라고 가정하면 분산이 \(\frac{N+1}{12}(\frac{N-n_i}{n_i} + \frac{N-n_j}{n_j}) = 5.33\)으로 나와 \(33\%\) 과소평가된다. 그러면 \(z\)가 부풀려져 검정이 지나치게 많이 기각하게 된다.