상관과 인과¶

import numpy as np
from scipy import stats

np.random.seed(42)
n = 100

# 선형 관계. 세 측도가 모두 크게 나온다.
x_lin = np.random.normal(0, 1, n)
y_lin = 2 * x_lin + np.random.normal(0, 1, n)

# 단조이지만 곡선인 관계. Pearson 은 떨어지지만 순위를 쓰는 두 측도는 버틴다.
x_mono = np.random.uniform(0, 3, n)
y_mono = np.exp(x_mono) + np.random.normal(0, 2, n)

# 이차함수 관계. 관계는 아주 강한데 세 측도 모두 0 근처로 나온다.
# 상관계수가 0 이라는 말이 "관계가 없다"는 뜻이 아님을 보여 주는 자리다.
x_quad = np.random.normal(0, 2, n)
y_quad = x_quad**2 + np.random.normal(0, 1, n)

datasets = [
    ("Linear", x_lin, y_lin),
    ("Monotonic Nonlinear", x_mono, y_mono),
    ("Quadratic (r ~ 0)", x_quad, y_quad),
]

for name, x, y in datasets:
    r_p, _ = stats.pearsonr(x, y)
    r_s, _ = stats.spearmanr(x, y)
    r_k, _ = stats.kendalltau(x, y)
    print(f"{name:<25} r={r_p:.4f}  rho_s={r_s:.4f}  tau={r_k:.4f}")

# 모집단 값과 견준다.
rho_lin = 2 / np.sqrt(5)
print(f"\n선형: 참 rho = 2/sqrt(5) = {rho_lin:.6f}")
print(f"  (x,y) 가 이변량 정규이므로 닫힌 꼴이 있다")
print(f"  rho_s = (6/pi) arcsin(rho/2) = {6 / np.pi * np.arcsin(rho_lin / 2):.6f}")
print(f"  tau   = (2/pi) arcsin(rho)   = {2 / np.pi * np.arcsin(rho_lin):.6f}")

rng = np.random.default_rng(3)
def pop(gen, m=200_000, reps=20):
    rs, ss, ks = [], [], []
    for _ in range(reps):
        a, b = gen(m)
        rs.append(stats.pearsonr(a, b).statistic)
        ss.append(stats.spearmanr(a, b).statistic)
        ks.append(stats.kendalltau(a[:40_000], b[:40_000]).statistic)
    return np.mean(rs), np.mean(ss), np.mean(ks)

gens = [
    ("Linear", lambda m: (z := rng.standard_normal(m), 2 * z + rng.standard_normal(m))),
    ("Monotonic Nonlinear",
     lambda m: (z := rng.uniform(0, 3, m), np.exp(z) + rng.normal(0, 2, m))),
    ("Quadratic (r ~ 0)",
     lambda m: (z := rng.normal(0, 2, m), z**2 + rng.standard_normal(m))),
]
print(f"\n{'관계':<22s} {'모집단 r':>10s} {'rho_s':>10s} {'tau':>10s}")
for name, g in gens:
    a, b, c = pop(g)
    print(f"{name:<22s} {a:+10.5f} {b:+10.5f} {c:+10.5f}")

Linear                    r=0.8724  rho_s=0.8686  tau=0.6853
Monotonic Nonlinear       r=0.8676  rho_s=0.9032  tau=0.7402
Quadratic (r ~ 0)         r=0.1014  rho_s=-0.0229  tau=-0.0376

선형: 참 rho = 2/sqrt(5) = 0.894427
  (x,y) 가 이변량 정규이므로 닫힌 꼴이 있다
  rho_s = (6/pi) arcsin(rho/2) = 0.885502
  tau   = (2/pi) arcsin(rho)   = 0.704833

관계                          모집단 r      rho_s        tau
Linear                   +0.89427   +0.88537   +0.70497
Monotonic Nonlinear      +0.87282   +0.88110   +0.70520
Quadratic (r ~ 0)        +0.00038   +0.00026   -0.00096

def fisher_z_ci(x, y, alpha=0.05):
    """Fisher z 변환으로 상관계수의 신뢰구간을 구한다.

    r 은 -1 과 1 사이에 갇혀 있어 분포가 치우친다. arctanh 를 씌우면
    그 눈금이 실수 전체로 펴지며 분포가 거의 정규가 되고, 표준오차도
    1/sqrt(n-3) 로 간단해진다. 구간을 만든 뒤 tanh 로 되돌린다.
    """
    n = len(x)
    r, p_val = stats.pearsonr(x, y)
    z = np.arctanh(r)
    se = 1 / np.sqrt(n - 3)
    z_crit = stats.norm.ppf(1 - alpha / 2)
    z_lo, z_hi = z - z_crit * se, z + z_crit * se
    rho_lo, rho_hi = np.tanh(z_lo), np.tanh(z_hi)
    print(f"r = {r:.4f}, 95% CI for rho: ({rho_lo:.4f}, {rho_hi:.4f})")
    return rho_lo, rho_hi

np.random.seed(42)
x = np.random.normal(0, 1, 80)
y = 0.6 * x + np.random.normal(0, 0.8, 80)
fisher_z_ci(x, y)

# 참 rho 와 구간의 포함률을 확인한다.
rho_true = 0.6 / np.sqrt(0.6**2 + 0.8**2)
print(f"\n참 rho = 0.6/sqrt(0.6^2+0.8^2) = {rho_true:.4f}")
r_obs = stats.pearsonr(x, y)[0]
z = np.arctanh(r_obs); se = 1 / np.sqrt(len(x) - 3)
lo, hi = np.tanh(z - 1.96 * se), np.tanh(z + 1.96 * se)
print(f"구간이 참값을 덮는가: {lo < rho_true < hi}")
print(f"양팔: 아래 {r_obs - lo:.4f}, 위 {hi - r_obs:.4f}  "
      f"(비 {(r_obs - lo) / (hi - r_obs):.3f})")

rng = np.random.default_rng(2024)
B, n = 50_000, 80
cov_fisher = cov_naive = 0
for _ in range(B):
    a = rng.standard_normal(n)
    b = 0.6 * a + rng.normal(0, 0.8, n)
    r = np.corrcoef(a, b)[0, 1]
    zz = np.arctanh(r)
    if np.tanh(zz - 1.96 * se) < rho_true < np.tanh(zz + 1.96 * se):
        cov_fisher += 1
    half = (1 - r**2) / np.sqrt(n)        # 변환 없이 r 눈금에서 바로 만든 구간
    if r - 1.96 * half < rho_true < r + 1.96 * half:
        cov_naive += 1
print(f"\n포함률 ({B:,} 번, n = {n}, 명목 95%)")
print(f"  Fisher z 구간      = {cov_fisher / B:.4f}")
print(f"  변환 없는 대칭 구간 = {cov_naive / B:.4f}")

r = 0.6110, 95% CI for rho: (0.4519, 0.7324)

참 rho = 0.6/sqrt(0.6^2+0.8^2) = 0.6000
구간이 참값을 덮는가: True
양팔: 아래 0.1591, 위 0.1214  (비 1.310)

포함률 (50,000 번, n = 80, 명목 95%)
  Fisher z 구간      = 0.9497
  변환 없는 대칭 구간 = 0.9367

np.random.seed(42)
# 집단마다 x 의 평균이 커질수록 y 의 기준선도 함께 올라간다. 이것이 교란이다.
# 집단 안의 기울기는 셋 다 -0.5 로 음인데, 합쳐 놓으면 양이 된다.
groups = {"Group A": (50, 0.2, 2, -0.5),
          "Group B": (50, 0.5, 5, -0.5),
          "Group C": (50, 0.8, 8, -0.5)}

all_x, all_y = [], []
for name, (n, xm, yb, slope) in groups.items():
    x = np.random.normal(xm, 0.15, n)
    y = yb + slope * x + np.random.normal(0, 0.3, n)
    all_x.extend(x)
    all_y.extend(y)

all_x, all_y = np.array(all_x), np.array(all_y)
m_all, b_all = np.polyfit(all_x, all_y, 1)
r_overall, _ = stats.pearsonr(all_x, all_y)

print(f"Overall slope: {m_all:.2f} (positive)")
print(f"Within-group slope: -0.5 (negative)")
print(f"Overall r = {r_overall:.4f}")

# 합친 기울기의 이론값: 집단 사이 기울기와 집단 안 기울기의 분산가중 평균
mx = np.array([0.2, 0.5, 0.8])
yb = np.array([2.0, 5.0, 8.0])
within_slope = -0.5
my = yb + within_slope * mx                      # 집단 중심의 y 좌표
b_between = (my[2] - my[0]) / (mx[2] - mx[0])    # 중심 세 점을 잇는 기울기
var_between = mx.var()
var_within = 0.15**2
slope_theory = ((var_between * b_between + var_within * within_slope)
                / (var_between + var_within))
print(f"\n집단 중심 = {list(zip(mx, my.round(2)))}")
print(f"집단 사이 기울기 b = {b_between:.4f},  집단 안 기울기 = {within_slope}")
print(f"Var(집단 사이) = {var_between:.4f},  Var(집단 안) = {var_within:.4f}")
print(f"합친 기울기 이론값 = {slope_theory:.4f}   (표본 {m_all:.2f})")

vy_within = within_slope**2 * var_within + 0.3**2
cov = var_between * b_between + var_within * within_slope
r_theory = cov / np.sqrt((var_between + var_within) * (my.var() + vy_within))
print(f"합친 상관 이론값   = {r_theory:.4f}   (표본 {r_overall:.4f})")

print(f"\n기준선을 뒤집으면 (8, 5, 2)")
my2 = np.array([8.0, 5.0, 2.0]) + within_slope * mx
b2 = (my2[2] - my2[0]) / (mx[2] - mx[0])
print(f"  집단 사이 기울기 = {b2:.4f}  ->  합친 기울기 = "
      f"{(var_between * b2 + var_within * within_slope) / (var_between + var_within):.4f}")
print(f"부호가 뒤집히는 문턱: b > {-within_slope * var_within / var_between:.4f}")

Overall slope: 6.50 (positive)
Within-group slope: -0.5 (negative)
Overall r = 0.8593

집단 중심 = [(0.2, 1.9), (0.5, 4.75), (0.8, 7.6)]
집단 사이 기울기 b = 9.5000,  집단 안 기울기 = -0.5
Var(집단 사이) = 0.0600,  Var(집단 안) = 0.0225
합친 기울기 이론값 = 6.7727   (표본 6.50)
합친 상관 이론값   = 0.8287   (표본 0.8593)

기준선을 뒤집으면 (8, 5, 2)
  집단 사이 기울기 = -10.5000  ->  합친 기울기 = -7.7727
부호가 뒤집히는 문턱: b > 0.1875

n = 200
np.random.seed(42)

# 앞과 같은 구조다. Z 가 X 와 Y 를 함께 끌어 상관을 만든다.
Z = np.random.normal(0, 1, n)
X = 0.7 * Z + np.random.normal(0, 0.5, n)
Y = 0.6 * Z + np.random.normal(0, 0.5, n)

r_xy, p_xy = stats.pearsonr(X, Y)
r_xz, _ = stats.pearsonr(X, Z)
r_yz, _ = stats.pearsonr(Y, Z)

r_xy_z = (r_xy - r_xz * r_yz) / np.sqrt((1 - r_xz**2) * (1 - r_yz**2))

print(f"r(X, Y)    = {r_xy:.4f}  (appears significant)")
print(f"r(X,Y | Z) = {r_xy_z:.4f}  (nearly vanishes)")

# 모집단 값을 닫힌 꼴로 적는다.
rho_xz = 0.7 / np.sqrt(0.7**2 + 0.5**2)
rho_yz = 0.6 / np.sqrt(0.6**2 + 0.5**2)
rho_xy = (0.7 * 0.6) / np.sqrt((0.7**2 + 0.5**2) * (0.6**2 + 0.5**2))
part = (rho_xy - rho_xz * rho_yz) / np.sqrt((1 - rho_xz**2) * (1 - rho_yz**2))
print(f"\n모집단: rho_XZ = {rho_xz:.6f},  rho_YZ = {rho_yz:.6f}")
print(f"        rho_XY = rho_XZ * rho_YZ = {rho_xz * rho_yz:.6f}  (= {rho_xy:.6f})")
print(f"        편상관 = {part:.2e}   (정확히 0)")
se_xy = (1 - rho_xy**2) / np.sqrt(n)
print(f"\n표본 r_XY = {r_xy:.4f}   SE = {se_xy:.4f}   z = {(r_xy - rho_xy) / se_xy:+.3f}")
se_p = 1 / np.sqrt(n - 3)
print(f"표본 편상관 = {r_xy_z:.4f}   SE = {se_p:.4f}   z = {r_xy_z / se_p:+.3f}")
print(f"\n잔차로 구해도 같은가: ", end="")
ex = X - np.polyfit(Z, X, 1)[0] * Z - np.polyfit(Z, X, 1)[1]
ey = Y - np.polyfit(Z, Y, 1)[0] * Z - np.polyfit(Z, Y, 1)[1]
print(f"{stats.pearsonr(ex, ey)[0]:.6f}  (닫힌 꼴 {r_xy_z:.6f})")

r(X, Y)    = 0.5699  (appears significant)
r(X,Y | Z) = -0.0201  (nearly vanishes)

모집단: rho_XZ = 0.813733,  rho_YZ = 0.768221
        rho_XY = rho_XZ * rho_YZ = 0.625127  (= 0.625127)
        편상관 = 0.00e+00   (정확히 0)

표본 r_XY = 0.5699   SE = 0.0431   z = -1.283
표본 편상관 = -0.0201   SE = 0.0712   z = -0.283

잔차로 구해도 같은가: -0.020130  (닫힌 꼴 -0.020130)

def spurious_correlations_demo(n_vars=100, n_obs=30):
    """서로 완전히 무관한 변수 100개에서 유의한 상관이 몇 쌍이나 나오는지 센다.

    쌍이 4950개이므로 유의수준 5%에서 247쌍쯤은 그냥 나온다. 자료를 훑다가
    찾아낸 상관 하나를 그대로 보고하면 안 되는 까닭이 여기 있다.
    """
    np.random.seed(42)
    data = np.random.normal(0, 1, (n_obs, n_vars))
    n_pairs = n_vars * (n_vars - 1) // 2
    p_values = []
    for i in range(n_vars):
        for j in range(i + 1, n_vars):
            _, p = stats.pearsonr(data[:, i], data[:, j])
            p_values.append(p)
    p_values = np.array(p_values)
    n_sig = np.sum(p_values < 0.05)
    print(f"Pairs tested: {n_pairs}")
    print(f"Significant at 0.05: {n_sig} ({100*n_sig/n_pairs:.1f}%)")
    print(f"Expected false positives: {0.05 * n_pairs:.0f}")

spurious_correlations_demo()

# 240 이라는 값이 247.5 에서 얼마나 벗어난 것인가
rng = np.random.default_rng(11)
n_vars, n_obs, B = 100, 30, 300
iu = np.triu_indices(n_vars, 1)
counts = []
for _ in range(B):
    d = rng.standard_normal((n_obs, n_vars))
    R = np.corrcoef(d, rowvar=False)[iu]
    t = R * np.sqrt((n_obs - 2) / (1 - R**2))
    counts.append((2 * stats.t.sf(np.abs(t), n_obs - 2) < 0.05).sum())
counts = np.array(counts)
print(f"\n같은 실험을 {B} 번 되풀이")
print(f"  유의 쌍 개수: 평균 {counts.mean():.1f}  표준편차 {counts.std(ddof=1):.1f}")
print(f"  이론 기댓값 0.05 * 4950 = {0.05 * 4950:.1f}")
print(f"  쌍이 서로 독립이라면 표준편차 = sqrt(4950*0.05*0.95) = "
      f"{np.sqrt(4950 * 0.05 * 0.95):.1f}")
print(f"  관측값 240 의 z = {(240 - 247.5) / counts.std(ddof=1):+.2f}")
print(f"  5/50/95 백분위 = {np.percentile(counts, [5, 50, 95]).round(0)}")

Pairs tested: 4950
Significant at 0.05: 240 (4.8%)
Expected false positives: 248

같은 실험을 300 번 되풀이
  유의 쌍 개수: 평균 246.0  표준편차 15.3
  이론 기댓값 0.05 * 4950 = 247.5
  쌍이 서로 독립이라면 표준편차 = sqrt(4950*0.05*0.95) = 15.3
  관측값 240 의 z = -0.49
  5/50/95 백분위 = [223. 245. 271.]

def compare_two_correlations(r1, n1, r2, n2, alpha=0.05):
    """서로 독립인 두 표본의 상관계수가 다른지 검정한다.

    각각을 z 로 옮기면 차이의 분포가 정규가 되므로 z 검정을 쓸 수 있다.
    두 표본이 겹치지 않을 때만 이 방법이 맞다.
    """
    z1, z2 = np.arctanh(r1), np.arctanh(r2)
    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"r1={r1:.4f} (n={n1}), r2={r2:.4f} (n={n2})")
    print(f"z = {z_stat:.4f}, p = {p_value:.4f}")

compare_two_correlations(r1=0.72, n1=100, r2=0.65, n2=120)

# 차이에 대한 신뢰구간을 z 눈금에서 만든다.
z1, z2 = np.arctanh(0.72), np.arctanh(0.65)
se = np.sqrt(1 / (100 - 3) + 1 / (120 - 3))
print(f"\nz1 = {z1:.6f},  z2 = {z2:.6f},  차이 = {z1 - z2:.6f}")
print(f"SE = sqrt(1/97 + 1/117) = {se:.6f}")
lo, hi = (z1 - z2) - 1.96 * se, (z1 - z2) + 1.96 * se
print(f"zeta1 - zeta2 의 95% 구간 = ({lo:+.4f}, {hi:+.4f})   0 을 품는다")
print(f"\nrho2 = 0.65 로 고정하면 rho1 이 들어갈 수 있는 범위")
print(f"  ({np.tanh(z2 + lo):.4f},  {np.tanh(z2 + hi):.4f})")

r1=0.7200 (n=100), r2=0.6500 (n=120)
z = 0.9638, p = 0.3351

z1 = 0.907645,  z2 = 0.775299,  차이 = 0.132346
SE = sqrt(1/97 + 1/117) = 0.137318
zeta1 - zeta2 의 95% 구간 = (-0.1368, +0.4015)   0 을 품는다

rho2 = 0.65 로 고정하면 rho1 이 들어갈 수 있는 범위
  (0.5639,  0.8264)

import numpy as np
from scipy import stats

np.random.seed(42)
n = 150
x = np.random.uniform(-np.pi, np.pi, n)
y = np.cos(x) + np.random.normal(0, 0.3, n)

r_p, _ = stats.pearsonr(x, y)
r_s, _ = stats.spearmanr(x, y)
r_k, _ = stats.kendalltau(x, y)
print(f"Pearson r  = {r_p:.4f}")   # 0.0580
print(f"Spearman   = {r_s:.4f}")   # 0.1024
print(f"Kendall    = {r_k:.4f}")   # 0.0851

Pearson r  = 0.0580
Spearman   = 0.1024
Kendall    = 0.0851

import numpy as np
from scipy import stats

r = 0.45
n = 50
z = np.arctanh(r)
se = 1 / np.sqrt(n - 3)
z_crit = stats.norm.ppf(0.995)
lo = np.tanh(z - z_crit * se)
hi = np.tanh(z + z_crit * se)
print(f"99% CI: ({lo:.4f}, {hi:.4f})")

99% CI: (0.1085, 0.6965)

import numpy as np

np.random.seed(0)
# 1집단: 기준선도 낮고 x 도 작다
x1 = np.random.normal(1, 0.3, 60)
y1 = 2 - 1.0 * x1 + np.random.normal(0, 0.5, 60)

# 2집단: 기준선도 높고 x 도 크다
x2 = np.random.normal(4, 0.3, 60)
y2 = 10 - 1.0 * x2 + np.random.normal(0, 0.5, 60)

# 집단 안에서의 기울기
m1, _ = np.polyfit(x1, y1, 1)
m2, _ = np.polyfit(x2, y2, 1)

# 합쳐 놓았을 때의 기울기 — 부호가 뒤집힌다
x_all = np.concatenate([x1, x2])
y_all = np.concatenate([y1, y2])
m_all, _ = np.polyfit(x_all, y_all, 1)

print(f"Group 1 slope: {m1:.3f}")
print(f"Group 2 slope: {m2:.3f}")
print(f"Overall slope: {m_all:.3f}")

Group 1 slope: -1.130
Group 2 slope: -1.119
Overall slope: 1.531

import numpy as np
from scipy import stats

np.random.seed(42)
n_vars, n_obs = 200, 25
data = np.random.normal(0, 1, (n_obs, n_vars))

p_values = []
for i in range(n_vars):
    for j in range(i + 1, n_vars):
        _, p = stats.pearsonr(data[:, i], data[:, j])
        p_values.append(p)

p_values = np.array(p_values)
n_pairs = len(p_values)
n_sig = np.sum(p_values < 0.05)
bonf_threshold = 0.05 / n_pairs
n_sig_bonf = np.sum(p_values < bonf_threshold)

print(f"Total pairs: {n_pairs}")
print(f"Significant at 0.05: {n_sig}")
print(f"Bonferroni threshold: {bonf_threshold:.2e}")
print(f"Significant after Bonferroni: {n_sig_bonf}")

Total pairs: 19900
Significant at 0.05: 961
Bonferroni threshold: 2.51e-06
Significant after Bonferroni: 0