Stata：一行代码实现安慰剂检验-permute

理论模型可以简洁、凝练地抽离出经济现象的本质，使我们能够进行更深层次的思考和分析。然而，建立理论模型并非易事，若能将 理论和实证有机结合，那更加难能可贵了。

为此，我们邀请到了中山大学岭南学院郭凯明副教授，与大家一同学习理论模型的构建。郭老师一直专注于经济转型与中国经济方面的研究，发表论文近 40 篇，其中《经济研究》7 篇。

郭老师将从模型设定初衷、最基本的假设条件入手，通过讨论各种可能的建模思路和弯路，让学生不自觉中已经建立起理论分析的思维模式。最终的目标是：让学生不仅能「读懂模型」，还能「折腾模型」—— 可以自己修改甚至新设模型。

扫码直达课程主页：

作者：陈波 ( 深圳大学 )
邮箱：1900123011@email.szu.edu.cn

目录

  • 1. 引言

  • 2. 理论逻辑

  • 3. permute 命令

  • 4. 实例

    • 4.1 代码演示

    • 4.2 绘图

  • 5. 拓展

  • 6. 参考文献

  • 7. 相关推文

温馨提示： 文中链接在微信中无法生效。请点击底部「阅读原文」。或直接长按/扫描如下二维码，直达原文：

1. 引言

在现今广为流行的双重差分模型中，除了必不可少的平行趋势外，安慰剂检验也是重要的一环。其中应用最为普遍的安慰剂检验，无疑是检验不可观测的特征的影响。其基本思路是随机设定处理组，查看估计系数分布是否依然显著。道理虽然简单，但是在操作上却涌现出各异的方法，且似乎都不是那么简单。因此，本文试图通过一行代码来解决这个问题。

2. 理论逻辑

我们先引入一个经典的 DID 模型：

其中 为是否受政策影响的虚拟变量， 是政策实施时间的虚拟变量， 是我们关注的系数。 是控制变量矩阵， 和 是个体固定效应和时间固定效应。

在理想情况下，如果我们的政策是外生的，不受不可观测的因素的影响，此时可以直接通过 OLS 估计得到系数 的一致估计量 。

但是，在现实情况下，我们的政策会受到各种可观测因素与不可观测影响的影响，而我们无法穷尽所有控制变量，此时得到的估计结果可能是：

其中， 为非观测因素对被解释变量的影响。只有当 时， 非观测因素才不会影响到估计结果，即 是无偏的。但是，这一点却无法直接验证，因为它本身就是不可观测的。我们只能通过间接手段来验证其是否为 0 。

当前，许多中文文章的安慰剂检验思路，多依循周茂、陆毅老师 2018 年发表在《中国工业经济》上的《开发区设立与地区制造业升级》一文。其逻辑是找到一个理论上不会对结果变量产生影响的错误变量 替代真实的 ：

由于 是随机产生的，实际政策效应 。在此前提下，如果估计出的 ，则可以逆推出 。如果 ，则说明 ，文章的估计结果是有偏的，未观测的特征确实会影响估计结果。

3. permute 命令

现存的教程，多使用 forvalue 循环，随机抽取样本进行一定次数的回归。这种方法在操作的过程中，至少会生成好几个文件。命令也少则几行，多则几十行。而现在我们引入 permute 命令，一行代码即可实现安慰剂检验。

permute 的基本语法如下：

permute permvar exp_list [, options] : command
  • permvar : 需要进行随机抽样的变量，即 DID 中的 ，或交互项
  • exp_list : 需要提取的统计量，一般是回归系数
  • options 有以下设定：
    • reps(#) : 抽样次数
    • enumerate : 计算所有可能的不同排列
    • rseed(#) : 设定抽样种子
    • strara(varlist) : 分层抽样
    • saving(file) : 保存抽样值
  • command : 回归命令

4. 实例

4.1 代码演示

我们使用 permute 命令自带的数据进行演示：

. webuse permute2, clear

. list group y

     +------------+     | group    y |     |------------|  1. |     1    6 |  2. |     1   11 |  3. |     1   20 |  4. |     1    2 |  5. |     1    9 |     |------------|  6. |     1    5 |  7. |     0    2 |  8. |     0    1 |  9. |     0    6 | 10. |     0    0 |     |------------| 11. |     0    2 | 12. |     0    3 | 13. |     0    3 | 14. |     0   12 | 15. |     0    4 |     |------------| 16. |     0    1 | 17. |     0    5 |     +------------+

该数据共有 17 个样本，按照 group 分为两组。我们假设 group = 1 是处理组，group = 0 是控制组。

用 group 对 y 进行回归， 可以发现两组之间存在显著差异。

. reg y group

------------------------------------------------------------------------
y | Coefficient Std. err. t P>|t| [95% conf. interval]
-------+----------------------------------------------------------------
group | 5.288 2.306 2.294 0.037 0.374 10.202
_cons | 3.545 1.370 2.588 0.021 0.626 6.465
------------------------------------------------------------------------

接下来，我们需要进行安慰剂检验，以检验该结论是否收到不可观测因素的影响。

我们对 group 进行随机抽样，重复 500 次：

permute group _b[group], reps(500) rseed(123): reg y group

其基本结果如下：

Monte Carlo permutation results Number of observations = 17
Permutation variable: group Number of permutations = 500

Command: regress y group
_pm_1: _b[group]

--------------------------------------------------------------------------
| Monte Carlo error
| -------------------
T | T(obs) Test c n p SE(p) [95% CI(p)]
--------+-----------------------------------------------------------------
_pm_1 | 5.287879 lower 491 500 .9820 .0059 .9661 .9917
| upper 10 500 .0200 .0063 .0096 .0365
| two-sided .0400 .0088 .0228 .0572
--------------------------------------------------------------------------
Notes: For lower one-sided test, c = #{T <= T(obs)} and p = p_lower = c/n.
For upper one-sided test, c = #{T >= T(obs)} and p = p_upper = c/n.
For two-sided test, p = 2*min(p_lower, p_upper); SE and CI approximate.

_pm_1 是我们的实际估计系数，即之前 OLS 的估计结果。Test 是指单侧检验还是双侧检验。可以看到，500 次抽样中，有 491 次的抽样估计结果小于 5.288，仅有 10 次大于 5.288 (有时候实际抽样次数相加不一定等于设定的抽样次数，但会十分接近)。

该结果说明在随机抽样的情况下，估计值大于 5.288 的概率仅为 2%，无疑是一个小概率事件。双侧检验的结果也显示，估计值的绝对值大于 5.288 的概率为 4%，同样是一个小概率事件。

通过上述操作，我们即可以推断出，不可观测的因素 ，OLS 估计结果受到不可观测因素影响的可能性较小。

4.2 绘图

作为安慰剂检验的标志性动作，我们一般都会绘制一个估计系数的核密度分布图。我们同样也可以这么做：

. permute group beta = _b[group],  ///          reps(500) rseed(123) saving('simulations.dta'):  ///          reg y group

引入 saving 命令，将抽样估计系数保存到 simulations.dta 文件中。随后使用 dpplot 命令进行绘图。

use 'simulations.dta', clear
#delimit ;
dpplot beta, xline(5.288, lc(black*0.5) lp(dash))
             xline(0, lc(black*0.5) lp(solid))
             xtitle('Estimator', size(*0.8)) 
             xlabel(-8(4)8, format(%4.1f) labsize(small))
             ytitle('Density', size(*0.8)) 
             ylabel(, nogrid format(%4.1f) labsize(small)) 
             note('') caption('') 
             graphregion(fcolor(white)) ;
#delimit cr
graph export '安慰剂检验.png', width(1000) replace

稍作修饰之后，下图就是很典型的安慰剂检验结果了。可以看到，估计系数分布在零的附近，且服从正态分布，符合安慰剂检验的预期。

5. 拓展

由于没有趁手的数据，本文只是用 Stata 自带的数据，简单演示了一下 permute 在进行安慰剂检验时的做法。在实际的 DID 应用中，大家可以将本文的 group 替换为 treat，或者是 treat 和 post 的交互项，这两种做法都是可行的。reg 命令也可以替换为 xtreg、reghdfe 等常用的估计命令，设置一系列固定效应等，permute 都可以完美兼容。

6. 参考文献

  • 周茂,陆毅,杜艳,姚星.开发区设立与地区制造业升级[J].中国工业经济,2018(03):62-79. -PDF-

  • 宋弘,孙雅洁,陈登科.政府空气污染治理效应评估——来自中国“低碳城市”建设的经验研究[J].管理世界,2019,35(06):95-108+195. -PDF-

7. 相关推文

Note：产生如下推文列表的 Stata 命令为
lianxh 稳健性 安慰剂

安装最新版 lianxh 命令：

ssc install lianxh, replace

  • 专题：回归分析
    • 稳健性检验！稳健性检验！
    • aoeplacebo：地理安慰剂检验
    • acreg：允许干扰项随意相关的稳健性标准误
  • 专题：倍分法DID
    • 多期DID之安慰剂检验、平行趋势检验
  • 专题：内生性-因果推断
    • 安慰剂检验！安慰剂检验！
    • Stata新命令：konfound - 因果推断的稳健性检验

New！ Stata 搜索神器：lianxh 和 songbl  GIF 动图介绍
搜： 推文、数据分享、期刊论文、重现代码 ……
👉 安装：
  . ssc install lianxh
  . ssc install songbl
👉  使用：
  . lianxh DID 倍分法
  . songbl all

赞 (0)

相关推荐