Pearson의 r 검정 (상관에 대한 t-검정)¶

import numpy as np
from scipy import stats

np.random.seed(42)
n = 25
x = np.random.normal(5, 2, n)
y = 0.5 * x + np.random.normal(0, 1, n)

# 귀무가설은 rho = 0 이다. 유의하다는 말은 "상관이 0 은 아니다"까지일 뿐,
# 상관이 크다는 뜻도 인과가 있다는 뜻도 아니다.
r, p_value = stats.pearsonr(x, y)
print(f"r = {r:.4f}")
print(f"p-value = {p_value:.4f}")

# 정의대로 손으로 구해 본다. 이 통계량은 자유도 n-2 인 t 를 따른다.
# 회귀에서 기울기를 검정할 때 쓰는 t 값과 정확히 같은 값이다.
t_stat = r * np.sqrt(n - 2) / np.sqrt(1 - r**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-value = {p_manual:.4f}")

# 검정은 rho = 0 만 다루지만 신뢰구간은 rho 가 놓일 범위를 말해 준다.
# 그래서 검정보다 구간 쪽이 대개 더 쓸모 있다.
z = np.arctanh(r)
se_z = 1 / np.sqrt(n - 3)
ci_z = (z - 1.96 * se_z, z + 1.96 * se_z)
ci_r = (np.tanh(ci_z[0]), np.tanh(ci_z[1]))
print(f"95% CI for rho: ({ci_r[0]:.3f}, {ci_r[1]:.3f})")

# (1) 0.0000 뒤에 무엇이 있는지 자릿수를 열어 본다.
print(f"\nscipy p = {p_value:.6e}")
print(f"손계산 p = {p_manual:.6e}")
print(f"두 값의 상대차 = {abs(p_value - p_manual) / p_value:.2e}")
print("1 - cdf 는 t 가 커지면 무너진다 (df = 23):")
for t in (6.5638, 10, 20, 30):
    print(f"  t = {t:6.2f}   2*(1-cdf) = {2 * (1 - stats.t.cdf(t, 23)):.6e}   "
          f"2*sf = {2 * stats.t.sf(t, 23):.6e}")

# (2) 피셔 z 눈금에서는 대칭이고 되돌리면 비대칭이다.
print(f"\nz = arctanh(r) = {z:.6f},  SE = 1/sqrt(n-3) = {se_z:.6f}")
print(f"z 눈금 구간 = ({ci_z[0]:.6f}, {ci_z[1]:.6f})   "
      f"양팔 {z - ci_z[0]:.6f} / {ci_z[1] - z:.6f}  (대칭)")
print(f"r 눈금 구간 = ({ci_r[0]:.6f}, {ci_r[1]:.6f})   "
      f"양팔 {r - ci_r[0]:.4f} / {ci_r[1] - r:.4f}  "
      f"(비 {(r - ci_r[0]) / (ci_r[1] - r):.3f})")
half = (1 - r**2) * 1.96 * se_z
print(f"델타법 (1-r^2)*1.96*SE = {half:.4f}  ->  ({r - half:.4f}, {r + half:.4f})")

r = 0.8074
p-value = 0.0000
t-statistic = 6.5638
Manual p-value = 0.0000
95% CI for rho: (0.605, 0.912)

scipy p = 1.069634e-06
손계산 p = 1.069634e-06
두 값의 상대차 = 3.45e-11
1 - cdf 는 t 가 커지면 무너진다 (df = 23):
  t =   6.56   2*(1-cdf) = 1.069527e-06   2*sf = 1.069527e-06
  t =  10.00   2*(1-cdf) = 7.644623e-10   2*sf = 7.644623e-10
  t =  20.00   2*(1-cdf) = 4.440892e-16   2*sf = 4.835701e-16
  t =  30.00   2*(1-cdf) = 0.000000e+00   2*sf = 6.045668e-20

z = arctanh(r) = 1.119617,  SE = 1/sqrt(n-3) = 0.213201
z 눈금 구간 = (0.701744, 1.537490)   양팔 0.417873 / 0.417873  (대칭)
r 눈금 구간 = (0.605473, 0.911698)   양팔 0.2020 / 0.1043  (비 1.937)
델타법 (1-r^2)*1.96*SE = 0.1454  ->  (0.6620, 0.9529)

import numpy as np
from scipy import stats

rng = np.random.default_rng(0)
n = 20
print(f"{'rho':>6}{'평균 r':>10}{'왜도':>10}{'2.5%':>9}{'97.5%':>9}")
for rho in (0.0, 0.5, 0.8, 0.95):
    X = rng.multivariate_normal([0, 0], [[1, rho], [rho, 1]], size=(200_000, n))
    xa, ya = X[:, :, 0], X[:, :, 1]
    xa = xa - xa.mean(1, keepdims=True); ya = ya - ya.mean(1, keepdims=True)
    r = (xa*ya).sum(1) / np.sqrt((xa**2).sum(1) * (ya**2).sum(1))
    q = np.percentile(r, [2.5, 97.5])
    print(f"{rho:>6.2f}{r.mean():>10.4f}{stats.skew(r):>10.4f}"
          f"{q[0]:>9.4f}{q[1]:>9.4f}")

   rho    평균 r      왜도     2.5%    97.5%
  0.00    0.0001    0.0051  -0.4419   0.4437
  0.50    0.4894   -0.6462   0.0867   0.7774
  0.80    0.7920   -1.1909   0.5711   0.9212
  0.95    0.9474   -1.5739   0.8830   0.9813

import numpy as np

rng = np.random.default_rng(1)
print(f"{'n':>5}{'rho':>7}{'E[r]':>10}{'편향':>10}{'보정 후':>10}")
for n in (10, 20, 50):
    for rho in (0.3, 0.8):
        X = rng.multivariate_normal([0, 0], [[1, rho], [rho, 1]],
                                    size=(200_000, n))
        xa, ya = X[:, :, 0], X[:, :, 1]
        xa = xa - xa.mean(1, keepdims=True); ya = ya - ya.mean(1, keepdims=True)
        r = (xa*ya).sum(1) / np.sqrt((xa**2).sum(1) * (ya**2).sum(1))
        adj = r * (1 + (1 - r**2) / (2*n))
        print(f"{n:>5}{rho:>7.1f}{r.mean():>10.4f}"
              f"{r.mean()-rho:>+10.4f}{adj.mean():>10.4f}")

    n    rho      E[r]      편향     보정 후
   10    0.3    0.2851   -0.0149    0.2948
   10    0.8    0.7820   -0.0180    0.7950
   20    0.3    0.2925   -0.0075    0.2983
   20    0.8    0.7921   -0.0079    0.7991
   50    0.3    0.2970   -0.0030    0.2996
   50    0.8    0.7969   -0.0031    0.7997

import numpy as np
from scipy import stats
from scipy.spatial.distance import pdist, squareform

def dcor(x, y):
    x, y = np.asarray(x, float), np.asarray(y, float)
    A = squareform(pdist(x[:, None])); B = squareform(pdist(y[:, None]))
    A = A - A.mean(0) - A.mean(1)[:, None] + A.mean()
    B = B - B.mean(0) - B.mean(1)[:, None] + B.mean()
    dcov2 = (A*B).mean(); dvx = (A*A).mean(); dvy = (B*B).mean()
    return np.sqrt(dcov2 / np.sqrt(dvx*dvy)) if dvx*dvy > 0 else 0.0

rng = np.random.default_rng(2)
n = 300
cases = {
    "독립":        (rng.normal(0,1,n), rng.normal(0,1,n)),
    "선형":        None,
    "이차(U자)":   None,
    "원":          None,
}
x = rng.normal(0, 1, n)
cases["선형"]      = (x, 0.7*x + rng.normal(0, 0.5, n))
cases["이차(U자)"] = (x, x**2 + rng.normal(0, 0.3, n))
th = rng.uniform(0, 2*np.pi, n)
cases["원"]        = (np.cos(th), np.sin(th))

print(f"{'관계':>12}{'Pearson r':>12}{'Spearman':>11}{'거리상관':>11}")
for k, (a, b) in cases.items():
    print(f"{k:>12}{stats.pearsonr(a,b).statistic:>12.4f}"
          f"{stats.spearmanr(a,b).statistic:>11.4f}{dcor(a,b):>11.4f}")

        관계   Pearson r   Spearman     거리상관
        독립     -0.0670    -0.0557     0.1261
        선형      0.7950     0.7784     0.7481
   이차(U자)     -0.0210    -0.0703     0.5048
         원      0.0489     0.0434     0.2068

import numpy as np
from scipy import stats
import statsmodels.api as sm

rng = np.random.default_rng(4)
x = rng.normal(0, 1, 25); y = 0.5*x + rng.normal(0, 1, 25)
r, p = stats.pearsonr(x, y)
fit = sm.OLS(y, sm.add_constant(x)).fit()
print(f"  상관 t = {r*np.sqrt((len(x)-2)/(1-r**2)):.6f}   p = {p:.6f}")
print(f"  회귀 t = {fit.tvalues[1]:.6f}   p = {fit.pvalues[1]:.6f}")

import numpy as np
from scipy import stats

rng = np.random.default_rng(6)
n, rho = 200_000, 0.6
X = rng.multivariate_normal([0,0], [[1,rho],[rho,1]], n)
x, y = X[:,0], X[:,1]

print(f"  전체 r = {np.corrcoef(x,y)[0,1]:.4f}  (참 rho={rho})")
print(f"{'x 범위 제한':>16}{'남은 비율':>11}{'r':>9}{'기울기':>9}")
for lim in (3.0, 1.5, 1.0, 0.5):
    m = np.abs(x) < lim
    r = np.corrcoef(x[m], y[m])[0,1]
    b = np.polyfit(x[m], y[m], 1)[0]
    print(f"{f'|x| < {lim}':>16}{m.mean():>11.3f}{r:>9.4f}{b:>9.4f}")

  전체 r = 0.5978  (참 rho=0.6)
         x 범위 제한      남은 비율        r      기울기
       |x| < 3.0      0.997   0.5927   0.5985
       |x| < 1.5      0.866   0.4845   0.5994
       |x| < 1.0      0.683   0.3741   0.6012
       |x| < 0.5      0.384   0.2085   0.6031

import numpy as np
from scipy import stats

for r in (0.1, 0.2, 0.3, 0.5, 0.7):
    # 이변량 정규에서 x 상위 절반과 하위 절반의 y 평균 차 (sigma_y 단위)
    d = 2 * r * stats.norm.pdf(0) / 0.5      # E[x | x>0] = 2*phi(0)
    print(f"  r={r:.1f}: r^2={r**2:.3f}   "
          f"상위 절반과 하위 절반의 y 평균 차 = {d:.3f} sigma")

  r=0.1: r^2=0.010   상위 절반과 하위 절반의 y 평균 차 = 0.160 sigma
  r=0.2: r^2=0.040   상위 절반과 하위 절반의 y 평균 차 = 0.319 sigma
  r=0.3: r^2=0.090   상위 절반과 하위 절반의 y 평균 차 = 0.479 sigma
  r=0.5: r^2=0.250   상위 절반과 하위 절반의 y 평균 차 = 0.798 sigma
  r=0.7: r^2=0.490   상위 절반과 하위 절반의 y 평균 차 = 1.117 sigma

1. 산점도를 그렸는가?                      -> 모든 것의 출발점
2. 관계가 선형인가?                        -> 아니면 스피어만/거리상관/변환
3. 이상치가 있는가?                        -> 있으면 순위 측도 병기
4. 범위가 선발로 잘렸는가?                 -> 잘렸으면 기울기를 보고
5. 관측이 독립인가?                        -> 군집/시계열이면 보정 필요
6. 두 변수가 같은 것의 두 측정은 아닌가?   -> 비율/부분 상관 함정
7. 집계 자료인가?                          -> 생태학적 오류 주의