두 상관의 비교¶

import numpy as np
from scipy import stats

# 서로 다른 두 표본에서 얻은 상관계수를 견준다. 표본이 겹치지 않아야
# 이 방법을 쓸 수 있다. 같은 사람에게서 잰 두 상관이라면 다른 검정이 필요하다.
r1, n1 = 0.65, 50
r2, n2 = 0.40, 60

# 두 r 을 각각 z 로 옮긴다. z 는 거의 정규이고 분산이 표본크기로만 정해진다.
z1 = np.arctanh(r1)
z2 = np.arctanh(r2)
# 두 z 가 독립이므로 차이의 분산은 각 분산의 합이다.
se = np.sqrt(1 / (n1 - 3) + 1 / (n2 - 3))
Z_stat = (z1 - z2) / se
p_value = 2 * (1 - stats.norm.cdf(abs(Z_stat)))

print(f"z1 = {z1:.4f}, z2 = {z2:.4f}")
print(f"Z statistic = {Z_stat:.4f}")
print(f"Two-sided p-value = {p_value:.4f}")

# (1) zeta1 - zeta2 의 95% 구간과, 그것을 tanh 로 되돌리면 안 되는 까닭
lo, hi = (z1 - z2) - 1.96 * se, (z1 - z2) + 1.96 * se
print(f"\nSE = sqrt(1/{n1-3} + 1/{n2-3}) = {se:.6f}")
print(f"zeta1 - zeta2 의 95% 구간 = ({lo:+.4f}, {hi:+.4f})   0 을 품는가: {lo < 0 < hi}")
print(f"끝점에 tanh 를 씌우면      = ({np.tanh(lo):+.4f}, {np.tanh(hi):+.4f})")
print(f"그러나 rho1 - rho2 의 점추정은 {r1 - r2:+.4f} 이고 "
      f"tanh(z1-z2) = {np.tanh(z1 - z2):+.4f} 로 서로 다르다")

# (2) 검정력 0.80 에 필요한 집단당 표본크기
zeta1, zeta2 = np.arctanh(0.65), np.arctanh(0.40)
delta = zeta1 - zeta2
need = stats.norm.ppf(0.975) + stats.norm.ppf(0.80)
n_req = 3 + 2 * (need / delta) ** 2
print(f"\ndelta = zeta1 - zeta2 = {delta:.6f}")
print(f"z_0.025 + z_0.20 = {need:.6f}")
print(f"필요한 n = 3 + 2*(need/delta)^2 = {n_req:.2f}  ->  "
      f"집단당 {int(np.ceil(n_req))} 명")

n = int(np.ceil(n_req))
ncp = delta / np.sqrt(2 / (n - 3))
print(f"이론 검정력 = {stats.norm.sf(1.96 - ncp) + stats.norm.cdf(-1.96 - ncp):.4f}")

rng = np.random.default_rng(7)
B = 20_000
def sample_r(rho, m, B):
    a = rng.standard_normal((B, m))
    b = rho * a + np.sqrt(1 - rho**2) * rng.standard_normal((B, m))
    a = a - a.mean(axis=1, keepdims=True)
    b = b - b.mean(axis=1, keepdims=True)
    return (a * b).sum(axis=1) / np.sqrt((a**2).sum(axis=1) * (b**2).sum(axis=1))

R1, R2 = sample_r(0.65, n, B), sample_r(0.40, n, B)
Zs = (np.arctanh(R1) - np.arctanh(R2)) / np.sqrt(2 / (n - 3))
print(f"모의 검정력 ({B:,}회, n = {n}) = {np.mean(np.abs(Zs) > 1.96):.4f}  "
      f"+- {np.std(np.abs(Zs) > 1.96) / np.sqrt(B):.4f}")

# 1종 오류도 함께 본다 (rho1 = rho2 = 0.65)
R1, R2 = sample_r(0.65, n, B), sample_r(0.65, n, B)
Zs = (np.arctanh(R1) - np.arctanh(R2)) / np.sqrt(2 / (n - 3))
print(f"모의 1종 오류 (rho1 = rho2 = 0.65) = {np.mean(np.abs(Zs) > 1.96):.4f}  (명목 0.05)")

z1 = 0.7753, z2 = 0.4236
Z statistic = 1.7848
Two-sided p-value = 0.0743

SE = sqrt(1/47 + 1/57) = 0.197029
zeta1 - zeta2 의 95% 구간 = (-0.0345, +0.7378)   0 을 품는가: True
끝점에 tanh 를 씌우면      = (-0.0345, +0.6278)
그러나 rho1 - rho2 의 점추정은 +0.2500 이고 tanh(z1-z2) = +0.3378 로 서로 다르다

delta = zeta1 - zeta2 = 0.351650
z_0.025 + z_0.20 = 2.801585
필요한 n = 3 + 2*(need/delta)^2 = 129.95  ->  집단당 130 명
이론 검정력 = 0.8002
모의 검정력 (20,000회, n = 130) = 0.8009  +- 0.0028
모의 1종 오류 (rho1 = rho2 = 0.65) = 0.0498  (명목 0.05)

import warnings
warnings.filterwarnings("ignore")

import numpy as np
from scipy import stats

def williams_t(r12, r13, r23, n):
    """겹치는 종속 상관의 비교 (공통 변수 1)."""
    detR = 1 - r12**2 - r13**2 - r23**2 + 2 * r12 * r13 * r23
    rbar = (r12 + r13) / 2
    t = (r12 - r13) * np.sqrt(
        (n - 1) * (1 + r23)
        / (2 * detR * (n - 1) / (n - 3) + rbar**2 * (1 - r23)**3))
    return t, 2 * stats.t.sf(abs(t), n - 3)

def independent_z(r1, n1, r2, n2):
    z = (np.arctanh(r1) - np.arctanh(r2)) / np.sqrt(1 / (n1 - 3) + 1 / (n2 - 3))
    return z, 2 * stats.norm.sf(abs(z))

rng = np.random.default_rng(29001)
B = 4_000
print("H0: ρ12 = ρ13 이 참인 자료에서의 1종 오류 (명목 0.05)")
print(f"{'n':>5s} {'r23':>6s} {'독립표본 z 검정':>15s} {'윌리엄스 t':>12s}")
for n in [50, 100]:
    for r23 in [0.0, 0.5, 0.9]:
        R = np.array([[1, 0.5, 0.5], [0.5, 1, r23], [0.5, r23, 1]])
        L = np.linalg.cholesky(R)
        a = b = 0
        for _ in range(B):
            d = rng.standard_normal((n, 3)) @ L.T
            r12 = np.corrcoef(d[:, 0], d[:, 1])[0, 1]
            r13 = np.corrcoef(d[:, 0], d[:, 2])[0, 1]
            r23h = np.corrcoef(d[:, 1], d[:, 2])[0, 1]
            a += independent_z(r12, n, r13, n)[1] < 0.05
            b += williams_t(r12, r13, r23h, n)[1] < 0.05
        print(f"{n:5d} {r23:6.1f} {a / B:15.4f} {b / B:12.4f}")

H0: ρ12 = ρ13 이 참인 자료에서의 1종 오류 (명목 0.05)
    n    r23       독립표본 z 검정       윌리엄스 t
   50    0.0          0.0673       0.0473
   50    0.5          0.0092       0.0460
   50    0.9          0.0000       0.0505
  100    0.0          0.0602       0.0490
  100    0.5          0.0152       0.0475
  100    0.9          0.0000       0.0500

import pingouin as pg
pg.corr_test?  →  다음이 표준적이다

# 겹치는 종속
pg.rm_corr(...)  또는 직접 윌리엄스 t 구현

# cocor (R) 이 가장 완전하다: 10가지 이상의 검정을 제공

import numpy as np
from scipy import stats

za, zb = stats.norm.ppf(0.975), stats.norm.ppf(0.80)
print("검정력 0.80, 양측 0.05")
print(f"{'ρ1':>6s} {'ρ2':>6s} {'차이':>6s} {'z1-z2':>8s} {'군당 n':>8s}")
for r1, r2 in [(0.5, 0.3), (0.5, 0.4), (0.7, 0.5),
               (0.3, 0.1), (0.9, 0.8), (0.2, 0.0)]:
    d = abs(np.arctanh(r1) - np.arctanh(r2))
    print(f"{r1:6.2f} {r2:6.2f} {r1 - r2:6.2f} {d:8.4f} "
          f"{int(np.ceil(2 * ((za + zb) / d)**2 + 3)):8d}")

검정력 0.80, 양측 0.05
    ρ1     ρ2     차이    z1-z2     군당 n
  0.50   0.30   0.20   0.2398      277
  0.50   0.40   0.10   0.1257      998
  0.70   0.50   0.20   0.3180      159
  0.30   0.10   0.20   0.2092      362
  0.90   0.80   0.10   0.3736      116
  0.20   0.00   0.20   0.2027      385

rng = np.random.default_rng(29002)
B = 5_000
print("\n모의실험으로 검정력 확인")
for r1, r2, n in [(0.5, 0.3, 277), (0.5, 0.4, 998),
                  (0.7, 0.5, 159), (0.9, 0.8, 116)]:
    a = 0
    for _ in range(B):
        z = rng.standard_normal((n, 2))
        R1 = np.corrcoef(z[:, 0], r1 * z[:, 0]
                         + np.sqrt(1 - r1**2) * z[:, 1])[0, 1]
        z = rng.standard_normal((n, 2))
        R2 = np.corrcoef(z[:, 0], r2 * z[:, 0]
                         + np.sqrt(1 - r2**2) * z[:, 1])[0, 1]
        zz = (np.arctanh(R1) - np.arctanh(R2)) / np.sqrt(2 / (n - 3))
        a += 2 * stats.norm.sf(abs(zz)) < 0.05
    print(f"  ρ1={r1}, ρ2={r2}, 군당 n={n}: 검정력 {a / B:.4f}")

모의실험으로 검정력 확인
  ρ1=0.5, ρ2=0.3, 군당 n=277: 검정력 0.7996
  ρ1=0.5, ρ2=0.4, 군당 n=998: 검정력 0.8040
  ρ1=0.7, ρ2=0.5, 군당 n=159: 검정력 0.7948
  ρ1=0.9, ρ2=0.8, 군당 n=116: 검정력 0.8040

import warnings
warnings.filterwarnings("ignore")

import numpy as np
from scipy import stats

def fisher_ci(r, n, a=0.05):
    q = stats.norm.ppf(1 - a / 2) / np.sqrt(n - 3)
    return np.tanh(np.arctanh(r) - q), np.tanh(np.arctanh(r) + q)

def zou_ci(r1, n1, r2, n2, a=0.05):
    """두 독립 상관의 차이에 대한 조우의 구간."""
    l1, u1 = fisher_ci(r1, n1, a)
    l2, u2 = fisher_ci(r2, n2, a)
    d = r1 - r2
    return (d - np.sqrt((r1 - l1)**2 + (u2 - r2)**2),
            d + np.sqrt((u1 - r1)**2 + (r2 - l2)**2))

rng = np.random.default_rng(29002)
B = 20_000
print("피복확률 (명목 95%)")
print(f"{'ρ1':>6s} {'ρ2':>6s} {'n':>6s} {'Zou':>8s} {'z 차이를 tanh':>14s}")
for r1, r2, n in [(0.5, 0.3, 50), (0.5, 0.3, 200),
                  (0.8, 0.3, 50), (0.9, 0.1, 30)]:
    a = b = 0
    true = r1 - r2
    for _ in range(B):
        z = rng.standard_normal((n, 2))
        R1 = np.corrcoef(z[:, 0], r1 * z[:, 0]
                         + np.sqrt(1 - r1**2) * z[:, 1])[0, 1]
        z = rng.standard_normal((n, 2))
        R2 = np.corrcoef(z[:, 0], r2 * z[:, 0]
                         + np.sqrt(1 - r2**2) * z[:, 1])[0, 1]
        lo, hi = zou_ci(R1, n, R2, n)
        a += lo <= true <= hi
        se = np.sqrt(2 / (n - 3))
        d = np.arctanh(R1) - np.arctanh(R2)
        b += np.tanh(d - 1.96 * se) <= true <= np.tanh(d + 1.96 * se)
    print(f"{r1:6.2f} {r2:6.2f} {n:6d} {a / B:8.4f} {b / B:14.4f}")

피복확률 (명목 95%)
    ρ1     ρ2      n      Zou     z 차이를 tanh
  0.50   0.30     50   0.9509         0.9465
  0.50   0.30    200   0.9486         0.9317
  0.80   0.30     50   0.9478         0.7733
  0.90   0.10     30   0.9490         0.8158

남성(n=180)과 여성(n=210)에서 스트레스와 수면의 상관

  남성  r = 0.52  [0.40, 0.62]
  여성  r = 0.31  [0.19, 0.43]

  차이  Δr = 0.21,  95% CI [0.04, 0.37]   (조우의 방법)
  피셔 z 검정: z = 2.36, p = 0.018

구간이 0 을 포함하지 않으나 상한이 0.37 로 넓어,
차이의 크기는 정밀하게 추정되지 않았다.

import warnings
warnings.filterwarnings("ignore")

import numpy as np
from scipy import stats

def homogeneity(rs, ns):
    """k 개 독립 상관의 동질성 검정."""
    z = np.arctanh(np.asarray(rs, float))
    w = np.asarray(ns, float) - 3
    zbar = (w * z).sum() / w.sum()
    Q = (w * (z - zbar)**2).sum()
    return Q, stats.chi2.sf(Q, len(rs) - 1), np.tanh(zbar)

print("예: 네 지역의 상관")
rs = [0.55, 0.48, 0.62, 0.30]
ns = [80, 95, 70, 110]
Q, p, pooled = homogeneity(rs, ns)
print(f"  r = {rs},  n = {ns}")
print(f"  Q = {Q:.4f}, df = {len(rs) - 1}, p = {p:.4f}")
print(f"  통합 추정 r = {pooled:.4f}")

rng = np.random.default_rng(29003)
B = 5_000
print("\n1종 오류와 검정력 (명목 0.05)")
print(f"{'설정':>28s} {'기각률':>8s}")
for lab, rhos in [("모두 0.5 (귀무)", [0.5] * 4),
                  ("하나만 다름 0.5,0.5,0.5,0.3", [0.5, 0.5, 0.5, 0.3]),
                  ("넓게 퍼짐 0.2~0.8", [0.2, 0.4, 0.6, 0.8])]:
    a = 0
    for _ in range(B):
        rs = []
        for rho, n in zip(rhos, [80, 95, 70, 110]):
            z = rng.standard_normal((n, 2))
            rs.append(np.corrcoef(z[:, 0], rho * z[:, 0]
                                  + np.sqrt(1 - rho**2) * z[:, 1])[0, 1])
        a += homogeneity(rs, [80, 95, 70, 110])[1] < 0.05
    print(f"{lab:>28s} {a / B:8.4f}")

예: 네 지역의 상관
  r = [0.55, 0.48, 0.62, 0.3],  n = [80, 95, 70, 110]
  Q = 8.2997, df = 3, p = 0.0402
  통합 추정 r = 0.4756

1종 오류와 검정력 (명목 0.05)
                          설정      기각률
                 모두 0.5 (귀무)   0.0542
      하나만 다름 0.5,0.5,0.5,0.3   0.3810
               넓게 퍼짐 0.2~0.8   1.0000

import warnings
warnings.filterwarnings("ignore")

import numpy as np
from scipy import stats

rng = np.random.default_rng(29004)
n = 100_000
print("두 집단: 기울기는 같고 x 의 퍼짐만 다르다")
print(f"{'집단':>8s} {'SD(x)':>8s} {'기울기':>8s} {'SD(잔차)':>9s} {'r':>8s}")
for lab, sx in [("A", 1.0), ("B", 3.0)]:
    x = rng.normal(0, sx, n)
    y = 2.0 * x + rng.normal(0, 4, n)
    print(f"{lab:>8s} {sx:8.1f} {np.polyfit(x, y, 1)[0]:8.4f} "
          f"{4.0:9.1f} {np.corrcoef(x, y)[0, 1]:8.4f}")

print("\n두 집단: 상관은 같고 기울기만 다르다")
print(f"{'집단':>8s} {'SD(x)':>8s} {'기울기':>8s} {'SD(잔차)':>9s} {'r':>8s}")
for lab, beta, se in [("A", 1.0, 1.333), ("B", 4.0, 5.333)]:
    x = rng.normal(0, 1, n)
    y = beta * x + rng.normal(0, se, n)
    print(f"{lab:>8s} {1.0:8.1f} {np.polyfit(x, y, 1)[0]:8.4f} "
          f"{se:9.3f} {np.corrcoef(x, y)[0, 1]:8.4f}")

두 집단: 기울기는 같고 x 의 퍼짐만 다르다
      집단    SD(x)      기울기    SD(잔차)        r
       A      1.0   2.0022       4.0   0.4473
       B      3.0   1.9966       4.0   0.8312

두 집단: 상관은 같고 기울기만 다르다
      집단    SD(x)      기울기    SD(잔차)        r
       A      1.0   1.0067     1.333   0.6028
       B      1.0   3.9775     5.333   0.5995

y ~ x + 집단 + x:집단

x:집단 계수의 유의성  =  두 기울기가 다른가

두 상관을 비교하려 한다
    │
    ├─ 같은 것을 비교하는 게 맞는가?
    │     └─ σ_x 가 집단마다 크게 다르면 → 기울기를 비교한다
    │
    ├─ 서로 다른 표본인가?
    │     └─ 예 → 피셔 z 검정 + 조우의 차이 구간
    │
    └─ 같은 표본인가?
          ├─ 공통 변수가 있다 (r12 vs r13) → 윌리엄스 t
          └─ 공통 변수가 없다 (r12 vs r34) → 스타이거 z*

□ 두 표본이 독립인가 종속인가 (검정이 달라진다)
□ 각 상관의 n 과 신뢰구간
□ 차이의 신뢰구간 (조우의 방법)
□ 각 집단의 SD(x), SD(y)  ← 범위 차이를 독자가 판단할 수 있게
□ 검정력 또는 검출 가능한 최소 차이
□ 여러 비교를 했다면 보정 방법