Stata连享会 主页 || 视频 || 推文 || 知乎 || Bilibili 站
温馨提示: 定期 清理浏览器缓存,可以获得最佳浏览体验。
New!
lianxh
命令发布了:
随时搜索推文、Stata 资源。安装:
. ssc install lianxh
详情参见帮助文件 (有惊喜):
. help lianxh
连享会新命令:cnssc
,ihelp
,rdbalance
,gitee
,installpkg
⛳ Stata 系列推文:
作者:李烨阳 (浙江大学)
邮箱:li_yeyang95@163.com
目录
本文介绍了基于 Bootstrap (自抽样 / 自举法) 的组间系数检验方法及其 Stata 实现。具体思路如下:
第一种思路:首先通过有放回的自抽样方法获得一系列经验样本 (Empirical Sample);然后在经验样本中根据其实际分组情况进行分组回归,从而获得分组回归系数差异统计量
第二种思路:首先通过有放回的自抽样方法获得一系列经验样本;然后按照真实分组的比例,但随机的为每个样本分配组别,并进行分组回归,从而构造出随机分组情况下,组间系数差异统计量的
bootstrap 命令语法如下:
bootstrap exp_list [, options eform_option] : command
其中,exp_list
是想要得到的统计量,也就是每次抽样之后计算并记录的统计指标,在本文中为 command
可以是估计命令 (例如 regress
)、非估计命令 (例如 summarize
)、或者用户自己编写的命令。在本文中,需要进行分组回归,分别记录两次回归的系数结果并作差,因此无法直接通过一次估计命令实现,需要用户自己编写程序。
bootstrap
命令的核心作用是自动执行自抽样过程,并在经验样本中运行 command
命令,记录 exp_list
统计指标结果。关于自抽样方法以及编写 bootstrap
程序方法的详细内容,可以参考连享会推文「Stata: Bootstrap-自抽样-自举法」。
第一种思路,我们借鉴了 Lu 等 (2019) 的做法,详见原文的 Table 10。
首先导入 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
我们将编写的程序命名为 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
在 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
------------------------------------------------------------------------------
需要说明的是,结果中得到的组间差异系数,是在抽样样本中分组回归得到的组间系数差
更一般地,我们不对
经验
为了计算经验 bootstrap
命令中增加 saving()
选项,从而存储全部的
. 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. }
计算经验
. 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 时,组间系数差对应的经验
第二种思路除了可以自己写程序并结合bootstrap
命令实现,也可以通过使用 bdiff
命令并加入 bsample
选项更便捷地实现。该思路的具体实现步骤,可以参考连享会推文「Stata: 如何检验分组回归后的组间系数差异?」。
第二种思路与第一种思路的不同之处在于,需要对自抽样得到的经验样本进行随机分组。为了保证随机分组中两个组样本的比例不变,首先需要计算样本中两个组的比例 (由于只有两个分组,且自抽样时抽取与总样本等量的样本,因此只需要计算一个组的个数即可)。具体的 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
与 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 在统计量经验分布中的相对位置,因此不可以直接用结果中的
由于这里的
. 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
第二种思路可以通过 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 程序基本一致,但由于程序运行中涉及到随机数的生成,虽然设置了相同的种子值,但具体实现步骤不同,结果仍有所差异。
Note:产生如下推文列表的 Stata 命令为:
lianxh 分组 bootstrap, m
安装最新版lianxh
命令:
ssc install lianxh, replace
免费公开课
最新课程-直播课
专题 | 嘉宾 | 直播/回看视频 |
---|---|---|
⭐ 最新专题 | 文本分析、机器学习、效率专题、生存分析等 | |
研究设计 | 连玉君 | 我的特斯拉-实证研究设计,-幻灯片- |
面板模型 | 连玉君 | 动态面板模型,-幻灯片- |
面板模型 | 连玉君 | 直击面板数据模型 [免费公开课,2小时] |
⛳ 课程主页
⛳ 课程主页
关于我们
课程, 直播, 视频, 客服, 模型设定, 研究设计, stata, plus, 绘图, 编程, 面板, 论文重现, 可视化, RDD, DID, PSM, 合成控制法
等
连享会小程序:扫一扫,看推文,看视频……
扫码加入连享会微信群,提问交流更方便
✏ 连享会-常见问题解答:
✨ https://gitee.com/lianxh/Course/wikis
New!
lianxh
和songbl
命令发布了:
随时搜索连享会推文、Stata 资源,安装命令如下:
. ssc install lianxh
使用详情参见帮助文件 (有惊喜):
. help lianxh