目录
LaLonde 数据集
经济学家长期以来假设培训项目可以改善参与者的劳动力市场前景。为了测试（或证明）这一点，国家支持性工作示范项目（National Supported Work，NSW）使用私有和联邦资金联合启动。该项目在1975年至1979年间在美国15个地点实施。该项目为面临经济和社会问题的个人提供了6至18个月的培训（例如接受儿童抚养援助的家庭中的女性、前吸毒者、前罪犯和前青少年犯罪者）。
参与者被随机分配到实验组（支持性工作项目）和对照组。然而，由于研究周期较长，项目初期加入的参与者与后期加入的人具有不同的特征。
因此，为了估计工作项目对未来就业的真实因果效应，需要调整这种协变量偏移。
多年来，这个数据集已成为因果分析的一个常见基准。
最初的研究分析由Robert LaLonde完成，并发表在他1986年的论文《使用实验数据评估培训项目的计量经济学评价》中。
然而，我们将遵循Dehejia和Wahba在其1999年的论文《非实验研究中的因果效应：重新评估培训项目的评价》中所做的基于倾向性的分析。
数据
首先，让我们从Rajeev Dehejia的网页下载数据集。
import pandas as pd
columns = ["training", # Treatment assignment indicator
"age", # Age of participant
"education", # Years of education
"black", # Indicate whether individual is black
"hispanic", # Indicate whether individual is hispanic
"married", # Indicate whether individual is married
"no_degree", # Indicate if individual has no high-school diploma
"re74", # Real earnings in 1974, prior to study participation
"re75", # Real earnings in 1975, prior to study participation
"re78"] # Real earnings in 1978, after study end
treated = pd.read_csv("http://www.nber.org/~rdehejia/data/nswre74_treated.txt",
delim_whitespace=True, header=None, names=columns)
control = pd.read_csv("http://www.nber.org/~rdehejia/data/nswre74_control.txt",
delim_whitespace=True, header=None, names=columns)
lalonde = pd.concat([treated, control], ignore_index=True)
lalonde = lalonde.sample(frac=1.0, random_state=42) # Shuffle
print(lalonde.shape)
lalonde.head()
(445, 10)
我们可以看到数据集中有445个个体，正如Dehejia网页上所描述的。
收入指示变量
按照Gelman等人在他们的arm R库中进行的分析，我们将创建两个指示变量，表示1974年和1975年没有收入。
lalonde = lalonde.join((lalonde[["re74", "re75"]] == 0).astype(int), rsuffix=("=0"))
lalonde.head()
教育年限的因子化
既然学校的年数不应以其数值形式来考虑，我们将将其因子化为指示变量。
lalonde = pd.get_dummies(lalonde, columns=["education"], drop_first=True)
print(lalonde.shape)
lalonde.head()
(445, 24) 
变量选择
最后，我们提取协变量、治疗变量和结果变量。
a = lalonde.pop("training")
y = lalonde.pop("re78")
X = lalonde
X.shape, a.shape, y.shape
((445, 22), (445,), (445,))模型
在定义设计矩阵X之后，我们可以继续定义因果模型。
遵循Dehejia和Wahba基于倾向性的分析精神，我们将使用逆治疗概率加权（IPTW，或IPW）因果模型。 简而言之，这个模型将模拟参与者被分配到就业培训计划的概率，并利用它来模拟两个等规模的人群：一个被分配到该计划的人群，另一个是没有被分配的人群。在这个合成的人群中，我们可以使用1978年的实际收入来估计如果每个人都参加该计划或根本没有人参加会发生什么。
在定义因果模型本身之前，我们需要使用机器学习模型来估计倾向得分——每个参与者被分配到就业培训的概率。 鉴于我们上面准备的设计矩阵，以及我们的干预措施是二元性质的，我们将选择逻辑回归来完成这项任务。
from sklearn.linear_model import LogisticRegression
learner = LogisticRegression(penalty='none', # No regularization, new in scikit-learn 0.21.*
solver='lbfgs',
max_iter=500) # Increaed to achieve convergence with 'lbfgs' solver一旦我们定义了学习器，我们就可以简单地将其插入到因果模型中。
from causallib.estimation import IPW
ipw = IPW(learner)估计因果效应
一旦我们定义了因果模型（是的，只需要这么简单），我们就可以继续估计就业培训对年收入的影响。
首先，我们将拟合我们的因果模型。
其次，我们将预测潜在结果：如果每个人都参加就业培训或根本不参加培训计划，1978年的收入会是多少。第三，我们将使用这两个潜在结果来估计效应：两个潜在结果之间的差值。
ipw.fit(X, a)
outcomes = ipw.estimate_population_outcome(X, a, y)
effect = ipw.estimate_effect(outcomes[1], outcomes[0])
检查潜在结果，我们可以看到如果每个人都参加了该计划（1），研究结束时（1978年）的平均收入为6141美元；而如果每个人都未参加该计划（0），平均收入为4598美元。
outcomes
0.0 4598.414070
1.0 6141.051856
dtype: float64因此，我们可以得出结论，就业培训对收入的平均加性效应（差值）为1543美元。
未经调整的估计
为了进行比较，我们可以问如果不控制混杂因素我们会得出什么样的结论。
如果我们没有调整治疗分配偏差，而是将数据当作来自随机对照试验的数据，结果会怎样？
from causallib.estimation import MarginalOutcomeEstimator
moe = MarginalOutcomeEstimator(None).fit(X, a, y)
outcomes = moe.estimate_population_outcome(X, a, y)
moe.estimate_effect(outcomes[1], outcomes[0])
diff 1794.342404
dtype: float64我们可以看到，当未能调整混杂因素时，我们对就业培训对收入的影响高估了252美元。
这种看似微不足道的偏差占到了我们估计值的大约16%，这很可能是因为研究本身的随机性质造成的。
为了了解不同估计方法之间的差异，我们可以检查在IP加权之前和之后治疗组的平衡性。
我们可以通过Love图来检查这一点，这是一种图形化展示治疗组和未治疗组之间分布距离的方法。它绘制了每个协变量边际上治疗组之间的标准化均值差异。
%matplotlib inline
from causallib.evaluation import evaluate
import matplotlib.pyplot as plt
evaluation_results = evaluate(ipw, X, a, y)
fig, ax = plt.subplots(figsize=(6, 6))
evaluation_results.plot_covariate_balance(kind="love", ax=ax);
我们可以看到，在应用IPW加权之前（橙色三角形），治疗组间的最大绝对平均差异约为0.2个标准差。这很容易导致简单的边缘比较出现偏差（即，简单地比较治疗组与未治疗组的结果），因为我们看到这些组并非边缘相似。
例如，有人可能会认为结果的差异归因于两个治疗组在高中文凭持有者比例上的不同（no_degree）。
由于治疗组在特征分布上有差异，因此在因果推断分析中需要对此进行调整，比如：IPW模型采用的重要性抽样方案。经过平衡处理后（蓝色点），我们发现各组协变量分布的边缘差异平均变得更小，表明IPW模型成功地平衡了数据集。