상관 분석 시연¶

import numpy as np
from scipy import stats

np.random.seed(42)

# 기울기 0.8 의 선형 관계에 표준편차 8 짜리 잡음을 얹는다.
n = 120
x = np.random.uniform(10, 60, n)
noise = np.random.normal(0, 8, n)
y = 0.8 * x + 5 + noise

sigma_x = np.sqrt(50 ** 2 / 12)      # 균등분포의 표준편차
signal = 0.8 * sigma_x               # 기울기가 만들어 내는 Y 의 변동 폭
print(f"sigma_X = {sigma_x:.4f},  0.8 sigma_X = {signal:.4f}")
for sigma in (8, 16, 100):
    print(f"  sigma = {sigma:3d}  ->  rho = {signal / np.sqrt(signal ** 2 + sigma ** 2):.4f}")

sigma_X = 14.4338,  0.8 sigma_X = 11.5470
  sigma =   8  ->  rho = 0.8220
  sigma =  16  ->  rho = 0.5852
  sigma = 100  ->  rho = 0.1147

# 세 측도를 함께 구한다. 관계가 선형이고 이상치가 없으면 셋이 비슷하게 나온다.
# 값이 크게 갈린다면 관계가 곡선이거나 이상치가 있다는 신호다.
r_pearson, p_pearson = stats.pearsonr(x, y)
r_spearman, p_spearman = stats.spearmanr(x, y)
r_kendall, p_kendall = stats.kendalltau(x, y)

print(f"Pearson  r = {r_pearson:.4f}  (p = {p_pearson:.2e})")
print(f"Spearman rho = {r_spearman:.4f}  (p = {p_spearman:.2e})")
print(f"Kendall  tau = {r_kendall:.4f}  (p = {p_kendall:.2e})")

# (1) 의 닫힌 꼴과 견준다.
var_x = 50 ** 2 / 12
rho = 0.8 * var_x / np.sqrt(var_x * (0.64 * var_x + 64))
se = (1 - rho ** 2) / np.sqrt(n)
print(f"\n이론 rho   = {rho:.4f}   SE = {se:.4f}   z = {(r_pearson - rho) / se:+.3f}")

# tau 와 rho_s 에는 닫힌 꼴이 없다. 같은 모형에서 큰 표본을 거듭 뽑아 모집단 값을 잰다.
rng = np.random.default_rng(1)
taus, rhos = [], []
for _ in range(40):
    bx = rng.uniform(10, 60, 50_000)
    by = 0.8 * bx + 5 + rng.normal(0, 8, 50_000)
    taus.append(stats.kendalltau(bx, by).statistic)
    rhos.append(stats.spearmanr(bx, by).statistic)
print(f"모의 tau   = {np.mean(taus):.4f} +- {np.std(taus) / np.sqrt(40):.4f}   표본 {r_kendall:.4f}")
print(f"모의 rho_s = {np.mean(rhos):.4f} +- {np.std(rhos) / np.sqrt(40):.4f}   표본 {r_spearman:.4f}")

# 이변량 정규에서만 성립하는 닫힌 꼴. 여기서는 x 가 균등이므로 맞지 않는다.
print(f"\n정규 가정 공식: tau = (2/pi)arcsin(rho)   = {2 / np.pi * np.arcsin(rho):.4f}")
print(f"                rho_s = (6/pi)arcsin(rho/2) = {6 / np.pi * np.arcsin(rho / 2):.4f}")

Pearson  r = 0.8221  (p = 1.21e-30)
Spearman rho = 0.8363  (p = 1.38e-32)
Kendall  tau = 0.6420  (p = 2.54e-25)

이론 rho   = 0.8220   SE = 0.0296   z = +0.003
모의 tau   = 0.6292 +- 0.0002   표본 0.6420
모의 rho_s = 0.8326 +- 0.0002   표본 0.8363

정규 가정 공식: tau = (2/pi)arcsin(rho)   = 0.6143
                rho_s = (6/pi)arcsin(rho/2) = 0.8089

# Spearman 은 "순위에 대한 Pearson"이라는 정의를 그대로 확인한다.
rx, ry = stats.rankdata(x), stats.rankdata(y)
r_rank = stats.pearsonr(rx, ry)[0]
print(f"Pearson r on ranks = {r_rank:.4f}")
print(f"Spearman rho       = {r_spearman:.4f}")
# 두 값이 같아야 한다.

# 간편 공식도 같은 값을 주는지 본다.
d2 = np.sum((rx - ry) ** 2)
shortcut = 1 - 6 * d2 / (n * (n ** 2 - 1))
print(f"\n동점 개수 = {n - len(np.unique(rx))} (x), {n - len(np.unique(ry))} (y)")
print(f"sum d^2 = {d2:.0f},  S = n(n^2-1)/12 = {n * (n ** 2 - 1) / 12:.0f}")
print(f"간편 공식 = {shortcut:.6f}   spearmanr = {r_spearman:.6f}")
print(f"세 값이 같은가: {np.allclose([r_rank, shortcut], r_spearman)}")

Pearson r on ranks = 0.8363
Spearman rho       = 0.8363

동점 개수 = 0 (x), 0 (y)
sum d^2 = 47138,  S = 143990
간편 공식 = 0.836315   spearmanr = 0.836315
세 값이 같은가: True

import matplotlib.pyplot as plt

# 상관계수는 숫자 하나일 뿐이므로 반드시 그림과 함께 본다.
slope, intercept, _, _, _ = stats.linregress(x, y)

fig, ax = plt.subplots(figsize=(8, 5))
ax.scatter(x, y, alpha=0.6, edgecolors='k', linewidths=0.3)
x_line = np.array([x.min(), x.max()])
ax.plot(x_line, intercept + slope * x_line, 'r-', linewidth=2,
        label=f'OLS: y = {slope:.2f}x + {intercept:.2f}')
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_title(f'Pearson r = {r_pearson:.3f}')
ax.legend()
plt.tight_layout()
plt.show()

# 그림에서 읽을 수치를 찍어 둔다.
sx, sy = x.std(ddof=1), y.std(ddof=1)
print(f"기울기 {slope:.4f}  =  r * s_y/s_x = {r_pearson:.4f} * {sy:.4f}/{sx:.4f} = {r_pearson * sy / sx:.4f}")
print(f"절편 {intercept:.4f}   (참값: 기울기 0.8, 절편 5)")
print(f"SE(beta) = sigma/(sigma_X sqrt(n)) = {8 / (np.sqrt(50 ** 2 / 12) * np.sqrt(n)):.4f}, "
      f"z = {(slope - 0.8) / (8 / (np.sqrt(50 ** 2 / 12) * np.sqrt(n))):+.3f}")
print(f"r^2 = {r_pearson ** 2:.4f}   잔차 표준편차 = {np.std(y - (intercept + slope * x), ddof=2):.4f}  (참값 8)")

기울기 0.7622  =  r * s_y/s_x = 0.8221 * 13.8062/14.8914 = 0.7622
절편 6.7816   (참값: 기울기 0.8, 절편 5)
SE(beta) = sigma/(sigma_X sqrt(n)) = 0.0506, z = -0.748
r^2 = 0.6758   잔차 표준편차 = 7.8939  (참값 8)

import numpy as np
from scipy import stats

np.random.seed(0)
n = 200
x = np.random.uniform(0, 20, n)
y = 3 * x + 2 + np.random.normal(0, 5, n)

r_p, p_p = stats.pearsonr(x, y)
r_s, p_s = stats.spearmanr(x, y)
r_k, p_k = stats.kendalltau(x, y)

print(f"Pearson  r = {r_p:.4f}, p = {p_p:.2e}")
print(f"Spearman rho = {r_s:.4f}, p = {p_s:.2e}")
print(f"Kendall  tau = {r_k:.4f}, p = {p_k:.2e}")

Pearson  r = 0.9601, p = 1.42e-111
Spearman rho = 0.9593, p = 1.14e-110
Kendall  tau = 0.8227, p = 4.63e-67

import numpy as np
from scipy import stats

np.random.seed(42)
x = np.random.uniform(0, 5, 100)
y = np.exp(x) + np.random.normal(0, 1, 100)

r_p, _ = stats.pearsonr(x, y)
r_s, _ = stats.spearmanr(x, y)
print(f"Pearson r = {r_p:.4f}")
print(f"Spearman rho = {r_s:.4f}")

Pearson r = 0.8556
Spearman rho = 0.9821

import numpy as np
from scipy import stats

def all_correlations(x, y):
    r_p, _ = stats.pearsonr(x, y)
    r_s, _ = stats.spearmanr(x, y)
    r_k, _ = stats.kendalltau(x, y)
    return {"pearson": r_p, "spearman": r_s, "kendall": r_k}

np.random.seed(42)
n = 100

# (가) 강한 선형 관계
x_a = np.random.normal(0, 1, n)
y_a = 2 * x_a + np.random.normal(0, 0.5, n)
print("Linear:", all_correlations(x_a, y_a))

# (나) 약한 비선형 관계
x_b = np.random.uniform(-3, 3, n)
y_b = x_b ** 2 + np.random.normal(0, 1, n)
print("Quadratic:", all_correlations(x_b, y_b))

# (다) 이상치가 섞인 경우
x_c = np.random.normal(0, 1, n)
y_c = 0.8 * x_c + np.random.normal(0, 0.3, n)
x_c[:3] = [6, -6, 7]
y_c[:3] = [-6, 6, -7]
print("Outliers:", all_correlations(x_c, y_c))

Linear: {'pearson': 0.9654943669720492, 'spearman': 0.9640324032403239, 'kendall': 0.8408080808080809}
Quadratic: {'pearson': -0.28540752875587533, 'spearman': -0.15697569756975696, 'kendall': -0.09292929292929294}
Outliers: {'pearson': -0.1603686865830897, 'spearman': 0.7779657965796579, 'kendall': 0.6993939393939395}

import warnings
warnings.filterwarnings("ignore")

import numpy as np
import pandas as pd

rng = np.random.default_rng(26002)

def blk(r, n):
    z = rng.standard_normal((n, 2))
    return z[:, 0], r * z[:, 0] + np.sqrt(1 - r**2) * z[:, 1]

n = 100
rows = []
a, b = blk(0.9, n)
rows += [dict(A=a[i], B=b[i], C=np.nan) for i in range(n)]
a, c = blk(0.9, n)
rows += [dict(A=a[i], B=np.nan, C=c[i]) for i in range(n)]
b, c = blk(-0.5, n)
rows += [dict(A=np.nan, B=b[i], C=c[i]) for i in range(n)]
d = pd.DataFrame(rows)

print(f"  결측 패턴: 행마다 한 변수씩 결측 (총 {len(d)} 행)")
print(f"  listwise 삭제 후 남는 행: {d.dropna().shape[0]}")
P = d.corr()
print("\n  pairwise 상관행렬:")
print(P.round(4).to_string())
ev = np.linalg.eigvalsh(P.values)
print(f"\n  고윳값 = {ev.round(4).tolist()}")
print(f"  최소 고윳값 = {ev.min():+.4f}")
lo = (P.loc['A', 'B'] * P.loc['A', 'C']
      - np.sqrt((1 - P.loc['A', 'B']**2) * (1 - P.loc['A', 'C']**2)))
print(f"  r_AB={P.loc['A','B']:.3f}, r_AC={P.loc['A','C']:.3f} 이면 "
      f"r_BC ≥ {lo:.4f} 여야 하는데 관측값은 {P.loc['B','C']:.4f}")

  결측 패턴: 행마다 한 변수씩 결측 (총 300 행)
  listwise 삭제 후 남는 행: 0

  pairwise 상관행렬:
        A       B       C
A  1.0000  0.9412  0.9022
B  0.9412  1.0000 -0.4710
C  0.9022 -0.4710  1.0000

  고윳값 = [-0.5602, 1.4704, 2.0898]
  최소 고윳값 = -0.5602
  r_AB=0.941, r_AC=0.902 이면 r_BC ≥ 0.7033 여야 하는데 관측값은 -0.4710

import numpy as np
from scipy import stats

rng = np.random.default_rng(26002)
B = 2_000
print("모든 변수가 독립인 자료 (n=30, 명목 0.05)")
print(f"{'p':>4s} {'쌍의 수':>8s} {'하나라도 유의':>13s} "
      f"{'이론 1-0.95^m':>15s} {'평균 유의 개수':>13s}")
for p in [5, 10, 20, 50]:
    m = p * (p - 1) // 2
    any_sig = cnt = 0
    for _ in range(B):
        X = rng.standard_normal((30, p))
        r = np.corrcoef(X, rowvar=False)[np.triu_indices(p, 1)]
        t = r * np.sqrt(28 / (1 - r**2))
        s = (2 * stats.t.sf(np.abs(t), 28) < 0.05).sum()
        cnt += s
        any_sig += s > 0
    print(f"{p:4d} {m:8d} {any_sig / B:13.4f} {1 - 0.95**m:15.4f} "
          f"{cnt / B:13.2f}")

모든 변수가 독립인 자료 (n=30, 명목 0.05)
   p     쌍의 수       하나라도 유의     이론 1-0.95^m      평균 유의 개수
   5       10        0.4095          0.4013          0.50
  10       45        0.9065          0.9006          2.27
  20      190        1.0000          0.9999          9.41
  50     1225        1.0000          1.0000         61.24

설문 문항 30 개를 다 넣고 상관행렬을 본다
  → 435 쌍
  → 우연히 22 개가 p < 0.05
  → 그중 "말이 되는" 것을 골라 이야기를 만든다
  → 재현되지 않는다

from statsmodels.stats.multitest import multipletests
rej, p_adj, _, _ = multipletests(pvals, alpha=0.05, method="fdr_bh")

import warnings
warnings.filterwarnings("ignore")

import numpy as np
from scipy import stats

rng = np.random.default_rng(26003)
B, n = 4_000, 100

def ar1(phi, n, rng):
    """AR(1) 계열을 정상 상태에서 시작해 생성한다."""
    e = rng.standard_normal(n)
    x = np.empty(n)
    x[0] = e[0] / np.sqrt(1 - phi**2)
    for i in range(1, n):
        x[i] = phi * x[i - 1] + e[i]
    return x

print("두 계열이 완전히 독립인데도 유의하게 나온다 (명목 0.05)")
print(f"{'구조':>22s} {'1종 오류':>9s} {'|r| 의 평균':>11s}")
gens = [
    ("독립 백색잡음", lambda: (rng.standard_normal(n), rng.standard_normal(n))),
    ("AR(1) φ=0.5", lambda: (ar1(0.5, n, rng), ar1(0.5, n, rng))),
    ("AR(1) φ=0.9", lambda: (ar1(0.9, n, rng), ar1(0.9, n, rng))),
    ("AR(1) φ=0.99", lambda: (ar1(0.99, n, rng), ar1(0.99, n, rng))),
    ("확률보행", lambda: (np.cumsum(rng.standard_normal(n)),
                       np.cumsum(rng.standard_normal(n)))),
    ("선형 추세 + 잡음", lambda: (np.arange(n) * 0.05 + rng.standard_normal(n),
                            np.arange(n) * 0.05 + rng.standard_normal(n))),
]
for lab, gen in gens:
    a, rs = 0, []
    for _ in range(B):
        x, y = gen()
        res = stats.pearsonr(x, y)
        a += res.pvalue < 0.05
        rs.append(abs(res.statistic))
    print(f"{lab:>22s} {a / B:9.4f} {np.mean(rs):11.4f}")

두 계열이 완전히 독립인데도 유의하게 나온다 (명목 0.05)
                    구조     1종 오류    |r| 의 평균
               독립 백색잡음    0.0505      0.0805
           AR(1) φ=0.5    0.1368      0.1043
           AR(1) φ=0.9    0.5065      0.2248
          AR(1) φ=0.99    0.7410      0.3871
                  확률보행    0.7542      0.4181
            선형 추세 + 잡음    1.0000      0.6791

def eff_n(x, y, n):
    """베틀렛 근사 유효 표본수."""
    r1 = np.corrcoef(x[:-1], x[1:])[0, 1]
    s1 = np.corrcoef(y[:-1], y[1:])[0, 1]
    return n * (1 - r1 * s1) / (1 + r1 * s1)

a = b = c = 0
for _ in range(B):
    x, y = ar1(0.9, n, rng), ar1(0.9, n, rng)
    a += stats.pearsonr(x, y).pvalue < 0.05
    b += stats.pearsonr(np.diff(x), np.diff(y)).pvalue < 0.05
    r = np.corrcoef(x, y)[0, 1]
    ne = max(eff_n(x, y, n), 4)
    t = r * np.sqrt((ne - 2) / (1 - r**2))
    c += 2 * stats.t.sf(abs(t), ne - 2) < 0.05
print("\nAR(1) φ=0.9 인 두 독립 계열 (n=100)")
print(f"  보정 없음:        {a / B:.4f}")
print(f"  1차 차분 후:      {b / B:.4f}")
print(f"  유효표본수 보정:  {c / B:.4f}")

AR(1) φ=0.9 인 두 독립 계열 (n=100)
  보정 없음:        0.5195
  1차 차분 후:      0.0575
  유효표본수 보정:  0.0375

import warnings
warnings.filterwarnings("ignore")

import numpy as np
from scipy import stats

rng = np.random.default_rng(26004)

def perm_p(x, y, stat, B=20_000, rng=None):
    """y 를 재배열해 귀무분포를 만든다."""
    obs = abs(stat(x, y))
    cnt = sum(abs(stat(x, rng.permutation(y))) >= obs for _ in range(B))
    return (cnt + 1) / (B + 1)

n = 12
x = rng.standard_normal(n)
y = 0.6 * x + np.sqrt(1 - 0.36) * rng.standard_normal(n)
print(f"정규 자료 (n={n})")
for lab, f in [("피어슨", lambda a, b: stats.pearsonr(a, b)),
               ("스피어만", lambda a, b: stats.spearmanr(a, b)),
               ("켄들", lambda a, b: stats.kendalltau(a, b))]:
    r = f(x, y)
    pp = perm_p(x, y, lambda a, b: f(a, b).statistic, rng=rng)
    print(f"  {lab:>6s}: 통계량={r.statistic:+.4f}  "
          f"이론 p={r.pvalue:.4f}  순열 p={pp:.4f}")

y2 = y.copy()
y2[0] += 12                                  # 이상점 하나
print(f"\n이상점을 하나 넣으면")
for lab, f in [("피어슨", lambda a, b: stats.pearsonr(a, b)),
               ("스피어만", lambda a, b: stats.spearmanr(a, b)),
               ("켄들", lambda a, b: stats.kendalltau(a, b))]:
    r = f(x, y2)
    pp = perm_p(x, y2, lambda a, b: f(a, b).statistic, rng=rng)
    print(f"  {lab:>6s}: 통계량={r.statistic:+.4f}  "
          f"이론 p={r.pvalue:.4f}  순열 p={pp:.4f}")

정규 자료 (n=12)
     피어슨: 통계량=+0.6471  이론 p=0.0229  순열 p=0.0243
    스피어만: 통계량=+0.6713  이론 p=0.0168  순열 p=0.0208
      켄들: 통계량=+0.4545  이론 p=0.0447  순열 p=0.0457

이상점을 하나 넣으면
     피어슨: 통계량=+0.4333  이론 p=0.1594  순열 p=0.1544
    스피어만: 통계량=+0.7343  이론 p=0.0065  순열 p=0.0082
      켄들: 통계량=+0.5152  이론 p=0.0210  순열 p=0.0219

1. 자료를 본다        ─ 결측, 이상점, 자료형, 관측의 독립성
2. 산점도를 그린다     ─ 모양, 군집, 이상점
3. 계수를 고른다       ─ 피어슨 / 스피어만 / 켄들
4. 계산한다            ─ 유효 표본수를 함께 기록
5. 구간을 만든다       ─ 피셔 z 또는 부트스트랩
6. 진단한다            ─ 잔차, 쿡 D, 두 계수의 차이
7. 보고한다            ─ 인과를 주장하지 않는다

□ 산점도를 그렸는가
□ 결측을 어떻게 처리했는가 (쌍마다 n 이 다른가)
□ 관측이 서로 독립인가 (시계열·군집·반복측정)
□ 이상점이 결과를 바꾸는가 (빼고 다시 계산)
□ 피어슨과 스피어만이 크게 다른가
□ 여러 쌍을 검정했다면 보정했는가
□ 신뢰구간을 보고했는가
□ 인과 해석을 피했는가

광고비와 매출의 관계 (월별, 2019-2024, n = 72)

  수준 자료:   r = 0.91  ← 보고하지 않음 (둘 다 시간 추세를 가짐)
  1차 차분:    r = 0.34,  95% CI [0.11, 0.53],  p = 0.004
               유효 표본수 71

  두 계열 모두 단위근 검정에서 비정상으로 판정되어 차분 후
  분석했다. 차분 자료의 자기상관은 유의하지 않았다.

  관찰자료이므로 "광고가 매출을 늘린다"고 결론지을 수 없다.
  경기·계절 등 공통 요인이 두 계열을 함께 움직였을 수 있다.