Stata：ritest-随机推断（Randomization Inference）

Stata 连享会   主页 || 视频 || 推文

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

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

✨ 课程主页： https://gitee.com/arlionn/Course

作者： 王如玉（中国人民银行金融研究所）
E-mail: 15242968@qq.com

Source: 本文部分内容摘译自以下文章：

• Mckenzie, David. Finally, a way to do easy randomization inference in Stata, October 02, 2017. URL
• Kerwin, Jason. Randomization inference vs. bootstrapping for p-values, September 05, 2017. URL

目录

1. 背景

「随机推断（Randomization Inference）」被认为是分析随机实验数据的一种较好的方法。尤其在样本量较少、具有聚类随机性、以及高杠杆点 (潜在离群值) 等情形下，其优势更加明显。

遗憾的是，Stata 中一直没有相应的程序实现该方法，进而限制了其在经济学中的应用。直到 2017 年，一位法兰克福大学的博士生 Simon Heß 编写了 ritest 命令，使得「随机推断」具有了可操作性。

关于 ritest 命令更为详细介绍，请参考 Heß (2017)³。接下来，本文将对「随机推断」和 ritest 命令进行简单介绍。

2. 什么是随机推断？

我们通常会面临这样一个难题，在对给定地区所有学校进行一项实验，或美国各州进行一项政策的实验结果评估时，「标注误差」和「$p$」的实际意义是什么。看似简答，实则不然，毕竟，这里不存在采样误差，而我们通常在回归分析中使用的推断方法都是基于采样误差的。为此，我们引入「随机推断」对上述问题予以解答。

2.1 随机推断的定义

「随机推断」考虑了在随机实验中，随机本身引起的任何变化。当研究者控制实验分配时，差异来自对实验组的随机分配，而不是抽样策略。此时，「随机推断」不仅要考虑随机分配下会发生什么，还包括在所有可能的随机分配下会发生什么，即实验结果是否都会成立？

2.2 随机推断的步骤

1. 保留样本的原始实验分配；
2. 根据原始实验组分配方式，重新随机分配实验组；
3. 用 "假的"，或者说 "安慰剂" 实验组作为一个附加项，重新估计原始回归方程；
4. 重复以上 2-3 过程；
5. 随机推断 [$p$] 值是安慰剂处理效果大于估计实验效果的次数的比例。

3. Stata 实操：ritest 命令

我们可以通过 ssc install ritest, replace 安装该命令。

3.1 ritest 具体语法格式

ritest resampvar exp_list [, options]: command

• resampvar 为需要重新抽样的变量名，例如，处理组指示变量 treatment；
• exp_list 为需要比较的表达式，例如，回归系数 _b[treatment]；
• options 定义了重新抽样过程的细节、对 $p$ 值的计算、结果输出等。其中部分参数具体说明如下：

  • reps(#) 执行随机抽样的次数，默认为 100 次；
  • left|right 计算单边 $p$ 值，默认为双边；
  • level(#) 置信区间，默认为 95；
  • fixlevels() 将重新随机过程限制于部分实验变量，可被用于多重实验组成对测试的情况；
  • seed(#) 设定随机数种子。
• command 为每次重复时执行的回归程序，例如 regress testscore treatment age。

3.2 ritest 三种不同的重新抽样方法

方法1：自动重新采样

ritest resampvar exp_list, reps(#) strata(varlist) cluster(varlist): command

主要参数说明：

• strata(varlist) 在分层内重新排列 resampvar；
• cluster(varlist) 在聚类内保持 resampvar 为常数。

该语句根据分层 strata() 和聚类 cluster()，随机排列resampvar $n$ 次，并且每次都执行语句并收集 exp_list 中变量的实现值。不设置分层即假设所有观测值在同一个单一分层中，不设置聚类即假设每个观测值为单一聚类。该方法最简单但最受限的抽样方法，在连续变量时可能会存在很大的局限，建议使用如下两种方法之一。

方法2：基于文件的重新抽样

ritest resampvar exp_list, reps(#) samplingsourcefile(filename) samplingmatchvar(varlist): command

主要参数说明：

• samplingsourcefile(filename) 从 Stata 数据文件 filename 中获取 resampvar 并重新排列，文件包含命名为resampvar1, resampvar2, resampvar3, ... 的变量；

• samplingmatchvar(varlist) 将 samplingsourcefile() 中的排列与现有数据合并，用 varlist 中的变量 (1:1 或 m:1)。

  该语句用 samplingsourcefile() 中的文件，根据 samplingmatchvar() 中列明的解释变量，将数据合并 $n$ 次 (1:1，如果解释变量非唯一，则 m:1)。每次都执行语句，并用 resampvar1, resampvar2, resampvar3, ... 逐次代换 resampvar。samplingsourcefile() 必须在执行语句前手动创建。

方法3：基于程序的重新抽样

. ritest resampvar exp_list, reps(#) samplingprogram(progname) samplingprogramoptions(string): command

主要参数说明：

• samplingprogram(programname) 根据用户编写的程序 programname 对 resampvar 重新排列；
• samplingprogramoptions(string) 可选，将 string 作为选项传递给 programname。

该语句通过执行 samplingprogram() 并将 samplingprogramoptions()作为其可选项，重新排列 resampvar。这是最多样化，应用最广的方法。如果实验的原始随机过程是由一个 do 文件生成的，原始程序可以用于 samplingprogram()。通过编写程序可以一次重新抽样多个变量。

4. 应用案例

4.1 对分层内的个体随机分组

Mckenzie (2017) ¹ 中提供了「Identifying and Spurring High-Growth Entrepreneurship: Experimental Evidence from a Business Plan Competition⁴」案例的「完整数据和代码」，因此，我们将以该案例为例予以演示说明。

这个实验对尼日利亚「一个青年创业竞赛对就业和公司利润影响」进行了跟踪调研。实验随机选取该竞赛的获奖者，分配为实验组与对照组，奖金根据他们的商业计划的实施情况，分四个阶段支付，并进行了三轮跟踪问卷调查。

例如，为了研究第二轮调查时还在运行的公司的阶段利润，在 Stata 中使用如下命令：

use "Nigeriablogdata.dta", clear
areg s_prof_trunc assigntreat, a(strata) robust
ritest assigntreat _b[assigntreat], reps(5000) strata(strata) seed(125): areg s_prof_trunc assigntreat, a(strata) robust

首先进行控制随机分层的 OLS 回归，其中 s_prof_trunc 为阶段利润，assigntreat 为实验组控制变量，得到该回归的 $p$ 值为 0.051。

然后使用 ritest，告诉它控制变量是什么，想要什么系数的 $p$ 值。之后，在随机分层中进行随机推断，并重复 5000 次，输出结果如下：

..................................................  5000

      command:  areg s_prof_trunc assigntreat, a(strata) robust
        _pm_1:  _b[assigntreat]
  res. var(s):  assigntreat
   Resampling:  Permuting assigntreat
Clust. var(s):  __000000
     Clusters:  541
Strata var(s):  strata
       Strata:  12

------------------------------------------------------------------------------
T            |     T(obs)       c       n   p=c/n   SE(p) [95% Conf. Interval]
-------------+----------------------------------------------------------------
       _pm_1 |   69.23354     256    5000  0.0512  0.0031  .0452543   .0576763
------------------------------------------------------------------------------
Note: Confidence interval is with respect to p=c/n.
Note: c = #{|T| >= |T(obs)|}

此处 $p$ 值为 0.0512 ，与原始回归的结果非常相近。这是个体级别的随机化，样本量为 497 ，因此我们可以认为随机推断的结果与回归结果相似。

注意： 此处原作者使用的原版 ritest 命令可能会忽略随机种子，即每次随机得出的结果不同。原作者使用 1000 次重复和种子 125 ，获得的 $p$ 值为 0.47 。更新版本修正了这一问题，因此笔者根据原数据及 do 文件运行获得的 $p$ 值有变化，为 0.37 。使用 5000 次重复获得的结果更接近于回归所得的 $p$ 值。

4.2 对有多个时间段和聚类标准差的个体的随机分组

此处同样使用以上数据，但对其进行了重新调整，以便估计第二轮和第三轮随访数据的合并效应。OLS 回归得出的 $p$ 值为 0.096。然后，增加 cluster 选项 (以 uid 为个体标识符)，以便在进行重新排列时将聚类性质考虑在内：

//合并第二、三轮数据
gen prof_trunc2=s_prof_trunc
gen prof_trunc3=t_prof_trunc
gen uid=_n
keep uid prof_trunc2 prof_trunc3 group existing assigntreat strata
reshape long prof_trunc ,i(uid) j(time)

gen time3=time==3
areg prof_trunc assigntreat time3 if group<=2 & existing==1, a(strata) robust cluster(uid) //回归
ritest assigntreat _b[assigntreat], reps(1000) strata(strata) cluster(uid) seed(124): areg prof_trunc assigntreat time3 if group<=2 & existing==1, a(strata) robust cluster(uid)

可以看出，$p$ 值为 0.118 (注：原作者在此处得到的结果为 0.100 ，原因同上)，与估计的回归值非常相似。如果没有指定 cluster() 选项，$p$ 值则为 0.068。

4.3 对有分层的聚类进行随机分组

本例来自 Mckenzie (2017)¹ 的一篇关于缩短供应链的工作论文，并提供了「完整数据和代码」。在该例中，公司所在市场被划分为 63 个区域，这些区域又被划分为 31 个层次 (30 对和一个三元组)，并在市场区域层级中产生集群随机化。我们想要知道，一周中有多少天公司在中央市场购物。在具有聚类标准差的回归中获得的 $p$ 值是 0.024 。由于随机推断可以使聚类相对较少的情况下聚类随机产生更大的变化，此处又能得到什么不同呢？使用如下命令：

use "clusterColombia.dta", clear
areg dayscorab b_treat b_dayscorab miss_b_dayscorab round2 round3, cluster(b_block) a(b_pair)
ritest b_treat _b[b_treat], reps(5000) cluster(b_block) strata(b_pair) seed(546): areg dayscorab b_treat b_dayscorab miss_b_dayscorab round2 round3, cluster(b_block) a(b_pair)

获得结果如下：

..................................................  5000

      command:  areg dayscorab b_treat b_dayscorab miss_b_dayscorab round2 round3, cluster(b_block) a(b_pair)
        _pm_1:  _b[b_treat]/_se[b_treat]
  res. var(s):  b_treat
   Resampling:  Permuting b_treat
Clust. var(s):  b_block
     Clusters:  63
Strata var(s):  b_pair
       Strata:  31

------------------------------------------------------------------------------
T            |     T(obs)       c       n   p=c/n   SE(p) [95% Conf. Interval]
-------------+----------------------------------------------------------------
       _pm_1 |  -2.312006     494    5000  0.0988  0.0042  .0906651   .1074066
------------------------------------------------------------------------------
Note: Confidence interval is with respect to p=c/n.
Note: c = #{|T| >= |T(obs)|}

可见，随机推断确实会产生很大的不同：将 $p$ 值从 0.024 增加到 0.0988 (注：原作者在此处得到的结果为 0.106，原因同上)。因此，如果要进行聚类随机化，那么使用随机推断就显得特别重要。

5. 注意事项

1. 抽样多少次为宜？默认值为 100 次，但在许多情况下，这可能太低了。Alwyn Young 使用 10,000 次重复，并发现一旦超过 2000 次，结果差别就很小；
2. 花费的时间。在本文例子的简单情况下，即从数据集中删除了所有其他变量后，使用笔记本电脑运行 5000 次大约需要 3 分钟。如果使用完整的数据集，时间将会大大加长。如果要加速运行，则需从数据集中删除回归所需的所有其他数据；
3. 该命令不允许在回归中使用 xi。因此，如果回归已使用 xi，则删除它。Jason Kerwin (2017)² 指出，可以只使用 i.groupvariable 这种形式来添加组变量固定效果，而无需使用 xi；
4. 关于更复杂的随机过程，例如两级聚类随机，则需要做更多的工作。以Mckenzie (2017) ¹ 中提供的「Identifying and Spurring High-Growth Entrepreneurship: Experimental Evidence from a Business Plan Competition⁴」案例为例，要实现双重随机化，首先在市场分层层面，将市场随机分为实验组或对照组市场，然后个体分层层面，在实验组市场内，将公司分别随机分配是否接受培训。这时就需要通过调用自行编写的指定随机分配方式的程序，或使用自行创建的包括重新分配方式的文件来处理此问题；
5. 关于多重实验组。标准教科书中关于随机推断的解释通常只涉及一个二进制实验，一个单元不是属于实验组就是属于控制组。在严格零假设下，即所有实验效果均为零，那么无论将其分配给实验组还是对照组，每个单元的结果都应相同，因此我们可以仅估算所有可能的随机排列下的实验效果，并计算有多少种排列得到的系数至少与实际数据所得的系数一样大。但是，在「Teaching personal initiative beats traditional training in boosting small business in West Africa⁵」实验中，作者将公司分为对照组 (T = 0)，传统商业培训(T = 1) 和个人主动性培训 (T = 2)。此时，我们至少要检验三个假设：传统培训没有效果，个人主动性培训没有效果，任何一种培训都没有效果，以上可以通过更新版本的 fixlevels() 选项来实现。

参考文献

相关课程

连享会-直播课 上线了！
http://lianxh.duanshu.com

免费公开课：

• 直击面板数据模型 - 连玉君，时长：1小时40分钟
• Stata 33 讲 - 连玉君, 每讲 15 分钟.
• 部分直播课 课程资料下载 (PPT，dofiles等)

课程一览

支持回看，所有课程可以随时购买观看。

Note: 部分课程的资料，PPT 等可以前往 连享会-直播课 主页查看，下载。

关于我们

• Stata连享会 由中山大学连玉君老师团队创办，定期分享实证分析经验。直播间 有很多视频课程，可以随时观看。
• 连享会-主页 和 知乎专栏，300+ 推文，实证分析不再抓狂。
• 公众号推文分类： 计量专题 | 分类推文 | 资源工具。推文分成 内生性 | 空间计量 | 时序面板 | 结果输出 | 交乘调节 五类，主流方法介绍一目了然：DID, RDD, IV, GMM, FE, Probit 等。
• 公众号关键词搜索/回复 功能已经上线。大家可以在公众号左下角点击键盘图标，输入简要关键词，以便快速呈现历史推文，获取工具软件和数据下载。常见关键词：课程, 直播, 视频, 客服, 模型设定, 研究设计, stata, plus, 绘图, 编程, 面板, 论文重现, 可视化, RDD, DID, PSM, 合成控制法 等

连享会小程序：扫一扫，看推文，看视频……

扫码加入连享会微信群，提问交流更方便

✏ 连享会学习群-常见问题解答汇总：
✨ https://gitee.com/arlionn/WD