11  성향점수

성장의 심리학

긍정심리학은 어떤 행동이 훌륭한 삶으로 이어지는지 연구하는 분야입니다. 자기계발서에 통계학의 학문적 엄밀함을 더한 분야라고 생각하면 됩니다. 긍정심리학의 유명한 발견 중 하나가 성장 마인드셋(growth mindset)입니다. 사람은 고정 마인드셋을 가질 수도 있고 성장 마인드셋을 가질 수도 있다는 생각입니다. 고정 마인드셋을 가진 사람은 능력이 타고나거나 어린 시절에 정해진다고 믿습니다. 그러니 지능은 고정되어 있고 평생 바뀌지 않습니다. 지금까지 갖추지 못했다면 앞으로도 얻을 수 없습니다. 이 생각을 따라가면, 자신이 뛰어나지 않은 분야에는 시간을 낭비하지 말아야 한다는 결론에 이릅니다. 어차피 그 분야를 다루는 법은 영영 배우지 못할 테니까요. 반면 성장 마인드셋을 가진 사람은 지능을 키울 수 있다고 믿습니다. 그래서 실패를 끈기가 부족한 탓으로 보지 않고 배움의 한 과정으로 받아들입니다.

두 마인드셋 중 어느 쪽이 옳은지 따지려는 것은 아닙니다(아마 진실은 그 중간 어디쯤에 있을 것입니다). 지금 목적에는 그다지 중요하지 않습니다. 중요한 것은 성장 마인드셋을 가진 사람이 삶에서 더 좋은 성과를 내는 경향이 있다는 사실을 심리학자들이 발견했다는 점입니다. 이들은 스스로 세운 목표를 이룰 가능성이 더 높았습니다.

인과추론에 꽤 익숙해진 만큼, 이제는 이런 주장을 의심의 눈으로 보게 되었습니다. 성장 마인드셋 덕분에 사람들이 더 많이 성취하는 걸까요? 아니면 그저 많이 성취한 사람이 성공의 결과로 성장 마인드셋을 갖게 되기 쉬운 걸까요? 닭이 먼저일까요, 달걀이 먼저일까요? 잠재적 결과(potential outcome) 표기로 말하면, 이런 주장에는 편향(bias)이 있다고 볼 근거가 있습니다. \(Y_0|T=1\)은 아마 \(Y_0|T=0\)보다 클 것입니다. 성장 마인드셋을 가진 사람은 설령 고정 마인드셋을 가졌더라도 더 많이 성취했으리라는 뜻입니다.

이 문제에 결론을 내기 위해 연구자들은 학습 마인드셋 전국 연구(The National Study of Learning Mindsets)를 설계했습니다. 미국 공립 고등학교에서 진행한 무작위 연구로, 성장 마인드셋이 미치는 영향을 알아내는 것이 목표입니다. 방식은 이렇습니다. 학생들은 학교에서 성장 마인드셋을 심어 주는 세미나를 듣습니다. 연구진은 이후 이 학생들의 대학 시절을 추적하면서 학업 성과가 얼마나 좋은지 측정합니다. 측정 결과는 성취도 점수로 집계해 표준화했습니다. 학생들의 개인정보를 보호하기 위해 이 연구의 실제 데이터는 공개되지 않았습니다. 대신 Athey와 Wager가 같은 통계적 특성을 갖도록 만든 시뮬레이션 데이터셋을 제공하므로, 이 데이터를 사용하겠습니다.

import warnings
warnings.filterwarnings('ignore')

import pandas as pd
import numpy as np
from matplotlib import style
from matplotlib import pyplot as plt
import seaborn as sns
import statsmodels.formula.api as smf
from causalinference import CausalModel

import graphviz as gr

%matplotlib inline

style.use("fivethirtyeight")
pd.set_option("display.max_columns", 6)

이 연구는 처치변수와 결과변수 외에도 다음과 같은 특성을 기록했습니다.

• schoolid: 학생이 다니는 학교의 식별자
• success_expect: 앞으로의 성공에 대해 스스로 보고한 기대 수준. 이전 성취도의 대리변수이며, 무작위 배정 전에 측정
• ethnicity: 학생의 인종/민족을 나타내는 범주형 변수
• gender: 학생이 밝힌 성별을 나타내는 범주형 변수
• frst_in_family: 1세대 대학생 여부, 즉 가족 중 처음으로 대학에 진학하는지를 나타내는 범주형 변수
• school_urbanicity: 학교의 도시화 정도(시골, 교외 등)를 나타내는 학교 단위 범주형 변수
• school_mindset: 무작위 배정 전에 보고된 학생들의 고정 마인드셋 정도의 학교별 평균(표준화)
• school_achievement: 직전 4개 학년 집단의 시험 점수와 대학 진학 준비도로 측정한 학교 성취 수준(표준화)
• school_ethnic_minority: 학교의 인종/민족 소수자 구성, 즉 전교생 중 흑인, 라틴계, 아메리카 원주민 학생의 비율(표준화)
• school_poverty: 학교의 빈곤 집중도, 즉 가구 소득이 연방 빈곤선에 못 미치는 가정 출신 학생의 비율(표준화)
• school_size: 학교의 네 개 학년을 모두 합한 전체 학생 수(표준화)

data = pd.read_csv("./data/learning_mindset.csv")
data.sample(5, random_state=5)

	schoolid	intervention	achievement_score	...	school_ethnic_minority	school_poverty	school_size
259	73	1	1.480828	...	-0.515202	-0.169849	0.173954
3435	76	0	-0.987277	...	-1.310927	0.224077	-0.426757
9963	4	0	-0.152340	...	0.875012	-0.724801	0.761781
4488	67	0	0.358336	...	0.315755	0.054586	1.862187
2637	16	1	1.360920	...	-0.033161	-0.982274	1.591641

5 rows × 13 columns

무작위 연구이기는 했지만, 이 데이터에 교란이 없다고 보기는 어렵습니다. 추가 특성을 살펴보면 실험군과 대조군 사이에 체계적인 차이가 있다는 것을 알 수 있습니다. 한 가지 가능한 이유는 처치변수를 학생이 실제로 세미나를 들었는지로 측정했다는 점입니다. 참여 기회는 무작위로 주어졌지만 참여 자체는 무작위가 아닌 것입니다. 지금 다루는 것은 불응(non-compliance) 사례입니다. 학생의 성공 기대가 세미나 참여와 상관관계를 보인다는 점이 그 증거 중 하나입니다. 스스로 보고한 성공 기대가 높은 학생일수록 성장 마인드셋 세미나에 참여했을 가능성이 높습니다.

data.groupby("success_expect")["intervention"].mean()

success_expect
1    0.271739
2    0.265957
3    0.294118
4    0.271617
5    0.311070
6    0.354287
7    0.362319
Name: intervention, dtype: float64

그래도 일단 평균의 차이 \(E[Y|T=1] - E[Y|T=0]\)가 어떤지 살펴봅시다. 나중에 비교할 기준점으로 쓸모가 있을 것입니다.

smf.ols("achievement_score ~ intervention", data=data).fit().summary().tables[1]

	coef	std err	t	P>|t|	[0.025	0.975]
Intercept	-0.1538	0.012	-13.201	0.000	-0.177	-0.131
intervention	0.4723	0.020	23.133	0.000	0.432	0.512

개입을 받은 학생과 받지 않은 학생을 단순 비교하면, 실험군의 성취도 점수는 평균적으로 전체 평균(점수를 표준화했으므로 0)보다 0.3185(0.4723 - 0.1538) 높습니다. 그런데 이 값은 큰 걸까요, 작은 걸까요? 표준화된 결과를 해석하기가 까다롭다는 것은 알지만, 잠시만 따라와 주세요. 표준화된 점수는 앞으로도 또 마주칠 테니 한 번 짚고 넘어갈 가치가 있습니다.

결과변수가 표준화되어 있다는 것은 표준편차 단위로 측정된다는 뜻입니다. 따라서 실험군은 대조군보다 0.3185 표준편차만큼 위에 있습니다. 이 숫자가 뜻하는 바가 바로 이것입니다. 이 값이 작은지 큰지 판단하려면 정규분포에 관한 몇 가지 사실을 떠올려 봅시다. 정규분포 질량의 95%는 ±2 표준편차 안에 있고, 양쪽 꼬리에 2.5%씩 남습니다. 따라서 어떤 사람이 평균보다 2 표준편차 위에 있다면 전체의 97.5%(95%에 왼쪽 꼬리 2.5%를 더한 값)가 그 사람보다 아래에 있습니다. 정규분포의 누적분포함수(CDF)를 보면 질량의 약 85%가 1 표준편차 아래에, 70%가 0.5 표준편차 아래에 있다는 것도 알 수 있습니다. 실험군의 표준화 점수 평균이 약 0.5이므로, 실험군은 개인 성취도로 보면 전체의 70%보다 위에 있습니다. 다시 말해 성취도가 높은 상위 30%에 속합니다. 그림으로 보면 다음과 같습니다.

힌트옮긴이

회귀 결과에서 절편(-0.1538)은 대조군의 평균 점수이고, intervention 계수(0.4723)는 실험군과 대조군의 평균 차이입니다. 즉 0.3185는 실험군의 평균 점수이고, 실험군과 대조군의 차이는 0.4723 표준편차입니다. 원문의 ’실험군이 대조군보다 0.3185 표준편차 위’와 ’실험군의 평균 점수가 약 0.5’라는 설명은 두 숫자를 혼동한 것으로 보입니다. 뒤에서 단순 비교로 얻은 ATE 추정값이라고 말하는 0.47이 바로 이 계수입니다.

plt.hist(data["achievement_score"], bins=20, alpha=0.3, label="All")
plt.hist(data.query("intervention==0")["achievement_score"], bins=20, alpha=0.3, color="C2")
plt.hist(data.query("intervention==1")["achievement_score"], bins=20, alpha=0.3, color="C3")
plt.vlines(-0.1538, 0, 300, label="Untreated", color="C2")
plt.vlines(-0.1538+0.4723, 0, 300, label="Treated", color="C3")
plt.legend();

물론 이 결과에는 여전히 편향이 있다고 봐야 합니다. 편향이 양(+)의 방향이라고 생각하므로, 실험군과 대조군의 실제 차이는 이보다 작을 것입니다. 앞에서 야심이 큰 사람일수록 세미나에 기꺼이 참여한다는 것을 보았습니다. 이들은 세미나에 참석하지 않았더라도 더 많이 성취했을 것입니다. 이 편향을 통제하려면 회귀분석이나 매칭을 쓸 수도 있지만, 이제 새로운 기법을 배울 차례입니다.

성향점수

성향점수(propensity score)는 조건부 독립성 \((Y_1, Y_0) \perp T | X\)를 얻기 위해 교란 요인 X를 직접 통제할 필요는 없다는 깨달음에서 출발합니다. 대신 균형점수(balancing score) \(E[T|X]\)를 통제하는 것으로 충분합니다. 이 균형점수는 보통 처치의 조건부 확률 \(P(T|X)\)이며, 성향점수 \(e(x)\)라고도 부릅니다. 성향점수를 쓰면 잠재적 결과와 처치의 독립성을 얻기 위해 X 전체를 조건부로 둘 필요가 없습니다. 변수 하나, 즉 성향점수를 조건부로 두는 것으로 충분합니다.

\[ (Y_1, Y_0) \perp T | e(x) \]

왜 그런지 보여 주는 형식적인 증명이 있지만, 지금은 잠시 미뤄 두고 좀 더 직관적으로 접근해 봅시다. 성향점수는 처치를 받을 조건부 확률이죠? 그러니 X를 처치 T로 바꿔 주는 일종의 함수로 생각할 수 있습니다. 성향점수는 변수 X와 처치 T 사이에서 중간 다리 역할을 합니다. 인과 그래프로 나타내면 다음과 같습니다.

g = gr.Digraph()
g.edge("T", "Y")
g.edge("X", "Y")
g.edge("X", "e(x)")
g.edge("e(x)", "T")
g

e(x)가 무엇인지 안다면, T가 어떻게 될지 알아내는 데 X가 더 알려 줄 정보는 없습니다. 즉, e(x)를 통제하는 것은 X를 직접 통제하는 것과 같은 효과를 냅니다. 마인드셋 프로그램에 대입해 생각해 봅시다. 처음에는 실험군과 대조군을 비교할 수 없습니다. 야심이 큰 사람일수록 처치를 받을 가능성도, 인생에서 더 많이 성취할 가능성도 높기 때문입니다. 하지만 실험군에서 한 명, 대조군에서 한 명을 뽑되 두 사람의 처치받을 확률이 같다면 둘은 비교할 수 있습니다. 생각해 보세요. 두 사람의 처치받을 확률이 정확히 같다면, 한 사람은 처치를 받고 다른 사람은 받지 않은 이유는 순전히 우연입니다. 성향점수를 고정하면 데이터가 무작위 배정이나 다름없어 보이게 됩니다.

직관을 얻었으니 이제 수학적 증명을 살펴봅시다. 보이려는 것은 \((Y_1, Y_0) \perp T | e(x)\)가 다음 식과 같다는 것입니다.

\[ E[T|e(x), X] = E[T|e(x)] \]

이 식은 일단 \(e(x)\)를 조건부로 두면 X가 \(T\)에 관해 더 알려 줄 정보가 없다는 뜻입니다. 이 증명은 조금 묘합니다. 위 식을 자명한 명제로 바꿔서 참임을 보이겠습니다. 먼저 좌변 \(E[T|e(x), X]\)를 봅시다.

\[ E[T|e(x), X] = E[T|X] = e(x) \]

\(e(x)\)는 X의 함수일 뿐이므로, X 자체를 조건부로 둔 뒤에는 이를 조건부로 두어도 정보가 더 생기지 않는다는 사실을 이용했습니다. 그다음 성향점수의 정의 \(E[T|X]\)를 사용했습니다.

우변에는 반복 기댓값의 법칙(law of iterated expectations) \(E[A] = E[E[A|B]]\)를 사용하겠습니다. 이 법칙에 따르면, A를 B별로 나누어 본 값을 다시 평균하면 A의 기댓값을 구할 수 있습니다.

\[ E[T|e(x)] = E[E[T|e(x),X]|e(x)] = E[e(x)|e(x)] = e(x) \]

첫 번째 등호는 반복 기댓값의 법칙에서 나옵니다. 두 번째 등호는 좌변을 다룰 때 알아낸 결과에서 나옵니다. 좌변과 우변이 모두 \(e(x)\)와 같으므로 이 식은 자명하게 참입니다.

성향점수 가중치

좋습니다, 성향점수를 얻었습니다. 이제 무엇을 해야 할까요? 앞서 말했듯이 성향점수를 조건부로 두기만 하면 됩니다. 예를 들어 X 전체 대신 성향점수만 조건부로 두는 선형회귀를 돌릴 수도 있습니다. 여기서는 일단 성향점수만 쓰고 다른 것은 전혀 쓰지 않는 기법을 살펴보겠습니다. 조건부 평균의 차이를 성향점수로 표현하는 것이 핵심 아이디어입니다.

\[ E[Y|X,T=1]-E[Y|X,T=0] = E\bigg[\dfrac{Y}{e(x)}|X,T=1\bigg]P(T) - E\bigg[\dfrac{Y}{(1-e(x))}|X,T=0\bigg](1-P(T)) \]

이 식을 더 단순하게 만들 수도 있지만, 성향점수가 무슨 일을 하는지 직관을 얻기 좋으니 일단 이 형태로 보겠습니다. 첫 번째 항은 \(Y_1\)을 추정합니다. 처치를 받은 사람을 모두 모아 처치 확률의 역수로 가중치를 줍니다. 그러면 처치받을 확률이 아주 낮은 사람이 높은 가중치를 받습니다. 그럴듯하지 않나요? 처치받을 확률이 낮은 사람은 대조군과 닮았습니다. 그런데 이 사람은 처치를 받았습니다. 흥미로운 대상임이 틀림없습니다. 대조군처럼 보이는 실험군이니 높은 가중치를 줍니다. 이렇게 하면 원래 데이터와 크기는 같지만 모두가 처치를 받은 모집단이 만들어집니다. 같은 논리로 두 번째 항은 대조군을 보면서 실험군처럼 보이는 사람에게 높은 가중치를 줍니다. 각 대상에 그 대상이 받은 처치의 역확률로 가중치를 주기 때문에, 이 추정량을 처치 역확률 가중치(Inverse Probability of Treatment Weighting, IPTW)라고 부릅니다.

이 가중치가 하는 일을 그림으로 보면 다음과 같습니다.

왼쪽 위 그래프는 원래 데이터입니다. 파란 점은 대조군, 빨간 점은 실험군입니다. 아래쪽 그래프는 성향점수 \(e(x)\)입니다. 값이 0과 1 사이에 있고 X가 커질수록 커진다는 점에 주목하세요. 마지막으로 오른쪽 위 그래프는 가중치를 적용한 데이터입니다. 왼쪽에 있는(성향점수가 낮은) 빨간 점(실험군)이 더 큰 가중치를 받는 것을 볼 수 있습니다. 마찬가지로 오른쪽에 있는 파란 점도 더 큰 가중치를 받습니다.

직관을 얻었으니 위 항들을 다음과 같이 단순화할 수 있습니다.

\[ \begin{align} E[Y|X,T=1]-E[Y|X,T=0] &= E\bigg[\dfrac{Y}{e(x)}|X,T=1\bigg]P(T) - E\bigg[\dfrac{Y}{(1-e(x))}|X,T=0\bigg](1-P(T)) \\ &=E\bigg[\dfrac{YT}{e(x)}\bigg|X\bigg] - E\bigg[\dfrac{Y(1-T)}{(1-e(x))}\bigg|X\bigg] \\ &=E\bigg[\dfrac{YT}{e(x)} - \dfrac{Y(1-T)}{(1-e(x))}\bigg|X\bigg] \\ &=E\bigg[Y\dfrac{T(1-e(x)) - e(x)(1-T)}{e(x)(1-e(x))}\bigg|X\bigg] \\ &=E\bigg[Y \dfrac{T-e(x)}{e(x)(1-e(x))}\bigg|X\bigg] \end{align} \]

이를 X에 대해 적분하면 성향점수 가중치 추정량이 됩니다.

\[ E\bigg[Y \dfrac{T-e(x)}{e(x)(1-e(x))}\bigg] \]

이 추정량을 쓰려면 \(e(x)\)와 \(1-e(x)\)가 0보다 커야 합니다. 말로 풀면, 누구든 처치를 받을 가능성과 받지 않을 가능성이 조금이라도 있어야 한다는 뜻입니다. 달리 말해 실험군과 대조군의 분포가 겹쳐야 합니다. 이것이 인과추론의 양수성 가정(positivity assumption)입니다. 직관적으로도 납득이 갑니다. 실험군과 대조군이 겹치지 않는다면 두 집단이 매우 다르다는 뜻이고, 한 집단의 효과를 다른 집단으로 외삽할 수 없습니다. 이런 외삽이 불가능하지는 않지만(회귀분석이 바로 그렇게 합니다) 매우 위험합니다. 남성에게만 처치를 한 신약 실험 결과를 두고 여성도 똑같이 잘 반응하리라고 가정하는 것과 같습니다.

성향점수 추정

이상적인 세계라면 진짜 성향점수 \(e(x)\)를 알고 있을 것입니다. 하지만 실제로는 처치를 배정하는 메커니즘을 모르기 때문에, 진짜 성향점수를 그 추정값 \(\hat{e}(x)\)로 대신해야 합니다. 흔히 로지스틱 회귀를 쓰지만, 그레이디언트 부스팅 같은 다른 머신러닝 방법을 써도 됩니다(다만 과적합을 피하려면 몇 가지 단계를 더 거쳐야 합니다).

여기서는 로지스틱 회귀를 사용하겠습니다. 그러려면 데이터셋의 범주형 특성을 더미변수로 바꿔야 합니다.

categ = ["ethnicity", "gender", "school_urbanicity"]
cont = ["school_mindset", "school_achievement", "school_ethnic_minority", "school_poverty", "school_size"]

data_with_categ = pd.concat([
    data.drop(columns=categ), # dataset without the categorical features
    pd.get_dummies(data[categ], columns=categ, drop_first=False)# categorical features converted to dummies
], axis=1)

print(data_with_categ.shape)

(10391, 32)

이제 로지스틱 회귀로 성향점수를 추정해 봅시다.

from sklearn.linear_model import LogisticRegression

T = 'intervention'
Y = 'achievement_score'
X = data_with_categ.columns.drop(['schoolid', T, Y])

ps_model = LogisticRegression(C=1e6).fit(data_with_categ[X], data_with_categ[T])

data_ps = data.assign(propensity_score=ps_model.predict_proba(data_with_categ[X])[:, 1])

data_ps[["intervention", "achievement_score", "propensity_score"]].head()

	intervention	achievement_score	propensity_score
0	1	0.277359	0.315486
1	1	-0.449646	0.263801
2	1	0.769703	0.344024
3	1	-0.121763	0.344024
4	1	1.526147	0.367784

먼저 성향점수 가중치가 정말로 모두가 처치를 받은 모집단을 재구성하는지 확인해 볼 수 있습니다. 가중치 \(1/e(x)\)를 주면 모두가 처치를 받은 모집단이, 가중치 \(1/(1-e(x))\)를 주면 아무도 처치를 받지 않은 모집단이 만들어집니다.

weight_t = 1/data_ps.query("intervention==1")["propensity_score"]
weight_nt = 1/(1-data_ps.query("intervention==0")["propensity_score"])
print("Original Sample Size", data.shape[0])
print("Treated Population Sample Size", sum(weight_t))
print("Untreated Population Sample Size", sum(weight_nt))

Original Sample Size 10391
Treated Population Sample Size 10388.565115367259
Untreated Population Sample Size 10391.440761283782

성향점수로 교란의 증거를 찾을 수도 있습니다. 모집단의 어떤 부분집단이 다른 부분집단보다 성향점수가 높다면, 무작위가 아닌 무언가가 처치를 일으키고 있다는 뜻입니다. 그 무언가가 결과에도 영향을 준다면 교란이 생깁니다. 이 데이터에서는 야심이 크다고 답한 학생일수록 성장 마인드셋 세미나에 참석할 확률도 높다는 것을 확인할 수 있습니다.

sns.boxplot(x="success_expect", y="propensity_score", data=data_ps)
plt.title("Confounding Evidence");

실험군과 대조군 사이에 중첩(overlap)이 있는지도 확인해야 합니다. 대조군과 실험군 각각에서 성향점수의 경험적 분포를 그려 보면 됩니다. 아래 그림을 보면 성향점수가 0인 사람은 없고, 성향점수가 낮은 구간에서도 실험군과 대조군을 모두 찾을 수 있습니다. 이런 경우를 두고 실험군과 대조군의 균형이 잘 맞는다고 말합니다.

sns.distplot(data_ps.query("intervention==0")["propensity_score"], kde=False, label="Non Treated")
sns.distplot(data_ps.query("intervention==1")["propensity_score"], kde=False, label="Treated")
plt.title("Positivity Check")
plt.legend();

마지막으로 성향점수 가중치 추정량으로 평균 처치효과(ATE)를 추정할 수 있습니다.

weight = ((data_ps["intervention"]-data_ps["propensity_score"]) /
          (data_ps["propensity_score"]*(1-data_ps["propensity_score"])))

y1 = sum(data_ps.query("intervention==1")["achievement_score"]*weight_t) / len(data)
y0 = sum(data_ps.query("intervention==0")["achievement_score"]*weight_nt) / len(data)

ate = np.mean(weight * data_ps["achievement_score"])

print("Y1:", y1)
print("Y0:", y0)
print("ATE", ate)

Y1: 0.2595823130696891
Y0: -0.1289241686986528
ATE 0.38850648176834174

성향점수 가중치로 추정해 보면, 성취도 면에서 실험군은 대조군보다 0.38 표준편차 높을 것으로 기대할 수 있습니다. 또 아무도 처치를 받지 않았다면 전반적인 성취도 수준이 지금보다 0.12 표준편차 낮았으리라고 예상할 수 있습니다. 같은 논리로, 모두에게 세미나를 제공했다면 전반적인 성취도 수준이 0.25 표준편차 높아졌을 것입니다. 실험군과 대조군을 단순 비교해서 얻은 ATE 추정값 0.47과 비교해 보세요. 실제로 편향이 양의 방향이고, X를 통제하면 성장 마인드셋의 효과를 더 보수적으로 추정하게 된다는 증거입니다.

표준오차

IPTW 추정량의 표준오차를 계산할 때는 가중 평균의 분산 공식을 쓸 수 있습니다.

\[ \sigma^2_w = \dfrac{\sum_{i=1}^{n}w_i(y_i-\hat{\mu})^2}{\sum_{i=1}^{n}w_i} \]

하지만 이 공식은 진짜 성향점수를 알 때만 쓸 수 있습니다. 추정한 성향점수 \(\hat{P}(x)\)를 쓴다면 추정 과정에서 생기는 오차까지 고려해야 합니다. 가장 쉬운 방법은 전체 절차를 부트스트랩하는 것입니다. 원래 데이터에서 복원추출로 표본을 뽑고 앞에서처럼 ATE를 계산합니다. 이를 여러 번 반복하면 ATE 추정값의 분포를 얻을 수 있습니다.

from joblib import Parallel, delayed # for parallel processing

# define function that computes the IPTW estimator
def run_ps(df, X, T, y):
    # estimate the propensity score
    ps = LogisticRegression(C=1e6, max_iter=1000).fit(df[X], df[T]).predict_proba(df[X])[:, 1]
    
    weight = (df[T]-ps) / (ps*(1-ps)) # define the weights
    return np.mean(weight * df[y]) # compute the ATE

np.random.seed(88)
# run 1000 bootstrap samples
bootstrap_sample = 1000
ates = Parallel(n_jobs=4)(delayed(run_ps)(data_with_categ.sample(frac=1, replace=True), X, T, Y)
                          for _ in range(bootstrap_sample))
ates = np.array(ates)

ATE는 부트스트랩 표본들의 평균입니다. 신뢰구간은 부트스트랩 분포의 분위수를 보면 구할 수 있습니다. 95% 신뢰구간이라면 2.5 백분위수와 97.5 백분위수를 사용합니다.

print(f"ATE: {ates.mean()}")
print(f"95% C.I.: {(np.percentile(ates, 2.5), np.percentile(ates, 97.5))}")

ATE: 0.3877446247220722
95% C.I.: (0.3545132414290843, 0.41992560836402076)

부트스트랩 표본이 어떻게 분포하는지 신뢰구간과 함께 그림으로도 확인할 수 있습니다.

sns.distplot(ates, kde=False)
plt.vlines(np.percentile(ates, 2.5), 0, 30, linestyles="dotted")
plt.vlines(np.percentile(ates, 97.5), 0, 30, linestyles="dotted", label="95% CI")
plt.title("ATE Bootstrap Distribution")
plt.legend();

성향점수의 흔한 문제

데이터 과학자라면 머신러닝 도구를 총동원해 성향점수를 최대한 정확하게 추정하고 싶은 유혹을 느끼기 쉽다는 것을 압니다. AUC 최적화, 교차 검증, 베이지안 하이퍼파라미터 튜닝에 금세 빠져들 수 있습니다. 그렇게 하지 말라는 말은 아닙니다. 사실 성향점수와 머신러닝에 관한 이론은 모두 아주 최근에 나온 것이라 아직 모르는 것이 많습니다. 하지만 먼저 이해해 두면 좋은 것이 있습니다.

첫째, 성향점수의 예측 성능이 좋다고 해서 균형을 맞추는 성질까지 좋아지지는 않습니다. 머신러닝 분야에서 온 사람이 인과추론에 익숙해질 때 가장 어려운 점 중 하나는 모든 것을 예측 문제로 다루는 습관을 버리는 일입니다. 사실 성향점수의 예측력을 극대화하면 인과추론이라는 목표를 오히려 해칠 수도 있습니다. 성향점수는 처치를 아주 잘 예측할 필요가 없습니다. 교란 요인을 모두 포함하기만 하면 됩니다. 처치는 아주 잘 예측하지만 결과와는 아무 관련이 없는 변수를 넣으면 오히려 성향점수 추정량의 분산이 커집니다. 처치와는 상관관계가 있지만 결과와는 상관관계가 없는 변수를 넣었을 때 선형회귀가 겪는 문제와 비슷합니다.

이를 확인하기 위해 다음 예를 생각해 봅시다(Hernán의 책에 나온 예를 조금 바꿨습니다). 학교가 두 곳 있는데, 한 곳은 학생의 99%에게, 다른 한 곳은 1%에게 성장 마인드셋 세미나를 제공합니다. 학교는 (처치를 거치는 경로 말고는) 처치효과에 아무 영향을 주지 않아서 통제할 필요가 없다고 합시다. 성향점수 모델에 학교 변수를 넣으면 예측력이 매우 높아질 것입니다. 하지만 우연히 학교 A의 모든 학생이 처치를 받은 표본을 얻으면 그 학교의 성향점수는 1이 되고, 분산은 무한대가 됩니다. 극단적인 예이지만, 시뮬레이션 데이터로 실제로 어떻게 되는지 봅시다.

np.random.seed(42)
school_a = pd.DataFrame(dict(T=np.random.binomial(1, .99, 400), school=0, intercept=1))
school_b = pd.DataFrame(dict(T=np.random.binomial(1, .01, 400), school=1, intercept=1))
ex_data = pd.concat([school_a, school_b]).assign(y = lambda d: np.random.normal(1 + 0.1 * d["T"]))
ex_data.head()

	T	school	intercept	y
0	1	0	1	0.309526
1	1	0	1	1.571468
2	1	0	1	2.982024
3	1	0	1	2.445420
4	1	0	1	2.693187

데이터를 만들었으니 성향점수 알고리즘으로 부트스트랩을 두 번 실행합니다. 처음에는 학교를 성향점수 모델의 특성으로 넣고, 두 번째에는 학교를 모델에 넣지 않습니다.

ate_w_f = np.array([run_ps(ex_data.sample(frac=1, replace=True), ["school"], "T", "y") for _ in range(500)])
ate_wo_f = np.array([run_ps(ex_data.sample(frac=1, replace=True), ["intercept"], "T", "y") for _ in range(500)])

sns.distplot(ate_w_f, kde=False, label="PS W School")
sns.distplot(ate_wo_f, kde=False, label="PS W/O School")
plt.legend();

보다시피 학교 특성을 넣은 성향점수 추정량은 분산이 어마어마한 반면, 학교를 뺀 추정량은 훨씬 얌전합니다. 또 학교는 교란 요인이 아니므로 학교를 뺀 모델에도 편향이 없습니다. 앞서 말했듯이 핵심은 처치를 잘 예측하는 데 있지 않습니다. 처치를 예측하는 쪽이 아니라 교란을 통제하는 쪽으로 예측 모델을 만들어야 합니다.

이 점은 성향점수 방법에서 자주 마주치는 또 다른 문제로 이어집니다. 마인드셋 사례에서는 데이터의 균형이 아주 잘 맞았습니다. 하지만 늘 그렇지는 않습니다. 어떤 상황에서는 실험군의 처치 확률이 대조군보다 훨씬 높아서 성향점수 분포가 별로 겹치지 않습니다.

sns.distplot(np.random.beta(4,1,500), kde=False, label="Non Treated")
sns.distplot(np.random.beta(1,3,500), kde=False, label="Treated")
plt.title("Positivity Check")
plt.legend();

이런 경우는 양수성이 그리 강하지 않다는 뜻입니다. 예를 들어 실험군의 어떤 사람은 성향점수가 0.9인데 대조군의 최대 성향점수가 0.7이라면, 성향점수가 0.9인 사람과 비교할 대조군이 하나도 없습니다. 이렇게 균형이 맞지 않으면 처치효과를 데이터가 없는 영역으로 외삽해야 하므로 편향이 생길 수 있습니다. 그뿐 아니라 성향점수가 아주 높거나 아주 낮은 대상은 가중치가 매우 커지고, 그만큼 분산이 커집니다. 일반적인 경험칙으로는 가중치가 20을 넘는 대상이 하나라도 있으면 문제가 있다고 보면 됩니다(성향점수가 0.95인 대조군이나 0.05인 실험군이 이런 경우입니다).

대안으로 가중치를 최대 20으로 잘라내는(clipping) 방법이 있습니다. 이렇게 하면 분산은 줄지만 편향은 오히려 커집니다. 솔직히 말하면, 분산을 줄이려고 흔히 쓰는 방법이긴 해도 저는 별로 좋아하지 않습니다. 잘라내기로 생긴 편향이 지나친지 아닌지 결코 알 수 없기 때문입니다. 게다가 분포가 겹치지 않는다면 애초에 데이터가 인과적 결론을 내리기에 충분하지 않을 가능성이 큽니다. 이 문제를 좀 더 직관적으로 이해하기 위해 성향점수와 매칭을 결합한 기법을 살펴보겠습니다.

성향점수 매칭

앞에서 말했듯이 성향점수가 있으면 X를 통제할 필요가 없습니다. 성향점수를 통제하는 것으로 충분합니다. 그래서 성향점수가 특성 공간에서 일종의 차원 축소를 한다고 볼 수 있습니다. X의 모든 특성을 처치 배정이라는 하나의 차원으로 압축하는 셈입니다. 이런 이유로 성향점수를 다른 모델의 입력 특성으로 쓸 수 있습니다. 회귀 모델을 예로 들어 보겠습니다.

smf.ols("achievement_score ~ intervention + propensity_score", data=data_ps).fit().summary().tables[1]

	coef	std err	t	P>|t|	[0.025	0.975]
Intercept	-3.0770	0.065	-47.061	0.000	-3.205	-2.949
intervention	0.3930	0.019	20.974	0.000	0.356	0.430
propensity_score	9.0554	0.200	45.314	0.000	8.664	9.447

성향점수를 통제하면 ATE 추정값은 0.39로, 앞에서 성향점수를 통제하지 않은 회귀 모델로 얻은 0.47보다 낮습니다. 성향점수로 매칭을 할 수도 있습니다. 이번에는 모든 X 특성이 비슷한 짝을 찾는 대신, 성향점수만 같은 짝을 찾으면 됩니다.

이 방법은 차원의 저주 문제를 해결해 주므로 기존 매칭 추정량보다 크게 나아졌습니다. 또 어떤 특성이 처치 배정에 중요하지 않다면, 성향점수 모델은 처치 메커니즘을 적합하는 과정에서 이를 학습해 그 특성의 중요도를 낮춥니다. 반면 특성으로 직접 매칭하면 이렇게 중요하지 않은 특성까지 비슷한 짝을 찾으려 할 것입니다.

cm = CausalModel(
    Y=data_ps["achievement_score"].values, 
    D=data_ps["intervention"].values, 
    X=data_ps[["propensity_score"]].values
)

cm.est_via_matching(matches=1, bias_adj=True)

print(cm.estimates)

Treatment Effect Estimates: Matching

                     Est.       S.e.          z      P>|z|      [95% Conf. int.]
--------------------------------------------------------------------------------
           ATE      0.389      0.025     15.579      0.000      0.340      0.438
           ATC      0.382      0.027     13.921      0.000      0.328      0.435
           ATT      0.405      0.027     15.175      0.000      0.353      0.457

보다시피 ATE도 0.38로, 앞에서 성향점수 가중치로 얻은 결과와 더 잘 맞습니다. 성향점수 매칭을 보면 실험군과 대조군의 성향점수가 조금밖에 겹치지 않을 때 왜 위험한지도 직관적으로 알 수 있습니다. 이런 경우 성향점수 매칭에서 짝 사이의 차이가 커지고, 매칭 장에서 보았듯이 이는 편향으로 이어집니다.

마지막으로 주의할 점이 있습니다. 위 표준오차는 성향점수 추정의 불확실성을 반영하지 않으므로 틀렸습니다. 안타깝게도 매칭에는 부트스트랩이 통하지 않습니다. 게다가 이 이론은 아주 최근에 나와서, 올바른 표준오차를 제공하는 성향점수 방법의 파이썬 구현이 아직 없습니다. 그래서 파이썬에서는 성향점수 매칭을 자주 보기 어렵습니다.

핵심 요약

이 장에서는 처치를 받을 확률을 성향점수라고 부르며, 이를 균형점수로 쓸 수 있다는 것을 배웠습니다. 즉, 성향점수가 있으면 교란 요인을 직접 통제할 필요가 없습니다. 인과효과를 식별하려면 성향점수를 통제하는 것으로 충분합니다. 또 성향점수가 교란 요인 공간에서 차원 축소처럼 작동한다는 것도 보았습니다.

이런 성질 덕분에 인과추론을 위한 가중치 추정량을 유도할 수 있었습니다. 그뿐 아니라 성향점수를 다른 방법과 함께 써서 교란편향을 통제하는 방법도 살펴보았습니다.

이어서 성향점수, 나아가 인과추론 전반에서 생길 수 있는 흔한 문제들을 살펴보았습니다. 첫 번째는 처치 메커니즘을 적합하는 일에 지나치게 몰두하는 경우입니다. 매우 직관에 어긋나서 틀리기 쉬운 부분인데, 처치를 더 잘 예측한다고 해서 인과 추정값이 좋아지지 않으며 오히려 분산이 커질 수 있다는 것을 보았습니다.

마지막으로 실험군과 대조군의 성향점수 분포가 잘 겹치지 않을 때 마주칠 수 있는 외삽 문제를 살펴보았습니다.

참고 자료

저는 이 책 전체를 Joshua Angrist, Alberto Abadie, Christopher Walters의 훌륭한 계량경제학 강의에 바치는 헌사라고 생각합니다. 이 책에 담긴 아이디어 대부분은 이들이 미국경제학회(American Economic Association)에서 한 강의에서 가져왔습니다. 힘든 2020년을 버틸 수 있는 것도 이 강의들을 보는 덕분입니다. * Cross-Section Econometrics * Mastering Mostly Harmless Econometrics

Angrist의 멋진 책들도 참고 자료로 소개하고 싶습니다. 이 책들 덕분에 계량경제학, 이들 표현으로는 ’Metrics가 아주 유용할 뿐 아니라 무척 재미있다는 것을 알게 되었습니다.

• Mostly Harmless Econometrics
• Mastering ’Metrics

마지막 참고 자료는 Miguel Hernan과 Jamie Robins의 책입니다. 까다로운 인과적 질문에 답해야 할 때마다 믿음직한 동반자가 되어 주었습니다.

• Causal Inference Book

이 장에서 사용한 데이터는 Susan Athey와 Stefan Wager의 논문 Estimating Treatment Effects with Causal Forests: An Application에서 가져왔습니다.

기여하기

Causal Inference for the Brave and True는 과학의 통계학이라 할 수 있는 인과추론을 다루는 오픈소스 자료입니다. 파이썬 기반의 무료 소프트웨어만 사용하며, 돈이 없어도, 배경지식이 많지 않아도 누구나 볼 수 있는 자료를 목표로 합니다. 이 책이 도움이 되었고 후원하고 싶다면 Patreon을 방문해 주세요. 금전적으로 기여하기 어렵다면 오타를 고치거나, 수정할 부분을 제안하거나, 이해하기 어려웠던 부분에 의견을 주는 것으로도 도울 수 있습니다. 책의 저장소에서 이슈를 열어 주세요. 마지막으로 이 내용이 마음에 들었다면 도움이 될 만한 사람들과 공유하고 GitHub에서 스타를 눌러 주세요.

힌트옮긴이

한국어판의 오역이나 오탈자는 한국어판 이슈 페이지에 알려 주세요. 한국어판이 도움이 되었다면 한국어판 저장소에도 스타를 부탁드립니다.