import matplotlib
if not hasattr(matplotlib.RcParams, "_get"):
matplotlib.RcParams._get = dict.get
Regression Discontinuity Design (RDD)#
멕시코는 1990년대 후반 극심한 빈곤 문제에 직면해 있었고, 1997년 조건부 현금지원 프로그램인 Progresa를 도입했습니다.
pov_index 라는 빈곤 인덱스가 0 이상이면 무조건 현금지원을 받고, 0 미만이면 받지 못합니다. (Sharp RDD 조건 만족)
Set up#
!pip install rdd rdrobust lightgbm
import pandas as pd
from sklearn.linear_model import LinearRegression, LassoCV
from sklearn.ensemble import RandomForestRegressor
from sklearn.model_selection import KFold
from lightgbm import LGBMRegressor
import matplotlib.pyplot as plt
import numpy as np
import patsy
from rdd.rdd import optimal_bandwidth
from rdrobust import rdrobust
Requirement already satisfied: rdd in /opt/homebrew/Caskroom/miniforge/base/lib/python3.12/site-packages (0.0.3)
Requirement already satisfied: rdrobust in /opt/homebrew/Caskroom/miniforge/base/lib/python3.12/site-packages (1.3.0)
Requirement already satisfied: pandas in /opt/homebrew/Caskroom/miniforge/base/lib/python3.12/site-packages (from rdd) (2.2.3)
Requirement already satisfied: numpy in /opt/homebrew/Caskroom/miniforge/base/lib/python3.12/site-packages (from rdd) (2.1.3)
Requirement already satisfied: statsmodels in /opt/homebrew/Caskroom/miniforge/base/lib/python3.12/site-packages (from rdd) (0.14.5)
Requirement already satisfied: scipy in /opt/homebrew/Caskroom/miniforge/base/lib/python3.12/site-packages (from rdrobust) (1.15.2)
Requirement already satisfied: scikit-learn>=1.2.0 in /opt/homebrew/Caskroom/miniforge/base/lib/python3.12/site-packages (from rdrobust) (1.5.2)
Requirement already satisfied: plotnine>=0.13.6 in /opt/homebrew/Caskroom/miniforge/base/lib/python3.12/site-packages (from rdrobust) (0.15.0)
Requirement already satisfied: matplotlib>=3.9.2 in /opt/homebrew/Caskroom/miniforge/base/lib/python3.12/site-packages (from rdrobust) (3.10.1)
Requirement already satisfied: contourpy>=1.0.1 in /opt/homebrew/Caskroom/miniforge/base/lib/python3.12/site-packages (from matplotlib>=3.9.2->rdrobust) (1.3.1)
Requirement already satisfied: cycler>=0.10 in /opt/homebrew/Caskroom/miniforge/base/lib/python3.12/site-packages (from matplotlib>=3.9.2->rdrobust) (0.12.1)
Requirement already satisfied: fonttools>=4.22.0 in /opt/homebrew/Caskroom/miniforge/base/lib/python3.12/site-packages (from matplotlib>=3.9.2->rdrobust) (4.56.0)
Requirement already satisfied: kiwisolver>=1.3.1 in /opt/homebrew/Caskroom/miniforge/base/lib/python3.12/site-packages (from matplotlib>=3.9.2->rdrobust) (1.4.8)
Requirement already satisfied: packaging>=20.0 in /opt/homebrew/Caskroom/miniforge/base/lib/python3.12/site-packages (from matplotlib>=3.9.2->rdrobust) (24.1)
Requirement already satisfied: pillow>=8 in /opt/homebrew/Caskroom/miniforge/base/lib/python3.12/site-packages (from matplotlib>=3.9.2->rdrobust) (11.1.0)
Requirement already satisfied: pyparsing>=2.3.1 in /opt/homebrew/Caskroom/miniforge/base/lib/python3.12/site-packages (from matplotlib>=3.9.2->rdrobust) (3.2.1)
Requirement already satisfied: python-dateutil>=2.7 in /opt/homebrew/Caskroom/miniforge/base/lib/python3.12/site-packages (from matplotlib>=3.9.2->rdrobust) (2.9.0.post0)
Requirement already satisfied: pytz>=2020.1 in /opt/homebrew/Caskroom/miniforge/base/lib/python3.12/site-packages (from pandas->rdd) (2025.1)
Requirement already satisfied: tzdata>=2022.7 in /opt/homebrew/Caskroom/miniforge/base/lib/python3.12/site-packages (from pandas->rdd) (2025.1)
Requirement already satisfied: mizani~=0.14.0 in /opt/homebrew/Caskroom/miniforge/base/lib/python3.12/site-packages (from plotnine>=0.13.6->rdrobust) (0.14.2)
Requirement already satisfied: joblib>=1.2.0 in /opt/homebrew/Caskroom/miniforge/base/lib/python3.12/site-packages (from scikit-learn>=1.2.0->rdrobust) (1.4.2)
Requirement already satisfied: threadpoolctl>=3.1.0 in /opt/homebrew/Caskroom/miniforge/base/lib/python3.12/site-packages (from scikit-learn>=1.2.0->rdrobust) (3.5.0)
Requirement already satisfied: patsy>=0.5.6 in /opt/homebrew/Caskroom/miniforge/base/lib/python3.12/site-packages (from statsmodels->rdd) (1.0.1)
Requirement already satisfied: six>=1.5 in /opt/homebrew/Caskroom/miniforge/base/lib/python3.12/site-packages (from python-dateutil>=2.7->matplotlib>=3.9.2->rdrobust) (1.17.0)
df = pd.read_csv("https://raw.githubusercontent.com/CausalAIBook/MetricsMLNotebooks/main/data/progresa.csv",
index_col=0)
df = df.dropna()
df.rename(columns={"index": "pov_index"}, inplace=True)
print("Shape of Data:")
print(df.shape)
print("Variable Names:")
print(df.columns)
df.head()
Shape of Data:
(1944, 27)
Variable Names:
Index(['hhpiso', 'hhrooms', 'hhwater', 'hhwaterin', 'hhbano', 'hhownhouse',
'hhsize', 'hhelect', 'clus', 'headmale', 'headage', 'heademp',
'wifeage', 'wifeeduc', 'headeduc', 'child_0to5', 'boy_0to5',
'pov_index', 'conspcfood_t0', 'conspcfood_t1', 'conspcfood_t2',
'conspcnonfood_t0', 'conspcnonfood_t1', 'conspcnonfood_t2', 'conspc_t0',
'conspc_t1', 'conspc_t2'],
dtype='object')
| hhpiso | hhrooms | hhwater | hhwaterin | hhbano | hhownhouse | hhsize | hhelect | clus | headmale | ... | pov_index | conspcfood_t0 | conspcfood_t1 | conspcfood_t2 | conspcnonfood_t0 | conspcnonfood_t1 | conspcnonfood_t2 | conspc_t0 | conspc_t1 | conspc_t2 | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 1 | 1.0 | 1.0 | 1.0 | 1.0 | 1.0 | 0.0 | 5 | 1 | 1 | 1 | ... | 1.5410 | 260.223999 | 374.071991 | 379.208008 | 106.186668 | 183.149323 | 57.127998 | 366.410675 | 557.221313 | 436.335999 |
| 2 | 1.0 | 1.0 | 1.0 | 0.0 | 0.0 | 1.0 | 5 | 1 | 1 | 1 | ... | -0.6615 | 265.359985 | 574.375977 | 195.167999 | 87.266663 | 160.013336 | 112.136002 | 352.626648 | 734.389282 | 307.304016 |
| 3 | 1.0 | 1.0 | 0.0 | 0.0 | 1.0 | 1.0 | 4 | 1 | 1 | 1 | ... | -0.3240 | 604.549988 | 433.350006 | 367.010010 | 348.813324 | 263.083344 | 123.323334 | 953.363281 | 696.433350 | 490.333344 |
| 4 | 0.0 | 1.0 | 0.0 | 0.0 | 1.0 | 0.0 | 6 | 1 | 1 | 1 | ... | 1.1500 | 233.259995 | 281.766663 | 260.468567 | 100.877777 | 140.233337 | 36.325714 | 334.137756 | 422.000000 | 296.794281 |
| 6 | 1.0 | 1.0 | 1.0 | 0.0 | 1.0 | 1.0 | 4 | 1 | 1 | 1 | ... | -0.8260 | 280.339996 | 185.110001 | 364.869995 | 158.416656 | 201.183334 | 68.000000 | 438.756653 | 386.293335 | 432.869995 |
5 rows × 27 columns
처치 변수: pov_index (빈곤 지수: 프로그램 참여 자격 결정 기준, RDD의 running variable)
결과 변수: 소비 (Consumption Outcomes)
conspcfood_t0/t1/t2: 식료품 소비 (baseline, 1년 후, 2년 후)
conspcnonfood_t0/t1/t2: 비식료품 소비 (baseline, 1년 후, 2년 후)
conspc_t0/t1/t2: 총 소비 (baseline, 1년 후, 2년 후)
# 처치 변수 생성 (pov_index >= 0인 경우 처치)
df['treatment'] = (df['pov_index'] >= 0).astype(int)
print("처치군 vs 통제군:")
print(df['treatment'].value_counts())
처치군 vs 통제군:
treatment
1 1395
0 549
Name: count, dtype: int64
# 변화량(delta) 변수들 생성
df['conspcfood_delta_t1'] = df['conspcfood_t1'] - df['conspcfood_t0']
df['conspcfood_delta_t2'] = df['conspcfood_t2'] - df['conspcfood_t0']
df['conspcnonfood_delta_t1'] = df['conspcnonfood_t1'] - df['conspcnonfood_t0']
df['conspcnonfood_delta_t2'] = df['conspcnonfood_t2'] - df['conspcnonfood_t0']
df['conspc_delta_t1'] = df['conspc_t1'] - df['conspc_t0']
df['conspc_delta_t2'] = df['conspc_t2'] - df['conspc_t0']
Basic RDD#
result = []
outcomes = [
('conspcfood_delta_t1', 'Food Consumption (t1-t0)'),
('conspcfood_delta_t2', 'Food Consumption (t2-t0)'),
('conspcnonfood_delta_t1', 'Non-Food Consumption (t1-t0)'),
('conspcnonfood_delta_t2', 'Non-Food Consumption (t2-t0)'),
('conspc_delta_t1', 'Total Consumption (t1-t0)'),
('conspc_delta_t2', 'Total Consumption (t2-t0)')
]
for (outcome_var, label) in outcomes:
rdd_result = rdrobust(y=df[outcome_var], x=df['pov_index'], c=0, masspoints="off")
# 일관되게 Robust estimates 사용 (index=2)
coef = rdd_result.coef.iloc[2].values[0] # Robust coefficient
p_value = rdd_result.pv.iloc[2].values[0] # Robust p-value
ci_bc_lower = rdd_result.ci.loc["Robust", "CI Lower"]
ci_bc_upper = rdd_result.ci.loc["Robust", "CI Upper"]
result.append([coef, ci_bc_lower, ci_bc_upper, p_value])
res_dataframe = pd.DataFrame(result,
columns=["LATE", "CI Lower", "CI Upper", "p-value"],
index=["Food T_1", "Food T_2", "Non-Food T_1",
"Non-Food T_2", "Total T_1", "Total T_2"])
# 결과 포맷팅
res_dataframe = res_dataframe.round(4)
res_dataframe
| LATE | CI Lower | CI Upper | p-value | |
|---|---|---|---|---|
| Food T_1 | -24.6361 | -71.3140 | 22.0417 | 0.3009 |
| Food T_2 | 57.3015 | -36.8738 | 151.4767 | 0.2330 |
| Non-Food T_1 | 10.4383 | -34.5067 | 55.3833 | 0.6490 |
| Non-Food T_2 | 71.7462 | 6.4437 | 137.0487 | 0.0313 |
| Total T_1 | -12.2731 | -89.2205 | 64.6743 | 0.7546 |
| Total T_2 | 129.8931 | -19.0838 | 278.8701 | 0.0875 |
RDD with DML#
여러 공변량에 의한 교란 효과를 추가로 제거하여 Robustness Check를 수행할 수 있습니다.
아래 코드는 CausalAIBook/MetricsMLNotebooks 에서 가져오고, second stage 코드만 수정하였습니다. (rdrobust 라이브러리의 내부 함수 확인을 통해 일관성 있고 사용하기 편리하도록 조정하였음)
def first_stage(df_ml, b_covs, h_fs, Z_lasso, Kf=5, random_seed=123):
'''
df_ml: dataframe
b_covs: which columns of the data frame will be used as baseline covariates
h_fs: bandwidth around discontinuity for training points
Z_lasso: dataframe with extra baseline covariates in high-dimensional specification
'''
# Set up the cross-fitting
n = df_ml.shape[0]
# Matrix to store eta predictions
eta_fit = np.empty((n, 5))
# Create vector of observations to be considered in the first stage model
weights = np.abs(df_ml.X) < h_fs
for train, test in KFold(shuffle=True, n_splits=Kf, random_state=random_seed).split(df_ml.X, df_ml.Y):
df_train = df_ml.iloc[train]
treated_train = (df_train.X > 0) & (weights.iloc[train] > 0)
control_train = (df_train.X < 0) & (weights.iloc[train] > 0)
data_treated = df_train[treated_train]
data_control = df_train[control_train]
data_fold = df_ml.iloc[test]
rf1 = RandomForestRegressor(max_features=4, n_estimators=1000, random_state=random_seed)
rf1.fit(data_treated[b_covs], data_treated.Y)
rf0 = RandomForestRegressor(max_features=4, n_estimators=1000, random_state=random_seed)
rf0.fit(data_control[b_covs], data_control.Y)
eta_fit[test, 0] = (rf1.predict(data_fold[b_covs]) + rf0.predict(data_fold[b_covs])) / 2
lgbm1 = LGBMRegressor(verbosity=-1, random_state=random_seed)
lgbm1.fit(data_treated[b_covs], data_treated.Y)
lgbm0 = LGBMRegressor(verbosity=-1, random_state=random_seed)
lgbm0.fit(data_control[b_covs], data_control.Y)
eta_fit[test, 1] = (lgbm1.predict(data_fold[b_covs]) + lgbm0.predict(data_fold[b_covs])) / 2
lm1 = LinearRegression()
lm1.fit(data_treated[b_covs], data_treated.Y)
lm0 = LinearRegression()
lm0.fit(y=data_control.Y, X=data_control[b_covs])
eta_fit[test, 2] = (lm1.predict(data_fold[b_covs]) + lm0.predict(data_fold[b_covs])) / 2
las_base1 = LassoCV(random_state=random_seed)
las_base1.fit(data_treated[b_covs], data_treated.Y)
las_base0 = LassoCV(random_state=random_seed)
las_base0.fit(data_control[b_covs], data_control.Y)
eta_fit[test, 3] = (las_base1.predict(data_fold[b_covs]) + las_base0.predict(data_fold[b_covs])) / 2
X_flex_treated = pd.concat([Z_lasso.loc[data_treated.index], data_treated[b_covs]], axis=1)
X_flex_control = pd.concat([Z_lasso.loc[data_control.index], data_control[b_covs]], axis=1)
X_flex_fold = pd.concat([Z_lasso.loc[data_fold.index], data_fold[b_covs]], axis=1)
X_flex_treated.columns = X_flex_treated.columns.astype(str)
X_flex_control.columns = X_flex_control.columns.astype(str)
X_flex_fold.columns = X_flex_fold.columns.astype(str)
las_flex1 = LassoCV(random_state=random_seed)
las_flex1.fit(X_flex_treated, data_treated.Y)
las_flex0 = LassoCV(random_state=random_seed)
las_flex0.fit(X_flex_control, data_control.Y)
eta_fit[test, 4] = (las_flex1.predict(X_flex_fold) + las_flex0.predict(X_flex_fold)) / 2
return eta_fit
methods = ["Random Forest", "Gradient Boosting", "Linear Regression",
"Lasso Baseline", "Lasso Flexible"]
def second_stage(df_ml, eta_fit, methods, level=95):
"""
DML second stage with confidence intervals and p-values
"""
import scipy.stats as stats
adj_results = []
quant = stats.norm.ppf(1 - (1-level/100)/2) # Critical value for CI
for i in range(len(methods)):
M_Y = df_ml.Y - eta_fit[:, i] # Residualized outcome
rd_call = rdrobust(y=M_Y, x=df_ml.X, masspoints="off")
# Robust estimates 추출 (index=2)
coef = rd_call.coef.iloc[2].values[0] # Robust coefficient
se_robust = rd_call.se.iloc[2].values[0] # Robust standard error
# p-value 계산
t_stat = coef / se_robust
p_value = 2 * (1 - stats.norm.cdf(abs(t_stat)))
# Confidence interval
ci_lower = coef - quant * se_robust
ci_upper = coef + quant * se_robust
adj_results.append([coef, ci_lower, ci_upper, p_value])
return adj_results
investigated_outcome = "conspcnonfood_delta_t2"
df_ml = df.rename(columns={"pov_index": "X", investigated_outcome: "Y"})
b_covs = df_ml.columns[[0, 1, 2, 3, 4, 5, 6, 7, 10, 11, 12, 13, 14, 15, 16, 18, 21]]
h_fs = rdrobust(y=df_ml.Y, x=df_ml.X, masspoints="off").bws.values[1, 0]
# Fixed effects for localities
i_fe = pd.get_dummies(df_ml['clus'], drop_first=True)
# Flexible covariates including localities indicators
f_covs = patsy.dmatrix('~ (' + ' + '.join(b_covs) + ')**2', data=df_ml, return_type='dataframe')
# Dropping the intercept column that is automatically added by patsy
f_covs = f_covs.iloc[:, 1:]
Z_lasso = pd.concat([i_fe, f_covs], axis=1)
eta_fit = first_stage(df_ml, b_covs, h_fs, Z_lasso)
adj_frame = pd.DataFrame(second_stage(df_ml, eta_fit, methods),
columns=["LATE", "CI Lower", "CI Upper", "p-value"],
index=methods)
# 포맷팅
adj_frame
| LATE | CI Lower | CI Upper | p-value | |
|---|---|---|---|---|
| Random Forest | 61.363726 | 1.352838 | 121.374613 | 0.045054 |
| Gradient Boosting | 51.080203 | -11.528300 | 113.688707 | 0.109805 |
| Linear Regression | 46.685732 | -12.175647 | 105.547111 | 0.120056 |
| Lasso Baseline | 50.078328 | -9.696978 | 109.853634 | 0.100588 |
| Lasso Flexible | 50.553800 | -9.454359 | 110.561960 | 0.098704 |