Kendall의 tau (재론)¶

검정:

\(\alpha = 0.05\)에서 \(H_0\)을 기각하지 못한다. 관측값이 6개뿐이라 중간 정도의 일치도를 탐지할 검정력이 부족하다.

from scipy import stats
R = [1, 2, 3, 4, 5, 6]; S = [2, 1, 4, 3, 6, 5]
print(stats.kendalltau(R, S))
# SignificanceResult(statistic=0.6, pvalue=0.13611)

SignificanceResult(statistic=0.6, pvalue=0.1361111111111111)

(a) \(n = 8\)이므로 쌍이 \(\binom{8}{2} = 28\)개이다. 수학 순위가 이미 \(1, 2, \ldots, 8\)이므로, \(i < j\)인 쌍은 학생 \(j\)의 영어 순위가 학생 \(i\)의 것보다 클 때 부합이다.

수학 순위 순서로 나열한 영어 순위: 3, 1, 2, 5, 4, 8, 6, 7.

각 학생에 대해 뒤에 나오는 영어 순위 중 더 큰 것의 개수를 센다.

• 학생 1 (영어 3): 뒤의 1, 2, 5, 4, 8, 6, 7 중 큰 것은 5, 4, 8, 6, 7 → \(C_1 = 5\), \(D_1 = 2\)
• 학생 2 (영어 1): 뒤의 2, 5, 4, 8, 6, 7 중 큰 것은 전부 6개 → \(C_2 = 6\), \(D_2 = 0\)
• 학생 3 (영어 2): 뒤의 5, 4, 8, 6, 7 중 큰 것은 전부 5개 → \(C_3 = 5\), \(D_3 = 0\)
• 학생 4 (영어 5): 뒤의 4, 8, 6, 7 중 큰 것은 8, 6, 7 → \(C_4 = 3\), \(D_4 = 1\)
• 학생 5 (영어 4): 뒤의 8, 6, 7 중 큰 것은 전부 3개 → \(C_5 = 3\), \(D_5 = 0\)
• 학생 6 (영어 8): 뒤의 6, 7 중 큰 것은 없음 → \(C_6 = 0\), \(D_6 = 2\)
• 학생 7 (영어 6): 뒤의 7 중 큰 것은 7 → \(C_7 = 1\), \(D_7 = 0\)

합계: \(C = 5 + 6 + 5 + 3 + 3 + 0 + 1 = 23\), \(D = 2 + 0 + 0 + 1 + 0 + 2 + 0 = 5\).

검산: \(C + D = 28 = \binom{8}{2}\). \(\checkmark\)

(b)

import numpy as np
from scipy import stats
math = np.arange(1, 9); eng = np.array([3, 1, 2, 5, 4, 8, 6, 7])
print(stats.kendalltau(math, eng))
# SignificanceResult(statistic=0.64286, pvalue=0.031151)

출력:

SignificanceResult(statistic=0.6428571428571428, pvalue=0.03115079365079365)

(c) 이 자료에서 \(\tau = 0.643\)인 반면 \(r_s = 0.833\)이다. 같은 자료에 대해 Kendall tau가 Spearman rho보다 작은 것은 두 계수가 다른 척도를 쓰기 때문이다.

• \(r_s\)는 순위차이의 제곱을 쓰므로 큰 불일치를 증폭한다.
• \(\tau\)는 쌍별 부합 여부라는 이진 판정을 세므로, 순위가 1칸 어긋나든 7칸 어긋나든 비부합 하나로 똑같이 센다.

이변량 정규분포에서 Pearson 상관 \(\rho\)와의 정확한 관계는

이다. 이 둘에서 \(\rho\)를 소거하면 \(\tau\)와 \(r_s\)의 관계가 나오며, \(|\rho|\)가 작을 때는 \(\tau \approx \frac{2}{3} r_s\)로 근사된다.

크기는 다르지만 두 측도 모두 강한 양의 단조 연관성을 가리킨다. \(p\)값도 \(\tau\)가 \(0.0312\)(정확), \(r_s\)가 \(0.0154\)(정확)로 같은 결론을 준다.

존재한다. \(n = 6\)에서 \(720\)개 순열을 모두 열거하여 확인한다.

import numpy as np, itertools
from scipy import stats

n = 6
R = np.arange(1, n + 1)
rows = [(p, stats.kendalltau(R, p).statistic, stats.spearmanr(R, p).statistic)
        for p in itertools.permutations(R)]

cnt, first = 0, None
for (pa, ta, sa), (pb, tb, sb) in itertools.combinations(rows, 2):
    if (ta - tb) * (sa - sb) < 0:
        cnt += 1
        if first is None:
            first = (pa, round(ta, 3), round(sa, 3), pb, round(tb, 3), round(sb, 3))
print(first)
print(cnt, len(rows) * (len(rows) - 1) // 2)

출력:

((1, 2, 4, 5, 6, 3), 0.6, 0.657, (1, 3, 2, 6, 5, 4), 0.467, 0.714)
3062 258840

첫 번째로 발견되는 역전 사례:

첫 번째 자료가 \(\tau\)로는 더 강한 연관성을 보이지만 \(r_s\)로는 더 약하다. 순서가 완전히 뒤집힌다.

\(720\)개 순열에서 만들 수 있는 \(258{,}840\)개 쌍 중 \(3{,}062\)개(\(1.18\%\))에서 이런 역전이 일어난다.

왜 이런 일이 생기는가. 두 자료의 구조를 보자.

• \((1, 2, 4, 5, 6, 3)\): \(6\)이 마지막에서 세 번째로 밀려 있다. 즉 한 원소가 크게 어긋나 있다. 비부합쌍은 \((4,3)\), \((5,3)\), \((6,3)\)의 3개뿐이지만 순위차이 \(d = (0,0,1,1,1,-3)\)의 제곱합이 \(12\)로 크다.
• \((1, 3, 2, 6, 5, 4)\): 인접한 원소들이 조금씩 어긋나 있다. 비부합쌍이 4개로 더 많지만 \(d = (0,-1,1,-2,0,2)\)의 제곱합이 \(10\)으로 작다.

\(r_s\)는 어긋남의 크기를 제곱해서 반영하므로 두 번째 자료를 더 좋게 본다. \(\tau\)는 어긋남의 개수만 세므로 첫 번째 자료를 더 좋게 본다.

실무적 함의: "\(\tau\)와 \(r_s\) 중 무엇을 쓸지"는 취향의 문제가 아니라 무엇을 재고 싶은지의 문제이다. 큰 불일치에 더 큰 벌점을 주고 싶으면 \(r_s\), 불일치의 개수만 세고 싶으면 \(\tau\)를 쓴다. 그리고 두 계수를 모두 계산한 뒤 큰 쪽을 보고하는 것은 명백한 오용이다.

5점 Likert 자료를 만들자.

import numpy as np
from scipy import stats
rng = np.random.default_rng(0)
n = 200
latent = rng.normal(0, 1, n)
x = np.digitize(latent, [-1, -0.3, 0.3, 1]) + 1              # 1~5
y = np.digitize(latent + rng.normal(0, 0.7, n),
                [-1, -0.3, 0.3, 1]) + 1                       # 1~5

C = D = 0
for i in range(n):
    for j in range(i + 1, n):
        s = np.sign(x[i] - x[j]) * np.sign(y[i] - y[j])
        C += s > 0; D += s < 0
n0 = n * (n - 1) / 2
print("tau-a =", (C - D) / n0)
print("tau-b =", stats.kendalltau(x, y, variant='b').statistic)
print("tau-c =", stats.kendalltau(x, y, variant='c').statistic)
print("r_s   =", stats.spearmanr(x, y).statistic)

출력:

tau-a = 0.5389949748743719
tau-b = 0.6724352234242746
tau-c = 0.670375
r_s   = 0.7745703393918136

tau-a가 심하게 축소된다. \(\binom{200}{2} = 19{,}900\)개 쌍 중 약 20%가 \(x\)에서 동점인데, tau-a는 이 동점쌍을 분모에 그대로 두면서 분자에는 기여시키지 않는다. 결과적으로 도달 가능한 최댓값이 1보다 훨씬 작아진다.

실제로 \(x\)와 \(y\)가 완전히 같아도 tau-a가 1이 되지 않는다.

C = D = 0
for i in range(n):
    for j in range(i + 1, n):
        s = np.sign(x[i] - x[j]) ** 2
        C += s > 0
print("완전 일치 시 tau-a =", C / n0)   # 0.8003  ← 1이 아니다

출력:

완전 일치 시 tau-a = 0.800251256281407

tau-b는 분모에서 동점쌍을 제거하여(\(\sqrt{(n_0-n_1)(n_0-n_2)}\)) 이 문제를 고친다. 완전 일치일 때 정확히 1이 된다.

tau-c는 표가 정사각형이 아닐 때(행 범주 수 \(\ne\) 열 범주 수) 쓰는 변형으로, 여기서는 둘 다 5범주라 tau-b와 비슷하다.

권고: 동점이 있으면 언제나 tau-b를 쓴다. scipy.stats.kendalltau의 기본값이 variant='b'이므로 별도 지정이 필요 없다. tau-a는 동점이 전혀 없는 순위 자료에서만 의미가 있다.

import numpy as np, itertools
from scipy import stats

n = 6
R = np.arange(1, n + 1)
S_vals = []
for p in itertools.permutations(R):
    tau = stats.kendalltau(R, p).statistic
    S_vals.append(tau * n * (n - 1) / 2)
S_vals = np.array(S_vals)
print(len(S_vals), S_vals.mean().round(10), S_vals.var().round(4))
print(n * (n - 1) * (2 * n + 5) / 18)

출력:

720 0.0 28.3333
28.333333333333332

출력: \(720\)개 순열, 평균 \(0.0\), 분산 \(28.3333\). 공식값 \(6 \times 5 \times 17 / 18 = 28.3333\)과 정확히 일치한다.

유도의 개요. \(S = \sum_{i<j} a_{ij}\)로 쓰자. 여기서 \(a_{ij} = \operatorname{sign}(R_i - R_j)\operatorname{sign}(S_i - S_j) \in \{-1, +1\}\)이다.

\(H_0\) 아래에서 각 \(a_{ij}\)는 \(\pm 1\)을 확률 \(1/2\)씩 가지므로 \(E[a_{ij}] = 0\), \(\text{Var}(a_{ij}) = 1\)이다. 따라서

이다. 첫 항은 \(\binom{n}{2} = n(n-1)/2\)이다.

공분산 항이 핵심이다. 두 쌍이 관측값을 공유하지 않으면 독립이라 공분산이 0이다. 하나를 공유하면(예: \(\{i,j\}\)와 \(\{i,k\}\)) 공분산이 \(1/3\)이 된다. 세 관측값의 순서 배열 6가지 중 \(a_{ij}\)와 \(a_{ik}\)의 부호가 같은 경우가 4가지이므로 \(E[a_{ij}a_{ik}] = (4 - 2)/6 = 1/3\)이다.

한 관측값을 공유하는 순서쌍의 개수는 \(n(n-1)(n-2)\)이므로

이 되어야 할 것 같지만, 실제 공식은 \(\frac{n(n-1)(2n+5)}{18}\)이다. 차이가 나는 것은 위 계산이 \(R_i\)를 고정하고 \(S_i\)만 순열하는 조건부 구조를 제대로 반영하지 못했기 때문이다. 정확한 유도는 \(E[a_{ij}a_{ik}]\)를 조건부 순열 아래에서 다시 계산해야 하며, 그 결과 \(\frac{n(n-1)(2n+5)}{18}\)이 나온다.

\(n = 6\)에서 두 값을 비교하면 \(\frac{6 \cdot 5 \cdot 11}{6} = 55\)와 \(\frac{6 \cdot 5 \cdot 17}{18} = 28.33\)이다. 열거 결과 \(28.33\)이 맞으므로 후자가 옳다.

실무적 요점: 공식을 외우기보다 \(\text{Var}(S)\)가 \(n^3\) 차수로 커지고 \(\tau = S/\binom{n}{2}\)의 분산이 \(\text{Var}(\tau) \approx \frac{2(2n+5)}{9n(n-1)} = O(1/n)\)임을 기억하는 편이 유용하다. 표준오차가 \(\sqrt{2(2n+5)/(9n(n-1))}\)이므로 \(n = 100\)이면 약 \(0.067\)이다.