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

发布时间：2021-09-22 阅读 30564

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下载 - 推文合集

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

1. 引言

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

2. 理论逻辑

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

Y_{i t} = α + β T r e a t \times P o s t + ϕ X + η_{i} + λ_{t} + ε_{i t}

其中 $T r e a t$ 为是否受政策影响的虚拟变量， $P o s t$ 是政策实施时间的虚拟变量， $β$ 是我们关注的系数。 $X$ 是控制变量矩阵， $η_{i}$ 和 $λ_{t}$ 是个体固定效应和时间固定效应。

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

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

\hat{β} = β + γ \times \frac{c o v (T r e a t \times P o s t, ε_{i t} | W)}{v a r (T r e a t \times P o s t | W)}

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

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

Y_{i t} = α + β T r e a t_{f a k e} \times P o s t + + ϕ X + η_{i} + λ_{t} + ε_{i t}

由于 $T r e a t_{f a k e}$ 是随机产生的，实际政策效应 $β = 0$ 。在此前提下，如果估计出的 $\hat{β} = 0$ ，则可以逆推出 $γ = 0$ 。如果 $\hat{β} \neq 0$ ，则说明 $γ \neq 0$ ，文章的估计结果是有偏的，未观测的特征确实会影响估计结果。

3. permute 命令

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

permute 的基本语法如下：

permute permvar exp_list [, options] : command

permvar : 需要进行随机抽样的变量，即 DID 中的 $T r e a t$ ，或交互项 $T r e a t \times P o s t$
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%，同样是一个小概率事件。

通过上述操作，我们即可以推断出，不可观测的因素 $γ = 0$ ，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

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

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

1. 引言

2. 理论逻辑

3. permute 命令

4. 实例

4.1 代码演示

4.2 绘图

5. 拓展

6. 参考文献

7. 相关推文

相关课程

免费公开课

最新课程-直播课

⛳ 课程主页

⛳ 课程主页

关于我们

热门资讯