Stata：自己动手做组间系数差异检验-bootstrap-bdiff

发布时间：2022-05-18 阅读 6999

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

作者：李烨阳 (浙江大学)
邮箱：li_yeyang95@163.com

1. 引言

本文介绍了基于 Bootstrap (自抽样 / 自举法) 的组间系数检验方法及其 Stata 实现。具体思路如下：

第一种思路：首先通过有放回的自抽样方法获得一系列经验样本 (Empirical Sample)；然后在经验样本中根据其实际分组情况进行分组回归，从而获得分组回归系数差异统计量 $d$ 的经验分布；最后通过检验 0 在 $d$ 分布中的相对位置来检验 $H_{0} : d = 0$ 。

第二种思路：首先通过有放回的自抽样方法获得一系列经验样本；然后按照真实分组的比例，但随机的为每个样本分配组别，并进行分组回归，从而构造出随机分组情况下，组间系数差异统计量的 $d$ 经验分布。其背后的逻辑是，如果不存在组间系数差异，那么无论样本属于哪个组别， $x$ 对 $y$ 的边际影响都是相同的；最后检验实际观察到的组间系数差异 $\hat{d_{0}}$ 在 $d$ 分布中的相对位置来判断我们实际观察到 $\hat{d_{0}}$ 的概率，若小于给定置信水平，则拒绝 $H_{0} : d = 0$ 。

2. bootstrap 命令

bootstrap 命令语法如下：

bootstrap exp_list [, options eform_option] : command

其中，exp_list 是想要得到的统计量，也就是每次抽样之后计算并记录的统计指标，在本文中为 $x$ 在两组样本之间的系数差。command 可以是估计命令 (例如 regress)、非估计命令 (例如 summarize)、或者用户自己编写的命令。在本文中，需要进行分组回归，分别记录两次回归的系数结果并作差，因此无法直接通过一次估计命令实现，需要用户自己编写程序。

bootstrap 命令的核心作用是自动执行自抽样过程，并在经验样本中运行 command 命令，记录 exp_list 统计指标结果。关于自抽样方法以及编写 bootstrap 程序方法的详细内容，可以参考连享会推文「Stata: Bootstrap-自抽样-自举法」。

3. 第一种思路

第一种思路，我们借鉴了 Lu 等 (2019) 的做法，详见原文的 Table 10。

3.1 导入数据和变量设定

首先导入 Stata 自带的 nlsw88.dta 数据，并根据种族 (race) 生成分组变量。

. * 调用数据
. sysuse "nlsw88.dta", clear
. gen agesq = age*age    // 年龄平方项
. gen lnwage = log(wage) // wage 取对数

. * 设置分组变量
. drop if race == 3
. gen black = 2.race 

. * 设置自变量
. global xx "c.ttl_exp married south c.hours c.tenure c.age c.agesq i.industry"

由于白人组中没有处于 Mining (industry=2) 的观察值，可以根据研究设计选择或保留这部分样本。无论是否删除这部分样本，程序运行都不会报错，这里选择删除了这部分样本。

. drop if industry == 2

3.2 编写计算统计量的程序

我们将编写的程序命名为 bse。首先进行分组回归，需要注意的是，这里与第二种思路不同，我们是按照样本实际的组别进行分组回归。在每次回归之后，存储我们关注的自变量系数 (ttl_exp)。然后计算两个组别回归系数差，将其存储在矩阵中并返回至 e(b)。具体的 Stata 程序如下：

capture program drop bse
program bse, eclass
    *-分组回归
    reg dep $xx if black == 0
    scalar b1= _b[ttl_exp]
    reg dep $xx if black == 1
    scalar b2= _b[ttl_exp]

    *-计算组间系数差异
    scalar diff= b1- b2

    *-将组间系数差存储在矩阵中，设置列名方便调取
    matrix b = diff
    matrix colnames b = diff

    *-将组间系数差矩阵返回 e() 中
    ereturn post b
    ereturn display
end

3.3 完成自抽样并返回统计量

在 bse 程序中，因变量被要求命名为 dep，在执行 bootstrap 命令之前需要将因变量重命名为 dep。这里以 wage 和 lnwage 作为因变量为例，设置了循环方法。由于数据集中不能有重名变量，因此加入了 preserve-restore 命令。设置重复自抽样次数 (本例中为 1000)，为了结果可复现设置种子值 (本例中为 1234)。具体的 Stata 代码如下：

. local D "wage lnwage"
. foreach dep of local D{
  2.     preserve
  3.     rename `dep' dep
  4.     bse
  5.     bootstrap _b[diff], reps(1000) seed(1234) nowarn : bse
  6.     restore
  7. }
 
. * 变量 wage
Bootstrap results                                        Number of obs = 2,216
                                                         Replications  = 1,000
      Command: bse
        _bs_1: _b[diff]
------------------------------------------------------------------------------
             |   Observed   Bootstrap                         Normal-based
             | coefficient  std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
       _bs_1 |     -0.017      0.063    -0.27   0.790       -0.141       0.107
------------------------------------------------------------------------------
  
. * 变量 lnwage 
Bootstrap results                                        Number of obs = 2,216
                                                         Replications  = 1,000
      Command: bse
        _bs_1: _b[diff]
------------------------------------------------------------------------------
             |   Observed   Bootstrap                         Normal-based
             | coefficient  std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
       _bs_1 |     -0.003      0.007    -0.48   0.628       -0.016       0.010
------------------------------------------------------------------------------

3.4 统计量检验

需要说明的是，结果中得到的组间差异系数，是在抽样样本中分组回归得到的组间系数差 $\hat{d_{0}}$ ；标准误是在自抽样获得的经验样本中，计算的组间差异系数差标准误。因此，这里的 $z$ 值和 $p$ 值实际上是在预先假设组间系数差 $d$ 符合正态分布的基础上得到的结果。

更一般地，我们不对 $d$ 的分布进行预先假设，而是根据自抽样 $K$ 次得到的系数差统计量的经验分布 ${d^{S_{1}}, d^{S_{2}}, . . ., d^{S_{k}}}$ 。然后通过分析 0 在 $d$ 分布中的相对位置来检验 $H_{0} : d = 0$ 。

经验 $p$ 值的计算方法为：若 $\hat{d_{0}} > 0$ ，则 $\hat{p} = \frac{c o u n t i f (d^{S_{k}} < 0)}{K}$ 。其中 $c o u n t i f (d^{S_{k}} < 0)$ 代表自抽样获得的经验样本的组间系数差小于 0 的个数， $K$ 代表自抽样次数。若 $\hat{d_{0}} < 0$ ，则 $\hat{p} = \frac{c o u n t i f (d^{S_{k}} > 0)}{K}$ 。

为了计算经验 $p$ 值，我们需要对 3.3 部分进行修改，即在 bootstrap 命令中增加 saving() 选项，从而存储全部的 $d^{S_{k}} (k = 1, 2, . . . K)$ 结果。修改后的 Stata 程序如下：

. local D "wage lnwage"
. foreach dep of local D{
  2.     preserve
  3.     rename `dep' dep
  4.     bse
  5.     bootstrap _b[diff], reps(1000) seed(1234) saving(diff_`dep') nowarn : bse
  6.     restore
  7. }

计算经验 $p$ 值的 Stata 程序为：

. local D "wage lnwage"
. foreach dep of local D{
  2.     qui{
  3.         use diff_`dep' ,clear
  4.         count if _bs_1>0
  5.         local num = r(N)
  6.         local p = `num'/_N
  7.         if `p'>0.5{
  8.             local p = 1-`p' 
  9.         }
 10.     }
 11.     dis "`dep':`p'"
 12. }

wage:.379
lnwage:.271

本例中，因变量为 wage 时，组间系数差对应的经验 $p$ 值为 0.379。

4. 第二种思路

第二种思路除了可以自己写程序并结合bootstrap命令实现，也可以通过使用 bdiff 命令并加入 bsample 选项更便捷地实现。该思路的具体实现步骤，可以参考连享会推文「Stata: 如何检验分组回归后的组间系数差异？」。

4.1 编写计算统计量的程序

第二种思路与第一种思路的不同之处在于，需要对自抽样得到的经验样本进行随机分组。为了保证随机分组中两个组样本的比例不变，首先需要计算样本中两个组的比例 (由于只有两个分组，且自抽样时抽取与总样本等量的样本，因此只需要计算一个组的个数即可)。具体的 Stata 程序如下：

. * 调用数据
. sysuse "nlsw88.dta", clear
. gen agesq = age*age    // 年龄平方项
. gen lnwage = log(wage) // wage取对数

. * 设置分组变量
. drop if race==3
. gen black = 2.race 

. * 设置自变量
. global xx "c.ttl_exp married south c.hours c.tenure c.age c.agesq i.industry"  
. drop if industry == 2  

. * 获得实际分组中第一组的个数，用于保证随机分组中分组比例不变
. tab black, matcell(freq)
. global n1 = el(freq, 1, 1)

随机分组的具体实现方法：设置一个变量，该变量为每一个样本赋值一个 (0, 1) 之间的随机数。然后根据这个变量排序，将前 n1 (组别 1 的个数) 个样本分配至组别 1，其余分配至组别 2。剩余步骤与第一种思路相同。我们将编写的程序命名为 bse2。具体的 Stata 程序如下：

capture program drop bse2
program bse2, eclass
    *-生成随机变量，用于为样本随机分组
    tempvar random 
    qui gen `random' = runiform()
    sort `random'           

    *-分组回归
    reg dep $xx in 1/$n1
    scalar b1= _b[ttl_exp]
    reg dep $xx if ~e(sample)
    scalar b2= _b[ttl_exp]

    *-计算组间系数差异
    scalar diff= b1- b2

    *-将组间系数差存储在矩阵中，设置列名方便调取
    matrix b= diff
    matrix colnames b= diff

    *-将组间系数差矩阵返回e()中
    ereturn post b
    ereturn display
end

4.2 完成自抽样并返回统计量

与 3.3 部分内容类似，不同之处是将其中的 bse 命令更换为 bse2 命令，并修改存储统计量结果的数据集名称。具体的 Stata 程序如下：


. local D "wage lnwage"

. foreach dep of local D{
  2.     preserve
  3.     rename `dep' dep
  4.     bse
  5.     bootstrap _b[diff], reps(1000) seed(1234) saving(diff_`dep'_2) nowarn : bse2
  6.     restore
  7. }

. * 变量 wage 
Bootstrap results                                        Number of obs = 2,216
                                                         Replications  = 1,000
      Command: bse2
        _bs_1: _b[diff]
------------------------------------------------------------------------------
             |   Observed   Bootstrap                         Normal-based
             | coefficient  std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
       _bs_1 |      0.000      0.070     0.00   0.997       -0.137       0.137
------------------------------------------------------------------------------

. * 变量 lnwage 
Bootstrap results                                        Number of obs = 2,216
                                                         Replications  = 1,000
      Command: bse2
        _bs_1: _b[diff]
------------------------------------------------------------------------------
             |   Observed   Bootstrap                         Normal-based
             | coefficient  std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
       _bs_1 |     -0.001      0.007    -0.18   0.855       -0.015       0.013
------------------------------------------------------------------------------

需要注意，与第一种思路不同，这里不再是检验 0 在统计量经验分布中的相对位置，因此不可以直接用结果中的 $z$ 值和其对应的 $p$ 值。

4.3 统计量检验

由于这里的 $d^{S_{k}}$ 是根据随机分配的组别进行分组回归得到的系数差，而不是真实分组的系数差，因此计算经验 $p$ 值的方法有所不同， $\hat{p} = \frac{c o u n t i f (d^{S_{k}} > \hat{d_{0}})}{K}$ ，且由于 $d^{S_{k}}$ 分布未必是对称的， $\hat{p} = 0.03$ 与 $\hat{p} = 0.97$ 都可以视为在 5% 水平上拒绝原假设 $H_{0} : d = 0$ 的证据。计算经验 $p$ 值的 Stata 程序为：

. local D "wage lnwage"
. foreach dep of local D{
  2.     preserve
  3.     qui{
  4.         rename `dep' dep
  5.         reg dep $xx if black==0
  6.         scalar b1= _b[ttl_exp]
  7.         reg dep $xx if black==1
  8.         scalar b2= _b[ttl_exp]
  9.         global d0= b1- b2
 10.         use diff_`dep'_2 ,clear
 11.         count if _bs_1>$d0
 12.         local num = r(N)
 13.         local p = `num'/_N
 14.         if `p'>0.5{
 15.             local p = 1-`p' 
 16.         }
 17.     }
 18.     dis "`dep':`p'"
 19.     restore
 20. }
 
wage:.395
lnwage:.33

4.4 与 bdiff 命令的对比

第二种思路可以通过 bdiff 命令便捷的实现，仍然才用于前面相同的数据和模型设定一致，具体的 Stata 程序如下：

. * 导入数据和变量设定
. sysuse "nlsw88.dta", clear
. drop if race==3
. gen black = 2.race
. gen agesq = age*age
. gen lnwage = log(wage)
. // drop if industry==2  // 白人组中没有处于 Mining (industry=2) 的观察值
. tab industry, gen(d)    // 手动生成行业虚拟变量
. local dumind "d2 d3 d4 d5 d6 d7 d8 d9 d10 d11" // 行业虚拟变量
. global xx "c.ttl_exp married south c.hours c.tenure c.age c.agesq `dumind'"  

. * 使用 bdiff 命令实现思路二 
. // ssc install bdiff, replace
. bdiff, group(black) model(reg wage $xx) reps(1000) seed(1234) bsample

   Variables |      b0-b1    Freq     p-value
-------------+-------------------------------
     ttl_exp |     -0.028     656       0.344
     married |     -0.873     936       0.064
       south |      1.272      11       0.011
       hours |      0.016     247       0.247
      tenure |      0.031     289       0.289
         age |     -1.202     695       0.305
       agesq |      0.017     284       0.284
          d2 |      5.635     355       0.355
          d3 |     -1.159     704       0.296
          d4 |      0.258     364       0.364
          d5 |     -1.075     764       0.236
          d6 |      0.379     345       0.345
          d7 |      1.721     116       0.116
          d8 |      1.087     260       0.260
          d9 |      0.275     372       0.372
         d10 |      1.906     160       0.160
         d11 |      0.519     232       0.232
       _cons |     20.633     333       0.333
---------------------------------------------

. bdiff, group(black) model(reg lnwage $xx) reps(5000) seed(1234) bsample

   Variables |      b0-b1    Freq     p-value
-------------+-------------------------------
     ttl_exp |     -0.006    3968       0.206
     married |     -0.045    4033       0.193
       south |      0.221       0       0.000
       hours |     -0.000    2750       0.450
      tenure |      0.004    1069       0.214
         age |     -0.100    3419       0.316
       agesq |      0.002    1436       0.287
          d2 |      0.256    2130       0.426
          d3 |     -0.127    3529       0.294
          d4 |      0.089     826       0.165
          d5 |     -0.078    3654       0.269
          d6 |      0.089     870       0.174
          d7 |      0.063    1491       0.298
          d8 |      0.329     111       0.022
          d9 |      0.122     820       0.164
         d10 |      0.553     173       0.035
         d11 |      0.152     149       0.030
       _cons |      1.558    1766       0.353
---------------------------------------------

结果与编写的 Bootstrap 程序基本一致，但由于程序运行中涉及到随机数的生成，虽然设置了相同的种子值，但具体实现步骤不同，结果仍有所差异。

5. 参考资料

连玉君, 廖俊平. 如何检验分组回归后的组间系数差异?[J]. 郑州航空工业管理学院学报, 2017, 35(6): 97-109. -Link-
Lu Y, Wang J, Zhu L. Place-based policies, creation, and agglomeration economies: Evidence from China's economic zone program[J]. American Economic Journal: Economic Policy, 2019, 11(3): 325-60. -Link-

6. 相关推文

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

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