Simpson의 역설¶

dept = ['A','B','C','D','E','F']
m_app = [825,560,325,417,191,373]; m_adm = [512,353,120,138,53,22]
f_app = [108, 25,593,375,393,341]; f_adm = [ 89, 17,202,131, 94,24]

print(f"{'학과':>5}{'남 지원':>8}{'남 합격률':>10}{'여 지원':>8}"
      f"{'여 합격률':>10}{'차(여-남)':>11}")
for i, d in enumerate(dept):
    mr, fr = m_adm[i]/m_app[i], f_adm[i]/f_app[i]
    print(f"{d:>5}{m_app[i]:>8}{mr:>10.4f}{f_app[i]:>8}{fr:>10.4f}{fr-mr:>+11.4f}")
M, F = sum(m_adm)/sum(m_app), sum(f_adm)/sum(f_app)
print(f"\n  전체: 남 {M:.4f}  여 {F:.4f}   차 {F-M:+.4f}")

   학과    남 지원     남 합격률    여 지원     여 합격률     차(여-남)
    A     825    0.6206     108    0.8241    +0.2035
    B     560    0.6304      25    0.6800    +0.0496
    C     325    0.3692     593    0.3406    -0.0286
    D     417    0.3309     375    0.3493    +0.0184
    E     191    0.2775     393    0.2392    -0.0383
    F     373    0.0590     341    0.0704    +0.0114

  전체: 남 0.4452  여 0.3035   차 -0.1416

import numpy as np

rate   = [(m_adm[i]+f_adm[i])/(m_app[i]+f_app[i]) for i in range(6)]
fshare = [f_app[i]/(m_app[i]+f_app[i]) for i in range(6)]
print(f"{'학과':>5}{'전체 합격률':>12}{'여성 지원 비율':>14}")
for i, d in enumerate(dept):
    print(f"{d:>5}{rate[i]:>12.4f}{fshare[i]:>14.4f}")
print(f"\n  둘의 상관 = {np.corrcoef(rate, fshare)[0,1]:+.4f}")

   학과     전체 합격률      여성 지원 비율
    A      0.6442        0.1158
    B      0.6325        0.0427
    C      0.3508        0.6460
    D      0.3396        0.4735
    E      0.2517        0.6729
    F      0.0644        0.4776

  둘의 상관 = -0.7858

import numpy as np

a1, a2 = 0.80, 0.20      # 집단 A 의 층별 성공률
b1, b2 = 0.75, 0.15      # 집단 B — 두 층 모두 A 보다 낮다
print(f"{'w(A의 층1 비중)':>16}{'v(B의 층1 비중)':>16}"
      f"{'A 전체':>9}{'B 전체':>9}{'역설?':>7}")
for w, v in [(0.5,0.5), (0.2,0.8), (0.1,0.9), (0.05,0.95)]:
    A = w*a1 + (1-w)*a2
    B = v*b1 + (1-v)*b2
    print(f"{w:>16.2f}{v:>16.2f}{A:>9.4f}{B:>9.4f}"
          f"{'예' if A < B else '아니오':>7}")

  w(A의 층1 비중)   v(B의 층1 비중)    A 전체    B 전체   역설?
            0.50            0.50   0.5000   0.4500  아니오
            0.20            0.80   0.3200   0.6300      예
            0.10            0.90   0.2600   0.6900      예
            0.05            0.95   0.2300   0.7200      예

import numpy as np

rng = np.random.default_rng(3)
G, m = 5, 60
gx = np.array([1., 2., 3., 4., 5.])       # 집단 평균 x 가 오른다
gy = np.array([5., 4., 3., 2., 1.])       # 집단 평균 y 는 내린다
X, Y, lab = [], [], []
for i in range(G):
    xi = rng.normal(gx[i], 0.4, m)
    yi = gy[i] + 0.8*(xi - gx[i]) + rng.normal(0, 0.2, m)  # 내부는 +0.8
    X.append(xi); Y.append(yi); lab += [i]*m
X = np.concatenate(X); Y = np.concatenate(Y); lab = np.array(lab)

print(f"  전체 기울기 = {np.polyfit(X, Y, 1)[0]:+.4f}")
for i in range(G):
    print(f"    집단 {i+1} 기울기 = "
          f"{np.polyfit(X[lab==i], Y[lab==i], 1)[0]:+.4f}")

  전체 기울기 = -0.8429
    집단 1 기울기 = +0.7264
    집단 2 기울기 = +0.7380
    집단 3 기울기 = +0.8250
    집단 4 기울기 = +0.8940
    집단 5 기울기 = +0.7943

import numpy as np

rng = np.random.default_rng(4)
n = 200_000
A = rng.normal(0, 1, n)          # 질병 A
B = rng.normal(0, 1, n)          # 질병 B — A 와 독립
hosp = (A + B) > 2.0             # 둘 중 하나라도 심하면 입원
print(f"  전체       corr(A,B) = {np.corrcoef(A, B)[0,1]:+.4f}")
print(f"  입원자만   corr(A,B) = {np.corrcoef(A[hosp], B[hosp])[0,1]:+.4f}")
print(f"  (입원 비율 {hosp.mean():.4f})")

  전체       corr(A,B) = -0.0025
  입원자만   corr(A,B) = -0.7331

1. 층 변수가 처치보다 시간적으로 앞서는가?
     앞선다 -> 교란 가능성 -> 보정
     뒤따른다 -> 매개 또는 충돌 -> 보정 금지
2. 층 변수가 처치와 결과의 공통 결과인가?
     그렇다 -> 충돌자 -> 절대 보정 금지
3. 확신이 없으면 양쪽 결과를 모두 보고한다

import numpy as np

rng = np.random.default_rng(5)
REP = 500_000
cnt = 0
for _ in range(REP):
    a1, a2, b1, b2 = rng.uniform(0, 1, 4)     # 층별 성공률
    w, v = rng.uniform(0, 1, 2)               # 층 1 의 비중
    A = w*a1 + (1-w)*a2
    B = v*b1 + (1-v)*b2
    if (a1 > b1 and a2 > b2 and A < B) or (a1 < b1 and a2 < b2 and A > B):
        cnt += 1
print(f"  무작위 2x2x2 표에서 심슨 역설이 일어나는 비율: {cnt/REP:.4f}")

  무작위 2x2x2 표에서 심슨 역설이 일어나는 비율: 0.0243

import numpy as np
import pandas as pd
import statsmodels.api as sm

rows = []
for i, d in enumerate(dept):
    rows += [(d, 1, 1)]*m_adm[i] + [(d, 1, 0)]*(m_app[i]-m_adm[i])
    rows += [(d, 0, 1)]*f_adm[i] + [(d, 0, 0)]*(f_app[i]-f_adm[i])
df = pd.DataFrame(rows, columns=['dept', 'male', 'admit'])

m1 = sm.Logit(df.admit, sm.add_constant(df[['male']])).fit(disp=0)
D = pd.get_dummies(df.dept, drop_first=True).astype(float)
m2 = sm.Logit(df.admit,
              sm.add_constant(pd.concat([df[['male']], D], axis=1))).fit(disp=0)

print(f"  성별만       : male 계수 = {m1.params['male']:+.4f}  "
      f"오즈비 = {np.exp(m1.params['male']):.4f}")
print(f"  성별 + 학과  : male 계수 = {m2.params['male']:+.4f}  "
      f"오즈비 = {np.exp(m2.params['male']):.4f}")

  성별만       : male 계수 = +0.6104  오즈비 = 1.8411
  성별 + 학과  : male 계수 = -0.0999  오즈비 = 0.9050

ci = m2.conf_int().loc['male']
print(f"  95% CI(오즈비) = [{np.exp(ci[0]):.4f}, {np.exp(ci[1]):.4f}]")

1. 결과와 관련될 만한 변수를 나열한다 (도메인 지식)
2. 그중 처치/집단과도 관련된 것을 고른다 -> 교란 후보
3. 각 후보로 층화해 결과를 본다
     방향이 바뀌는가?        -> 심슨 역설
     크기가 크게 바뀌는가?   -> 교호작용 또는 교란
     거의 같은가?            -> 그 변수는 문제 아님
4. 시간 순서를 확인한다 (처치 전 변수만 보정 대상)
5. 충돌자를 걸러 낸다 (처치와 결과의 공통 결과는 보정 금지)

전체: 여성 합격률 30.4%, 남성 44.5% (차 -14.2%p)

학과별로 층화하면:
  6개 학과 중 4개에서 여성이 높음 (+20.4%p ~ -3.8%p)
  여성은 합격률이 낮은 학과에 더 많이 지원
  (학과 합격률과 여성 지원비율의 상관 -0.79)

학과 보정 후 오즈비 0.90 (95% CI 포함)
-> 전체 격차는 주로 학과 선택 차이로 설명된다.
   단, 학과 선택 자체가 무엇의 결과인지는 이 자료로 답할 수 없다.