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

文摘   教育   2024-11-05 22:00   中国  

👇 连享会 · 推文导航 | www.lianxh.cn

🍓 课程推荐:2024 机器学习与因果推断专题
主讲老师:司继春 (上海对外经贸大学) ;张宏亮(浙江大学)
课程时间:2024 年 11 月 9-10 日 ;16-17日
课程咨询:王老师 18903405450(微信)

 课程特色 · 2024机器学习与因果推断

  • 懂原理、会应用。本次课程邀请了两位老师合作讲授,目的在于最大限度地实现理论与应用的有机结合。为期四天的课程,分成两个部分:第一部分讲解常用的机器学习算法和适用条件,以及文本分析和大语言模型;第二部分通过精讲 4-6 篇发表于 Top 期刊的论文,帮助大家理解各类机器学习算法的应用场景,以及它们与传统因果推断方法的巧妙结合。
  • 以 Top 期刊论文为范例。目前多数人的困惑是不清楚如何将传统因果推断方法与机器学习结合起来。事实上,即便是 MIT 和 Harvard 的大牛们也都在「摸着石头过河」。为此,通过论文精讲和复现来学习这部分内容或许是目前最有效的方式了。张宏亮老师此前在浙江大学按照这一模式教授了「因果推断和机器学习」课程,效果甚佳:学生们能够逐渐建立起研究设计的理念,并在构造识别策略时适当地嵌入机器学习方法。

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

温馨提示: 文中链接在微信中无法生效。请点击底部「阅读原文」。或直接长按/扫描如下二维码,直达原文:

   


目录

  • 1. 引言

  • 2. bootstrap 命令

  • 3. 第一种思路

    • 3.1 导入数据和变量设定

    • 3.2 编写计算统计量的程序

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

    • 3.4 统计量检验

  • 4. 第二种思路

    • 4.1 编写计算统计量的程序

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

    • 4.3 统计量检验

    • 4.4 与 bdiff 命令的对比

  • 5. 参考资料

  • 6. 相关推文



1. 引言

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

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

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

2. bootstrap 命令

bootstrap 命令语法如下:

bootstrap exp_list [, options eform_option] : command

其中,exp_list 是想要得到的统计量,也就是每次抽样之后计算并记录的统计指标,在本文中为 在两组样本之间的系数差。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 统计量检验

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

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

经验 值的计算方法为:若 ,则 。其中 代表自抽样获得的经验样本的组间系数差小于 0 的个数, 代表自抽样次数。若 ,则

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

计算经验 值的 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 时,组间系数差对应的经验 值为 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 在统计量经验分布中的相对位置,因此不可以直接用结果中的 值和其对应的 值。

4.3 统计量检验

由于这里的 是根据随机分配的组别进行分组回归得到的系数差,而不是真实分组的系数差,因此计算经验 值的方法有所不同,,且由于 分布未必是对称的, 都可以视为在 5% 水平上拒绝原假设 的证据。计算经验 值的 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

  • 专题:数据处理
    • Stata:原始聚类自助法(wild cluster bootstrap)-boottest
    • Stata数据处理:用-astile-快速创建分组
  • 专题:Stata绘图
    • forest-森林图:分组回归系数可视化
    • Stata绘图:用-bytwoway-实现快速分组绘图
  • 专题:Stata程序
    • Stata: Bootstrap-自抽样-自举法
    • Stata:runby - 一切皆可分组计算!
  • 专题:结果输出
    • sumup:快速呈现分组统计量
  • 专题:回归分析
    • Stata:分组回归系数比较的新思路
    • Stata: 如何检验分组回归后的组间系数差异?
    • Stata: 获取分组回归系数的三种方式
  • 专题:倍分法DID
    • 倍分法:DID是否需要随机分组?

🍓 课程推荐:2024 机器学习与因果推断专题
主讲老师:司继春 (上海对外经贸大学) ;张宏亮(浙江大学)
课程时间:2024 年 11 月 9-10 日 ;16-17日
课程咨询:王老师 18903405450(微信)

尊敬的老师 / 亲爱的同学们:

连享会致力于不断优化和丰富课程内容,以确保每位学员都能获得最有价值的学习体验。为了更精准地满足您的学习需求,我们诚挚地邀请您参与到我们的课程规划中来。请您在下面的问卷中,分享您 感兴趣的学习主题或您希望深入了解的知识领域 。您的每一条建议都是我们宝贵的资源,将直接影响到我们课程的改进和创新。我们期待您的反馈,因为您的参与和支持是我们不断前进的动力。感谢您抽出宝贵时间,与我们共同塑造更加精彩的学习旅程!https://www.wjx.cn/vm/YgPfdsJ.aspx# 再次感谢大家宝贵的意见!

New! Stata 搜索神器:lianxh 和 songbl  GIF 动图介绍
搜: 推文、数据分享、期刊论文、重现代码 ……
👉 安装:
  . ssc install lianxh
  . ssc install songbl
👉  使用:
  . lianxh DID 倍分法
  . songbl all

🍏 关于我们

  • 连享会 ( www.lianxh.cn,推文列表) 由中山大学连玉君老师团队创办,定期分享实证分析经验。
  • 直通车: 👉【百度一下: 连享会】即可直达连享会主页。亦可进一步添加 「知乎」,「b 站」,「面板数据」,「公开课」 等关键词细化搜索。

连享会
连玉君老师团队分享,主页:lianxh.cn。白话计量,代码实操;学术路上,与君同行。
 最新文章