ebalance：基于熵平衡法的协变量平衡性检验
👇 连享会 · 推文导航 | www.lianxh.cn
🍎 Stata：Stata基础 | Stata绘图 | Stata程序 | Stata新命令 📘 论文：数据处理 | 结果输出 | 论文写作 | 数据分享 💹 计量：回归分析 | 交乘项-调节 | IV-GMM | 时间序列 | 面板数据 | 空间计量 | Probit-Logit | 分位数回归 ⛳ 专题：SFA-DEA | 生存分析 | 爬虫 | 机器学习 | 文本分析 🔃 因果：DID | RDD | 因果推断 | 合成控制法 | PSM-Matching 🔨 工具：工具软件 | Markdown | Python-R-Stata 🎧 课程：公开课-直播 | 计量专题 | 关于连享会
连享会 · 2022 空间计量专题
连享会 · 2022 面板数据因果推断专题
作者：王文韬 (山东大学)
邮箱：190138445@qq.com
目录
1. 引言
2. 熵平衡
2.1 基本含义
2.2 总体平均处理效应与熵平衡
3. Stata 实现
3.1 命令安装和数据介绍
3.2 ebalance 的基本语法
3.3 平衡性检验
3.4 变量交互项的匹配
3.5 NSW 实验效果的无偏估计
4. 结论
5. 参考文献
6. 相关推文
温馨提示： 文中链接在微信中无法生效。请点击底部「阅读原文」。或直接长按/扫描如下二维码，直达原文：
1. 引言
在进行随机实验时，我们需要对处理组和对照组的协变量进行平衡性检验，以保证实验的随机性或外生性，增强回归结果的有效性。Hainmueller (2012) 构建了「熵平衡方法」 (Entropy Balancing)，使研究者能够同时控制处理组与对照组样本协变量多维平衡性，如同时考虑协变量的一阶矩、二阶矩、三阶矩和交叉矩等，进而最大程度上使两组样本在实现精确匹配。
需要强调的是，熵平衡方法并不能完全克服内生性问题。在实际工作中，研究者仍需搭配其他计量工具进行系统检验，以保证回归结果的稳健性和合理性。
本文的主要目的是介绍熵平衡方法理论、及 Stata 中的实现命令 ebalance。
2. 熵平衡
2.1 基本含义
在高中物理课堂上，我们曾学习过 “熵” 的概念。它是热力学中表征物质状态的参数之一，能够度量一个体系的混乱程度。在物理热力学中，熵越大表示该体系中的热量越高。在计量经济学中，“熵平衡” 则是一个基于最大化熵的权重计算方法。
所谓最大化熵，指的是研究者进行随机性实验研究时，在对数据进行预处理的过程中，最大化处理组与对照组的协变量匹配效率。
2.2 总体平均处理效应与熵平衡
在进行随机性实验时，研究者主要关注的结果是总体平均处理效应 (the Population Average Treament Effect on the Treated)，可以由以下公式表示：
其中， 表示处理组实验发生后的效应大小， 表示处理组在没有发生实验后的效应大小。
显然， 是一个反事实指标，在现实中不可能存在对应的数据。为此，学者尽可能寻找与处理组协变量分布相似的对照组样本，来近似估计这一项。
目前，我们熟知的协变量匹配方法包括最近邻匹配 (NNM)、粗化精确匹配 (CEM)、倾向得分匹配 (PSM) 等。以倾向得分匹配为例，可以通过以下表达式计算得到：
其中，， 为倾向得分值。
不过，这些匹配方法面临的一个共同问题是：研究者无法保证联合并平衡所有协变量，难以避免倾向得分模型被错误指定的可能。为了解决这一问题，学者通常利 Logit 回归挑选可以通过平衡性检验的协变量组合。但是，这样做费时费力，同时无法完全保证最终得到高水平的协变量组合。
在上述基础上，Hainmueller (2012) 推导出熵平衡方法，以熵权重 代替 ，得到反事实指标 的表达式如下所示：
权重 通过以下公式决定：
其中，， 表示对照组样本数量。 具有一定约束条件，其核心思想是，给予对照组协变量增加一组矩约束，使其与处理组的协变量相平衡。矩约束包括均值 (一阶矩)、方差 (二阶矩) 和偏度 (三阶矩)。有关最小化熵距离和 的相关数理验证这里就不展开描述了，对此感兴趣的读者可以自行阅读 Hainmueller 的相关文献。
综上所述，Hainmueller (2012) 创造的熵平衡方法为研究者在平衡处理组与对照组的协变量分布上提供了新的科学方法。接下来，我们使用 Hainmueller 和 Xu (2013) 提供的案例及其数据来说明 ebalance 命令在 Stata 中的具体应用。
3. Stata 实现
3.1 命令安装和数据介绍
首先，安装 ebalance 命令，以及下载所需要的数据，具体命令如下：
. net describe ebalance, from(http://fmwww.bc.edu/RePEc/bocode/e)
. ssc install ebalance, all replace
. net get ebalance.pkg, replace // 下载数据到当前工作路径
我们使用的案例数据是 Dehejia 和 Wahba (1999) 提供的「cps1re74.dta」。它包括 185 名参与 NSW (the National Supported Work Demonstration，一个关于补贴资助工作的随机评估) 项目的成员和 15992 名未参与项目的成员 (来自美国社会保障管理档案 CPS-1)。其主要变量如下所示：
. use "cps1re74.dta", clear // 载入数据
. des
Contains data from cps1re74.dta
obs: 16,177 DW Subset of LaLonde Data
vars: 12 6 Oct 2011 14:57
-------------------------------------------------------------
storage display
variable name type format variable label
-------------------------------------------------------------
re78 double %9.0g real earnings 78
treat long %9.0g 1 if treated, 0 control
age long %9.0g age in years
educ long %9.0g years of schooling
black long %9.0g indicator for black
hispan long %9.0g indicator for hispanic
married long %9.0g indicator for married
nodegree long %9.0g indicator for no HS degree
re74 double %9.0g real earnings 74
re75 double %9.0g real earnings 75
u74 double %9.0g indicator for unemployed 74
u75 double %9.0g indicator for unemployed 75
-------------------------------------------------------------
本文衡量政策效果的变量为 re78，根据 Dehejia 和 Wahba (1999) 的研究分析，NSW 这一补贴资助计划平均提高了 1794 美元的收入水平，95% 的置信区间为 [551, 3038]。如果我们直接使用 OLS 回归进行检验，结果如下所示：
. use cps1re74.dta, clear
(DW Subset of LaLonde Data)
. reg re78 treat age-u75
Source | SS df MS Number of obs = 16,177
----------+---------------------------------- F(11, 16165) = 1343.88
Model | 7.2418e+11 11 6.5835e+10 Prob > F = 0.0000
Residual | 7.9190e+11 16,165 48988567.3 R-squared = 0.4777
----------+---------------------------------- Adj R-squared = 0.4773
Total | 1.5161e+12 16,176 93724175.2 Root MSE = 6999.2
---------------------------------------------------------------------------
re78 | Coef. Std. Err. t P>|t| [95% Conf. Interval]
----------+----------------------------------------------------------------
treat | 1067.546 554.060 1.93 0.054 -18.472 2153.564
age | -94.541 6.000 -15.76 0.000 -106.302 -82.780
educ | 175.225 28.697 6.11 0.000 118.977 231.474
black | -811.089 212.849 -3.81 0.000 -1228.296 -393.881
hispan | -230.535 218.610 -1.05 0.292 -659.034 197.965
married | 153.228 142.775 1.07 0.283 -126.626 433.083
nodegree | 342.927 177.878 1.93 0.054 -5.734 691.587
re74 | 0.291 0.013 22.89 0.000 0.266 0.316
re75 | 0.443 0.013 34.35 0.000 0.417 0.468
u74 | 355.556 231.600 1.54 0.125 -98.406 809.519
u75 | -1612.758 239.803 -6.73 0.000 -2082.798 -1142.717
_cons | 5762.180 445.614 12.93 0.000 4888.726 6635.634
---------------------------------------------------------------------------
可以发现，OLS 估计的处理效应为 1067.546，这大大低于真正的政策实施效果。于是我们考虑利用熵平衡对数据进行预处理，看看是否能使我们更准确地评估实际情况。
3.2 ebalance 的基本语法
ebalance 的基本语法如下所示：
ebalance [treat] [covar], targets(numlist)
其中，treat 代表是否为处理组变量，covar 代表需要在熵平衡中进行调整的协变量组合，targets(numlist) 是最为关键的选择指标，numlist 对应协变量组合需要调整的矩阶数。具体来看：
. ebalance treat age black educ, targets(3) //表示对 age、black 和 educ 三个协变量的一阶、二阶和三阶矩进行调整
Data Setup
Treatment variable: treat
Covariate adjustment: age black educ (1st order). age black educ (2nd order). age black educ (3rd order).
Optimizing...
Iteration 1: Max Difference = 580799.347
Iteration 2: Max Difference = 213665.688
Iteration 3: Max Difference = 78604.7628
...
Iteration 15: Max Difference = .420310244
Iteration 16: Max Difference = .037151116
Iteration 17: Max Difference = .008791339
maximum difference smaller than the tolerance level; convergence achieved
Treated units: 185 total of weights: 185
Control units: 15992 total of weights: 185
Before: without weighting
| Treat | Control
| mean variance skewness | mean variance skewness
------+---------------------------+----------------------------
age | 25.82 51.19 1.115 | 33.23 122 .3478
black | .8432 .1329 -1.888 | .07354 .06813 3.268
educ | 10.35 4.043 -.7212 | 12.03 8.242 -.4233
After: _webal as the weighting variable
| Treat | Control
| mean variance skewness | mean variance skewness
------+---------------------------+---------------------------
age | 25.82 51.19 1.115 | 25.75 51.22 1.14
black | .8432 .1329 -1.888 | .8423 .1328 -1.879
educ | 10.35 4.043 -.7212 | 10.35 4.039 -.7224
可以看出，在调整后，处理组和加权调整后的对照组的年龄均值、方差和偏度都非常相似。需要指出的是，ebalance treat age, targets(2) 与 ebalance treat age age2, targets(1) 表示的执行指令相同，这里 age2 代表年龄变量 age 的平方项。
3.3 平衡性检验
为了验证某一个协变量匹配后是否平衡，可以使用 tabstat 进行验证，以年龄 age 为例：
. tabstat age [aweight=_webal], by(treat) s(N me v) nototal //检验是否平衡
Summary for variables: age
by categories of: treat (1 if treated, 0 control)
treat | N mean variance
---------+------------------------------
0 | 15992 25.75072 51.22003
1 | 185 25.81622 51.1943
----------------------------------------
默认情况下，ebalance 命令执行后，Stata 会自动将匹配权重储存为变量 _webal。从上面的检验结果看，经过熵平衡调整之后，age 变量在对照组的分布几乎与处理组一致。
3.4 变量交互项的匹配
在实际研究中，我们还需要考虑协变量之间的交互项。例如，为了调整协变量组，我们应保证黑人和非黑人的年龄均值相似，这里便要引入 black 和 age 的交互项。
. gen ageXblack = age*black //设置 age 和 black 的交互项 ageXblack
. ebalance treat age educ black ageXblack, targets(3 2 1 1) //进行熵平衡
Data Setup
Treatment variable: treat
Covariate adjustment: age educ black ageXblack (1st order). age educ (2nd order). age (3rd order).
Optimizing...
Iteration 1: Max Difference = 573647.601
Iteration 2: Max Difference = 211032.079
Iteration 3: Max Difference = 77633.2843
...
Iteration 15: Max Difference = .277671087
Iteration 16: Max Difference = .015938697
Iteration 17: Max Difference = .000055717
maximum difference smaller than the tolerance level; convergence achieved
Treated units: 185 total of weights: 185
Control units: 15992 total of weights: 185
Before: without weighting
| Treat | Control
| mean variance skewness | mean variance skewness
----------+---------------------------+----------------------------
age | 25.82 51.19 1.115 | 33.23 122 .3478
educ | 10.35 4.043 -.7212 | 12.03 8.242 -.4233
black | .8432 .1329 -1.888 | .07354 .06813 3.268
ageXblack | 21.91 134.6 -.4435 | 2.402 81.55 3.893
After: _webal as the weighting variable
| Treat | Control
| mean variance skewness | mean variance skewness
----------+---------------------------+---------------------------
age | 25.82 51.19 1.115 | 25.82 51.19 1.115
educ | 10.35 4.043 -.7212 | 10.35 4.043 -.7231
black | .8432 .1329 -1.888 | .8432 .1322 -1.888
ageXblack | 21.91 134.6 -.4435 | 21.91 133.2 -.4572
. bysort black: tabstat age [aweight=_webal], by(treat) s(N me v) nototal //根据 _webal 验证年龄均值在黑人和非黑人之间是平衡的
--------------------------------------------------------
-> black = 0
Summary for variables: age
by categories of: treat (1 if treated, 0 control)
treat | N mean variance
---------+------------------------------
0 | 14816 24.93109 45.28087
1 | 29 24.93103 40.49507
----------------------------------------
--------------------------------------------------------
-> black = 1
Summary for variables: age
by categories of: treat (1 if treated, 0 control)
treat | N mean variance
---------+------------------------------
0 | 1176 25.98077 52.16224
1 | 156 25.98077 53.2835
----------------------------------------
直播课程：实证研究设计 (2.4小时)
3.5 NSW 实验效果的无偏估计
为最大程度的平衡样本，我们使用所有协变量的一阶、二阶、三阶以及所有一阶交互项进行 ebalance 匹配。
. use cps1re74.dta, clear
(DW Subset of LaLonde Data)
. *协变量组合生成二阶矩及一阶交互项
. foreach v in age educ black hispan married nodegree re74 re75 u74 u75 {
foreach m in age educ black hispan married nodegree re74 re75 u74 u75 {
gen `v'X`m'=`v'*`m'
}
}
. *age、educ、re74 和 re75 作为连续变量，还需设置他们的三阶矩
. foreach v in age educ re74 re75 {
gen `v'X`v'X`v' = `v'^3
}
. *进行熵平衡，并将结果保存到 baltable.dta
. ebalance treat-re75Xre75Xre75, keep(baltable) replace
. reg re78 treat [pweight=_webal]
Linear regression Number of obs = 16,177
F(1, 16175) = 5.58
Prob > F = 0.0182
R-squared = 0.0161
Root MSE = 6889.6
---------------------------------------------------------------
| Robust
re78 | Coef. Std. Err. t P>|t| [95% Conf. Interval]
-------+-------------------------------------------------------
treat | 1761.344 745.533 2.36 0.018 300.017 3222.671
_cons | 4587.800 472.243 9.71 0.000 3662.151 5513.449
---------------------------------------------------------------
可以发现，上述回归结果显示的处理效果为 1761 美元，这表明熵平衡预处理步骤能够让我们非常接近随机实验的实际评估结果，95% 置信区间为 [300, 3223]。
4. 结论
熵平衡可以与回归方法相结合。研究者首先通过调整协变量对协变量进行重新加权，然后将加权后的数据应用于回归模型。这一过程可应用于回归前对处理组和对照组数据的预处理工作以及后续的稳健性检验。
在实证回归中，ebalance 命令的关键点是如何设置阶数 (numlist)：
对于虚拟变量及包含虚拟变量的交互项，只需进行一阶矩加权调整即可实现平衡； 对于年龄等连续变量，通常选择其一阶、二阶、三阶矩以及一阶与其他协变量一阶矩的交互项进行加权调整； 在进行熵平衡之后，权重值自动储存为变量 _webal。
熵平衡不是万能的。在实际研究中，为保证检验结果的稳健性，研究者最好搭配 NNM、PSM、CEM 等其他匹配方法一起进行实证分析。Hainmueller 和 Xu (2013) 也提出，未来会考虑如何将熵平衡与 Stata 中其他匹配方法相结合，进一步提高匹配效率。
5. 参考文献
Hainmueller J. Entropy balancing for causal effects: A multivariate reweighting method to produce balanced samples in observational studies[J]. Political analysis, 2012: 25-46. -PDF- Hainmueller J, Xu Y. Ebalance: A Stata package for entropy balancing[J]. Journal of Statistical Software, 2013, 54(7). -PDF- Dehejia R H, Wahba S. Causal effects in nonexperimental studies: Reevaluating the evaluation of training programs[J]. Journal of the American statistical Association, 1999, 94(448): 1053-1062. -PDF- 张海峰, 梁若冰, 林细细. 子女数量对农村家庭经济决策的影响——兼谈对 “二孩政策” 的启示[J]. 中国经济问题, 2019 (3): 68-80. -Link-
6. 相关推文
Note：产生如下推文列表的 Stata 命令为：
lianxh 匹配, m
安装最新版lianxh命令：
ssc install lianxh, replace
温馨提示： 文中链接在微信中无法生效。请点击底部「阅读原文」。或直接长按/扫描如下二维码，直达原文：
专题：Stata命令 Stata-Matching：肾脏交换匹配问题 Stata：iematch-近邻贪婪匹配 专题：数据处理 Stata：数据合并与匹配-merge-reclink 专题：倍分法DID 面板PSM+DID如何做匹配？ 专题：PSM-Matching Stata：终极匹配 ultimatch Stata 手动：各类匹配方法大全 A——理论篇 Stata：psestimate-倾向得分匹配(PSM)中协变量的筛选 Stata：广义精确匹配-Coarsened-Exact-Matching-(CEM) Stata：psestimate-倾向得分匹配(PSM)中匹配变量的筛选 Stata+PSM：倾向得分匹配分析简介 Stata-从匹配到回归：精确匹配、模糊匹配和PSM Stata：PSM-倾向得分匹配分析的误区 Stata：模糊匹配之matchit
课程推荐：面板数据因果推断
主讲老师：徐轶青 (斯坦福大学)
🍓 课程主页：https://gitee.com/arlionn/Course
New！ Stata 搜索神器：
lianxh和songblGIF 动图介绍
搜： 推文、数据分享、期刊论文、重现代码 ……
👉 安装：
. ssc install lianxh
. ssc install songbl
👉 使用：
. lianxh DID 倍分法
. songbl all
🍏 关于我们
连享会 ( www.lianxh.cn，推文列表) 由中山大学连玉君老师团队创办，定期分享实证分析经验。 直通车： 👉【百度一下： 连享会】即可直达连享会主页。亦可进一步添加 「知乎」,「b 站」,「面板数据」,「公开课」 等关键词细化搜索。