부분상관¶

import numpy as np
from scipy import stats

np.random.seed(42)
n = 100

# 기온이 아이스크림 판매와 익사 사고를 함께 끌어올린다. 둘 사이에 직접
# 연결은 없다. 그런데도 상관은 크게 나온다.
z = np.random.normal(0, 1, n)          # 기온 — 교란변수
x = 2 * z + np.random.normal(0, 1, n)  # 아이스크림 판매량
y = 1.5 * z + np.random.normal(0, 1, n)  # 익사 사고 건수

# 먼저 쌍마다의 상관을 구한다.
r_xy = stats.pearsonr(x, y)[0]
r_xz = stats.pearsonr(x, z)[0]
r_yz = stats.pearsonr(y, z)[0]

# 부분상관 공식. Z 로 설명되는 몫을 걷어 낸 뒤 남는 상관이다.
# 위 자료에서는 거의 0 으로 떨어져야 한다.
r_xy_z = (r_xy - r_xz * r_yz) / (
    np.sqrt(1 - r_xz**2) * np.sqrt(1 - r_yz**2)
)

print(f"Marginal r(X,Y)     = {r_xy:.4f}")
print(f"Partial r(X,Y | Z)  = {r_xy_z:.4f}")

# (1)(2) 이론값과 표본값을 나란히
t_xz, t_yz = 2 / np.sqrt(5), 1.5 / np.sqrt(3.25)
t_xy = 3 / np.sqrt(16.25)
print(f"\n이론: rho_xz = {t_xz:.6f}  rho_yz = {t_yz:.6f}  rho_xy = {t_xy:.6f}")
print(f"      rho_xz * rho_yz = {t_xz * t_yz:.6f}  <- rho_xy 와 같다"
      f"   차 {t_xy - t_xz * t_yz:.1e}")
for name, s, th in [("r_xy", r_xy, t_xy), ("r_xz", r_xz, t_xz), ("r_yz", r_yz, t_yz)]:
    se = (1 - th**2) / np.sqrt(n)
    print(f"  {name}: 표본 {s:.4f}  이론 {th:.4f}  차 {s - th:+.4f}"
          f"  SE {se:.4f}  ({(s - th) / se:+.2f} SE)")
se_p = 1 / np.sqrt(n - 3)
print(f"  부분상관: 표본 {r_xy_z:+.4f}  이론 0  SE {se_p:.4f}"
      f"  ({r_xy_z / se_p:+.3f} SE)")
zf = np.arctanh(r_xy_z)
print(f"  피셔 z 구간 [{np.tanh(zf - 1.96 * se_p):+.4f},"
      f" {np.tanh(zf + 1.96 * se_p):+.4f}]")

# (3) 부호가 뒤집히는 문턱
print(f"\n부분상관을 0 으로 만드는 r_yz = r_xy / r_xz = {r_xy / r_xz:.6f}")
print(f"실제 r_yz = {r_yz:.6f}   문턱을 {r_yz - r_xy / r_xz:+.6f} 넘었다")

# (4) 잔차 경로로도 같은 값이 나오는가
e_x = x - np.polyval(np.polyfit(z, x, 1), z)
e_y = y - np.polyval(np.polyfit(z, y, 1), z)
print(f"\n잔차의 상관 {stats.pearsonr(e_x, e_y)[0]:.6f}"
      f"   공식 {r_xy_z:.6f}")

Marginal r(X,Y)     = 0.7190
Partial r(X,Y | Z)  = -0.0109

이론: rho_xz = 0.894427  rho_yz = 0.832050  rho_xy = 0.744208
      rho_xz * rho_yz = 0.744208  <- rho_xy 와 같다   차 0.0e+00
  r_xy: 표본 0.7190  이론 0.7442  차 -0.0252  SE 0.0446  (-0.57 SE)
  r_xz: 표본 0.8724  이론 0.8944  차 -0.0220  SE 0.0200  (-1.10 SE)
  r_yz: 표본 0.8276  이론 0.8321  차 -0.0045  SE 0.0308  (-0.15 SE)
  부분상관: 표본 -0.0109  이론 0  SE 0.1015  (-0.107 SE)
  피셔 z 구간 [-0.2069, +0.1859]

부분상관을 0 으로 만드는 r_yz = r_xy / r_xz = 0.824157
실제 r_yz = 0.827584   문턱을 +0.003427 넘었다

잔차의 상관 -0.010897   공식 -0.010897

import warnings
warnings.filterwarnings("ignore")

import numpy as np

def pcorr(x, y, z):
    """z 를 통제한 x, y 의 부분상관 (잔차의 상관)."""
    M = np.column_stack([np.ones(len(z)), z])
    bx = np.linalg.lstsq(M, x, rcond=None)[0]
    by = np.linalg.lstsq(M, y, rcond=None)[0]
    return np.corrcoef(x - M @ bx, y - M @ by)[0, 1]

rng = np.random.default_rng(22001)
n = 5_000
z = rng.standard_normal(n)
x = 0.8 * z + rng.normal(0, 0.6, n)
y = 0.7 * z + rng.normal(0, 0.714, n)

rxy = np.corrcoef(x, y)[0, 1]
rxz = np.corrcoef(x, z)[0, 1]
ryz = np.corrcoef(y, z)[0, 1]
form = (rxy - rxz * ryz) / np.sqrt((1 - rxz**2) * (1 - ryz**2))
print(f"  r_xy = {rxy:.6f}, r_xz = {rxz:.6f}, r_yz = {ryz:.6f}")
print(f"  공식:      {form:.6f}")
print(f"  잔차 상관: {pcorr(x, y, z):.6f}")

  r_xy = 0.562896, r_xz = 0.806170, r_yz = 0.697538
  공식:      0.001325
  잔차 상관: 0.001325

z 를 만들고
x = 0.8z + 잡음      (x 와 y 를 잇는 직접 경로는 없다)
y = 0.7z + 잡음

M = [1, z, z², z³] 로 회귀한 뒤 잔차의 상관
  → "z 의 비선형 효과까지 제거한" 부분상관

import warnings
warnings.filterwarnings("ignore")

import numpy as np

def pcorr(x, y, z):
    M = np.column_stack([np.ones(len(z)), z])
    bx = np.linalg.lstsq(M, x, rcond=None)[0]
    by = np.linalg.lstsq(M, y, rcond=None)[0]
    return np.corrcoef(x - M @ bx, y - M @ by)[0, 1]

rng = np.random.default_rng(22001)
n = 200_000
cases = {}

z = rng.standard_normal(n)                       # (1) 교란
x = 0.8 * z + rng.normal(0, 0.6, n)
y = 0.8 * z + rng.normal(0, 0.6, n)
cases["교란 Z→X, Z→Y"] = (x, y, z)

x = rng.standard_normal(n)                       # (2) 매개
z = 0.8 * x + rng.normal(0, 0.6, n)
y = 0.8 * z + rng.normal(0, 0.6, n)
cases["매개 X→Z→Y"] = (x, y, z)

x = rng.standard_normal(n)                       # (3) 충돌
y = rng.standard_normal(n)
z = 0.7 * x + 0.7 * y + rng.normal(0, 0.5, n)
cases["충돌 X→Z←Y"] = (x, y, z)

z = rng.standard_normal(n)                       # (4) 억제
x = 0.8 * z + rng.normal(0, 0.6, n)
y = 0.6 * x - 0.8 * z + rng.normal(0, 0.5, n)
cases["억제 (부호 상쇄)"] = (x, y, z)

interp = {"교란 Z→X, Z→Y": "가짜 상관이 사라진다",
          "매개 X→Z→Y": "매개를 막으면 사라진다",
          "충돌 X→Z←Y": "없던 상관이 생긴다",
          "억제 (부호 상쇄)": "숨은 관계가 드러난다"}
print(f"{'구조':>20s} {'주변 r_xy':>10s} {'부분 r_xy·z':>12s} {'해석':>24s}")
for lab, (x, y, z) in cases.items():
    print(f"{lab:>20s} {np.corrcoef(x, y)[0, 1]:10.4f} "
          f"{pcorr(x, y, z):12.4f} {interp[lab]:>24s}")

                  구조    주변 r_xy    부분 r_xy·z                       해석
         교란 Z→X, Z→Y     0.6416       0.0052              가짜 상관이 사라진다
            매개 X→Z→Y     0.6417       0.0039             매개를 막으면 사라진다
            충돌 X→Z←Y     0.0013      -0.6617               없던 상관이 생긴다
          억제 (부호 상쇄)    -0.0589       0.5849              숨은 관계가 드러난다

교란:  X 와 Y 는 실제로 무관하다.  Z 를 통제한 0.00 이 참이다.
매개:  X 는 Z 를 통해 Y 에 영향을 준다.  총효과 0.64 가 참이고,
       Z 를 통제하면 그 효과를 지워 버린다.

병원 입원 환자만 보면 (Z = 입원)
  당뇨(X)와 골절(Y)이 음의 상관을 보인다
  → 실제로는 무관하지만, 둘 중 하나만 있어도 입원하므로
    입원자 중에서는 "둘 다 있는 사람"이 드물다

직접:   X → Y   (+)
간접:   X ← Z → Y  (Z→X 는 +, Z→Y 는 −)

주변 r = 직접 + 간접 ≈ 0
부분 r = 직접만 = +0.585

import warnings
warnings.filterwarnings("ignore")

import numpy as np

def pcorr(x, y, Z):
    Z = Z.reshape(len(x), -1)
    M = np.column_stack([np.ones(len(x)), Z])
    rx = x - M @ np.linalg.lstsq(M, x, rcond=None)[0]
    ry = y - M @ np.linalg.lstsq(M, y, rcond=None)[0]
    return np.corrcoef(rx, ry)[0, 1]

rng = np.random.default_rng(22002)
B = 4_000
print("참 부분상관 0.4 를 덮는 비율 (명목 95%)")
print(f"{'k':>4s} {'n':>5s} {'df=n-3-k 사용':>14s} {'df=n-3 사용(틀림)':>17s}")
for k in [1, 3, 5]:
    for n in [20, 50]:
        a = b = 0
        for _ in range(B):
            Z = rng.standard_normal((n, k))
            e1, e2 = rng.standard_normal(n), rng.standard_normal(n)
            u = Z.sum(1)
            x = u + e1
            y = u + 0.4 * e1 + np.sqrt(1 - 0.16) * e2
            r = pcorr(x, y, Z)
            for use, which in [(n - 3 - k, "a"), (n - 3, "b")]:
                se = 1 / np.sqrt(max(use, 1))
                lo = np.tanh(np.arctanh(r) - 1.96 * se)
                hi = np.tanh(np.arctanh(r) + 1.96 * se)
                if which == "a":
                    a += lo <= 0.4 <= hi
                else:
                    b += lo <= 0.4 <= hi
        print(f"{k:4d} {n:5d} {a / B:14.4f} {b / B:17.4f}")

참 부분상관 0.4 를 덮는 비율 (명목 95%)
   k     n    df=n-3-k 사용     df=n-3 사용(틀림)
   1    20         0.9433            0.9370
   1    50         0.9507            0.9487
   3    20         0.9467            0.9213
   3    50         0.9477            0.9403
   5    20         0.9495            0.9025
   5    50         0.9483            0.9325

from scipy import stats

rng = np.random.default_rng(22002)
B = 8_000
n = 30
print("\n참 부분상관 0 인데 통제변수를 k 개 넣는다 (n=30, 명목 0.05)")
print(f"{'k':>4s} {'오류율':>8s} {'자유도':>8s}")
for k in [0, 1, 3, 5, 10, 20]:
    a = 0
    for _ in range(B):
        Z = rng.standard_normal((n, max(k, 1)))
        x, y = rng.standard_normal(n), rng.standard_normal(n)
        if k == 0:
            r, df = np.corrcoef(x, y)[0, 1], n - 2
        else:
            r, df = pcorr(x, y, Z[:, :k]), n - 2 - k
        t = r * np.sqrt(df / (1 - r**2))
        a += 2 * stats.t.sf(abs(t), df) < 0.05
    print(f"{k:4d} {a / B:8.4f} {n - 2 - k:8d}")

참 부분상관 0 인데 통제변수를 k 개 넣는다 (n=30, 명목 0.05)
   k      오류율      자유도
   0   0.0495       28
   1   0.0529       27
   3   0.0498       25
   5   0.0527       23
  10   0.0509       18
  20   0.0519        8

import warnings
warnings.filterwarnings("ignore")

import numpy as np

def resid(a, Z):
    M = np.column_stack([np.ones(len(a)), Z.reshape(len(a), -1)])
    return a - M @ np.linalg.lstsq(M, a, rcond=None)[0]

rng = np.random.default_rng(22003)
n = 200_000
z = rng.standard_normal(n)
x = 0.7 * z + rng.normal(0, 0.714, n)
y = 0.5 * x + 0.5 * z + rng.normal(0, 0.6, n)

rxy = np.corrcoef(x, y)[0, 1]
part = np.corrcoef(resid(x, z), resid(y, z))[0, 1]
semi = np.corrcoef(resid(x, z), y)[0, 1]
print(f"주변 상관      r_xy      = {rxy:.4f}")
print(f"부분 상관      r_xy·z    = {part:.4f}")
print(f"준부분 상관    r_y(x·z)  = {semi:.4f}")
print(f"\n관계: 준부분 = 부분 × √(1 - r_yz²)")
ryz = np.corrcoef(y, z)[0, 1]
print(f"  r_yz = {ryz:.4f}")
print(f"  부분 × √(1-r_yz²) = {part * np.sqrt(1 - ryz**2):.4f}")

print(f"\n회귀와의 연결 (y ~ x + z)")
M = np.column_stack([np.ones(n), x, z])
b = np.linalg.lstsq(M, y, rcond=None)[0]
yhat = M @ b
R2_full = 1 - ((y - yhat)**2).sum() / ((y - y.mean())**2).sum()
M0 = np.column_stack([np.ones(n), z])
b0 = np.linalg.lstsq(M0, y, rcond=None)[0]
R2_z = 1 - ((y - M0 @ b0)**2).sum() / ((y - y.mean())**2).sum()
print(f"  R²(y ~ x + z) = {R2_full:.4f}")
print(f"  R²(y ~ z)     = {R2_z:.4f}")
print(f"  차이 ΔR²      = {R2_full - R2_z:.4f}")
print(f"  준부분상관의 제곱 = {semi**2:.4f}   ← 같다")

주변 상관      r_xy      = 0.7735
부분 상관      r_xy·z    = 0.5154
준부분 상관    r_y(x·z)  = 0.3279

관계: 준부분 = 부분 × √(1 - r_yz²)
  r_yz = 0.7716
  부분 × √(1-r_yz²) = 0.3279

회귀와의 연결 (y ~ x + z)
  R²(y ~ x + z) = 0.7028
  R²(y ~ z)     = 0.5953
  차이 ΔR²      = 0.1075
  준부분상관의 제곱 = 0.1075   ← 같다

pingouin.partial_corr(data, x, y, covar=...)      → 부분상관
pingouin.partial_corr(data, x, y, x_covar=...)    → 준부분상관
SPSS 의 "Part correlation" 열                     → 준부분상관

import warnings
warnings.filterwarnings("ignore")

import numpy as np
from scipy import stats

def pcorr(x, y, z):
    M = np.column_stack([np.ones(len(z)), z])
    rx = x - M @ np.linalg.lstsq(M, x, rcond=None)[0]
    ry = y - M @ np.linalg.lstsq(M, y, rcond=None)[0]
    return np.corrcoef(rx, ry)[0, 1]

rng = np.random.default_rng(22004)
n = 200_000
print("부분상관이 0 이지만 조건부 독립이 아닌 예")

z = rng.standard_normal(n)
x = rng.standard_normal(n)
y = x**2 + z + rng.normal(0, 0.5, n)      # y 는 x 에 의존하지만 선형이 아니다
print(f"  (가) y = x² + z:  r_xy·z = {pcorr(x, y, z):+.4f}")
print(f"       |x| 와 y 의 부분상관 = {pcorr(np.abs(x), y, z):+.4f}  ← 의존이 드러난다")

x = rng.standard_normal(n)
s = np.abs(x) + 0.2
y = rng.normal(0, 1, n) * s + z            # y 의 분산이 x 에 의존
print(f"\n  (나) y 의 분산이 |x| 에 비례:  r_xy·z = {pcorr(x, y, z):+.4f}")
print(f"       |x| 와 |y| 의 부분상관 = {pcorr(np.abs(x), np.abs(y), z):+.4f}  ← 드러난다")

print("\n선형 통제의 한계: z 의 효과가 비선형일 때")
z = rng.standard_normal(n)
x = z**2 + rng.normal(0, 0.5, n)
y = z**2 + rng.normal(0, 0.5, n)           # x, y 는 z² 를 통해서만 연결
print(f"  주변 r_xy          = {np.corrcoef(x, y)[0, 1]:+.4f}")
print(f"  z 로 선형 통제      = {pcorr(x, y, z):+.4f}   ← 남아 있다")
Z2 = np.column_stack([z, z**2])
M = np.column_stack([np.ones(n), Z2])
rx = x - M @ np.linalg.lstsq(M, x, rcond=None)[0]
ry = y - M @ np.linalg.lstsq(M, y, rcond=None)[0]
print(f"  z, z² 로 통제       = {np.corrcoef(rx, ry)[0, 1]:+.4f}   ← 사라진다")

부분상관이 0 이지만 조건부 독립이 아닌 예
  (가) y = x² + z:  r_xy·z = +0.0043
       |x| 와 y 의 부분상관 = +0.8827  ← 의존이 드러난다

  (나) y 의 분산이 |x| 에 비례:  r_xy·z = -0.0010
       |x| 와 |y| 의 부분상관 = +0.3550  ← 드러난다

선형 통제의 한계: z 의 효과가 비선형일 때
  주변 r_xy          = +0.8888
  z 로 선형 통제      = +0.8888   ← 남아 있다
  z, z² 로 통제       = +0.0019   ← 사라진다

변수 Z 를 통제할 것인가?
    │
    ├─ Z 가 X 와 Y 의 공통 원인 (교란)   ──→ 통제한다
    ├─ Z 가 X 의 결과이면서 Y 의 원인 (매개)
    │     ├─ 직접효과가 궁금하다  ──→ 통제한다
    │     └─ 총효과가 궁금하다    ──→ 통제하지 않는다
    ├─ Z 가 X 와 Y 의 공통 결과 (충돌부) ──→ 절대 통제하지 않는다
    ├─ Z 가 X 의 원인일 뿐 (도구 성격)   ──→ 통제하면 검정력만 잃는다
    └─ 모르겠다                          ──→ 인과 그림을 먼저 그린다

import pingouin as pg
pg.partial_corr(data=df, x="X", y="Y", covar=["Z1","Z2"])
  → r, CI95%, p-val  (부분상관)

pg.partial_corr(data=df, x="X", y="Y", x_covar=["Z"])
  → 준부분상관

method="spearman" 으로 순위 기반 부분상관도 가능

운동 시간과 수면의 질의 관계 (나이·성별·BMI 통제)

  주변 상관    r = 0.34  [0.25, 0.42]
  부분 상관    r = 0.21  [0.11, 0.30]   (k = 3, df = 195)

통제변수는 선행 문헌에서 두 변수 모두의 원인으로 알려진 것만
선택했다. 나이는 비선형 효과가 예상되어 자연 스플라인(df=3)으로
넣었다.

준부분상관 = 0.18 → 운동 시간이 수면의 질 변동의 3.2% 를
추가로 설명한다.