20  플러그 앤 플레이 추정량

지금까지는 처치가 무작위로 배정되지 않아 교란편향(confounding bias)이 생기는 경우에 데이터에서 편향을 제거하는 방법을 살펴봤습니다. 이 방법들은 인과추론의 식별(identification) 문제를 푸는 데 도움이 됩니다. 다시 말해, 분석단위가 교환 가능(exchangeable)하다면, 즉 \(Y(0), Y(1) \perp T|X\)이면 처치효과를 학습할 수 있습니다. 하지만 아직 갈 길이 멉니다.

식별할 수 있다는 것은 평균 처치효과(ATE)를 구할 수 있다는 뜻입니다. 즉, 처치가 평균적으로 얼마나 효과적인지 알 수 있습니다. 물론 이 정보는 처치를 전면 도입할지 말지 결정하는 데 도움이 되니 유용합니다. 하지만 그 이상을 원합니다. 처치에 더 잘 반응하거나 덜 반응하는 하위 집단이 있는지 알고 싶습니다. 그러면 훨씬 나은 정책, 즉 처치로 이득을 볼 대상에게만 처치하는 정책을 세울 수 있습니다.

문제 설정

관심 있는 설정을 다시 떠올려 봅시다. 잠재적 결과(potential outcome)가 주어지면 개별 처치효과(ITE)를 두 잠재적 결과의 차이로 정의할 수 있습니다.

\[ \tau_i = Y_i(1) − Y_i(0), \]

처치가 연속형이라면 \(t\)가 처치변수일 때 \(\tau_i = \partial Y(t)\)로 정의합니다. 물론 개별 처치효과는 결코 관측할 수 없습니다. 잠재적 결과 중 하나만 볼 수 있기 때문입니다.

\[ Y^{obs}_i(t)= \begin{cases} Y_i(1), & \text{if } t=1\\ Y_i(0), & \text{if } t=0 \end{cases} \]

평균 처치효과(ATE)는 다음과 같이 정의할 수 있고

\[ \tau = E[Y_i(1) − Y_i(0)] = E[\tau_i] \]

조건부 평균 처치효과(CATE)는 다음과 같이 정의할 수 있습니다.

\[ \tau(x) = E[Y_i(1) − Y_i(0)|X] = E[\tau_i|X] \]

이 책의 1부에서는 주로 ATE에 초점을 맞췄습니다. 이제부터는 CATE를 다룹니다. CATE는 의사결정 과정을 개인화하는 데 유용합니다. 예를 들어 약물이 처치 \(t\)라면, 어떤 유형의 환자가 약물에 더 잘 반응하는지(CATE가 높은지), 그리고 오히려 역효과가 나는 유형의 환자(CATE < 0)가 있는지 알고 싶을 것입니다.

앞에서 처치와 특성(feature) 사이의 상호작용을 넣은 선형회귀로 CATE를 추정하는 방법을 살펴봤습니다.

\[ y_i = \beta_0 + \beta_1 t_i + \beta_2 X_i + \beta_3 t_i X_i + e_i. \]

이 모델을 추정하면 \(\tau(x)\)의 추정값을 얻을 수 있습니다.

\[ \hat{\tau}(x) = \hat{\beta}_1 + \hat{\beta}_3 X_i \]

그래도 선형 모델에는 몇 가지 단점이 있습니다. 가장 큰 단점은 \(X\)에 대한 선형성 가정입니다. 이 모델에서 \(\beta_2\)에는 관심조차 없다는 점에 주목하세요. 하지만 특성 \(X\)가 결과와 선형 관계가 아니라면 인과 매개변수 \(\beta_1\)과 \(\beta_3\)의 추정값이 빗나갑니다.

선형 모델을 더 유연한 머신러닝 모델로 바꿀 수 있다면 좋을 것입니다. 부스팅 트리나 신경망 같은 머신러닝 모델에 처치를 특성으로 넣을 수도 있습니다.

\[ y_i = M(X_i, T_i) + e_i \]

하지만 그다음에 처치효과의 추정값을 어떻게 얻을지가 분명하지 않습니다. 이 모델은 \(\hat{\tau(x)}\)가 아니라 \(\hat{y}\)를 예측하기 때문입니다. 이상적으로는 머신러닝 회귀 모델이 다음과 같은 결과의 평균제곱오차(MSE)를 최소화하는 대신

\[ E[(Y_i - \hat{Y}_i)^2] \]

다음과 같은 처치효과의 MSE를 최소화하면 좋겠습니다.

\[ E[(\tau(x)_i - \hat{\tau}(x)_i)^2] = E[(Y_i(1) - Y_i(0) - \hat{\tau}(x)_i)^2]. \]

하지만 이 기준은 이른바 실현 불가능한(infeasible) 기준입니다. 이번에도 문제는 \(\tau(x)_i\)를 관측할 수 없어서 직접 최적화할 수 없다는 점입니다. 난처한 상황입니다… 문제를 조금 단순하게 만들어 보면 뭔가 떠오를지도 모릅니다.

목표 변환

처치가 이진이라고 해 봅시다. 투자 회사에서 일하고 있고, 금융 교육 이메일을 보내는 것이 효과가 있는지 실험한다고 가정하겠습니다. 이메일 덕분에 사람들이 투자를 더 많이 하기를 바랍니다. 그리고 고객의 50%는 이메일을 받고 나머지 50%는 받지 않는 무작위 실험을 했다고 하겠습니다.

여기서 정신 나간 아이디어를 하나 내 보겠습니다. 결과변수에 처치를 곱해서 결과변수를 변환하는 것입니다.

\[ Y^*_i = 2 Y_i * T_i - 2 Y_i*(1-T_i) \]

즉, 처치를 받은 분석단위라면 결과에 2를 곱하고, 처치를 받지 않았다면 -2를 곱합니다. 예를 들어 어떤 고객이 이메일을 받고 2,000헤알(BRL)을 투자했다면 변환된 목표는 4,000이 됩니다. 하지만 이메일을 받지 않았다면 -4,000이 됩니다.

이메일의 효과가 음수일 수도 있다고 말하는 셈이니 무척 이상해 보이겠지만, 조금만 참고 따라와 주세요. 약간의 수학을 거치면 이 변환된 목표가 평균적으로, 즉 기댓값으로 보면 처치효과와 같다는 것을 알 수 있습니다. 놀랍다는 말로도 부족합니다. 조금 엉뚱한 이 변환을 적용하기만 하면 관측할 수조차 없는 값을 추정할 수 있다는 말이니까요.

이를 이해하려면 수학이 조금 필요합니다. 무작위 배정 덕분에 \(T \perp Y(0), Y(1)\)이 성립합니다. 우리의 오랜 친구인 비교란성(unconfoundedness)입니다. 그리고 독립성에서 바로 \(E[T, Y(t)]=E[T]*E[Y(t)]\)가 따라 나옵니다.

또한 다음이 성립한다는 것도 알고 있습니다.

\[ Y_i*T_i = Y(1)_i*T_i \text{ and } Y_i*(1-T_i) = Y(0)_i*(1-T_i) \]

처치가 두 잠재적 결과 중 어느 쪽이 실현될지를 결정하기 때문입니다. 이 점을 염두에 두고 \(Y^*_i\)의 기댓값을 구해 보겠습니다.

\[ \begin{align} E[Y^*_i|X_i=x] &= E[2 Y(1)_i * T_i - 2 Y(0)_i*(1-T_i)|X_i=x] \\ &= 2E[Y(1)_i * T_i | X_i=x] - 2E[Y(0)_i*(1-T_i)|X_i=x]\\ &= 2E[Y(1)_i| X_i=x] * E[ T_i | X_i=x] - 2E[Y(0)_i| X_i=x]*E[(1-T_i)|X_i=x] \\ &= 2E[Y(1)_i| X_i=x] * 0.5 - 2E[Y(0)_i| X_i=x]*0.5 \\ &= E[Y(1)_i| X_i=x] - E[Y(0)_i| X_i=x] \\ &= \tau(x)_i \end{align} \]

결국 얼핏 정신 나간 듯했던 이 아이디어가 개별 처치효과 \(\tau(x)_i\)의 편향되지 않은 추정값이었던 셈입니다. 이제 실현 불가능했던 최적화 기준을 다음 식으로 바꿀 수 있습니다.

\[ E[(Y^*_i - \hat{\tau}(x)_i)^2] \]

쉽게 말해, 아무 머신러닝 회귀 모델로 \(Y^*_i\)를 예측하기만 하면 그 모델이 처치효과의 예측값을 내놓습니다.

간단한 경우는 해결했으니, 처치가 50 대 50으로 배정되지 않았거나 아예 무작위로 배정되지 않은 더 복잡한 경우는 어떨까요? 알고 보면 답이 조금 더 복잡해지지만, 아주 많이 복잡하지는 않습니다. 먼저 무작위 배정이 아니라면 최소한 조건부 독립성 \(T \perp Y(1), Y(0) | X\)가 필요합니다. 즉, \(X\)를 통제하면 \(T\)가 무작위 배정이나 다름없어야 합니다. 그러면 변환된 목표를 다음과 같이 일반화할 수 있습니다.

\[ Y^*_i = Y_i * \dfrac{T_i - e(X_i)}{e(X_i)(1-e(X_i))} \]

여기서 \(e(X_i)\)는 성향점수(propensity score)입니다. 따라서 처치가 50 대 50이 아니라 다른 확률 \(p\)로 무작위 배정됐다면 위 공식의 성향점수를 \(p\)로 바꾸기만 하면 됩니다. 처치가 무작위가 아니라면 저장해 둔 성향점수나 추정한 성향점수를 써야 합니다.

이 값의 기댓값을 구해 보면 역시 처치효과와 일치합니다. 증명은 독자 여러분의 연습 문제로 남겨 두겠습니다. 농담입니다, 여기 있습니다. 조금 번거로우니 건너뛰어도 괜찮습니다.

\[ \begin{align} E[Y^*_i|X_i=x] &= E\big[Y_i * \dfrac{T_i - e(X_i)}{e(X_i)(1-e(X_i))}|X_i=x\big] \\ &= E\big[Y_i T_i * \dfrac{T_i - e(X_i)}{e(X_i)(1-e(X_i))} + Y_i (1-T_i) * \dfrac{T_i - e(X_i)}{e(X_i)(1-e(X_i))}|X_i=x\big]\\ &= E\big[Y(1)_i * \dfrac{T_i(1 - e(X_i))}{e(X_i)(1-e(X_i))} | X_i=x\big] - E\big[Y(0)_i * \dfrac{(1-T_i)e(X_i)}{e(X_i)(1-e(X_i))}|X_i=x\big]\\ &= \dfrac{1}{e(X_i)} E[Y(1)_i * T_i|X_i=x] - \dfrac{1}{1-e(X_i)} E[Y(0)_i * (1-T_i)| X_i=x]\\ &= \dfrac{1}{e(X_i)} E[Y(1)_i|X_i=x] * E[T_i|X_i=x] - \dfrac{1}{1-e(X_i)} E[Y(0)_i|X_i=x] * E[(1-T_i)| X_i=x]\\ &= E[Y(1)_i|X_i=x] - E[Y(0)_i|X_i=x]\\ &= \tau(x)_i \end{align} \]

늘 그렇듯 예제를 보면 훨씬 구체적으로 와닿을 것입니다. 이번에도 사람들이 투자를 더 하도록 보낸 투자 이메일을 생각해 봅시다. 결과변수는 이진 변수(투자함 vs 투자하지 않음)인 converted입니다.

import pandas as pd
import numpy as np
from matplotlib import pyplot as plt
import seaborn as sns
from nb21 import cumulative_gain, elast

email = pd.read_csv("./data/invest_email_rnd.csv")
email.head()

	age	income	insurance	invested	em1	em2	em3	converted
0	44.1	5483.80	6155.29	14294.81	0	1	1	0
1	39.8	2737.92	50069.40	7468.15	1	0	0	0
2	49.0	2712.51	5707.08	5095.65	1	0	1	1
3	39.7	2326.37	15657.97	6345.20	1	1	1	0
4	35.3	2787.26	27074.44	14114.86	1	1	1	0

여기서 목표는 개인화(personalization)입니다. email-1에 집중해 보겠습니다. 이 이메일에 더 잘 반응할 고객에게만 이메일을 보내고 싶습니다. 다시 말해 email-1의 조건부 평균 처치효과를 추정하려고 합니다.

\[ E[Converted(1)_i - Converted(0)_i|X_i=x] = \tau(x)_i \]

그래야 이메일에 가장 잘 반응할(CATE가 높은) 고객을 골라서 보낼 수 있습니다.

하지만 그 전에 데이터셋을 훈련 셋과 검증 셋으로 나누겠습니다. 한쪽 셋에서 \(\tau(x)_i\)를 추정하고 다른 쪽 셋에서 그 추정값을 평가합니다.

from sklearn.model_selection import train_test_split

np.random.seed(123)
train, test = train_test_split(email, test_size=0.4)
print(train.shape, test.shape)

(9000, 8) (6000, 8)

이제 방금 배운 목표 변환(target transformation)을 적용해 보겠습니다. 이메일은 (50 대 50은 아니지만) 무작위로 배정됐으므로 성향점수를 신경 쓸 필요가 없습니다. 성향점수는 처치 확률과 같은 상수입니다.

y = "converted"
T = "em1"
X = ["age", "income", "insurance", "invested"]

ps = train[T].mean()

y_star_train = train[y] * (train[T] - ps)/(ps*(1-ps))

변환된 목표가 있으면 어떤 머신러닝 회귀 알고리즘으로든 이를 예측할 수 있습니다. 여기서는 부스팅 트리를 쓰겠습니다.

from lightgbm import LGBMRegressor

np.random.seed(123)
cate_learner = LGBMRegressor(max_depth=3, min_child_samples=300, num_leaves=5)
cate_learner.fit(train[X], y_star_train);

이제 이 모델로 \(\tau(x)_i\)를 추정할 수 있습니다. 즉, 모델이 출력하는 값은 \(\hat{\tau}(x)_i\)입니다. 예를 들어 테스트 셋에서 예측해 보면 어떤 분석단위는 다른 분석단위보다 CATE가 높게 나옵니다. 가령 고객 6958의 CATE는 0.1인데, 이 고객에게 이메일을 보내면 투자 상품을 구매할 확률이 0.1 오를 것으로 예측된다는 뜻입니다. 반면 고객 3903은 상품 구매 확률이 0.04만 오를 것으로 예측됩니다.

test_pred = test.assign(cate=cate_learner.predict(test[X]))
test_pred.head()

	age	income	insurance	invested	em1	em2	em3	converted	cate
6958	40.9	4486.14	37320.33	12559.25	0	0	1	0	0.105665
7534	42.6	6386.19	13270.47	29114.42	0	0	1	0	0.121922
2975	47.6	1900.26	25588.72	2420.39	0	0	1	0	0.034161
3903	41.0	5802.19	57087.37	20182.20	1	0	1	1	0.046805
8437	49.1	2202.96	5050.81	9245.88	1	0	1	0	-0.009099

이 모델이 얼마나 좋은지 평가하려면 훈련 셋과 테스트 셋에서 각각 누적 이득 곡선(cumulative gain curve)을 그려 보면 됩니다.

gain_curve_test = cumulative_gain(test_pred, "cate", y="converted", t="em1")
gain_curve_train = cumulative_gain(train.assign(cate=cate_learner.predict(train[X])), "cate", y="converted", t="em1")
plt.plot(gain_curve_test, color="C0", label="Test")
plt.plot(gain_curve_train, color="C1", label="Train")
plt.plot([0, 100], [0, elast(test, "converted", "em1")], linestyle="--", color="black", label="Baseline")
plt.legend();

보다시피 이 플러그 앤 플레이 추정량은 테스트 셋에서 무작위보다 낫습니다. 그래도 훈련 셋의 성능이 테스트 셋보다 훨씬 좋은 것으로 보아 과적합이 꽤 심한 듯합니다.

사실 이것이 목표 변환 기법의 가장 큰 단점 가운데 하나입니다. 목표 변환을 쓰면 목표만 변환한 뒤 아무 머신러닝 추정량으로 이질적 처치효과를 예측할 수 있으니 무척 간단합니다. 하지만 그 대가로 분산이 크게 늘어납니다. 변환된 목표는 개별 처치효과의 추정값으로서 잡음이 매우 심하고, 그 분산이 그대로 추정에 옮겨 가기 때문입니다. 데이터가 많지 않다면 큰 문제지만, 표본이 100만 개를 넘는 빅데이터 환경에서는 문제가 덜할 것입니다.

연속형 처치의 경우

목표 변환 방법의 또 다른 뚜렷한 단점은 이산형이나 이진 처치에서만 작동한다는 점입니다. 인과추론 문헌에서 흔히 볼 수 있는 상황입니다. 연구 대부분이 이진 처치를 다루고, 연속형 처치에 관한 내용은 찾기 어렵습니다. 저는 이 점이 무척 답답했습니다. 산업 현장에는 연속형 처치가 어디에나 있고, 대부분은 최적화해야 하는 가격의 형태이기 때문입니다. 결국 연속형 처치를 위한 목표 변환은 찾지 못했지만, 실무에서 통하는 방법을 직접 고안해 냈습니다. 다만 이 방법을 뒷받침하는 탄탄한 계량경제학 연구는 없다는 점을 염두에 두세요.

이 방법을 소개하기 위해 아이스크림 판매 예제로 돌아가 보겠습니다. 그 예제에서 맡은 과제는 가격에 대한 수요의 탄력성(elasticity)을 추정하는 것이었습니다. 그래야 아이스크림 가격을 더 잘 정해 수익을 최적화할 수 있습니다. 데이터셋의 각 표본은 하루이고, 사람들이 가격 인상에 덜 민감한 때가 언제인지 알고 싶었다는 점을 떠올려 보세요. 또한 이 데이터셋에서는 가격이 무작위로 배정됐으므로 교란편향을 걱정할 필요가 없다는 점도 기억하실 것입니다.

prices_rnd = pd.read_csv("./data/ice_cream_sales_rnd.csv")
prices_rnd.head()

	temp	weekday	cost	price	sales
0	25.8	1	0.3	7	230
1	22.7	3	0.5	4	190
2	33.7	7	1.0	5	237
3	23.0	4	0.5	5	193
4	24.4	1	1.0	3	252

앞에서처럼 먼저 데이터를 훈련 셋과 테스트 셋으로 나누겠습니다.

np.random.seed(123)
train, test = train_test_split(prices_rnd, test_size=0.3)
train.shape, test.shape

((3500, 5), (1500, 5))

이제 약간의 창의력이 필요한 대목입니다. 이산형 처치에서 조건부 평균 처치효과는 분석단위의 특성 \(X\)가 주어졌을 때 처치를 받지 않은 상태에서 처치를 받은 상태로 바뀌면 결과가 얼마나 변하는지로 정의됩니다.

\[ \tau(x) = E[Y_i(1) − Y_i(0)|X] = E[\tau_i|X] \]

쉽게 말해, 특성 \(X\)로 정의한 여러 분석단위 프로필에 처치가 미치는 영향을 추정하는 것입니다. 연속형 처치에는 이렇게 켜고 끄는 스위치가 없습니다. 분석단위는 처치를 받거나 안 받는 것이 아니라, 모두 처치를 받되 그 강도가 다릅니다. 그러므로 처치를 주는 효과는 말할 수 없고, 대신 처치를 늘리는 관점에서 이야기해야 합니다. 다시 말해, 처치를 어느 정도 늘리면 결과가 어떻게 바뀌는지 알고 싶습니다. 이는 결과 함수 \(Y\)를 처치 \(t\)에 대해 편미분한 값을 추정하는 것과 같습니다. 그리고 이를 그룹마다 알고 싶으므로(ATE가 아니라 CATE이므로) 특성 \(X\)를 조건으로 둡니다.

\[ \tau(x) = E[\partial Y_i(t)|X] = E[\tau_i|X] \]

이 값은 어떻게 추정할 수 있을까요? 먼저 결과가 처치에 대해 선형인 쉬운 경우부터 생각해 봅시다. 더운 날(노란색)과 추운 날(파란색), 두 종류의 날이 있다고 해 보겠습니다. 추운 날에는 사람들이 가격 인상에 더 민감합니다. 또한 가격이 오르면 수요는 선형으로 줄어듭니다.

이 경우 CATE는 각 수요 직선의 기울기입니다. 이 기울기를 보면 가격을 얼마만큼 올렸을 때 수요가 얼마나 줄어드는지 알 수 있습니다. 이 관계가 정말 선형이라면 더운 날과 추운 날 각각 따로 단순선형회귀를 추정해서, 그 계수로 탄력성을 추정할 수 있습니다.

\[ \hat{\tau(x)} = Cov(Y_i, T_i)/Var(T_i) = \dfrac{\sum(Y_i- \bar{Y})(T_i - \bar{T})}{\sum (T_i - \bar{T})^2} \]

이 추정량에서 영감을 얻어, 개별 분석단위라면 이 식이 어떤 모습일지 생각해 볼 수 있습니다. 즉, 위와 같은 값을 날마다 하나씩 정의한다면 어떨까요? 제 머릿속에서는 대략 이런 모습입니다.

\[ Y^*_i = (Y_i- \bar{Y})\dfrac{(T_i - \bar{T})}{\sigma^2_T} \]

말로 풀면, 원래 목표에서 평균을 빼고, 마찬가지로 평균을 뺀 처치를 곱한 다음, 마지막으로 처치의 분산으로 나누는 것입니다. 짜잔, 연속형 처치를 위한 목표 변환이 완성됐습니다.

이제 남은 질문은 이것이 정말 통하느냐입니다. 실제로 통합니다. 이진 처치에서 했던 것처럼 왜 통하는지 비슷하게 증명해 볼 수도 있습니다. 먼저 다음과 같이 두겠습니다.

\[ V_i = \dfrac{(T_i - \bar{T})}{\sigma^2_T} \]

무작위 배정에서는 \(E[T_i|X_i=x]=\bar{T}\)이므로 \(E[V_i|X_i=x]=0\)이라는 점에 주목하세요. 다시 말해, X의 모든 영역에서 \(E[T_i]=\bar{T}\)입니다. 또한 \(E[T_i(T_i - \bar{T})|X_i=x] = E[(T_i - \bar{T})^2|X_i=x]\)이고 이는 처치의 분산이므로 \(E[T_i V_i | X_i=x]=1\)입니다. 마지막으로, 조건부 독립성(무작위 처치 배정이라면 저절로 얻는 성질) 아래에서 \(E[T_i e_i | X_i=x] = E[T_i | X_i=x] E[e_i | X_i=x]\)입니다.

이 목표 변환이 통한다는 것을 보이려면, 지금 추정하는 것이 다음과 같은 국소 선형 모델의 매개변수라는 점을 기억해야 합니다.

\[ Y_i = \alpha + \beta T_i + e_i | X_i=x \]

이 예제에서는 더운 날과 추운 날 각각의 선형 모델이 여기에 해당합니다. 여기서 관심 있는 것은 \(\beta\) 매개변수이고, 이것이 조건부 탄력성, 즉 CATE입니다. 이 모든 것을 이용하면 다음을 증명할 수 있습니다.

\[ \begin{align} E[Y^*_i|X_i=X] &= E[(Y_i-\bar{Y})V_i | X_i=x] \\ &= E[(\alpha + \beta T_i + e_i - \bar{Y})V_i | X_i=x] \\ &= \alpha E[V_i | X_i=x] + \beta E[T_i V_i | X_i=x] + E[e_i V_i | X_i=x] \\ &= \beta + E[e_i V_i | X_i=x] \\ &= \beta = \tau(x) \end{align} \]

이 방법은 처치가 무작위로 배정됐을 때만 통한다는 점을 명심하세요. 무작위 배정이 아니라면 \(\bar{T}\)를 \(M(X_i)\)로 바꿔야 합니다. 여기서 \(M\)은 \(E[T_i|X_i=x]\)를 추정하는 모델입니다.

\[ Y^*_i = (Y_i- \bar{Y})\dfrac{(T_i - M(T_i))}{(T_i - M(T_i))^2} \]

이렇게 하면 셋째 줄의 \(\alpha E[V_i | X_i=x]\) 항이 0으로 사라지고 \(E[T_i V_i | X_i=x]\) 항이 1이 됩니다. 처치효과 순으로 분석단위를 정렬하기만 하면 된다면 \(E[T_i V_i | X_i=x]\)가 꼭 1이 될 필요는 없다는 점도 알아 두세요. 즉, 어느 날 수요가 가격 인상에 더 민감한지만 알면 되고 얼마나 더 민감한지는 몰라도 된다면, \(\beta\) 추정값이 일정한 배율로 커지거나 작아져도 상관없습니다. 그런 경우에는 분모를 생략할 수 있습니다.

\[ Y^*_i = (Y_i- \bar{Y})(T_i - M(T_i)) \]

수학이 지루하게 느껴지더라도 걱정하지 마세요. 코드는 정말 간단합니다. 이번에도 위의 공식으로 훈련 데이터의 목표를 변환합니다. 이 데이터는 처치가 무작위로 배정됐으므로 가격을 예측하는 모델을 만들 필요가 없습니다. 또 여기서는 처치효과의 순서만 중요하므로 분모도 생략하겠습니다.

y_star_cont = (train["price"] - train["price"].mean()) * (train["sales"] - train["sales"].mean())

그다음은 앞에서와 똑같이 머신러닝 회귀 모델을 학습시켜 이 목표를 예측합니다.

cate_learner = LGBMRegressor(max_depth=3, min_child_samples=300, num_leaves=5)

np.random.seed(123)
cate_learner.fit(train[["temp", "weekday", "cost"]], y_star_cont)

cate_test_transf_y = cate_learner.predict(test[["temp", "weekday", "cost"]])

test_pred = test.assign(cate=cate_test_transf_y)
test_pred.sample(5)

	temp	weekday	cost	price	sales	cate
2815	15.7	4	1.5	3	187	-1395.956278
257	29.4	3	1.0	3	209	-1607.400415
2585	24.6	6	1.0	10	197	-1497.197402
3260	20.2	1	0.5	4	246	-1629.798111
1999	10.0	4	0.5	10	139	-1333.690544

이번에는 CATE를 직관적으로 해석하기 어렵습니다. 목표 변환에서 분모를 뺐기 때문에 지금 보이는 CATE는 \(Var(X)\)만큼 배율이 달라진 값입니다. 그래도 이 예측값은 처치효과의 순서를 꽤 잘 정렬할 것입니다. 이를 확인하려면 앞에서처럼 누적 이득 곡선을 그려 보면 됩니다.

gain_curve_test = cumulative_gain(test.assign(cate=cate_test_transf_y),
                                "cate", y="sales", t="price")

gain_curve_train = cumulative_gain(train.assign(cate=cate_learner.predict(train[["temp", "weekday", "cost"]])),
                                   "cate", y="sales", t="price")


plt.plot(gain_curve_test, label="Test")
plt.plot(gain_curve_train, label="Train")
plt.plot([0, 100], [0, elast(test, "sales", "price")], linestyle="--", color="black", label="Baseline")
plt.legend();

이 데이터에서는 변환된 목표로 학습한 모델이 무작위보다 훨씬 나아 보입니다. 그뿐 아니라 훈련 셋과 테스트 셋의 결과가 꽤 비슷하므로 여기서는 분산이 문제가 되지 않습니다. 하지만 이는 이 데이터셋의 특성일 뿐입니다. 기억하시겠지만 이진 처치를 다룰 때는 그렇지 않았습니다. 그때는 모델 성능이 그리 좋지 않았습니다.

비선형 처치효과

연속형 처치까지 이야기했지만, 아직 외면하고 있던 큰 문제가 하나 남아 있습니다. 지금까지는 처치효과가 선형이라고 가정했습니다. 하지만 이것이 합리적인 가정인 경우는 매우 드뭅니다. 보통 처치효과는 어떤 식으로든 포화됩니다. 이 예제에서는 가격을 처음 몇 단위 올릴 때는 수요가 빠르게 줄다가 그 뒤로는 더 천천히 줄어든다고 보는 것이 합리적입니다.

문제는 탄력성 또는 처치효과가 처치 자체에 따라 변한다는 점입니다. 이 예제에서 처치효과는 곡선의 앞부분에서 더 강하고 가격이 높아질수록 작아집니다. 다시 더운 날(노란색)과 추운 날(파란색), 두 종류의 날이 있고 인과 모델로 둘을 구분하고 싶다고 해 봅시다. 그런데 인과 모델이 예측하는 것은 탄력성입니다. 비선형인 경우에는 곡선에서 서로 다른 가격 지점을 보면 더운 날과 추운 날의 탄력성이 같게 나올 수 있습니다(오른쪽 그림).

이 문제를 쉽게 빠져나갈 방법은 없고, 고백하자면 저도 무엇이 가장 잘 통하는지 아직 연구하고 있습니다. 지금으로서는 처치효과의 함수 형태를 생각해 보고 이를 어떻게든 선형화하는 방식을 씁니다. 예를 들어 수요는 보통 다음과 같은 함수 형태를 띠는데, \(\alpha\)가 클수록 가격을 올릴 때마다 수요가 더 빠르게 줄어듭니다.

\[ D_i = \dfrac{1}{P_i^{\alpha}} \]

따라서 수요 \(Y\)와 가격 \(T\)에 모두 로그 변환을 적용하면 선형인 형태를 얻을 수 있습니다.

\[ \begin{align} log(D)_i &= log\bigg(\dfrac{1}{P_i^{\alpha}}\bigg) \\ &= log(1) - log(P_i^{\alpha}) \\ &= log(1) - log(P_i^{\alpha}) \\ &= - \alpha * log(P_i) \\ \end{align} \]

선형화는 꽤 고민해야 하는 일이라 그리 쉽지 않습니다. 하지만 이것저것 시도해 보면서 무엇이 가장 잘 통하는지 확인해도 됩니다. 로그나 제곱근 같은 변환이 도움이 되는 경우가 많습니다.

핵심 요약

이제 머신러닝 모델로 조건부 평균 처치효과를 추정하는 방향으로 나아가고 있습니다. 이때 가장 큰 과제는 예측 모델을 인과효과를 추정하는 모델로 바꾸는 것입니다. 달리 생각해 보면, 예측 모델은 결과 Y를 특성 X와 (경우에 따라) 처치 T의 함수 \(Y = M(X, T)\)로 추정하는 데 초점을 두지만, 인과 모델은 이 결과 함수를 처치에 대해 편미분한 \(\partial Y = \partial M(X, T)\)를 추정해야 합니다. 이는 결코 간단하지 않습니다. 결과 Y는 관측할 수 있지만 \(\partial Y\)는, 적어도 개별 수준에서는 관측할 수 없기 때문입니다. 그래서 모델의 목적함수를 설계할 때 창의력을 발휘해야 합니다.

이 장에서는 목표 변환이라는 아주 간단한 기법을 살펴봤습니다. 원래 목표 Y와 처치 T를 결합해, 기댓값이 CATE와 같은 변환된 목표를 만드는 것이 핵심 아이디어입니다. 이 새 목표가 있으면 어떤 머신러닝 예측 모델이든 끼워 넣어 추정할 수 있고, 그 모델의 예측값이 곧 CATE의 추정값이 됩니다. 덧붙이면, 이 방법은 목표 변환 말고 F 러너라는 이름으로도 불립니다.

이렇게 간단한 만큼 치러야 할 대가도 있습니다. 변환된 목표는 개별 처치효과의 추정값으로서 잡음이 매우 심하고, 그 잡음은 분산의 형태로 모델 추정값에 옮겨 갑니다. 그래서 목표 변환은 표본 크기 자체가 커서 분산이 덜 문제 되는 빅데이터 환경에 더 잘 맞습니다. 목표 변환 방법의 또 다른 단점은 이진 처치나 범주형 처치에만 정의된다는 점입니다. 최선을 다해 연속형 버전을 고안했고 그럭저럭 통하는 듯한 방법까지 얻었지만, 아직까지는 이를 뒷받침할 탄탄한 이론적 틀이 없습니다.

마지막으로 비선형 처치효과와 그에 따르는 어려움을 논의했습니다. 처치효과가 처치 자체에 따라 변하면 착각하기 쉽습니다. 처치에 대한 반응도가 같다는 이유로 분석단위들의 처치 반응 곡선이 같다고 생각할 수 있는데, 실제로는 서로 다른 양의 처치를 받고 있을 뿐입니다.

참고 자료

여기 쓴 내용은 대부분 제 머릿속에서 나온 것입니다. 경험으로 배운 것들입니다. 즉, 좋은 과학이라면 보통 거치는 학계의 엄밀한 검증을 통과하지 않았다는 뜻입니다. 대신 실무에서 통하는 것들을 이야기하면서도 왜 그런지 설명하는 데는 시간을 많이 들이지 않았다는 점을 눈여겨보세요. 말하자면 길거리 과학인 셈입니다. 하지만 이 내용을 공개해 누구나 검토할 수 있게 해 두었으니, 터무니없는 내용을 발견하면 주저하지 말고 이슈를 열어 주세요. 최선을 다해 답하겠습니다.

이 장의 대부분은 Susan Athey와 Guido W. Imbens의 논문 Machine Learning Methods for Estimating Heterogeneous Causal Effects에서 가져왔습니다. 목표 변환에 관한 내용 일부는 Pierre Gutierrez와 Jean-Yves Gérardy의 논문 Causal Inference and Uplift Modeling: A review of the literature에서도 찾을 수 있습니다. 두 논문 모두 이진 처치만 다룬다는 점에 유의하세요. F 러너를 언급하는, CATE 추정용 인과 모델에 관한 또 다른 리뷰로는 Künzel et al.(2019)의 Meta-learners for Estimating Heterogeneous Treatment Effects using Machine Learning이 있습니다.

기여하기

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

힌트옮긴이

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