누락 변수 편향(OVB)을 통제하는 한 가지 방법은, 뭐, 빠뜨린 변수를 모델에 넣는 것입니다. 하지만 늘 그렇게 할 수 있는 것은 아닙니다. 대개는 누락된 변수에 대한 데이터가 아예 없기 때문입니다. 예를 들어 교육이 임금에 미치는 효과를 다룬 모델로 돌아가 보겠습니다.
\[ \log(\mathrm{wage})_i = \beta_0 + \kappa \ \mathrm{educ}_i + \pmb{\beta}\mathrm{Ability}_i + u_i \]
교육이 \(\log\mathrm{(wage)}\)에 미치는 인과효과 \(\kappa\)를 알아내려면 능력 요인 \(\mathrm{Ability}_i\)를 통제해야 합니다. 그러지 않으면 편향이 생길 가능성이 큽니다. 능력은 처치인 교육과 결과인 소득 모두에 영향을 주는 교란 요인(confounder)일 테니까요.
import warnings
warnings.filterwarnings('ignore')
import pandas as pd
import numpy as np
from scipy import stats
from matplotlib import style
import seaborn as sns
from matplotlib import pyplot as plt
import statsmodels.formula.api as smf
import graphviz as gr
from linearmodels.iv import IV2SLS
%matplotlib inline
pd.set_option("display.max_columns", 5)
style.use("fivethirtyeight")g = gr.Digraph()
g.edge("ability", "educ")
g.edge("ability", "wage")
g.edge("educ", "wage")
g
이 문제를 피하는 한 가지 방법은 교육이 임금에 미치는 효과를 측정할 때 능력을 일정한 수준으로 고정해 통제하는 것입니다. 선형회귀 모델에 능력을 포함하면 그렇게 할 수 있습니다. 하지만 능력을 제대로 측정한 데이터가 없습니다. 기껏해야 IQ처럼 상당히 미심쩍은 대리변수 정도가 있을 뿐입니다.
그렇다고 끝난 것은 아닙니다. 바로 여기서 도구변수(instrumental variable)가 등장합니다. 도구변수의 아이디어는 처치에 영향을 주면서 결과와는 처치를 거쳐서만 상관관계를 갖는 다른 변수를 찾는 것입니다. 달리 말하면, 이 도구변수 \(Z_i\)는 \(Y_0\), \(Y_1\)과는 상관관계가 없지만 \(T\)와는 상관관계가 있습니다. 이를 배제 제약(exclusion restriction)이라고 부르기도 합니다.
g = gr.Digraph()
g.edge("ability", "educ")
g.edge("ability", "wage")
g.edge("educ", "wage")
g.edge("instrument", "educ")
g
이런 변수가 있다면, 곧 살펴볼 도구변수 공식으로 인과효과 \(\kappa\)를 구할 수 있습니다. 이를 위해 먼저 이상적으로 돌리고 싶은 식을 생각해 보겠습니다. 처치를 \(T\), 교란 요인을 \(W\)처럼 좀 더 일반적인 기호로 쓰면, 원하는 식은 다음과 같습니다.
\[ Y_i = \beta_0 + \kappa \ T_i + \pmb{\beta}W_i + u_i \]
하지만 \(W\)에 대한 데이터가 없으므로 돌릴 수 있는 식은 이것뿐입니다.
\[ Y_i = \beta_0 + \kappa\ T_i + v_i \]
\[ v_i = \pmb{\beta}W_i + u_i \]
\(W\)는 교란 요인이므로 \(\mathrm{Cov}(T, v) \neq 0\)입니다. 지금 가진 것은 긴 식이 아니라 짧은 식입니다. 앞의 예로 말하면 능력이 교육과 상관관계가 있다는 뜻입니다. 이 경우 짧은 회귀식을 돌리면 누락 변수 때문에 \(\kappa\)에 대한 편향된 추정량을 얻게 됩니다.
이제 도구변수의 마법을 감상할 차례입니다! 도구변수 Z는 T를 거쳐서만 결과와 상관관계를 가지므로 \(\mathrm{Cov}(Z,v) = 0\)입니다. 그렇지 않다면 Z에서 W를 거쳐 Y로 가는 두 번째 경로가 있다는 뜻이 되니까요. 이 점을 염두에 두면 다음과 같이 쓸 수 있습니다.
\[ \mathrm{Cov}(Z,Y) = \mathrm{Cov}(Z,\beta_0 + \kappa\ T_i + v_i) = \kappa \mathrm{Cov}(Z,T) + \mathrm{Cov}(Z, v) = \kappa \mathrm{Cov}(Z,T) \]
양변을 \(V(Z_i)\)로 나누고 항을 정리하면 다음을 얻습니다.
\[ \kappa = \dfrac{\mathrm{Cov}(Y_i, Z_i)/V(Z_i)}{\mathrm{Cov}(T_i, Z_i)/V(Z_i)} = \dfrac{\text{Reduced Form}}{\text{1st Stage}} \]
분자와 분모가 모두 회귀계수(공분산을 분산으로 나눈 값)라는 점에 주목하세요. 분자는 Y를 Z에 회귀한 결과입니다. 다시 말해 Z가 Y에 미치는 ’영향’입니다. 그렇다고 Z가 Y의 원인이라는 말은 아닙니다. Z는 T를 통해서만 Y에 영향을 준다는 조건이 있으니까요. 결국 Z가 T를 거쳐 Y에 미치는 효과가 얼마나 큰지를 담고 있을 뿐입니다. 이 분자는 워낙 유명해서 축약형(reduced form) 계수라는 이름까지 붙어 있습니다.
분모 역시 회귀계수입니다. 이번에는 T를 Z에 회귀한 것입니다. 이 회귀는 Z가 T에 미치는 영향을 담고 있고, 이것 역시 유명해서 1단계(first stage) 계수라고 부릅니다.
이 식을 편미분으로 보는 것도 재미있습니다. T가 Y에 미치는 영향은 Z가 Y에 미치는 영향을 Z가 T에 미치는 영향으로 환산한 것과 같다는 것을 보일 수 있습니다.
\[ \kappa = \dfrac{\frac{\partial y}{\partial z}}{\frac{\partial T}{\partial z}} = \dfrac{\partial y}{\partial z} * \dfrac{\partial z}{\partial T} = \dfrac{\partial y}{\partial T} \]
이 식이 보여 주는 바는 대부분의 사람이 생각하는 것보다 미묘하고, 또 생각보다 훨씬 멋집니다. 도구변수를 이렇게 쓰면 이런 말을 하는 셈입니다. “봐, 교란 요인 때문에 T가 Y에 미치는 영향은 알아내기 어려워. 그런데 Z가 Y에 미치는 영향은 쉽게 구할 수 있지. Z와 Y에 함께 영향을 주는 건 아무것도 없으니까(배제 제약). 하지만 내가 궁금한 건 Z가 Y에 미치는 영향이 아니라 T가 Y에 미치는 영향이야. 그러니 구하기 쉬운 Z의 Y에 대한 효과를 먼저 추정하고, 그 값을 Z가 T에 미치는 효과로 환산해서 Z 단위가 아닌 T 단위의 효과로 바꾸면 돼.”
도구변수가 더미변수인 단순한 경우에서도 이를 확인할 수 있습니다. 이때 도구변수 추정량은 두 평균 차이의 비율로 더 간단해집니다.
\[ \kappa = \dfrac{E[Y|Z=1]-E[Y|Z=0]}{E[T|Z=1]-E[T|Z=0]} \]
이 비율을 왈드 추정량(Wald estimator)이라고 부르기도 합니다. 여기서도 같은 도구변수 이야기를 할 수 있습니다. 알고 싶은 것은 T가 Y에 미치는 효과인데, 이는 구하기 어렵습니다. 그래서 구하기 쉬운 Z의 Y에 대한 효과에 집중합니다. 정의상 Z는 T를 통해서만 Y에 영향을 주므로, Z가 Y에 미치는 영향을 T가 Y에 미치는 영향으로 바꿀 수 있습니다. Z의 Y에 대한 효과를 Z의 T에 대한 효과로 나누어 환산하면 됩니다.
지금까지는 도구변수를 처치를 통해서만 결과에 영향을 주는 기적 같은 성질을 지닌 마법의 변수 \(Z\)로 다뤄 왔습니다. 솔직히 말해 좋은 도구변수는 너무 찾기 어려워서 기적이라고 불러도 과언이 아닙니다. 마음 약한 사람이 할 일은 아니라고만 해 두겠습니다. 소문에 따르면 시카고학파의 잘나가는 친구들은 술집에서 자기가 이런저런 도구변수를 어떻게 생각해 냈는지 떠든다고 합니다.
그래도 이야기를 조금 더 구체적으로 만들어 줄 흥미로운 도구변수 사례가 몇 가지 있습니다. 다시 교육이 임금에 미치는 효과를 추정해 보겠습니다. 이번에는 그 사람이 태어난 분기를 도구변수 Z로 씁니다.
이 아이디어는 미국의 의무 취학법을 활용합니다. 이 법은 보통 아이가 입학하는 해의 1월 1일까지 만 6세가 되어야 한다고 정합니다. 그래서 연초에 태어난 아이는 더 늦은 나이에 학교에 들어갑니다. 의무 취학법은 또 학생이 만 16세가 될 때까지 학교에 다니도록 정하고 있으며, 그 나이가 되면 법적으로 중퇴할 수 있습니다. 그 결과 연말에 태어난 사람이 연초에 태어난 사람보다 평균적으로 교육 연수가 더 깁니다.
출생 분기가 능력 요인과 독립이라는 것, 즉 교육이 임금에 미치는 효과를 교란하지 않는다는 것을 받아들인다면 출생 분기를 도구변수로 쓸 수 있습니다. 달리 말해, 출생 분기는 교육에 미치는 영향을 거치지 않고는 임금에 아무런 영향을 주지 않는다고 믿어야 합니다. 별자리 운세를 믿지 않는다면 꽤 설득력 있는 주장입니다.
g = gr.Digraph()
g.edge("ability", "educ")
g.edge("ability", "wage")
g.edge("educ", "wage")
g.edge("qob", "educ")
g
분석에는 10년마다 실시하는 인구 총조사 세 차례분의 데이터를 사용합니다. Angrist and Krueger가 도구변수 논문에서 쓴 것과 같은 데이터입니다. 이 데이터셋에는 결과변수인 로그 임금과 처치변수인 교육 연수가 들어 있습니다. 도구변수인 출생 분기와, 출생 연도와 출생 주(state) 같은 추가 통제변수도 있습니다.
data = pd.read_csv("./data/ak91.csv")
data.head()| log_wage | years_of_schooling | year_of_birth | quarter_of_birth | state_of_birth | |
|---|---|---|---|---|---|
| 0 | 5.790019 | 12.0 | 30.0 | 1.0 | 45.0 |
| 1 | 5.952494 | 11.0 | 30.0 | 1.0 | 45.0 |
| 2 | 5.315949 | 12.0 | 30.0 | 1.0 | 45.0 |
| 3 | 5.595926 | 12.0 | 30.0 | 1.0 | 45.0 |
| 4 | 6.068915 | 12.0 | 30.0 | 1.0 | 37.0 |
출생 분기를 도구변수로 쓰기 전에 이것이 타당한 도구변수인지 확인해야 합니다. 즉, 도구변수의 두 가지 가정이 성립한다고 주장할 수 있어야 합니다.
다행히 첫 번째 가정은 검증할 수 있습니다. \(\mathrm{Cov}(Z, T)\)가 0이 아님을 데이터로 확인할 수 있으니까요. 이 예에서 출생 분기가 앞서 말한 대로 도구변수 역할을 한다면, 그해 마지막 분기에 태어난 사람이 연초에 태어난 사람보다 교육을 조금 더 오래 받았어야 합니다. 이를 확인하는 통계 검정을 돌리기 전에 먼저 데이터를 그려서 직접 눈으로 보겠습니다.
group_data = (data
.groupby(["year_of_birth", "quarter_of_birth"])
[["log_wage", "years_of_schooling"]]
.mean()
.reset_index()
.assign(time_of_birth = lambda d: d["year_of_birth"] + (d["quarter_of_birth"])/4))plt.figure(figsize=(15,6))
plt.plot(group_data["time_of_birth"], group_data["years_of_schooling"], zorder=-1)
for q in range(1, 5):
x = group_data.query(f"quarter_of_birth=={q}")["time_of_birth"]
y = group_data.query(f"quarter_of_birth=={q}")["years_of_schooling"]
plt.scatter(x, y, marker="s", s=200, c=f"C{q}")
plt.scatter(x, y, marker=f"${q}$", s=100, c=f"white")
plt.title("Years of Education by Quarter of Birth (first stage)")
plt.xlabel("Year of Birth")
plt.ylabel("Years of Schooling");
놀랍게도 교육 연수에 분기를 따르는 계절 패턴이 있습니다. 그래프를 보면 1분기에 태어난 사람은 거의 언제나 4분기에 태어난 사람보다 교육 연수가 짧습니다(출생 연도를 고정했을 때 그렇습니다. 일반적으로 나중에 태어난 세대일수록 교육을 더 많이 받으니까요).
좀 더 엄밀하게 확인하려면 1단계를 선형회귀로 돌려 볼 수 있습니다. 먼저 출생 분기를 더미변수로 바꿉니다.
factor_data = data.assign(**{f"q{int(q)}": (data["quarter_of_birth"] == q).astype(int)
for q in data["quarter_of_birth"].unique()})
factor_data.head()| log_wage | years_of_schooling | ... | q3 | q4 | |
|---|---|---|---|---|---|
| 0 | 5.790019 | 12.0 | ... | 0 | 0 |
| 1 | 5.952494 | 11.0 | ... | 0 | 0 |
| 2 | 5.315949 | 12.0 | ... | 0 | 0 |
| 3 | 5.595926 | 12.0 | ... | 0 | 0 |
| 4 | 6.068915 | 12.0 | ... | 0 | 0 |
5 rows × 9 columns
단순하게 하려고 일단 마지막 분기인 q4만 도구변수로 쓰겠습니다. 처치인 교육 연수를 도구변수인 출생 분기에 회귀합니다. 이 회귀를 보면 앞의 그래프에서 본 것처럼 출생 분기가 정말 교육 기간을 늘리는지 알 수 있습니다. 여기서는 출생 연도도 통제해야 하고, 출생 주도 추가 통제변수로 넣겠습니다.
first_stage = smf.ols("years_of_schooling ~ C(year_of_birth) + C(state_of_birth) + q4", data=factor_data).fit()
print("q4 parameter estimate:, ", first_stage.params["q4"])
print("q4 p-value:, ", first_stage.pvalues["q4"])q4 parameter estimate:, 0.10085809272786578
q4 p-value:, 5.4648294166163044e-15그해 마지막 분기에 태어난 사람은 다른 분기에 태어난 사람보다 평균적으로 교육 연수가 0.1년 더 긴 것으로 보입니다. p 값은 0에 가깝습니다. 출생 분기가 교육 연수를 늘리거나 줄이는지에 대한 수사는 이것으로 종결입니다.
안타깝게도 두 번째 도구변수 조건은 검증할 수 없습니다. 그 조건이 성립한다고 주장할 수 있을 뿐입니다. 출생 분기가 잠재적 소득에 영향을 주지 않는다는 믿음을 내세울 수는 있습니다. 다시 말해 태어난 시기는 개인의 능력이나, 교육에 미치는 효과 외에 소득 차이를 낳을 수 있는 다른 어떤 요인도 나타내지 않는다는 것입니다. 이를 설득하는 좋은 방법은 소득에 미치는 영향만 놓고 보면 출생 분기가 무작위 배정이나 다름없다고 말하는 것입니다. (무작위는 아닙니다. 사람들이 늦여름이나 명절 즈음에 아이를 갖는 경향이 있다는 증거가 있습니다. 하지만 이 패턴이 교육 말고 다른 경로로 소득에 영향을 줄 그럴듯한 이유는 떠오르지 않습니다.)
배제 제약을 옹호했으니 이제 축약형을 돌려 보겠습니다. 축약형의 목적은 도구변수가 결과에 어떤 영향을 주는지 알아내는 것입니다. 가정에 따라 이 영향은 모두 처치에 미치는 효과에서 비롯되므로, 축약형을 보면 처치가 결과에 어떻게 영향을 주는지 어느 정도 알 수 있습니다. 이번에도 본격적으로 회귀분석을 하기 전에 그래프로 먼저 살펴보겠습니다.
plt.figure(figsize=(15,6))
plt.plot(group_data["time_of_birth"], group_data["log_wage"], zorder=-1)
for q in range(1, 5):
x = group_data.query(f"quarter_of_birth=={q}")["time_of_birth"]
y = group_data.query(f"quarter_of_birth=={q}")["log_wage"]
plt.scatter(x, y, marker="s", s=200, c=f"C{q}")
plt.scatter(x, y, marker=f"${q}$", s=100, c=f"white")
plt.title("Average Weekly Wage by Quarter of Birth (reduced form)")
plt.xlabel("Year of Birth")
plt.ylabel("Log Weekly Earnings");
출생 분기에 따른 소득에서도 계절 패턴이 보입니다. 연말에 태어난 사람이 연초에 태어난 사람보다 소득이 조금 더 높습니다. 이 가설을 검정하기 위해 이번에는 로그 임금을 도구변수 q4에 회귀합니다. 1단계와 같은 추가 통제변수도 넣습니다.
reduced_form = smf.ols("log_wage ~ C(year_of_birth) + C(state_of_birth) + q4", data=factor_data).fit()
print("q4 parameter estimate:, ", reduced_form.params["q4"])
print("q4 p-value:, ", reduced_form.pvalues["q4"])q4 parameter estimate:, 0.008603484260140015
q4 p-value:, 0.0014949127183659584이번에도 유의한 결과가 나왔습니다. 그해 마지막 분기에 태어난 사람은 평균적으로 임금이 0.8% 더 높습니다. p 값이 이전만큼 0에 가깝지는 않지만 0.0015에 불과하니 여전히 꽤 유의합니다.
축약형과 1단계를 모두 구했으니, 이제 1단계의 효과로 축약형을 환산할 수 있습니다. 1단계 계수가 대략 0.1이었으므로, 축약형 계수에 거의 10을 곱하는 셈이 됩니다. 이렇게 하면 평균 인과효과에 대한 편향 없는 도구변수 추정값을 얻습니다.
\[ \mathrm{ATE}_{IV} = \dfrac{\text{Reduced Form}}{\text{1st Stage}} \]
원문은 ‘편향 없는’이라고 썼지만, 뒤의 ’약한 도구변수’ 절에서 보듯 도구변수 추정량은 일치성이 있을 뿐 유한 표본에서는 편향될 수 있습니다.
reduced_form.params["q4"] / first_stage.params["q4"]0.08530286492085315교육을 1년 더 받을 때마다 임금이 8% 오를 것으로 기대할 수 있다는 뜻입니다.
도구변수 추정값을 구하는 또 다른 방법은 2단계 최소제곱법(two-stage least squares, 2SLS)입니다. 이 방법에서는 앞에서처럼 1단계를 돌린 다음, 처치변수를 1단계의 적합값으로 바꾸어 2단계를 돌립니다.
\[ \mathrm{educ}_i = \gamma_0 + \gamma_1 \times \mathrm{q4}_i + \gamma_2\times \mathrm{yob}_i + \gamma_3\times \mathrm{sob}_i + v_i \]
\[ \log(\mathrm{wage})_i = \beta_0 + \beta_1\times \mathrm{educ}_i + \beta_2\times \mathrm{yob}_i + \beta_3\times \mathrm{sob}_i + u_i \]
\[ \log(\mathrm{wage})_i = \beta_0 + \beta_1 [\gamma_0 + \gamma_1 \times \mathrm{q4}_i + \gamma_2\times \mathrm{yob}_i + \gamma_3\times \mathrm{sob}_i + v_i ] + \beta_2\times \mathrm{yob}_i + \beta_3\times \mathrm{sob}_i + u_i \]
한 가지 주의할 점은 도구변수 분석에서 2단계에 넣는 추가 통제변수는 모두 1단계에도 넣어야 한다는 것입니다.
iv_by_hand = smf.ols("log_wage ~ C(year_of_birth) + C(state_of_birth) + years_of_schooling_fitted",
data=factor_data.assign(years_of_schooling_fitted=first_stage.fittedvalues)).fit()
iv_by_hand.params["years_of_schooling_fitted"]0.0853028649208674보다시피 매개변수가 정확히 같습니다. 도구변수를 이렇게 두 번째 방식으로 보면 직관을 얻는 데 도움이 됩니다. 2SLS에서 1단계는 누락 변수 편향을 걷어 낸 새로운 버전의 처치를 만듭니다. 그런 다음 이 정제된 처치, 즉 1단계의 적합값을 선형회귀에 사용합니다.
하지만 실무에서는 도구변수를 손으로 계산하지 않습니다. 번거로워서가 아니라, 이렇게 돌린 2단계에서 얻는 표준오차가 조금 어긋나기 때문입니다. 대신 늘 컴퓨터가 일을 하도록 맡겨야 합니다. 파이썬에서는 linearmodels 라이브러리를 쓰면 2SLS를 올바르게 돌릴 수 있습니다.
2SLS의 공식 표기는 조금 다릅니다. 공식 안에 1단계를 [ ] 사이에 넣어야 합니다. 이 예에서는 years_of_schooling ~ q4를 넣습니다. 추가 통제변수는 1단계에 넣을 필요가 없습니다. 2단계에 넣으면 컴퓨터가 알아서 1단계에도 넣어 주기 때문입니다. 그래서 year_of_birth와 state_of_birth는 1단계 공식 바깥에 넣습니다.
def parse(model, exog="years_of_schooling"):
param = model.params[exog]
se = model.std_errors[exog]
p_val = model.pvalues[exog]
print(f"Parameter: {param}")
print(f"SE: {se}")
print(f"95 CI: {(-1.96*se,1.96*se) + param}")
print(f"P-value: {p_val}")
formula = 'log_wage ~ 1 + C(year_of_birth) + C(state_of_birth) + [years_of_schooling ~ q4]'
iv2sls = IV2SLS.from_formula(formula, factor_data).fit()
parse(iv2sls)Parameter: 0.0853028649580665
SE: 0.02554081280759766
95 CI: [0.03524287 0.13536286]
P-value: 0.0008381914592261452이번에도 매개변수가 앞에서 구한 값과 정확히 같습니다. 덤으로 이제는 올바른 표준오차까지 얻었습니다. 이를 바탕으로 교육을 1년 더 받으면 임금이 평균 8.5% 오를 것으로 기대할 수 있다고 말할 수 있습니다.
2SLS를 컴퓨터로 돌릴 때의 또 다른 장점은 도구변수를 여러 개 넣기 쉽다는 것입니다. 이 예에서는 출생 분기 더미를 모두 교육 연수의 도구변수로 쓰겠습니다.
formula = 'log_wage ~ 1 + C(year_of_birth) + C(state_of_birth) + [years_of_schooling ~ q1+q2+q3]'
iv_many_zs = IV2SLS.from_formula(formula, factor_data).fit()
parse(iv_many_zs)Parameter: 0.10769370488924324
SE: 0.0195571490089523
95 CI: [0.06936169 0.14602572]
P-value: 3.657974678716869e-08더미 3개를 모두 쓰면 추정된 교육 수익률은 0.1이 됩니다. 교육을 1년 더 받을 때마다 소득이 평균 10% 늘어날 것으로 기대할 수 있다는 뜻입니다. 이를 전통적인 최소제곱법(OLS) 추정값과 비교해 보겠습니다. 2SLS를 다시 쓰되 이번에는 1단계 없이 돌리면 됩니다.
formula = "log_wage ~ years_of_schooling + C(state_of_birth) + C(year_of_birth) + C(quarter_of_birth)"
ols = IV2SLS.from_formula(formula, data=data).fit()
parse(ols)Parameter: 0.067325728176578
SE: 0.00038839984390487046
95 CI: [0.06656446 0.06808699]
P-value: 0.0OLS로 추정한 교육 수익률이 2SLS보다 낮습니다. 처음 생각했던 것만큼 누락 변수 편향이 크지 않을 수도 있다는 뜻입니다. 신뢰구간도 눈여겨보세요. 2SLS의 신뢰구간이 OLS 추정값의 신뢰구간보다 훨씬 넓습니다. 이 부분을 좀 더 살펴보겠습니다.
도구변수를 다룰 때는 평균 처치효과(ATE)를 간접적으로 추정하고 있다는 사실을 기억해야 합니다. 추정값은 1단계와 2단계 모두에 좌우됩니다. 처치가 결과에 미치는 영향이 정말 강하다면 2단계도 강할 것입니다. 하지만 1단계가 약하면 2단계가 아무리 강해도 소용이 없습니다. 1단계가 약하다는 것은 도구변수와 처치의 상관관계가 아주 작다는 뜻입니다. 그러니 도구변수로부터 처치에 대해 알아낼 수 있는 것이 별로 없습니다.
도구변수 표준오차 공식은 조금 복잡하고 직관적이지 않으므로, 이 문제를 이해하기 위해 다른 방법을 써 보겠습니다. 결과 Y에 2.0의 효과를 주는 처치 T, 관측되지 않은 교란 요인 U, 추가 통제변수 X가 있는 데이터를 시뮬레이션합니다. 또 1단계의 강도가 서로 다른 도구변수도 여러 개 만듭니다.
\[ \begin{align} X \sim & N(0, 2^2)\\ U \sim & N(0, 2^2)\\ T \sim & N(1+0.5U, 5^2)\\ Y \sim & N(2+ X - 0.5U + 2T, 5^2)\\ Z \sim & N(T, \sigma^2) \text{ for }\sigma^2 \text{ in 0.1 to 100} \end{align} \]
np.random.seed(12)
n = 10000
X = np.random.normal(0, 2, n) # observable variable
U = np.random.normal(0, 2, n) # unobservable (omitted) variable
T = np.random.normal(1 + 0.5*U, 5, n) # treatment
Y = np.random.normal(2 + X - 0.5*U + 2*T, 5, n) # outcome
stddevs = np.linspace(0.1, 100, 50)
Zs = {f"Z_{z}": np.random.normal(T, s, n) for z, s in enumerate(stddevs)} # instruments with decreasing \mathrm{Cov}(Z, T)
sim_data = pd.DataFrame(dict(U=U, T=T, Y=Y)).assign(**Zs)
sim_data.head()| U | T | ... | Z_48 | Z_49 | |
|---|---|---|---|---|---|
| 0 | 2.696148 | 8.056988 | ... | -117.798705 | -13.485292 |
| 1 | 2.570240 | 0.245067 | ... | -209.727577 | -70.792948 |
| 2 | 0.664741 | 5.597510 | ... | 60.562232 | 47.619414 |
| 3 | 1.037725 | 0.493532 | ... | 78.136513 | -108.322304 |
| 4 | -2.590591 | -6.263014 | ... | 78.776566 | -80.547214 |
5 rows × 53 columns
혹시 모르니 Z와 T의 상관관계가 실제로 줄어드는지 확인해 보겠습니다.
corr = (sim_data.corr()["T"]
[lambda d: d.index.str.startswith("Z")])
corr.head()Z_0 0.999807
Z_1 0.919713
Z_2 0.773434
Z_3 0.634614
Z_4 0.523719
Name: T, dtype: float64이제 도구변수마다 도구변수 모델을 하나씩 돌려서 ATE 추정값과 표준오차를 모읍니다.
se = []
ate = []
for z in range(len(Zs)):
formula = f'Y ~ 1 + X + [T ~ Z_{z}]'
iv = IV2SLS.from_formula(formula, sim_data).fit()
se.append(iv.std_errors["T"])
ate.append(iv.params["T"])plot_data = pd.DataFrame(dict(se=se, ate=ate, corr=corr)).sort_values(by="corr")
plt.scatter(plot_data["corr"], plot_data["se"])
plt.xlabel("Corr(Z, T)")
plt.ylabel("IV Standard Error");
plt.title("Variance of the IV Estimates by 1st Stage Strength");
plt.scatter(plot_data["corr"], plot_data["ate"])
plt.fill_between(plot_data["corr"],
plot_data["ate"]+1.96*plot_data["se"],
plot_data["ate"]-1.96*plot_data["se"], alpha=.5)
plt.xlabel("Corr(Z, T)")
plt.ylabel("$\hat{ATE}$");
plt.title("IV ATE Estimates by 1st Stage Strength");
위 그래프에서 보듯이 T와 Z의 상관관계가 약하면 추정값이 크게 요동칩니다. 상관관계가 낮을수록 표준오차도 크게 커지기 때문입니다.
또 하나 주목할 점은 2SLS가 편향되어 있다는 것입니다! 상관관계가 높아도 매개변수 추정값은 실제 ATE인 2.0에 이르지 못합니다. 심지어 2.0은 95% 신뢰구간 안에 들지도 않습니다! 2SLS는 일치성(consistency)만 갖습니다. 표본 크기가 충분히 크면 실제 매개변수 값에 가까워진다는 뜻입니다. 하지만 얼마나 커야 충분한지는 알 수 없습니다. 이 편향이 어떻게 움직이는지 이해하려면 몇 가지 경험칙에 기댈 수밖에 없습니다.
2SLS는 OLS 쪽으로 편향됩니다. OLS에 음의 편향이나 양의 편향이 있으면 2SLS에도 같은 편향이 있다는 뜻입니다. 누락 변수가 있을 때 OLS는 일치성이 없지만 2SLS는 적어도 일치성이 있다는 것이 2SLS의 장점입니다. 위 예에서 관측되지 않은 U는 결과에 음의 영향을 주지만 처치와는 양의 상관관계가 있으므로 음의 편향이 생깁니다. 그래서 ATE 추정값이 실제 값보다 낮게 나오는 것입니다(음의 편향).
도구변수를 많이 넣을수록 편향이 커집니다. 도구변수를 너무 많이 넣으면 2SLS는 점점 OLS와 비슷해집니다.
이 편향이 어떻게 움직이는지 아는 것과 더불어, 마지막으로 도구변수를 쓸 때 흔히 하는 실수를 피하라고 조언하고 싶습니다.
도구변수를 손으로 계산하는 것. 앞서 보았듯이 손으로 계산하면 매개변수 추정값은 맞더라도 표준오차가 틀립니다. 완전히 엉뚱한 값이 나오지는 않겠지만, 소프트웨어를 쓰면 올바른 표준오차를 얻을 수 있는데 굳이 손으로 할 이유가 있을까요?
1단계에 OLS가 아닌 다른 방법을 쓰는 것. 많은 데이터 과학자가 도구변수를 접하고는 더 잘할 수 있다고 생각합니다. 예를 들어 처치가 더미변수인 것을 보고 1단계를 로지스틱 회귀로 바꾸려고 합니다. 더미변수를 예측하는 거니까 그게 맞지 않느냐는 것입니다. 문제는 이것이 명백히 틀렸다는 것입니다. 도구변수의 일치성은 잔차의 직교성이라는, OLS만 줄 수 있는 성질에 기대고 있습니다. 그래서 1단계에 OLS가 아닌 방법을 쓰면 편향된 결과가 나옵니다. (참고로 도구변수에 머신러닝을 쓰는 최신 기법도 있지만, 그 결과는 좋게 봐줘도 미심쩍은 수준입니다.)
이 장에서는 도구변수가 있을 때 누락 변수 편향을 어떻게 우회할 수 있는지 시간을 들여 살펴보았습니다. 도구변수는 처치와 상관관계가 있으면서(1단계가 있음) 처치를 통해서만 결과에 영향을 주는(배제 제약) 변수입니다. 교육이 임금에 미치는 효과를 추정하기 위해 출생 분기를 도구변수로 쓰는 예를 보았습니다.
이어서 도구변수로 인과효과를 추정하는 방법, 즉 2SLS의 작동 방식을 자세히 들여다보았습니다. 도구변수가 만능 해결책이 아니라는 것도 배웠습니다. 1단계가 약하면 상당히 골치 아파질 수 있습니다. 또 2SLS는 일치성이 있긴 하지만 여전히 인과효과를 편향되게 추정하는 방법입니다.
저는 이 책 전체를 Joshua Angrist, Alberto Abadie, Christopher Walters의 훌륭한 계량경제학 강의에 바치는 헌사라고 생각합니다. 이 책에 담긴 아이디어 대부분은 이들이 미국경제학회(American Economic Association)에서 한 강의에서 가져왔습니다. 힘든 2020년을 버틸 수 있는 것도 이 강의들을 보는 덕분입니다. * Cross-Section Econometrics * Mastering Mostly Harmless Econometrics
Angrist의 멋진 책들도 참고 자료로 소개하고 싶습니다. 이 책들 덕분에 계량경제학, 이들 표현으로는 ’Metrics가 아주 유용할 뿐 아니라 무척 재미있다는 것을 알게 되었습니다.
마지막 참고 자료는 Miguel Hernan과 Jamie Robins의 책입니다. 까다로운 인과적 질문에 답해야 할 때마다 믿음직한 동반자가 되어 주었습니다.
Causal Inference for the Brave and True는 과학의 통계학이라 할 수 있는 인과추론을 다루는 오픈소스 자료입니다. 파이썬 기반의 무료 소프트웨어만 사용하며, 돈이 없어도, 배경지식이 많지 않아도 누구나 볼 수 있는 자료를 목표로 합니다. 이 책이 도움이 되었고 후원하고 싶다면 Patreon을 방문해 주세요. 금전적으로 기여하기 어렵다면 오타를 고치거나, 수정할 부분을 제안하거나, 이해하기 어려웠던 부분에 의견을 주는 것으로도 도울 수 있습니다. 책의 저장소에서 이슈를 열어 주세요. 마지막으로 이 내용이 마음에 들었다면 도움이 될 만한 사람들과 공유하고 GitHub에서 스타를 눌러 주세요.
한국어판의 오역이나 오탈자는 한국어판 이슈 페이지에 알려 주세요. 한국어판이 도움이 되었다면 한국어판 저장소에도 스타를 부탁드립니다.