인과추론 모의실험¶

import numpy as np
from scipy import stats

np.random.seed(21)
n = 300

# Z 가 X 와 Y 를 함께 움직인다. X 와 Y 사이에는 직접 연결이 전혀 없다.
# 그런데도 둘은 상관을 보인다 — 교란변수가 만드는 가짜 상관이다.
Z = np.random.randn(n)
X = 0.6 * Z + np.random.randn(n) * 0.5
Y = 0.8 * Z + np.random.randn(n) * 0.5

# 모집단 값을 닫힌 꼴로 적는다.
var_x, var_y, cov_xy = 0.6**2 + 0.5**2, 0.8**2 + 0.5**2, 0.6 * 0.8
rho_xy = cov_xy / np.sqrt(var_x * var_y)
rho_xz, rho_yz = 0.6 / np.sqrt(var_x), 0.8 / np.sqrt(var_y)
print(f"Var(X) = {var_x:.4f},  Var(Y) = {var_y:.4f},  Cov(X,Y) = {cov_xy:.4f}")
print(f"rho_XY = {rho_xy:.6f}")
print(f"rho_XZ * rho_YZ = {rho_xz:.6f} * {rho_yz:.6f} = {rho_xz * rho_yz:.6f}")
print(f"둘이 같은가: {abs(rho_xy - rho_xz * rho_yz) < 1e-15}")
print(f"\n표본 r(X,Y) = {stats.pearsonr(X, Y)[0]:.4f}   "
      f"SE = {(1 - rho_xy**2) / np.sqrt(n):.4f}   "
      f"z = {(stats.pearsonr(X, Y)[0] - rho_xy) / ((1 - rho_xy**2) / np.sqrt(n)):+.3f}")

# Z 가 없는 모형으로도 같은 (X, Y) 를 만들 수 있다.
beta = cov_xy / var_x
resid_var = var_y - beta * cov_xy
print(f"\nZ 없이 X -> Y 직접효과만으로 같은 분포를 만드는 모형:")
print(f"  X ~ N(0, {var_x:.4f}),  Y = {beta:.6f} X + N(0, {resid_var:.6f})")
X2 = np.random.normal(0, np.sqrt(var_x), 200_000)
Y2 = beta * X2 + np.random.normal(0, np.sqrt(resid_var), 200_000)
print(f"  그 모형의 상관 = {np.corrcoef(X2, Y2)[0, 1]:.4f}   (교란 모형 {rho_xy:.4f})")

Var(X) = 0.6100,  Var(Y) = 0.8900,  Cov(X,Y) = 0.4800
rho_XY = 0.651450
rho_XZ * rho_YZ = 0.768221 * 0.847998 = 0.651450
둘이 같은가: True

표본 r(X,Y) = 0.6512   SE = 0.0332   z = -0.007

Z 없이 X -> Y 직접효과만으로 같은 분포를 만드는 모형:
  X ~ N(0, 0.6100),  Y = 0.786885 X + N(0, 0.512295)
  그 모형의 상관 = 0.6516   (교란 모형 0.6515)

# 부분상관은 Z 로 설명되는 몫을 X 와 Y 에서 걷어 낸 뒤의 상관이다.
# 위 자료에서는 걷어 내고 나면 거의 0 만 남아야 한다.
r_xy, _ = stats.pearsonr(X, Y)
r_xz, _ = stats.pearsonr(X, Z)
r_yz, _ = stats.pearsonr(Y, Z)

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

print(f"Pearson r(X, Y)       = {r_xy:.3f}")
print(f"Partial r(X, Y | Z)   = {r_partial:.3f}")

# (1) 모집단 편상관은 정확히 0 이다.
var_x, var_y = 0.6**2 + 0.5**2, 0.8**2 + 0.5**2
rho_xy = (0.6 * 0.8) / np.sqrt(var_x * var_y)
rho_xz, rho_yz = 0.6 / np.sqrt(var_x), 0.8 / np.sqrt(var_y)
part_pop = (rho_xy - rho_xz * rho_yz) / np.sqrt((1 - rho_xz**2) * (1 - rho_yz**2))
se = 1 / np.sqrt(len(X) - 3)
print(f"\n모집단 편상관 = {part_pop:.2e}   (정확히 0)")
print(f"표본 편상관   = {r_partial:.4f}   SE = 1/sqrt(n-3) = {se:.4f}   "
      f"z = {r_partial / se:+.3f}")

# (2) 충돌부에서는 방향이 반대다. Z = X + Y + e 이면 X 와 Y 는 원래 독립인데
#     Z 를 통제하는 순간 음의 상관이 생긴다.
rho_xz_c = 1 / np.sqrt(3)          # Corr(X, X+Y+e)
part_collider = (0 - rho_xz_c**2) / (1 - rho_xz_c**2)
print(f"\n충돌부 Z = X + Y + e  (X, Y, e ~ N(0,1) 독립)")
print(f"  rho_XZ = rho_YZ = 1/sqrt(3) = {rho_xz_c:.6f},  rho_XY = 0")
print(f"  편상관 이론값 = (0 - 1/3)/(1 - 1/3) = {part_collider:.6f}")
rng = np.random.default_rng(7)
m = 400_000
Xc = rng.standard_normal(m)
Yc = rng.standard_normal(m)
Zc = Xc + Yc + rng.standard_normal(m)
a, b, c = (np.corrcoef(Xc, Yc)[0, 1], np.corrcoef(Xc, Zc)[0, 1],
           np.corrcoef(Yc, Zc)[0, 1])
print(f"  모의: r_XY = {a:+.4f}  ->  편상관 = "
      f"{(a - b * c) / np.sqrt((1 - b**2) * (1 - c**2)):+.4f}")

Pearson r(X, Y)       = 0.651
Partial r(X, Y | Z)   = 0.052

모집단 편상관 = 0.00e+00   (정확히 0)
표본 편상관   = 0.0518   SE = 1/sqrt(n-3) = 0.0580   z = +0.893

충돌부 Z = X + Y + e  (X, Y, e ~ N(0,1) 독립)
  rho_XZ = rho_YZ = 1/sqrt(3) = 0.577350,  rho_XY = 0
  편상관 이론값 = (0 - 1/3)/(1 - 1/3) = -0.500000
  모의: r_XY = +0.0018  ->  편상관 = -0.4996

rng = np.random.default_rng(42)

# 두 집단 모두 안에서는 기울기가 -0.4 로 음이다. 그런데 B 집단이 x 도 크고
# y 의 기준선도 높아, 둘을 합쳐 놓으면 전체 기울기가 양으로 뒤집힌다.
n_a, n_b = 100, 100
x_a = rng.uniform(10, 30, n_a)
y_a = -0.4 * x_a + 30 + rng.normal(0, 2, n_a)

x_b = rng.uniform(25, 50, n_b)
y_b = -0.4 * x_b + 45 + rng.normal(0, 2, n_b)

# 집단별 모집단 상관을 닫힌 꼴로 구한다.
print(f"\n{'집단':>4s} {'x 범위':>12s} {'Var(x)':>9s} {'sd(x)':>8s} {'sd(y)':>8s} "
      f"{'모집단 rho':>11s} {'표본 r':>9s}")
for lab, lo, hi, xs, ys in [("A", 10, 30, x_a, y_a), ("B", 25, 50, x_b, y_b)]:
    vx = (hi - lo) ** 2 / 12
    sy = np.sqrt(0.4**2 * vx + 2**2)
    rho = -0.4 * np.sqrt(vx) / sy
    print(f"{lab:>4s} {f'[{lo}, {hi}]':>12s} {vx:9.4f} {np.sqrt(vx):8.4f} {sy:8.4f} "
          f"{rho:+11.6f} {stats.pearsonr(xs, ys)[0]:+9.4f}")
print(f"\n집단 안 기울기는 둘 다 -0.4 로 같은데 상관은 다르다.")
print(f"  rho = -0.4 * sd(x) / sd(y) 이고 sd(x) 가 B 에서 더 크기 때문이다.")

  집단         x 범위    Var(x)    sd(x)    sd(y)     모집단 rho      표본 r
   A     [10, 30]   33.3333   5.7735   3.0551   -0.755929   -0.7371
   B     [25, 50]   52.0833   7.2169   3.5119   -0.821995   -0.8246

집단 안 기울기는 둘 다 -0.4 로 같은데 상관은 다르다.
  rho = -0.4 * sd(x) / sd(y) 이고 sd(x) 가 B 에서 더 크기 때문이다.

# 합친 상관과 집단별 상관의 부호가 갈리는 것을 확인한다. 이것이 심슨의 역설이다.
x_all = np.concatenate([x_a, x_b])
y_all = np.concatenate([y_a, y_b])

r_all, _ = stats.pearsonr(x_all, y_all)
r_a, _ = stats.pearsonr(x_a, y_a)
r_b, _ = stats.pearsonr(x_b, y_b)

print(f"Aggregate  r = {r_all:+.3f}")
print(f"Subgroup A r = {r_a:+.3f}")
print(f"Subgroup B r = {r_b:+.3f}")

# 합친 상관의 이론값: 집단 안 몫과 집단 사이 몫을 더한다.
mx = np.array([20.0, 37.5])                     # 집단별 x 평균
vx = np.array([(30 - 10) ** 2 / 12, (50 - 25) ** 2 / 12])
my = -0.4 * mx + np.array([30.0, 45.0])         # 집단별 y 평균
vy = 0.4**2 * vx + 2**2

var_within, var_between = vx.mean(), mx.var()
VX = var_within + var_between
VY = vy.mean() + my.var()
COV = -0.4 * var_within + np.mean((mx - mx.mean()) * (my - my.mean()))
print(f"\nVar(x) = 집단 안 {var_within:.4f} + 집단 사이 {var_between:.4f} = {VX:.4f}")
print(f"Var(y) = {VY:.4f}")
print(f"Cov    = 집단 안 {-0.4 * var_within:+.4f} + 집단 사이 "
      f"{np.mean((mx - mx.mean()) * (my - my.mean())):+.4f} = {COV:+.4f}")
print(f"합친 기울기 이론값 = {COV / VX:.6f}")
print(f"합친 상관   이론값 = {COV / np.sqrt(VX * VY):.6f}   (표본 {r_all:+.4f})")

b_between = (my[1] - my[0]) / (mx[1] - mx[0])
print(f"\n집단 사이 기울기 = {b_between:.6f}")
print(f"분산가중 평균    = {(var_between * b_between + var_within * (-0.4)) / VX:.6f}")
print(f"부호가 뒤집히는 문턱: b > 0.4 * {var_within:.4f} / {var_between:.4f} = "
      f"{0.4 * var_within / var_between:.6f}")

Aggregate  r = +0.321
Subgroup A r = -0.737
Subgroup B r = -0.825

Var(x) = 집단 안 42.7083 + 집단 사이 76.5625 = 119.2708
Var(y) = 26.8333
Cov    = 집단 안 -17.0833 + 집단 사이 +35.0000 = +17.9167
합친 기울기 이론값 = 0.150218
합친 상관   이론값 = 0.316703   (표본 +0.3208)

집단 사이 기울기 = 0.457143
분산가중 평균    = 0.150218
부호가 뒤집히는 문턱: b > 0.4 * 42.7083 / 76.5625 = 0.223129

import matplotlib.pyplot as plt

# 점을 집단별로 다른 표식으로 찍고, 그 위에 합친 자료의 회귀직선을 얹는다.
# 직선의 기울기가 각 무리의 기울기와 반대 방향인 것이 한눈에 보인다.
fig, ax = plt.subplots(figsize=(8, 5))
ax.scatter(x_a, y_a, label='Group A', alpha=0.6)
ax.scatter(x_b, y_b, label='Group B', alpha=0.6, marker='s')

slope, intercept = np.polyfit(x_all, y_all, 1)
xs = np.linspace(x_all.min(), x_all.max(), 100)
ax.plot(xs, slope * xs + intercept, 'k--', linewidth=2,
        label=f'Aggregate OLS (r={r_all:+.2f})')
ax.set_xlabel('X')
ax.set_ylabel('Y')
ax.set_title("Simpson's Paradox")
ax.legend()
plt.tight_layout()
plt.show()

# 그림에서 읽을 수치를 적어 둔다.
print(f"\n합친 회귀직선: 기울기 {slope:+.4f}, 절편 {intercept:.4f}, r = {r_all:+.4f}")
print(f"집단 안 기울기 (참값) = -0.4,  표본 A {np.polyfit(x_a, y_a, 1)[0]:+.4f}, "
      f"B {np.polyfit(x_b, y_b, 1)[0]:+.4f}")
print(f"\nx 범위: A [{x_a.min():.1f}, {x_a.max():.1f}], "
      f"B [{x_b.min():.1f}, {x_b.max():.1f}]")
ov_lo, ov_hi = max(x_a.min(), x_b.min()), min(x_a.max(), x_b.max())
print(f"겹치는 구간 [{ov_lo:.1f}, {ov_hi:.1f}] 에 든 점 = "
      f"A {((x_a >= ov_lo) & (x_a <= ov_hi)).sum()} 개, "
      f"B {((x_b >= ov_lo) & (x_b <= ov_hi)).sum()} 개  (전체 200 개 중)")
print(f"그 구간에서 y 평균: A {y_a[(x_a >= ov_lo) & (x_a <= ov_hi)].mean():.2f}, "
      f"B {y_b[(x_b >= ov_lo) & (x_b <= ov_hi)].mean():.2f}")
print(f"\n집단 꼬리표를 지우면: 합친 r = {r_all:+.4f} 하나만 남는다")
print(f"집단 꼬리표가 있으면: {r_a:+.4f}, {r_b:+.4f} 와 집단 사이 기울기 "
      f"{(30 - 22) / (37.5 - 20):+.4f} 를 갈라 볼 수 있다")

합친 회귀직선: 기울기 +0.1606, 절편 21.7781, r = +0.3208
집단 안 기울기 (참값) = -0.4,  표본 A -0.3923, B -0.3977

x 범위: A [10.1, 29.5], B [25.3, 49.5]
겹치는 구간 [25.3, 29.5] 에 든 점 = A 20 개, B 28 개  (전체 200 개 중)
그 구간에서 y 평균: A 19.04, B 33.68

집단 꼬리표를 지우면: 합친 r = +0.3208 하나만 남는다
집단 꼬리표가 있으면: -0.7371, -0.8246 와 집단 사이 기울기 +0.4571 를 갈라 볼 수 있다

import numpy as np
from scipy import stats

np.random.seed(10)
n = 500
Z = np.random.randn(n)
X = 0.9 * Z + np.random.normal(0, 0.3, n)
Y = 0.3 * Z + np.random.normal(0, 0.3, n)

r_xy, _ = stats.pearsonr(X, Y)
r_xz, _ = stats.pearsonr(X, Z)
r_yz, _ = stats.pearsonr(Y, Z)
r_partial = (r_xy - r_xz * r_yz) / np.sqrt((1 - r_xz**2) * (1 - r_yz**2))

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

r(X, Y)     = 0.6715
r(X, Y | Z) = -0.0115

import numpy as np
import matplotlib.pyplot as plt

np.random.seed(42)
groups = [(100, 0, 50), (100, 10, 30), (100, 20, 10)]
# (크기, x 중심, y 절편). 집단 안에서는 기울기가 양이다

fig, ax = plt.subplots()
all_x, all_y = [], []

for n, xc, yb in groups:
    x = np.random.normal(xc, 1.5, n)
    y = yb + 2 * (x - xc) + np.random.normal(0, 1, n)
    ax.scatter(x, y, alpha=0.5, s=15)
    all_x.extend(x)
    all_y.extend(y)

all_x, all_y = np.array(all_x), np.array(all_y)
m, b = np.polyfit(all_x, all_y, 1)
xs = np.linspace(all_x.min(), all_x.max(), 100)
ax.plot(xs, m * xs + b, 'k--', lw=2, label=f'Aggregate slope = {m:.2f}')
ax.legend()
ax.set_xlabel('X')
ax.set_ylabel('Y')
plt.tight_layout()
plt.show()

import numpy as np
from scipy import stats

np.random.seed(42)
n = 200
x_a = np.random.uniform(0, 20, n)
y_a = 10 + 0.5 * x_a + np.random.normal(0, 2, n)

x_b = np.random.uniform(0, 20, n)
y_b = 10 - 0.5 * x_b + np.random.normal(0, 2, n)

x_all = np.concatenate([x_a, x_b])
y_all = np.concatenate([y_a, y_b])

r_a, _ = stats.pearsonr(x_a, y_a)
r_b, _ = stats.pearsonr(x_b, y_b)
r_all, _ = stats.pearsonr(x_all, y_all)

print(f"Group A r = {r_a:+.3f}")
print(f"Group B r = {r_b:+.3f}")
print(f"Aggregate r = {r_all:+.3f}")

Group A r = +0.832
Group B r = -0.825
Aggregate r = -0.023