Stata：双层优化匹配估计量-blopmatch

发布时间：2022-12-19 阅读 954

Stata连享会主页 || 视频 || 推文 || 知乎 || Bilibili 站

温馨提示： 定期清理浏览器缓存，可以获得最佳浏览体验。

New！ lianxh 命令发布了：
随时搜索推文、Stata 资源。安装：
. ssc install lianxh
详情参见帮助文件 (有惊喜)：
. help lianxh
连享会新命令：cnssc, ihelp, rdbalance, gitee, installpkg

✌ 课程详情： https://gitee.com/lianxh/Course

⛳ 课程主页： https://gitee.com/lianxh/Course

⛳ Stata 系列推文：

全部 | Stata入门 | Stata教程 | Stata资源 | Stata命令
计量专题 | 论文写作 | 数据分享 | 专题课程
结果输出 | Stata绘图 | 数据处理 | Stata程序
回归分析 | 面板数据 | 交乘项-调节 | IV-GMM
内生性-因果推断 | 倍分法DID | 断点回归RDD | PSM-Matching | 合成控制法
Probit-Logit | 时间序列 | 空间计量 | 分位数回归 | 生存分析 | SFA-DEA
文本分析-爬虫 | Python-R-Matlab | 机器学习
Markdown | 工具软件 | 其它

☝ PDF下载 - 推文合集

作者：陈广宇 (中山大学)
邮箱：cgyyy_sysu@qq.com

1. 简介

在经济学中，我们常常希望评估某项公共政策实施后的效应。为了解决样本自选择问题，并得到政策的平均处理效应 (average treatment effect)，匹配估计量在评价文献中被广泛使用。

常用的匹配方法有最近邻特征匹配 (k-Nearest neighbors matching)、马氏距离匹配 (Mahalanobis matching)、半径匹配 (Radius matching)、倾向得分匹配 (Propensity score match, PSM) 等。其中，倾向得分匹配又需要基于前面提到的任意一种匹配的方法。

最流行匹配方法的莫过于基于最近邻匹配的倾向得分匹配。然而，在进行最邻近匹配的过程中，匹配的结果和效果会受到研究人员给定的近邻数量影响。双层优化匹配 (BLOP Matching) 通过将匹配和权重的选择与匹配后协变量平衡的改进联系起来，能够内生决定最佳近邻数量，且在协变量平衡方面表现出色。

接下来，本文主要介绍双层优化匹配 (BLOP Matching) 的理论内容和 Stata 操作。

2. 理论背景

2.1 基本定义

通常一个匹配问题包含三个关键变量：

结果变量 $Y_{i}$ ：即我们关心的实验结果，或者理解为回归中的被解释变量；
分组变量 $W_{i}$ ：即样本是否被处理，通常用哑变量表示；
协变量 $X_{i}$ ：即匹配时需要被控制，使得处理组和控制组尽可能接近的变量。当需要控制的协变量有 $K$ 个时，就构成一个 $K$ 维向量 $X_{i}$ 。

以下是关于双层优化匹配问题的一些重要定义：

2.1.1 反事实估计

{\hat{Y}}_{i} (0) = {\begin{cases} Y_{i} & if & W_{i} = 0 \\ {\hat{Y}}_{i}^{b} & if & W_{i} = 1 \end{cases}

{\hat{Y}}_{i} (1) = {\begin{cases} {\hat{Y}}_{i}^{b} & if & W_{i} = 0 \\ Y_{i} & if & W_{i} = 1 \end{cases}

${\hat{Y}}_{i}^{b}$ ：表示 $Y_{i}$ 的反事实估计，即在控制协变量、未被处理 (被处理) 的情况下，处理组 (对照组) 的估计结果

2.1.2 平均处理效应

{\hat{τ}}^{b} = \frac{1}{N} \sum_{i = 1}^{N} ({\hat{Y}}_{i} (1) - {\hat{Y}}_{i} (0))

{\hat{τ}}_{tre}^{b} = \frac{1}{N_{1}} \sum_{i = 1}^{N} W_{i} ({\hat{Y}}_{i} (1) - {\hat{Y}}_{i} (0))

${\hat{τ}}^{b}$ ：全样本的平均处理效应 (average treatment effect，ATE)。 ${\hat{τ}}_{tre}^{b}$ ：处理组的平均处理效应 (average treatment effect on the treated，ATET)。其中 $N_{1}$ 是处理组样本的个数， $N_{0}$ 是控制组样本的个数， $N = N_{1} + N_{0}$ 是总样本的个数。

2.1.3 凸组合和凸壳

凸组合 (Convex combination)： $λ_{1} X_{1} + . . . + λ_{N} X_{N}$ ，其中 $λ_{i}$ 非负且 $λ_{1} + . . . + λ_{N} = 1$ 。凸组合是多个协变量的一种特殊的线性组合。

凸壳 (Convex hull)： $co {X_{1}, \dots, X_{N}} = {\sum_{j = 1}^{N} λ_{j} X_{j}, (λ_{1}, \dots, λ_{N}) \in Δ_{N}}$ ，其中 $Δ_{N} = {(λ_{1}, \dots, λ_{N}) \in R_{+}^{N}, \sum_{j = 1}^{N} λ_{j} = 1}$ 。凸壳是所有凸组合的集合，形成了 K 维向量空间的一个壳。

2.1.4 距离和投影

距离：用 $| | \cdot | |$ 表示，可以是欧氏距离、马氏距离等。

投影：用 $P r o j (X)$ 表示，是凸壳 (Convex hull) 上到向量 $X$ 距离最短的向量。

2.2 双层优化匹配

双层优化匹配将匹配问题分解为两个问题：第一个问题是投影匹配问题，即在凸壳上求解待匹配协变量 $X$ 的投影 $P r o j (X)$ ，也就是在所有控制组向量组成的凸壳里寻找凸组合，使得这个凸组合距离向量 $X$ 最短。用数学公式表示为：

\sum_{j = 1}^{N_{0}} λ_{j} ∥ X_{i} - X_{j} ∥

第二个问题是优化问题。由于第一个问题给出的投影 $P r o j (X)$ 多数情况下不唯一，需要用一个描述协变量平衡好坏程度的指标 $\sum_{j = 1}^{N_{0}} λ_{j} ∥ X_{i} - X_{j} ∥$ 来确定最优的匹配估计量。用数学公式表示为：

\begin{aligned} min_{(λ_{1}, \dots, λ_{N_{0}})} \sum_{j = 1}^{N_{0}} λ_{j} ∥ X_{i} - X_{j} ∥ \\ s.t. \\ \sum_{j = 1}^{N_{0}} λ_{j} X_{j} = Proj (X_{i}), \\ \sum_{j = 1}^{N_{0}} λ_{j} = 1, λ_{j} \geq 0, j = 1, \dots, N_{0} \end{aligned}

解出上述两个线性优化问题能够对每个处理组样本求解出一组权重 $λ_{i} = (λ_{1}, . . ., λ_{N})$ ，使得凸组合 $λ_{1} X_{1} + . . . + λ_{N} X_{N}$ 到待匹配协变量 $X_{i}$ 距离最短，且使得协变量能够最好地进行平衡。这相当于对于处理组里面的每一个样本，都找到了最优的邻居数，并进行了匹配。

类似地，对于控制组的个体，将凸壳从 $Δ_{N_{0}}$ 换成 $Δ_{N_{1}}$ 即可。最终，反事实估计量可以由如下式子表示：

{\hat{Y}}_{i}^{b} = W_{i} \sum_{j = 1}^{N_{0}} λ_{j} (i) Y_{j} + (1 - W_{i}) \sum_{j = 1}^{N_{1}} λ_{j} (i) Y_{N_{0} + j}

从而通过下式计算出 ATE 和 ATT：

\begin{aligned} {\hat{τ}}^{b} = \frac{1}{N} \sum_{i = 1}^{N} ({\hat{Y}}_{i} (1) - {\hat{Y}}_{i} (0)) \\ {\hat{τ}}_{tre}^{b} = \frac{1}{N_{1}} \sum_{i = 1}^{N} W_{i} ({\hat{Y}}_{i} (1) - {\hat{Y}}_{i} (0)) \end{aligned}

值得一提的是，BLOP 估计量具有一致性和渐近正态性。

3. 命令介绍

blopmatch 命令是 Díaz et al. (2021) 根据以上估计方法所编写的 Stata 新命令。

命令安装：

net install st0632.pkg, replace

命令语法：

blopmatch (ovar omvarlist) (tvar) [if] [in] [weight] [, stat options]

其中，必填项包括：

ovar：结果变量 (outcome) ，可以是二进制、计数、连续、分数或非负结果；
omvarlist：指定结果模型中的协变量 (covariate)；
tvar：表示处理级别，必须是整数，只允许二值变量 (通常是哑变量)。

可选项 [stat] 包括：

ate：默认，估计全样本的平均处理效应 (average treatment effect)；
atet：估计处理组的平均处理效应 (average treatment effect on the treated)。

可选项 [options] 包括：

ematch(varlist)：对于某些变量采用精确匹配；
level：设置置信度；
dmvariables：输出匹配变量的名字；
control：tvar 控制组的值，默认是 tvar 较小的，比如哑变量中的 0；
tlevel：tvar 实验组的值，默认是 tvar 较大的，比如哑变量中的 1；
metric：指定多维度下距离的度量方式，默认为马氏距离。

4. Stata 实例

4.1 实例演示

接下来，我们使用 LaLonde (1986) 文章中美国国家劳动培训计划 (National Supported Work Demonstration) 的数据来演示双层优化匹配 (BLOP Matching) 的 Stata 操作。

. * 数据描述
. lxhuse nsw.dta, clear
. describe

 Observations:   722   Sample from the National Supported Work Demonstration
    Variables:    10   15 Nov 2019 12:24
                        (_dta has notes)
-----------------------------------------------------------------------------
Variable      Storage   Display    Value
    name         type    format    label      Variable label
-----------------------------------------------------------------------------
data_id         str14   %14s                  
treat           byte    %9.0g      treat      treatment indicator
age             byte    %8.0g                 age (in years)
education       byte    %8.0g                 schooling (in years)
black           byte    %9.0g      black      1(black)
hispanic        byte    %9.0g      hispanic   1(hispanic)
married         byte    %9.0g      married    1(married)
nodegree        byte    %21.0g     nodegree   1(no high school degree)
re75            float   %9.0g                 earnings in 1975 (in 1982 USD)
re78            float   %9.0g                 earnings in 1978 (in 1982 USD)
-----------------------------------------------------------------------------

其中，re78 是工人 1978 年的收入，也是我们感兴趣的结果。协变量包括年龄、教育、人种哑变量 (以白种人为基准组)、是否接受过中学教育、工人 1975 年的收入 re75。对照组与实验组用 treat 区分，1 为实验组，0 为对照组。

. bysort treat: fsum

-> treat = otherwise
  Variable |        N     Mean       SD      Min      Max                                                   
-----------+---------------------------------------------
     treat |      425     0.00     0.00     0.00     0.00  
       age |      425    24.45     6.59    17.00    55.00  
 education |      425    10.19     1.62     3.00    14.00  
     black |      425     0.80     0.40     0.00     1.00  
  hispanic |      425     0.11     0.32     0.00     1.00  
   married |      425     0.16     0.36     0.00     1.00  
  nodegree |      425     0.81     0.39     0.00     1.00  
      re75 |      425  3026.68  5201.25     0.00 36941.27  
      re78 |      425  5090.05  5718.09     0.00 39483.53  
---------------------------------------------------------

-> treat = treated
  Variable |        N     Mean       SD      Min      Max                                                     
-----------+---------------------------------------------
     treat |      297     1.00     0.00     1.00     1.00  
       age |      297    24.63     6.69    17.00    49.00  
 education |      297    10.38     1.82     4.00    16.00  
     black |      297     0.80     0.40     0.00     1.00  
  hispanic |      297     0.09     0.29     0.00     1.00  
   married |      297     0.17     0.37     0.00     1.00  
  nodegree |      297     0.73     0.44     0.00     1.00  
      re75 |      297  3066.10  4874.89     0.00 37431.66  
      re78 |      297  5976.35  6923.80     0.00 60307.93

可以看出，样本中属于控制组的有 425 人，属于处理组的有 297 人。控制组和处理组在年龄、教育、人种、是否结婚、是否是否接受过中学教育和 1975 年的工资均有较为相似的分布。

. egen stdre75 = std(re75) // 对 re75 进行标准化

. * 匹配并计算平均处理效应 (ATE)
. global Y = "re78"
. global W = "treat"
. global X = "age education black hispanic married nodegree stdre75"
. blopmatch ($Y $X) ($W), ate
                                                                    
Treatment-effects estimation                                        
Estimator      : blop matching                 Number of obs        =      722
Outcome model  : matching                      Control group size   =      425
Distance Metric: Mahalanobis                   Treatment group size =      297
------------------------------------------------------------------------------
        re78 | Coefficient  Std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
         ate |    975.625    503.853     1.94   0.053      -11.910    1963.159
------------------------------------------------------------------------------

根据 BLOP matching，所有样本都得到了匹配。ATE 显示，美国国家劳动培训计划 (National Supported Work Demonstration) 平均给工人 1978 年的工资带来了 975.625 美元的提升，该效应在 10% 的统计水平上显著。

4.2 结果对比

接下来我们对比 BLOP matching 和最近邻匹配的结果：


. * 首先是 k=1，即进行 1:1 匹配
. teffects nnmatch ($Y $X) ($W), ate nneighbor(1)  // ATE = 744.979

Treatment-effects estimation                   Number of obs      =        722
Estimator      : nearest-neighbor matching     Matches: requested =          1
Outcome model  : matching                                     min =          1
Distance metric: Mahalanobis                                  max =          8
-----------------------------------------------------------------------------------------
                        |              AI robust
                   re78 | Coefficient  std. err.      z    P>|z|     [95% conf. interval]
------------------------+----------------------------------------------------------------
ATE                     |
                  treat |
(treated vs otherwise)  |    744.979    565.032     1.32   0.187     -362.464    1852.422
-----------------------------------------------------------------------------------------

. * 更改 k=16
. teffects nnmatch ($Y $X) ($W), ate nneighbor(16) // ATE = 968.917

Treatment-effects estimation                   Number of obs      =        722
Estimator      : nearest-neighbor matching     Matches: requested =         16
Outcome model  : matching                                     min =         16
Distance metric: Mahalanobis                                  max =         23
-----------------------------------------------------------------------------------------
                        |              AI robust
                   re78 | Coefficient  std. err.      z    P>|z|     [95% conf. interval]
------------------------+----------------------------------------------------------------
ATE                     |
                  treat |
(treated vs otherwise)  |    968.917    517.388     1.87   0.061      -45.145    1982.978
-----------------------------------------------------------------------------------------

可以看出，随着 k 的改变，ATE 的计算结果也会改变。而怎样才能在确定最佳的 k 近邻数量，做好协变量平衡，最近邻匹配没有给出答案。

5. 参考资料

Díaz, J. D., Gutiérrez, I., & Rivera, J. (2021). Implementing blopmatching in Stata. The Stata Journal, 21(1), 180–194. -Link- -PDF-
Díaz, J. D., Rau, T. & Rivera, J. (2015). A Matching Estimator Based on a Bilevel Optimization Problem. The Review of Economics and Statistics, 97(4), 803-812. -Link- -PDF-

6. 相关推文

Note：产生如下推文列表的 Stata 命令为：
lianxh 匹配, m
安装最新版 lianxh 命令：
ssc install lianxh, replace

专题	嘉宾	直播/回看视频
⭐ 最新专题		文本分析、机器学习、效率专题、生存分析等
研究设计	连玉君	我的特斯拉-实证研究设计，-幻灯片-
面板模型	连玉君	动态面板模型，-幻灯片-
面板模型	连玉君	直击面板数据模型 [免费公开课，2小时]