安慰剂检验绘图：最强攻略

[toc]

一、 缘起

在目前流行的双重差分方法中，安慰剂检验已经成为和平行趋势一样必不可少的检验流程。

这本来是一种在思路上很直观，在具体操作上也并不复杂的方法。但是，目前网络上流行的教程，多数只讲怎么写代码。而且代码似乎也并不简洁，动辄十几行、乃至几十行代码。另外，作者也没有告诉我们其间的原理。这也是 po 主一直不赞成的，因为我们学的是计量经济学，不是 Stata 经济学。

因为，最终还是决定自己动手写一篇教程，记录一下自己的学习经历。

二、 理论逻辑

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

$Y_{it} = \alpha+\beta Treat \times Post + \phi \text{X} + \eta_i + \lambda_t + \varepsilon_{it}$

其中 $Treat$ 为是否受政策影响的虚拟变量，$Post$ 是政策实施时间的虚拟变量，$\beta$ 是我们关注的系数。$\text{X}$ 是控制变量矩阵，$\eta_i$ 和 $\lambda_t$ 是个体固定效应和时间固定效应。

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

但是，在现实情况下，我们的政策会受到各种可观测因素与不可观测影响的影响，而我们无法穷尽所有控制变量，此时得到的估计结果可能是：
$\hat{\beta} = \beta+\gamma \times \frac{cov(Treat \times Post, \varepsilon_{it}|W)}{var(Treat \times Post|W)}$
其中，$\gamma$ 为非观测因素对被解释变量的影响。只有当 $\gamma = 0$ 时， 非观测因素才不会影响到估计结果，即 $\hat{\beta}$ 是无偏的。但是，这一点却无法直接验证，因为它本身就是不可观测的。我们只能通过间接手段来验证其是否为 0 。

当前，许多中文文章的安慰剂检验思路，多依循周茂、陆毅老师 2018 年发表在《中国工业经济》上的《开发区设立与地区制造业升级》一文。其逻辑是找到一个理论上不会对结果变量产生影响的错误变量 $Treat_{fake}$ 替代真实的 $Treat$ ：
$Y_{it} = \alpha+\beta Treat_{fake} \times Post + +\phi \text{X} + \eta_i + \lambda_t + \varepsilon_{it}$
由于 $Treat_{fake}$ 是随机产生的，实际政策效应 $\beta = 0$。在此前提下，如果估计出的 $\hat{\beta} = 0$，则可以逆推出 $\gamma = 0$。如果 $\hat{\beta} \neq 0$，则说明 $\gamma \neq 0$ ，文章的估计结果是有偏的，未观测的特征确实会影响估计结果。

三、 permute 命令

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

permute 的基本语法如下：

permute permvar exp_list [, options] : command

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

四、数据

我们先调入数据。这份数据来自普林斯顿大学的 DID 教学课件，可以通过网站："http://dss.princeton.edu/training/Panel101.dta 下载。设定好处理组和处理时间之后，我们先进行简单的回归。

*- 安慰剂检验
cd "~/Desktop"
use "Panel101.dta", clear
gen time  = (year >= 1994) & !missing(year)
gen treat = (country > 4)  & !missing(country)
replace y = y / 1e+09
gen did = time * treat
reghdfe y did, absorb(country year) vce(robust)

回归结果如下。使用双固定效应之后，交互项 did 在 10% 的统计水平显著。

. reghdfe y did, absorb(country year) vce(robust)
(MWFE estimator converged in 2 iterations)

------------------------------------------------------------------------------
             |               Robust
           y | Coefficient  std. err.      t    P>|t|     [95% conf. interval]
-------------+----------------------------------------------------------------
         did |     -2.520      1.286   -1.960   0.055       -5.098       0.059
       _cons |      2.493      0.406    6.147   0.000        1.679       3.306
------------------------------------------------------------------------------

五、安慰剂检验：方法 1

接下来，我们进行安慰剂检验。对交互项随机抽取 500 次，查看系数是否与基准估计结果存在显著差异。

下面汇报了抽样估计结果。我们只需要关注 beta 行即可。500 次抽样中，仅有一次结果在基准回归系数的右侧，占总抽样次数的 0.2%，所以我们可以看到它的 p value = 0.0020。这是单侧检验的结果。对于双侧检验而言，p value 则等于 0.0040。

可以看到，无论是单侧检验还是双侧检验，都表明在随机抽样的情况下，基准回归系数 -2.50 是一个小概率事件。这说明我们的安慰剂检验是成立的。

cap erase "simulations.dta"
permute did beta = _b[did] se = _se[did] df = e(df_r), ///
	reps(500) rseed(123) saving("simulations.dta"): ///
		reghdfe y did, absorb(country year) vce(robust)

Monte Carlo permutation results                    Number of observations =  70
Permutation variable: did                          Number of permutations = 500

      Command: reghdfe y did, absorb(country year) vce(robust)
         beta: _b[did]
           se: _se[did]
           df: e(df_r)

-------------------------------------------------------------------------------
             |                                               Monte Carlo error
             |                                              -------------------
           T |    T(obs)       Test       c       n      p  SE(p)   [95% CI(p)]
-------------+-----------------------------------------------------------------
        beta | -2.519512      lower       1     500  .0020  .0020  .0001  .0111
             |                upper     499     500  .9980  .0020  .9889  .9999
             |            two-sided                  .0040  .0028  .0000  .0095
             |
          se |  1.285631      lower     499     500  .9980  .0020  .9889  .9999
             |                upper       1     500  .0020  .0020  .0001  .0111
             |            two-sided                  .0040  .0028  .0000  .0095
             |
          df |        53      lower     500     500 1.0000  .0000  .9926 1.0000
             |                upper     500     500 1.0000  .0000  .9926 1.0000
             |            two-sided                 1.0000  .0000      .      .
-------------------------------------------------------------------------------
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.

5.1 回归系数图

接下来，我们按照固定范式，绘制回归系数的分布图。具体代码如下：

		// 回归系数
use "simulations.dta", clear
gen t_value = beta / se
gen p_value = 2 * ttail(df, abs(beta/se))

#delimit ;
dpplot beta, 
	xline(-2.520, lc(black*0.5) lp(dash))
	xline(0, lc(black*0.5) lp(solid))
	xlabel(-3(1)3)
    xtitle("Estimator", size(*0.8)) xlabel(, format(%4.1f) labsize(small))
    ytitle("Density", size(*0.8)) ylabel(, nogrid format(%4.1f) labsize(small)) 
    note("") caption("") graphregion(fcolor(white)) ;
#delimit cr

绘图结果如下。通过系数分布图，我们也可以清洗的观察到，随机抽样系数以零为均值，呈正态分布。仅有一次抽样系数位于 -2.520 的右侧。

除了回归系数分布之外，我们也可以画 t 值的分布图：

5.2 t 值图

	// t 值
#delimit ;
dpplot t_value, 
	xline(-1.960, lc(black*0.5) lp(dash))
	xline(0, lc(black*0.5) lp(solid))
    xtitle("T Value", size(*0.8)) xlabel(, format(%4.1f) labsize(small))
    ytitle("Density", size(*0.8)) ylabel(, nogrid format(%4.1f) labsize(small)) 
    note("") caption("") graphregion(fcolor(white)) ;
#delimit cr

绘图结果如下。观察下图，我们可以得出相似的结论，大部分随机抽样结果的 t 值都位于零值附近，仅有少数估计结果的 t 值大于基准回归结果。

5.3 p 值图

考虑到部分文章也会绘制 p 值的分布结果，我们也同样可以实现：

	// p 值
#delimit ;
dpplot p_value, 
	xline(0.055, lc(black*0.5) lp(dash))
    xtitle("P Value", size(*0.8)) xlabel(, format(%4.1f) labsize(small))
    ytitle("Density", size(*0.8)) ylabel(, nogrid format(%4.1f) labsize(small)) 
    note("") caption("") graphregion(fcolor(white)) ;
#delimit cr

下图是 p 值的分布图。相对而言，p 值图没有前两者美观。

也有部分文章会同时绘制回归系数和 p 值的分布图，我们也尝试画一个：

5.4 系数与 p 值结合图

	// 系数和 p 值结合
#delimit ;
twoway (scatter p_value beta)(kdensity beta, yaxis(2)), 
	xline(-2.520, lc(black*0.5) lp(dash))
	xline(0, lc(black*0.5) lp(solid)) 
	yline(0.055, lc(black*0.5) lp(dash))
	xlabel(-3(1)3)
	xtitle("Estimator", size(*0.8)) xlabel(, format(%4.1f) labsize(small))
	ytitle("Density", size(*0.8)) ylabel(, nogrid format(%4.1f) labsize(small)) 
	ytitle("P Value", size(*0.8) axis(2)) ylabel(, nogrid format(%4.1f) labsize(small) axis(2))
	legend(r(1) order(1 "P Value" 2 "Estimator"))
	graphregion(color(white)) ;
#delimit cr

主坐标用来标识回归系数，副坐标用来标识 p 值。绘图结果如下。虽然这幅图将系数和 p 值统合在一起，但是却也牺牲了一定的美观度。因此，po 主还是倾向于仅汇报回归系数的分布图。

六、安慰剂检验：方法 2

与第一种方法直接抽取交互项不同的是，第二种方法是抽取处理组。只需要将permute 后面的 did 替换为 treat 即可：

*- 第二种抽样方法
cd "~/Desktop"
use "Panel101.dta", clear
gen time  = (year >= 1994) & !missing(year)
gen treat = (country > 4)  & !missing(country)
replace y = y / 1e+09
gen did = time * treat

cap erase "simulations.dta"
permute treat beta = _b[c.treat#c.time] se = _se[c.treat#c.time] df = e(df_r), ///
	reps(500) rseed(123) saving("simulations.dta"): ///
		reghdfe y c.treat#c.time, absorb(country year) vce(robust)

抽样结果如下。单侧检验 p 值为0.004，双侧检验结果为 0.008，都是小概率事件。安慰剂检验通过。

Monte Carlo permutation results                    Number of observations =  70
Permutation variable: treat                        Number of permutations = 500

      Command: reghdfe y c.treat#c.time, absorb(country year) vce(robust)
         beta: _b[c.treat#c.time]
           se: _se[c.treat#c.time]
           df: e(df_r)

-------------------------------------------------------------------------------
             |                                               Monte Carlo error
             |                                              -------------------
           T |    T(obs)       Test       c       n      p  SE(p)   [95% CI(p)]
-------------+-----------------------------------------------------------------
        beta | -2.519512      lower       2     500  .0040  .0028  .0005  .0144
             |                upper     498     500  .9960  .0028  .9856  .9995
             |            two-sided                  .0080  .0040  .0002  .0158
             |
          se |  1.285631      lower     467     500  .9340  .0111  .9086  .9541
             |                upper      33     500  .0660  .0111  .0459  .0914
             |            two-sided                  .1320  .0151  .1023  .1617
             |
          df |        53      lower     500     500 1.0000  .0000  .9926 1.0000
             |                upper     500     500 1.0000  .0000  .9926 1.0000
             |            two-sided                 1.0000  .0000      .      .
-------------------------------------------------------------------------------

6.1 回归系数图

我们同样可以绘制一个回归系数分布图：

		// 回归系数
use "simulations.dta", clear
gen t_value = beta / se
gen p_value = 2 * ttail(df, abs(beta/se))

#delimit ;
dpplot beta, 
	xline(-2.520, lc(black*0.5) lp(dash))
	xline(0, lc(black*0.5) lp(solid))
	xlabel(-3(1)3)
    xtitle("Estimator", size(*0.8)) xlabel(, format(%4.1f) labsize(small))
    ytitle("Density", size(*0.8)) ylabel(, nogrid format(%4.1f) labsize(small)) 
    note("") caption("") graphregion(fcolor(white)) ;
#delimit cr

所得结果与前面的基本一致：

6.2 t 值图

下面是 t 值分布图：

	// t 值
#delimit ;
dpplot t_value, 
	xline(-1.960, lc(black*0.5) lp(dash))
	xline(0, lc(black*0.5) lp(solid))
    xtitle("T Value", size(*0.8)) xlabel(, format(%4.1f) labsize(small))
    ytitle("Density", size(*0.8)) ylabel(, nogrid format(%4.1f) labsize(small)) 
    note("") caption("") graphregion(fcolor(white)) ;
#delimit cr

6.3 p 值图

下面是 p 值分布图：

	// p 值
#delimit ;
dpplot p_value, 
	xline(0.055, lc(black*0.5) lp(dash))
    xtitle("P Value", size(*0.8)) xlabel(, format(%4.1f) labsize(small))
    ytitle("Density", size(*0.8)) ylabel(, nogrid format(%4.1f) labsize(small)) 
    note("") caption("") graphregion(fcolor(white)) ;
#delimit cr

6.4 系数与 p 值结合图

绘图代码与前文一致：

	// 系数和 p 值结合
#delimit ;
twoway (scatter p_value beta)(kdensity beta, yaxis(2)), 
	xline(-2.520, lc(black*0.5) lp(dash))
	xline(0, lc(black*0.5) lp(solid)) 
	yline(0.055, lc(black*0.5) lp(dash))
	xlabel(-3(1)3)
	xtitle("Estimator", size(*0.8)) xlabel(, format(%4.1f) labsize(small))
	ytitle("Density", size(*0.8)) ylabel(, nogrid format(%4.1f) labsize(small)) 
	ytitle("P Value", size(*0.8) axis(2)) ylabel(, nogrid format(%4.1f) labsize(small) axis(2))
	legend(r(1) order(1 "P Value" 2 "Estimator"))
	graphregion(color(white)) ;
#delimit cr

七、总结

使用 permute 命令，我们基本可以一行代码搞定安慰剂检验。这个命令有两个明显的优点。其一，我们可以实时观看抽样进度，不用在电脑面前瞎等。同时，该命令也会提供单侧检验与双侧检验的 p 值，帮助我们直接判断安慰剂检验是否通过，相比于肉眼的直观判断更为客观。

以上介绍了两种抽样方法、八种安慰剂检验图的画法，大家可以视情况而定进行汇报。这也是计量经济学的魅力所在，只要我们熟悉特定的方法，就可以得到想要的结果 （bushi）。。。