Spearman의 rho 검정¶

import numpy as np
from scipy import stats

x = np.array([3, 1, 6, 4, 8, 2, 7, 5, 10, 9, 11, 12])
y = np.array([5, 2, 8, 3, 11, 1, 9, 6, 12, 7, 10, 4])

# 순위만 쓰므로 정규성이 필요 없고 이상치에도 강하다.
r_s, p_value = stats.spearmanr(x, y)
print(f"Spearman r_s = {r_s:.4f}")
print(f"p-value      = {p_value:.4f}")

# 피어슨과 같은 꼴의 t 통계량으로 손계산해 확인한다. 순위에 피어슨을
# 적용한 것이 스피어만이므로 검정 방식도 그대로 따라온다.
n = len(x)
t_stat = r_s * np.sqrt((n - 2) / (1 - r_s**2))
p_manual = 2 * (1 - stats.t.cdf(abs(t_stat), df=n - 2))
print(f"t-statistic  = {t_stat:.4f}")
print(f"Manual p     = {p_manual:.4f}")

# (1) 동점이 없으므로 간편 공식이 정확하다.
rx, ry = stats.rankdata(x), stats.rankdata(y)
print(f"\nx 가 1..{n} 의 치환인가: {sorted(x) == list(range(1, n + 1))},  "
      f"y 는: {sorted(y) == list(range(1, n + 1))}")
d2 = np.sum((rx - ry) ** 2)
print(f"sum d^2 = {d2:.0f},  n(n^2-1) = {n * (n**2 - 1)}")
print(f"간편 공식 1 - 6*{d2:.0f}/{n * (n**2 - 1)} = "
      f"{1 - 6 * d2 / (n * (n**2 - 1)):.6f}")
print(f"spearmanr                     = {r_s:.6f}")

# (2) 귀무분포를 순열로 직접 만든다. H0 에서 두 순위열은 서로 독립인 치환이다.
rng = np.random.default_rng(12)
B = 200_000
base = np.arange(1, n + 1)
perms = rng.permuted(np.tile(base, (B, 1)), axis=1)
vals = 1 - 6 * ((base - perms) ** 2).sum(axis=1) / (n * (n**2 - 1))
hit = (np.abs(vals) >= abs(r_s) - 1e-12).sum()
p_perm = (hit + 1) / (B + 1)
print(f"\n귀무분포: 평균 {vals.mean():+.5f}  표준편차 {vals.std(ddof=1):.5f}   "
      f"이론 0 과 1/sqrt(n-1) = {1 / np.sqrt(n - 1):.5f}")
print(f"순열 p   = {p_perm:.5f}  (B = {B:,}, 몬테카를로 SE = "
      f"{np.sqrt(p_perm * (1 - p_perm) / B):.5f})")
print(f"t 근사 p = {2 * stats.t.sf(abs(t_stat), n - 2):.5f}")
print(f"정규 근사 p = {2 * stats.norm.sf(abs(r_s) * np.sqrt(n - 1)):.5f}   "
      f"(z = r_s sqrt(n-1) = {r_s * np.sqrt(n - 1):.4f})")

Spearman r_s = 0.6573
p-value      = 0.0202
t-statistic  = 2.7584
Manual p     = 0.0202

x 가 1..12 의 치환인가: True,  y 는: True
sum d^2 = 98,  n(n^2-1) = 1716
간편 공식 1 - 6*98/1716 = 0.657343
spearmanr                     = 0.657343

귀무분포: 평균 -0.00028  표준편차 0.30186   이론 0 과 1/sqrt(n-1) = 0.30151
순열 p   = 0.02377  (B = 200,000, 몬테카를로 SE = 0.00034)
t 근사 p = 0.02019
정규 근사 p = 0.02925   (z = r_s sqrt(n-1) = 2.1802)

import numpy as np
from scipy import stats

rng = np.random.default_rng(0)
n = 30
rp, rs, rk = [], [], []
sp = ss = sk = 0
for _ in range(3000):
    X = rng.normal(0, 1, n)
    Y = 0.4*X + rng.normal(0, 1, n)
    X[0] += 10; Y[0] -= 10          # 이상치 하나
    a = stats.pearsonr(X, Y)
    b = stats.spearmanr(X, Y)
    c = stats.kendalltau(X, Y)
    rp.append(a.statistic); rs.append(b.statistic); rk.append(c.statistic)
    sp += (a.pvalue < 0.05) and a.statistic > 0
    ss += (b.pvalue < 0.05) and b.statistic > 0
    sk += (c.pvalue < 0.05) and c.statistic > 0

print(f"  평균 r: Pearson {np.mean(rp):+.4f}  "
      f"Spearman {np.mean(rs):+.4f}  Kendall {np.mean(rk):+.4f}")
print(f"  '양의 상관'을 옳게 잡은 비율:")
print(f"    Pearson {sp/3000:.4f}  Spearman {ss/3000:.4f}  Kendall {sk/3000:.4f}")

  평균 r: Pearson -0.6716  Spearman +0.2123  Kendall +0.1575
  '양의 상관'을 옳게 잡은 비율:
    Pearson 0.0000  Spearman 0.1640  Kendall 0.2033

import numpy as np
from scipy import stats

rng = np.random.default_rng(3)
x = rng.normal(0, 1, 12); y = 0.6*x + rng.normal(0, 1, 12)
R, S = stats.rankdata(x), stats.rankdata(y)
n = len(x); d = R - S
print(f"  간편식      {1 - 6*np.sum(d**2)/(n*(n**2-1)):.6f}")
print(f"  순위의 피어슨 {stats.pearsonr(R, S).statistic:.6f}")
print(f"  scipy       {stats.spearmanr(x, y).statistic:.6f}")

# 동점을 만들어 본다
xt = np.round(x); yt = np.round(y)
Rt, St = stats.rankdata(xt), stats.rankdata(yt)
dt = Rt - St
print(f"\n  동점 있음:")
print(f"  간편식      {1 - 6*np.sum(dt**2)/(n*(n**2-1)):.6f}")
print(f"  순위의 피어슨 {stats.pearsonr(Rt, St).statistic:.6f}")

  간편식      0.545455
  순위의 피어슨 0.545455
  scipy       0.545455

  동점 있음:
  간편식      0.692308
  순위의 피어슨 0.655729

import numpy as np
from scipy import stats

rng = np.random.default_rng(0)
REP, n = 3000, 30
print(f"{'모집단':>10}{'Pearson':>10}{'Spearman':>10}")
for name, gen in [("정규",  lambda m: rng.normal(0, 1, m)),
                  ("t(3)",  lambda m: rng.standard_t(3, m)),
                  ("지수",  lambda m: rng.exponential(1, m) - 1)]:
    rp = rs = 0
    for _ in range(REP):
        z, e = gen(n), gen(n)
        X, Y = z, 0.35*z + e
        rp += stats.pearsonr(X, Y).pvalue < 0.05
        rs += stats.spearmanr(X, Y).pvalue < 0.05
    print(f"{name:>10}{rp/REP:>10.4f}{rs/REP:>10.4f}")

       모집단   Pearson  Spearman
        정규    0.4390    0.3920
      t(3)    0.4800    0.4720
        지수    0.4503    0.5013

import numpy as np
from scipy import stats

rng = np.random.default_rng(2)
n, rho = 40, 0.6
REP = 20_000
cov_naive = cov_adj = 0
rho_s_true = 6/np.pi*np.arcsin(rho/2)

for _ in range(REP):
    X = rng.multivariate_normal([0,0], [[1,rho],[rho,1]], n)
    rs = stats.spearmanr(X[:,0], X[:,1]).statistic
    z = np.arctanh(rs)
    for se, box in ((1/np.sqrt(n-3), 'naive'),
                    (np.sqrt((1 + rs**2/2)/(n-3)), 'adj')):
        lo, hi = np.tanh(z - 1.96*se), np.tanh(z + 1.96*se)
        if lo <= rho_s_true <= hi:
            if box == 'naive': cov_naive += 1
            else: cov_adj += 1

print(f"  참 rho_s = {rho_s_true:.4f}")
print(f"  보정 없는 피셔 z : 포함률 {cov_naive/REP:.4f}")
print(f"  보노-웨어 보정   : 포함률 {cov_adj/REP:.4f}")

  참 rho_s = 0.5819
  보정 없는 피셔 z : 포함률 0.9395
  보노-웨어 보정   : 포함률 0.9573

import numpy as np
from scipy import stats

rng = np.random.default_rng(5)
n = 300
x = rng.uniform(-3, 3, n)
cases = {
    "단조 증가":   0.8*x + rng.normal(0, 0.5, n),
    "U자":        x**2 + rng.normal(0, 0.5, n),
    "주기(sin)":  np.sin(2*x) + rng.normal(0, 0.2, n),
    "단조 포화":   np.tanh(2*x) + rng.normal(0, 0.2, n),
}
print(f"{'관계':>12}{'Pearson':>10}{'Spearman':>11}")
for k, y in cases.items():
    print(f"{k:>12}{stats.pearsonr(x,y).statistic:>10.4f}"
          f"{stats.spearmanr(x,y).statistic:>11.4f}")

        관계   Pearson   Spearman
    단조 증가    0.9412     0.9485
        U자   -0.0057    -0.0249
  주기(sin)   -0.4010    -0.3729
    단조 포화    0.9037     0.8733

import numpy as np
from scipy import stats

def vdw(x, y):
    n = len(x)
    a = stats.norm.ppf(stats.rankdata(x) / (n + 1))
    b = stats.norm.ppf(stats.rankdata(y) / (n + 1))
    return stats.pearsonr(a, b)

rng = np.random.default_rng(7)
REP, n = 3000, 30
print(f"{'모집단':>10}{'Pearson':>10}{'Spearman':>10}{'정규점수':>10}")
for name, gen in [("정규", lambda m: rng.normal(0,1,m)),
                  ("t(3)", lambda m: rng.standard_t(3,m)),
                  ("지수", lambda m: rng.exponential(1,m)-1)]:
    rp = rs = rv = 0
    for _ in range(REP):
        z, e = gen(n), gen(n)
        X, Y = z, 0.35*z + e
        rp += stats.pearsonr(X,Y).pvalue < 0.05
        rs += stats.spearmanr(X,Y).pvalue < 0.05
        rv += vdw(X,Y).pvalue < 0.05
    print(f"{name:>10}{rp/REP:>10.4f}{rs/REP:>10.4f}{rv/REP:>10.4f}")

       모집단   Pearson  Spearman      정규점수
        정규    0.4390    0.3920    0.4297
      t(3)    0.4800    0.4720    0.5030
        지수    0.4503    0.5013    0.4733