Stata：pzms-RDD和RKD的最优模型选择

发布时间：2023-04-04 阅读 621

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

作者： 黄思佳 (湖南大学)
邮箱： ystone_17@163.com

编者按: 本文主要参考自下文，特此致谢！
Source： Kettlewell N. & Siminski P. (2022). Optimal Model Selection in RDD and Related Settings Using Placebo Zones. -PDF-

1. 简介

断点回归（RDD）和拐点回归（RKD）通常被用来评估政策干预的效果。以RRD为例，通过比较阈值(cut off)左右两侧的平均值来评估政策干预的效果。然而，该方法对模型的选择和参数估计非常敏感。在实证过程中，要面临的首要问题是如何在带宽、函数形式、核密度、协变量等多维度参数中进行选择。选择不合适的模型可能会导致错误的结论。因此，需要选择一个最优模型来准确地评估政策干预的效果。

为了解决上述问题，本文提出了一种新的方法——Placebo Zones。该方法有助于选择最优模型，以评估模型的准确性和稳健性。

具体而言，首先，选择阈值两侧的一定范围内的值来创建安慰剂区(Placebo Zones)，在范围足够大的情况下，安慰剂区不会受到政策冲击的影响。安慰剂区中的数据可以用来估计零效应 (null effect) ，如果估计模型是有效的，则该零效应的估计值为0（安慰剂区未受到政策冲击，其真实值为0）。

其次，通过将估计的处理效应与安慰剂区的零效应进行比较，可以评估RDD设定的有效性和所选模型的准确性。如果零效应估计为0，且估计的处理效应具有统计显著性，则说明具有因果效应。反之，如果零效应异于0，则说明模型设定无效，估计的治疗效应可能存在偏差。

综上，Placebo Zones是一种有助于选择最佳回归模型的新方法，能够帮助我们准确评估政策干预的效果。这种方法非常适用于RDD、RKD以及cohort-IV，因为它有助于排除其他因素的影响，从而获得准确和可靠的结果。

2. 推导

参考标准的精确 RDD 模型，考虑一个随机样本 $(Y_{i} (0), Y_{i} (1), X_{i})$ ， $i = 1, 2, \dots, n$ ，其中 $Y_{i} (0)$ 和 $Y_{i} (1)$ 分别为接受和未接受处理的潜在产出。处理 $(T)$ 由驱动变量 $X$ 决定， $X_{i} \geq 0$ 时， $T_{i} = 1$ 。感兴趣的参数是平均处理效应 $τ = E [Y_{i} (1) - Y_{i} (0) | X_{i} = 0]$ 。

本文的最优模型选择是渐进最优的，即随着安慰剂区重复次数 $m \to \infty$ ，最优模型最终收敛于在 $τ$ 的所有估计方式中选择具有最低 $M S E (\hat{τ})$ 的估计方式，连续安慰剂阈值之间的距离保持不变。

为了建立渐进最优性，要证明对于每一种候选估计模型，阈值估计的 $M S E$ 与每个相关安慰剂估计的 $M S E$ 相同，即证明：

M S E (\hat{τ}) = M S E (\hat{τ_{k}}), f o r a l l | k | > b (1)

其中，其中 $X = k$ 是每个安慰剂阈值的位置， $b$ 是用于此估计方式的带宽。对于每个候选估计方式，我们观察每个 $k$ 的 $\hat{τ_{k}}$ 。因为 $τ_{k} = 0$ ，所以 $M S E (\hat{τ_{k}}) = E (\hat{τ_{k}})$ 。如果等式 (1) 成立，那么 $E (\hat{τ_{k}}) E (^{τ} k)$ 与 $k$ 无关，且 $E (\hat{τ_{k}}) = M S E (\hat{τ})$ 。

本文证明上述等式 (1) 在以下假设条件下成立：

$X$ 均匀分布；
$Y_{i} (0)$ 和 $Y_{i} (1)$ 同方差；
$Y_{i} (0)$ 和 $Y_{i} (1)$ 关于 $X$ 的条件期望是连续的、平滑的，并且四阶导数为零；
平均处理效应关于 $X$ 的二阶导数为零。

考虑一组局部线性模型，每个估计量具有不同的对称带宽 $(b)$ 。 $\hat{τ_{b}}$ 表示带宽为 $b$ 的局部线性估计的处理效应。在该带宽内，在阈值的两侧均有 $b_{n}$ 个观测样本。估计的处理效应是 $\hat{τ_{b}} = \hat{α_{2 b}} - \hat{α_{1 b}}$ 。其中 $\hat{α_{2 b}}$ 和 $\hat{α_{1 b}}$ 可以使用阈值两侧的两个独立的线性回归来估计。简单起见，此处去掉下标 $i$ ，这些回归方程为：

y = α_{1} + β_{1} x + ε, - b < x < 0 (2)

y = α_{2} + β_{2} x + ε, 0 < x < b (3)

同理，将 $\hat{τ_{b k}}$ 定义为在安慰剂阈值 $X = K$ （真正的断点为零）处的断点估计， $\hat{τ_{k b}} = \hat{α_{k 2 b}} - \hat{α_{k 1 b}}$ ，其中 $\hat{α_{k 2 b}}$ 和 $\hat{α_{k 1 b}}$ 通过以下回归方程估计得出：

y = α_{k 1} + β_{1} (x - k) + ε, (k - b) < x < k (4)

y = α_{k 2} + β_{2} (x - k) + ε, k < x < (k + b) (5)

任意估计方式的处理效应的 $M S E$ 为：

M S E (\hat{τ}) = E (\hat{τ} - τ)^{2} = V a r (\hat{τ}) + B i a s (\hat{τ})^{2} (6)

对于每一个局部线性回归， $M S E (\hat{τ_{b}})$ 是两个回归（式2和式3）的方差及估计量偏差的函数：

M S E (\hat{τ_{b}}) = E (\hat{α_{2 b}} - \hat{α_{1 b}} - τ)^{2} = V a r (\hat{α_{1 b}}) + V a r (\hat{α_{2 b}}) + B i a s (\hat{τ_{b}})^{2} (7)

假设(3)和假设（4）意味着真正的DGP采用以下形式：

y = α + (τ + ϑ x) 1 (x > 0) + θ_{1} x + θ_{2} x^{2} + θ_{3} x^{3} + ε (8)

换言之，这些假设意味着CEF是全局立方的，在阈值处有潜在断点 $(τ)$ 和拐点 $(ϑ)$ 。

E ({\hat{α}}_{1 b} | x) = α + θ_{2} δ_{L 2 b} + θ_{3} δ_{L 3 b} (9)

其中， $δ_{L 2 b}$ 是 $x^{2}$ 在 $x$ 上映射的常数， $δ_{L 3 b}$ 是 $x^{3}$ 在 $x$ 上映射的常数，两者均满足 $- b ＜ x ＜ 0$ 。类似地，对于 ${\hat{α}}_{2 b}$ ，使用 $R H S$ 上的数据$（0

E ({\hat{α}}_{2 b} | x) = (α + τ) + θ_{2} δ_{R 2 b} + θ_{3} δ_{R 3 b} (10)

于是

E ({\hat{τ}}_{b} | x) = E ({\hat{α}}_{2 b} | x) - E ({\hat{α}}_{1 b} | x) = (α + τ) + θ_{2} δ_{R 2 b} + θ_{3} δ_{R 3 b} - (α + θ_{2} δ_{L 2 b} + θ_{3} δ_{L 3 b}) (11)

然而， $δ_{L 2 b}$ = $δ_{R 2 b}$ , $δ_{L 3 b}$ = $- δ_{R 3 b}$ 因此

E ({\hat{τ}}_{b} | x) = τ + 2 θ_{3} δ_{R 3 b} (12)

所以

B i a s ({\hat{τ}}_{b} | x) = τ + 2 θ_{3} δ_{R 3 b} (13)

这一偏差与 $θ_{3}$ 成正比，与真实DGP的任何其他参数无关。因此， $B i a s ({\hat{τ}}_{b} | x)$ 与DGP的CEF的三阶导数成正比。

现在证明，对于 $| k | \geq b$ ， $B i a s ({\hat{τ}}_{b} | x) = B i a s ({\hat{τ}}_{k b} | x)$ ，代入 $x_{k} = x - k$ ，方程（4）和（5）等价于：

y = α_{k 1} + β_{1} x_{k} + ε, - b ＜ x_{k} ＜ 0 (14)

y = α_{k 2} + β_{2} x_{k} + ε, 0 ＜ x_{k} ＜ b (15)

并且方程（A8）中的DGP可以表示为:

y = α + τ + (θ_{1} + ϑ) (x_{k} + k) + θ_{2} (x + k)^{2} + θ_{3} (x_{k} + k)^{3} + ε (16)

如果 $k \geq b$ ，

y = α + θ_{1} (x_{k} + k) + θ_{2} (x_{k} + k)^{2} + θ_{3} (x_{k} + k)^{3} + ε (17)

如果 $k \leq - b$ ，

y = π_{0} + π_{1} x_{k} + π_{2} x_{k}^{2} + θ_{3} x_{k}^{3} + ε (18)

其中， $π_{0} = α + τ + k ϑ + k θ_{1} + k^{2} θ_{2} + k^{3} θ_{3}$ ， $π_{1} = θ_{1} + ϑ + 2 k θ_{2} + 3 k^{2} θ_{3}$ ， $π_{2} = θ_{2} + 3 k θ_{3}$ 。

式（17）也可以做类似的展开。

式（14）、（15）、（18）分别等价于式（2）、（3）、（8）。其中， $x > 0$ 。同理，如果 $x < 0$ ，则阈值 $x_{k} = 0$ 。如前面所示，、RDD估计的偏差仅与真实DGP的CEF的三阶导数成正比。式（18）的三阶导数（ $6 θ_{3}$ ）与式（8）相同，因此，因此对于 $| k | > b ， B i a s (\hat{τ_{k b}}) = B i a s (\hat{τ_{b}})$ 。

在给定假设（1）和（2）的情况下，很容易证明 $V a r (\hat{τ_{k b}}) = V a r (\hat{τ_{b}})$ 。因此， $M S E (\hat{τ_{k b}}) = M S E (\hat{τ_{b}})$ ，在假设（1）-（4）下，本文的方法满足渐进最优性。

3.命令安装与介绍

3.1 命令安装

ssc install pzms, replace 
net get pzms.pkg, replace //下载作者提供的范例数据

3.2 命令介绍

该命令的基本语法为

pzms depvar runvar [if] [in] , maxbw(real) [options]

其中，

maxbw(real)表示用于估计的最大带宽

选择项的含义：

maxbw(#) : 在整个安慰剂区进行模型比较时考虑的最大带宽,没有默认值。
c(#): cut point, 阈值，默认值为0。
minbw(#):在整个安慰剂区进行模型比较时考虑的最大小带宽, 若未指定，则设置为 0.1*maxbw。
pzrange(##):安慰剂区的范围。在使用maxbw(#)进行范围划分后，c(#)须包含在其中。默认值是驱动变量的整个范围。
p(#):多项式阶数，默认为1，即一阶线性。若指定为2，则同时考虑一阶线性和二次项。
deriv(#)：导数阶数，采用精确RDD估计时默认为0，采用精确RKD估计时默认为1。
pzstepnum(#)：安慰剂区域内的迭代次数。若驱动变量观测值过少或pzstepnum设定值过大时，实际迭代次数可能少于指定次数。因此，迭代次数为安慰剂区域内的最大迭代次数。若已设定pzstepsize，则不能再设定pzstepnum。若pzstepnum和pzstepsize均未设定，则默认值为50。
pzstepsize(#)：相邻阈值之间的距离。若已设定pzstepnum，则不能再设定pzstepsize。
bwstepnum(#):介于最小带宽 minbw(#) 与最大带宽 maxbw(#) 之间的带宽数量，若pzstepnum和pzstepsize均未设定，则默认值为20。
bwstepsize(#): 候选估计模型之间的带宽增加（以驱动变量为单位）。
nolog:不显示迭代进度。
vce(vce type):设定标准误类型，默认为同方差标准误差。
covs（varlist）：设定候选模型中的协变量。若设定此项，则自动为每一个多项式p(#)进行模型估计。例如，若设定p(2)，则估计以下四个模型并进行比较：无协变量的线性模型、有协变量的非线性模型、无协变量二次模型和有协变量二次模型。
kernel(kernel type)：设定核密度类型，支持uniform和triangular，默认为uniform。
weight(varname)：指定分析权重aweights,除了mcustom#指定的模型之外，适用于所考虑的每个模型。
collapse(string):在算法运行前将数据折叠至驱动变量水平。若设定了 collapse(weight),将根据与折叠仓相对应的观测值数量对估计值进行加权。如果不希望在折叠数据时使用观测频率的权重，则选择collapse(noweight)。当驱动变量是离散值，且观测值较多，迭代次数较多时，设定collapse(weight)可以在不影响估计的情况下大大提高运行速度。若要引用其他分析权重，即使设定了collapse(weight)，也可以将collapse与weight组合使用。
bwlfix（#）：固定阈值左侧的带宽长度。若因数据限制使阈值左侧带宽较短，右侧带宽较长，或带宽不对称的情况下，则使用该选项。
bwrfix（#）：固定阈值右侧的带宽长度。若因数据限制使阈值右侧带宽较短，左侧带宽较长，或带宽不对称的情况下，则使用该选项。
pzplot：生成核密度图，描述最优模型中安慰剂估计值的分布。该图使用Stata kdensity命令和默认设置。
donut(# #): 删除指定范围内的观测值。第一个#指定在临界值左侧的驱动变量的该值内的观测值将被删除。第二个#指定在临界值右侧的驱动变量的该值内的观测值将被删除。
{opt mcustom#(# #, {it{mcustom_options})}:自定义模型。可以使用mcustom1()、mcustom2()至10个自定义模型。mcustom()中第一个参数是阈值值左侧的多项式阶数，第二个参数是阈值右侧的多项式阶数。例如，mcustom(1 2)表示左侧线性，右侧二项式的估计模型。多项式阶数可为0-10阶。其中，{mcustom_options} 包括：covs(varlist)、kernel(kernal type)、weight(varname)。

4. 示例

4.1 使用默认设置进行模型选择和估计

 . use "pzms_example_data.dta", clear
 . pzms y x, maxbw(0.99)

结果如下

. pzms y x, maxbw(0.99)

Results of Trial
Optimal bandwidth:      0.662
Polynomial order:       1
RMSE:                 .01358279
Effective sample size of placebo estimates from optimal model: 7.753

Estimated Treatment Effect Using Optimal Model   Observations: 1330
                                                 Threshold:       0
                                                 vce:              
-------------------------------------------------------------------
Inference method |  Coef.  Std. Err.   t     P>|t|      [95% CI]
-----------------+-------------------------------------------------
   Homoscedastic | .29204   .01644   17.768  0.0000  0.2598  0.3243
 KS rand. infer. | .29204     .013   22.466  0.0000  0.2611  0.3230
-------------------------------------------------------------------

4.2 使用二阶多项式和聚类标准误

 pzms y x, maxbw(0.99) p(2) vce(cluster x)

结果如下

Results of Trial
Optimal bandwidth:      0.662
Polynomial order:       1
RMSE:                 .01358279
Effective sample size of placebo estimates from optimal model: 7.753

Estimated Treatment Effect Using Optimal Model Observations:   1330
                                               Threshold:         0
                                               vce:       cluster x
-------------------------------------------------------------------
Inference method |  Coef.  Std. Err.   t     P>|t|  [95% CI]       
-----------------+-------------------------------------------------
          Robust | .29204   .01587     18.4  0.0000  0.2606  0.3234
 KS rand. infer. | .29204     .013   22.466  0.0000  0.2611  0.3230
-------------------------------------------------------------------

4.3 增加带宽数量、安慰剂断点之间的距离

pzms y x, maxbw(0.99) p(2) vce(cluster x) bwstepnum(25) pzstepsize(0.02)

结果如下

Results of Trial
Optimal bandwidth:      0.656
Polynomial order:       1
RMSE:                 .01343429
Effective sample size of placebo estimates from optimal model: 14.035

Estimated Treatment Effect with Optimal Model  Observations:   1310
                                               Threshold:         0
                                               vce:       cluster x
-------------------------------------------------------------------
Inference method |  Coef.  Std. Err.   t     P>|t| [95% CI]        
-----------------+-------------------------------------------------
          Robust | .29196   .01614   18.089  0.0000  0.2600  0.3239
 KS rand. infer. | .29196   .01235   23.645  0.0000  0.2653  0.3186
-------------------------------------------------------------------

4.4 考虑一个三角核的自定义线性模型，并绘制最优模型的安慰剂估计

. pzms y x, maxbw(0.99) p(2) vce(cluster x)  ///
          bwstepnum(25) pzstepsize(0.02)   ///
          mcustom(1 1, kernel(triangular)) pzplot

结果如下

. pzms y x, maxbw(0.99) p(2) vce(cluster x)  ///
          bwstepnum(25) pzstepsize(0.02)   ///
          mcustom(1 1, kernel(triangular)) pzplot

Results of Trial
Optimal bandwidth:      0.804
Optimal specification:  mcustom1
RMSE:                 .01307901
Effective sample size of placebo estimates from optimal model: 17.725


Estimated Treatment Effect wiht Optimal Model   Observations:  1610
                                                Threshold:        0
                                                vce:      cluster x
-------------------------------------------------------------------
Inference method |  Coef.  Std. Err.   t     P>|t|  [95% CI]       
-----------------+-------------------------------------------------
          Robust | .29787   .01592   18.712  0.0000  0.2664  0.3293
 KS rand. infer. | .29787    .0122    24.41  0.0000  0.2721  0.3236
-------------------------------------------------------------------

5. 参考资料

Kettlewell N. & Siminski P. (2022). Optimal Model Selection in RDD and Related Settings Using Placebo Zones. -PDF-

6. 相关推文

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

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