Spearman 순위상관¶

import numpy as np
from scipy import stats

# y 는 x 의 지수함수다. 곡선이지만 x 가 커지면 y 도 반드시 커지는 단조 관계다.
x = np.array([1, 2, 3, 4, 5, 6, 7, 8])
y = np.exp(x)

# 스피어만은 값 대신 순위만 보므로 단조이기만 하면 정확히 1 이 된다.
r_s, p_value = stats.spearmanr(x, y)
print(f"Spearman r_s = {r_s:.4f}, p-value = {p_value:.6f}")

# 피어슨은 직선에서 얼마나 벗어났는지를 재므로 1 에 못 미친다.
# 관계가 곡선일 때 두 측도가 갈리는 전형적인 모습이다.
r_p, p_p = stats.pearsonr(x, y)
print(f"Pearson  r   = {r_p:.4f}, p-value = {p_p:.6f}")

# (1) 로그를 씌우면 피어슨만 움직인다
print(f"\nlog y 로 바꾸면  Pearson {stats.pearsonr(x, np.log(y)).statistic:.6f}"
      f"   Spearman {stats.spearmanr(x, np.log(y)).statistic:.6f}")

# (2) 마지막 점이 차지하는 몫
Sxx = ((x - x.mean())**2).sum()
Syy = ((y - y.mean())**2).sum()
Sxy = ((x - x.mean()) * (y - y.mean())).sum()
print(f"\nybar = {y.mean():.3f}   e^8 = {y[-1]:.3f}   배수 {y[-1] / y.mean():.2f}")
print(f"마지막 점의 Syy 몫 {(y[-1] - y.mean())**2 / Syy:.4f}"
      f"   Sxy 몫 {(x[-1] - x.mean()) * (y[-1] - y.mean()) / Sxy:.4f}"
      f"   Sxx 몫 {(x[-1] - x.mean())**2 / Sxx:.4f}")
print(f"Sxx = {Sxx:.0f}   Syy = {Syy:.1f}   Sxy = {Sxy:.3f}"
      f"   r = {Sxy / np.sqrt(Sxx * Syy):.6f}")

# (3) p-값이 정말 0 인가. n = 8 의 순열을 전부 센다.
from itertools import permutations
from math import factorial

ranks = np.arange(1, 9)
exact = sum(1 for q in permutations(ranks)
            if abs(stats.spearmanr(ranks, q).statistic) >= 1 - 1e-12)
print(f"\nr_s = ±1 이 되는 배열 {exact} / {factorial(8)}"
      f"   정확 p = {exact / factorial(8):.6e}")
print(f"scipy 가 돌려준 p = {p_value}")

Spearman r_s = 1.0000, p-value = 0.000000
Pearson  r   = 0.7758, p-value = 0.023636

log y 로 바꾸면  Pearson 1.000000   Spearman 1.000000

ybar = 589.278   e^8 = 2980.958   배수 5.06
마지막 점의 Syy 몫 0.7628   Sxy 몫 0.6080   Sxx 몫 0.2917
Sxx = 42   Syy = 7498953.8   Sxy = 13768.872   r = 0.775842

r_s = ±1 이 되는 배열 2 / 40320   정확 p = 4.960317e-05
scipy 가 돌려준 p = 0.0

import warnings
warnings.filterwarnings("ignore")

import numpy as np
from scipy import stats

rng = np.random.default_rng(21001)
print(f"{'ρ (피어슨)':>11s} {'스피어만 이론':>13s} {'모의':>8s} "
      f"{'켄들 이론':>10s} {'모의':>8s}")
for rho in [0.0, 0.2, 0.4, 0.6, 0.8, 0.95]:
    z = rng.standard_normal((200_000, 2))
    x = z[:, 0]
    y = rho * z[:, 0] + np.sqrt(1 - rho**2) * z[:, 1]
    sp = stats.spearmanr(x, y).statistic
    idx = rng.choice(200_000, 4_000, replace=False)   # 켄들은 계산이 무겁다
    kt = stats.kendalltau(x[idx], y[idx]).statistic
    print(f"{rho:11.2f} {6 / np.pi * np.arcsin(rho / 2):13.4f} {sp:8.4f} "
          f"{2 / np.pi * np.arcsin(rho):10.4f} {kt:8.4f}")

    ρ (피어슨)       스피어만 이론       모의      켄들 이론       모의
       0.00        0.0000   0.0006     0.0000   0.0058
       0.20        0.1913   0.1926     0.1282   0.1054
       0.40        0.3846   0.3839     0.2620   0.2579
       0.60        0.5819   0.5831     0.4097   0.4139
       0.80        0.7859   0.7850     0.5903   0.5919
       0.95        0.9453   0.9453     0.7978   0.8050

import warnings
warnings.filterwarnings("ignore")

import numpy as np
from scipy import stats

rng = np.random.default_rng(20004)
B = 6_000
print("n=30, 명목 0.05")
print(f"{'관계':>20s} {'피어슨':>8s} {'스피어만':>9s} {'켄들':>8s}")
for lab in ["선형 ρ=0.4", "선형 ρ=0 (귀무)", "단조 비선형 (지수)",
            "이상점 오염 5%", "두꺼운 꼬리 t(3)"]:
    a = b = c = 0
    for _ in range(B):
        z = rng.standard_normal((30, 2))
        if lab == "선형 ρ=0.4":
            x, y = z[:, 0], 0.4 * z[:, 0] + np.sqrt(0.84) * z[:, 1]
        elif lab == "선형 ρ=0 (귀무)":
            x, y = z[:, 0], z[:, 1]
        elif lab == "단조 비선형 (지수)":
            x, y = z[:, 0], np.exp(0.6 * z[:, 0] + 0.8 * z[:, 1])
        elif lab == "이상점 오염 5%":
            x = z[:, 0].copy()
            y = 0.4 * z[:, 0] + np.sqrt(0.84) * z[:, 1]
            m = rng.random(30) < 0.05
            x[m] += rng.normal(0, 6, m.sum())
            y[m] += rng.normal(0, 6, m.sum())
        else:                                   # 두꺼운 꼬리
            s = np.sqrt(rng.chisquare(3, 30) / 3)
            x = z[:, 0] / s
            y = (0.4 * z[:, 0] + np.sqrt(0.84) * z[:, 1]) / s
        a += stats.pearsonr(x, y).pvalue < 0.05
        b += stats.spearmanr(x, y).pvalue < 0.05
        c += stats.kendalltau(x, y).pvalue < 0.05
    print(f"{lab:>20s} {a / B:8.4f} {b / B:9.4f} {c / B:8.4f}")

n=30, 명목 0.05
                  관계      피어슨      스피어만       켄들
            선형 ρ=0.4   0.6065    0.5435   0.5397
         선형 ρ=0 (귀무)   0.0523    0.0512   0.0485
         단조 비선형 (지수)   0.8657    0.9263   0.9232
           이상점 오염 5%   0.5350    0.4507   0.4677
         두꺼운 꼬리 t(3)   0.5993    0.5162   0.5437

import warnings
warnings.filterwarnings("ignore")

import numpy as np
from scipy import stats

rng = np.random.default_rng(21002)
z = rng.standard_normal((3_000, 2))
xc = z[:, 0]
yc = 0.6 * z[:, 0] + np.sqrt(1 - 0.36) * z[:, 1]

print("참 상관 0.6 인 연속 자료를 k 단계로 이산화")
print(f"{'척도 수준 수':>12s} {'스피어만':>9s} {'켄들 τ-b':>10s} {'피어슨':>8s}")
for k in [2, 3, 5, 7, 10, 100]:
    qs = np.quantile(xc, np.linspace(0, 1, k + 1)[1:-1])
    xd = np.digitize(xc, qs).astype(float)
    qs = np.quantile(yc, np.linspace(0, 1, k + 1)[1:-1])
    yd = np.digitize(yc, qs).astype(float)
    print(f"{k:12d} {stats.spearmanr(xd, yd).statistic:9.4f} "
          f"{stats.kendalltau(xd, yd, variant='b').statistic:10.4f} "
          f"{np.corrcoef(xd, yd)[0, 1]:8.4f}")
print(f"{'연속 (원본)':>12s} {stats.spearmanr(xc, yc).statistic:9.4f} "
      f"{stats.kendalltau(xc, yc).statistic:10.4f} "
      f"{np.corrcoef(xc, yc)[0, 1]:8.4f}")

참 상관 0.6 인 연속 자료를 k 단계로 이산화
     척도 수준 수      스피어만     켄들 τ-b      피어슨
           2    0.4053     0.4053   0.4053
           3    0.4935     0.4432   0.4935
           5    0.5470     0.4513   0.5470
           7    0.5638     0.4472   0.5638
          10    0.5724     0.4398   0.5724
         100    0.5848     0.4167   0.5848
     연속 (원본)    0.5849     0.4127   0.6112

import warnings
warnings.filterwarnings("ignore")

import itertools
import numpy as np
from scipy import stats

for n in [5, 6, 7]:
    xs = np.arange(n, dtype=float)
    sp, kt = [], []
    for p in itertools.permutations(range(n)):
        y = np.array(p, dtype=float)
        sp.append(stats.spearmanr(xs, y).statistic)
        kt.append(stats.kendalltau(xs, y).statistic)
    sp, kt = np.array(sp), np.array(kt)
    print(f"  n={n}: 순열 {len(sp):6,d}개,  "
          f"P(r_s ≥ 0.9) = {np.mean(sp >= 0.9):.5f},  "
          f"P(τ ≥ 0.6) = {np.mean(kt >= 0.6):.5f}")
    print(f"        r_s 의 SD = {sp.std():.5f} "
          f"(이론 1/√(n-1) = {1 / np.sqrt(n - 1):.5f}),  "
          f"τ 의 SD = {kt.std():.5f} (이론 = "
          f"{np.sqrt(2 * (2 * n + 5) / (9 * n * (n - 1))):.5f})")

print("\n작은 n 에서 scipy 의 p 와 t 근사")
rng = np.random.default_rng(7)
for n in [6, 8, 10, 15]:
    x = np.arange(n, dtype=float)
    y = rng.permutation(n).astype(float)
    ex = stats.spearmanr(x, y)
    r = ex.statistic
    t = r * np.sqrt((n - 2) / (1 - r**2)) if abs(r) < 1 else np.inf
    print(f"  n={n:3d}: r_s = {r:+.4f},  scipy p = {ex.pvalue:.5f},  "
          f"t 근사 p = {2 * stats.t.sf(abs(t), n - 2):.5f}")

  n=5: 순열    120개,  P(r_s ≥ 0.9) = 0.00833,  P(τ ≥ 0.6) = 0.11667
        r_s 의 SD = 0.50000 (이론 1/√(n-1) = 0.50000),  τ 의 SD = 0.40825 (이론 = 0.40825)
  n=6: 순열    720개,  P(r_s ≥ 0.9) = 0.00833,  P(τ ≥ 0.6) = 0.06806
        r_s 의 SD = 0.44721 (이론 1/√(n-1) = 0.44721),  τ 의 SD = 0.35486 (이론 = 0.35486)
  n=7: 순열  5,040개,  P(r_s ≥ 0.9) = 0.00337,  P(τ ≥ 0.6) = 0.03452
        r_s 의 SD = 0.40825 (이론 1/√(n-1) = 0.40825),  τ 의 SD = 0.31706 (이론 = 0.31706)

작은 n 에서 scipy 의 p 와 t 근사
  n=  6: r_s = -0.2571,  scipy p = 0.62279,  t 근사 p = 0.62279
  n=  8: r_s = -0.0476,  scipy p = 0.91085,  t 근사 p = 0.91085
  n= 10: r_s = -0.0061,  scipy p = 0.98674,  t 근사 p = 0.98674
  n= 15: r_s = +0.3607,  scipy p = 0.18655,  t 근사 p = 0.18655

scipy.stats.permutation_test 로 직접 순열검정을 하거나
n ≤ 10 이면 표를 참조한다

import warnings
warnings.filterwarnings("ignore")

import numpy as np
from scipy import stats

rng = np.random.default_rng(21005)
B = 4_000
print("참 ρ_s, τ 를 덮는 비율 (명목 95%)")
print(f"{'ρ':>6s} {'n':>5s} {'스피어만':>10s} {'켄들':>9s}")
for rho in [0.0, 0.3, 0.6, 0.85]:
    for n in [15, 30, 60]:
        sp_th = 6 / np.pi * np.arcsin(rho / 2)
        kt_th = 2 / np.pi * np.arcsin(rho)
        a = b = 0
        for _ in range(B):
            z = rng.standard_normal((n, 2))
            x = z[:, 0]
            y = rho * z[:, 0] + np.sqrt(1 - rho**2) * z[:, 1]
            rs = stats.spearmanr(x, y).statistic
            se = np.sqrt((1 + rs**2 / 2) / (n - 3))
            lo, hi = np.arctanh(rs) - 1.96 * se, np.arctanh(rs) + 1.96 * se
            a += np.tanh(lo) <= sp_th <= np.tanh(hi)
            tk = stats.kendalltau(x, y).statistic
            se2 = np.sqrt(0.437 / (n - 4))
            lo, hi = np.arctanh(tk) - 1.96 * se2, np.arctanh(tk) + 1.96 * se2
            b += np.tanh(lo) <= kt_th <= np.tanh(hi)
        print(f"{rho:6.2f} {n:5d} {a / B:10.4f} {b / B:9.4f}")

참 ρ_s, τ 를 덮는 비율 (명목 95%)
     ρ     n       스피어만        켄들
  0.00    15     0.9590    0.9475
  0.00    30     0.9565    0.9470
  0.00    60     0.9507    0.9453
  0.30    15     0.9620    0.9547
  0.30    30     0.9543    0.9465
  0.30    60     0.9537    0.9453
  0.60    15     0.9573    0.9483
  0.60    30     0.9595    0.9523
  0.60    60     0.9525    0.9463
  0.85    15     0.9550    0.9323
  0.85    30     0.9507    0.9485
  0.85    60     0.9583    0.9543

자료 쌍 (x_i, y_i) 를 쌍째로 재표집한다
각 부트스트랩 표본에서 r_s 를 계산한다
2.5% 와 97.5% 백분위수를 구간으로 쓴다

stats.spearmanr(x, y)                 → statistic, pvalue
stats.spearmanr(X)                    → 행렬 전체의 상관행렬
stats.spearmanr(x, y, nan_policy=..)  → 결측 처리

주의: 기본 p 는 t 근사이며 정확 검정이 아니다

직무 만족도(7점 척도)와 이직 의사(5점 척도)의 관계

  스피어만 r_s = -0.42,  95% CI [-0.55, -0.27],  n = 180,  p < 0.001
  (보너트-라이트 구간)
  켄들 τ_b = -0.32 로 방향과 크기가 일치한다
  두 변수 모두 서열 척도이므로 순위상관이 적절하다