교란과 인과 시연¶

import numpy as np
from scipy import stats

np.random.seed(42)


def simulate_confounded_data(n=500, rho_tc=0.8):
    mean = [0, 0]
    cov = [[1, rho_tc], [rho_tc, 1]]
    tc = np.random.multivariate_normal(mean, cov, n)
    t = tc[:, 0]
    c = tc[:, 1]
    # Y는 오직 C에만 의존한다. T가 Y에 미치는 참 효과는 정확히 0이다.
    y = c + np.random.normal(0, 1, n)
    return t, c, y


t, c, y = simulate_confounded_data()
print(f"corr(T, C) = {np.corrcoef(t, c)[0, 1]:.4f}")
print(f"corr(T, Y) = {np.corrcoef(t, y)[0, 1]:.4f}   (참 인과효과는 0)")

# 닫힌 꼴과 견준다.
n, rho_tc, sigma = len(t), 0.8, 1.0
for lab, pop, samp in [("corr(T,C)", rho_tc, np.corrcoef(t, c)[0, 1]),
                       ("corr(T,Y)", rho_tc / np.sqrt(1 + sigma**2),
                        np.corrcoef(t, y)[0, 1])]:
    se = (1 - pop**2) / np.sqrt(n)
    print(f"\n{lab}: 모집단 {pop:.6f}  표본 {samp:.4f}  "
          f"SE {se:.4f}  z {(samp - pop) / se:+.3f}")

print(f"\n교란만으로 만들 수 있는 상관의 상한 = 1/sqrt(1+sigma^2)")
for s in (0.5, 1.0, 2.0, 3.0):
    print(f"  sigma = {s:.1f}  ->  {1 / np.sqrt(1 + s**2):.4f}")

corr(T, C) = 0.7908
corr(T, Y) = 0.5282   (참 인과효과는 0)

corr(T,C): 모집단 0.800000  표본 0.7908  SE 0.0161  z -0.574

corr(T,Y): 모집단 0.565685  표본 0.5282  SE 0.0304  z -1.231

교란만으로 만들 수 있는 상관의 상한 = 1/sqrt(1+sigma^2)
  sigma = 0.5  ->  0.8944
  sigma = 1.0  ->  0.7071
  sigma = 2.0  ->  0.4472
  sigma = 3.0  ->  0.3162

def compute_regressions(t, c, y):
    """교란변수를 빼고 넣은 두 회귀를 나란히 돌린다.

    짧은 회귀는 T 만 넣고, 긴 회귀는 C 까지 넣는다. 참 효과가 0 인데도
    짧은 회귀의 계수가 크게 나오는 것이 누락변수 편향이다.
    """
    # 짧은 회귀: Y ~ T
    slope_short, _, r_short, p_short, _ = stats.linregress(t, y)

    # 긴 회귀: Y ~ T + C. 절편을 위해 1 로 된 열을 앞에 붙인다.
    X = np.column_stack([np.ones(len(t)), t, c])
    beta = np.linalg.lstsq(X, y, rcond=None)[0]

    return {
        "short_slope": slope_short,
        "short_p": p_short,
        "long_beta_T": beta[1],
        "long_beta_C": beta[2],
    }


res = compute_regressions(t, c, y)
print(f"short_slope  = {res['short_slope']:>8.4f}   (p = {res['short_p']:.1e})")
print(f"long_beta_T  = {res['long_beta_T']:>8.4f}")
print(f"long_beta_C  = {res['long_beta_C']:>8.4f}")

# 누락변수 편향 항등식을 확인한다.
gamma = stats.linregress(t, c).slope      # 빠뜨린 C 를 T 에 회귀한 기울기
delta = res["long_beta_C"]                # 긴 회귀의 C 계수
rhs = res["long_beta_T"] + delta * gamma
print(f"\n모집단 값: short = 0.8,  long = 0,  delta = 1,  gamma = 0.8")
print(f"표본  gamma = {gamma:.6f},  delta = {delta:.6f}")
print(f"  long + delta*gamma = {res['long_beta_T']:.6f} + "
      f"{delta:.6f}*{gamma:.6f} = {rhs:.6f}")
print(f"  short              = {res['short_slope']:.6f}")
print(f"  두 변의 차         = {abs(res['short_slope'] - rhs):.2e}")

short_slope  =   0.7772   (p = 2.8e-37)
long_beta_T  =  -0.1000
long_beta_C  =   1.0915

모집단 값: short = 0.8,  long = 0,  delta = 1,  gamma = 0.8
표본  gamma = 0.803622,  delta = 1.091524
  long + delta*gamma = -0.099974 + 1.091524*0.803622 = 0.777199
  short              = 0.777199
  두 변의 차         = 0.00e+00

# Frisch-Waugh-Lovell 정리: T 와 Y 에서 각각 C 로 설명되는 몫을 걷어 낸 뒤
# 남은 잔차끼리 회귀하면, 긴 회귀의 T 계수와 똑같은 값이 나온다.
# "C 를 통제한다"는 말이 실제로 무엇을 하는 일인지 보여 주는 계산이다.
t_resid = t - stats.linregress(c, t).slope * c
y_resid = y - stats.linregress(c, y).slope * c
slope_partial = stats.linregress(t_resid, y_resid).slope
print(f"slope_partial = {slope_partial:.4f}")

# (1) 직교성과 (2) Y 를 씻어 내지 않아도 되는지 확인한다.
print(f"\ncorr(e_T, C) = {np.corrcoef(t_resid, c)[0, 1]:.2e}   (0 이어야 한다)")
print(f"e_Y 를 e_T 에 회귀  = {stats.linregress(t_resid, y_resid).slope:.6f}")
print(f"원래 Y 를 e_T 에 회귀 = {stats.linregress(t_resid, y).slope:.6f}")
print(f"긴 회귀의 beta_T    = {res['long_beta_T']:.6f}")
print(f"셋이 모두 같은가: "
      f"{np.allclose([stats.linregress(t_resid, y_resid).slope, stats.linregress(t_resid, y).slope], res['long_beta_T'])}")

slope_partial = -0.1000

corr(e_T, C) = 1.24e-17   (0 이어야 한다)
e_Y 를 e_T 에 회귀  = -0.099974
원래 Y 를 e_T 에 회귀 = -0.099974
긴 회귀의 beta_T    = -0.099974
셋이 모두 같은가: True

def simpson_paradox_demo(n=1000):
    """중증도가 치료 배정과 결과를 함께 좌우할 때 생기는 역설을 보인다.

    중증 환자가 치료를 더 많이 받고(0.7 대 0.3), 중증 자체가 결과를 크게
    낮춘다. 그래서 치료가 실제로는 결과를 5 만큼 올리는데도, 층을 나누지
    않고 보면 치료군의 결과가 더 나빠 보인다.
    """
    severity = np.random.binomial(1, 0.5, n)
    p_treat = np.where(severity == 1, 0.7, 0.3)
    treatment = np.random.binomial(1, p_treat)

    y = (50 - 20 * severity + 5 * treatment
         + np.random.normal(0, 5, n))

    # 층을 나누지 않은 순진한 평균처치효과
    ate_naive = y[treatment == 1].mean() - y[treatment == 0].mean()

    # 중증도로 층을 나눠 각 층에서 효과를 구한 뒤, 층의 크기로 가중해 합친다.
    ate_mild = (y[(treatment == 1) & (severity == 0)].mean()
                - y[(treatment == 0) & (severity == 0)].mean())
    ate_severe = (y[(treatment == 1) & (severity == 1)].mean()
                  - y[(treatment == 0) & (severity == 1)].mean())
    p_severe = severity.mean()
    ate_adjusted = (1 - p_severe) * ate_mild + p_severe * ate_severe

    return ate_naive, ate_mild, ate_severe, ate_adjusted


# 모의 한 번 돌린 결과
ate_naive, ate_mild, ate_severe, ate_adjusted = simpson_paradox_demo(1000)
print(f"ate_naive    = {ate_naive:>6.2f}   (이론 -3)")
print(f"ate_mild     = {ate_mild:>6.2f}   (이론 +5)")
print(f"ate_severe   = {ate_severe:>6.2f}   (이론 +5)")
print(f"ate_adjusted = {ate_adjusted:>6.2f}   (이론 +5)")

# n = 1000 에서 네 추정량이 얼마나 흔들리는지 재어, 위 한 번의 값을 판정한다.
reps = np.array([simpson_paradox_demo(1000) for _ in range(5000)])
theory = [-3.0, 5.0, 5.0, 5.0]
names = ["ate_naive", "ate_mild", "ate_severe", "ate_adjusted"]
print(f"\n5000 번 되풀이 (n = 1000)")
print(f"{'추정량':>14s} {'이론':>6s} {'모의 평균':>10s} {'표준편차':>9s} "
      f"{'위 한 번':>9s} {'z':>7s}")
once = [ate_naive, ate_mild, ate_severe, ate_adjusted]
for j, nm in enumerate(names):
    sd = reps[:, j].std(ddof=1)
    print(f"{nm:>14s} {theory[j]:6.1f} {reps[:, j].mean():+10.4f} {sd:9.4f} "
          f"{once[j]:+9.4f} {(once[j] - theory[j]) / sd:+7.3f}")

ate_naive    =  -2.62   (이론 -3)
ate_mild     =   5.22   (이론 +5)
ate_severe   =   6.07   (이론 +5)
ate_adjusted =   5.64   (이론 +5)

5000 번 되풀이 (n = 1000)
           추정량     이론      모의 평균      표준편차     위 한 번       z
     ate_naive   -3.0    -3.0109    0.6591   -2.6184  +0.579
      ate_mild    5.0    +4.9957    0.4848   +5.2238  +0.462
    ate_severe    5.0    +4.9927    0.4909   +6.0716  +2.183
  ate_adjusted    5.0    +4.9945    0.3385   +5.6400  +1.891

import numpy as np
from scipy import stats

rhos = [-0.9, -0.5, 0, 0.5, 0.9]
short_slopes = []

for rho in rhos:
    np.random.seed(42)
    n = 500
    tc = np.random.multivariate_normal([0, 0],
                                       [[1, rho], [rho, 1]], n)
    t, c = tc[:, 0], tc[:, 1]
    y = c + np.random.normal(0, 1, n)

    short_slope = stats.linregress(t, y).slope
    X = np.column_stack([np.ones(n), t, c])
    beta = np.linalg.lstsq(X, y, rcond=None)[0]

    short_slopes.append(short_slope)
    print(f"rho={rho:+.1f}: short={short_slope:.3f}, "
          f"long_T={beta[1]:.3f}")

import matplotlib.pyplot as plt
plt.plot(rhos, short_slopes, 'o-')
plt.xlabel('rho(T, C)')
plt.ylabel('Short regression slope')
plt.axhline(0, color='gray', linestyle='--')
plt.title('Omitted Variable Bias')
plt.show()

rho=-0.9: short=-0.877, long_T=0.131
rho=-0.5: short=-0.444, long_T=0.056
rho=+0.0: short=-0.021, long_T=0.008
rho=+0.5: short=0.454, long_T=-0.065
rho=+0.9: short=0.885, long_T=-0.140

import numpy as np

np.random.seed(42)
n = 1000
severity = np.random.binomial(1, 0.5, n)
treatment = np.random.binomial(1, 0.5, n)  # independent of severity

y = 50 - 20 * severity + 5 * treatment + np.random.normal(0, 5, n)

ate_naive = y[treatment == 1].mean() - y[treatment == 0].mean()
print(f"Naive ATE = {ate_naive:.2f} (true = 5)")   # 4.86

Naive ATE = 4.86 (true = 5)

import numpy as np

np.random.seed(42)
n = 1500
severity = np.random.choice([0, 1, 2], size=n, p=[1/3, 1/3, 1/3])
p_treat = np.where(severity == 0, 0.2,
                   np.where(severity == 1, 0.5, 0.8))
treatment = np.random.binomial(1, p_treat)

y = 50 - 15 * severity + 5 * treatment + np.random.normal(0, 5, n)

ate_naive = y[treatment == 1].mean() - y[treatment == 0].mean()

ates = []
for s in [0, 1, 2]:
    mask_t = (treatment == 1) & (severity == s)
    mask_c = (treatment == 0) & (severity == s)
    ate_s = y[mask_t].mean() - y[mask_c].mean()
    ates.append(ate_s)
    print(f"ATE (severity={s}): {ate_s:.2f}")

p_sev = [np.mean(severity == s) for s in [0, 1, 2]]
ate_adjusted = sum(p * a for p, a in zip(p_sev, ates))

print(f"Naive ATE:    {ate_naive:.2f}")
print(f"Adjusted ATE: {ate_adjusted:.2f}")
print(f"True effect:  5.00")

ATE (severity=0): 4.51
ATE (severity=1): 4.17
ATE (severity=2): 5.60
Naive ATE:    -7.28
Adjusted ATE: 4.76
True effect:  5.00