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 系列推文:

PDF下载 - 推文合集

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


目录


1. 引言

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

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

第二种思路:首先通过有放回的自抽样方法获得一系列经验样本;然后按照真实分组的比例,但随机的为每个样本分配组别,并进行分组回归,从而构造出随机分组情况下,组间系数差异统计量的 d 经验分布。其背后的逻辑是,如果不存在组间系数差异,那么无论样本属于哪个组别,x 对 y 的边际影响都是相同的;最后检验实际观察到的组间系数差异 d0^ 在 d 分布中的相对位置来判断我们实际观察到 d0^ 的概率,若小于给定置信水平,则拒绝 H0: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。这里以 wagelnwage 作为因变量为例,设置了循环方法。由于数据集中不能有重名变量,因此加入了 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 统计量检验

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

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

经验 p 值的计算方法为:若 d0^>0,则 p^=countif(dSk<0)K。其中 countif(dSk<0) 代表自抽样获得的经验样本的组间系数差小于 0 的个数,K 代表自抽样次数。若 d0^<0,则 p^=countif(dSk>0)K

为了计算经验 p 值,我们需要对 3.3 部分进行修改,即在 bootstrap 命令中增加 saving() 选项,从而存储全部的 dSk(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 统计量检验

由于这里的 dSk 是根据随机分配的组别进行分组回归得到的系数差,而不是真实分组的系数差,因此计算经验 p 值的方法有所不同,p^=countif(dSk>d0^)K,且由于 dSk 分布未必是对称的,p^=0.03 与 p^=0.97 都可以视为在 5% 水平上拒绝原假设 H0: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小时]
  • Note: 部分课程的资料,PPT 等可以前往 连享会-直播课 主页查看,下载。

课程主页

课程主页

关于我们

  • Stata连享会 由中山大学连玉君老师团队创办,定期分享实证分析经验。
  • 连享会-主页知乎专栏,700+ 推文,实证分析不再抓狂。直播间 有很多视频课程,可以随时观看。
  • 公众号关键词搜索/回复 功能已经上线。大家可以在公众号左下角点击键盘图标,输入简要关键词,以便快速呈现历史推文,获取工具软件和数据下载。常见关键词:课程, 直播, 视频, 客服, 模型设定, 研究设计, stata, plus, 绘图, 编程, 面板, 论文重现, 可视化, RDD, DID, PSM, 合成控制法

连享会小程序:扫一扫,看推文,看视频……

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

✏ 连享会-常见问题解答:
https://gitee.com/lianxh/Course/wikis

New! lianxhsongbl 命令发布了:
随时搜索连享会推文、Stata 资源,安装命令如下:
. ssc install lianxh
使用详情参见帮助文件 (有惊喜):
. help lianxh