作者： 吴雄(湘潭大学)，童天天(中南财经政法大学)
连享会
Source: The Bootstrap in Stata
原文链接： 连享会-Bootstrap简介
文章目录
- 1. Bootstrap 简介
- 1.1 有放回抽样
- 1.2 标准差与标准误
- 2. 编写 bootstrap程序
- 2.1 Stata 范例 1：OLS 回归的 RMSE 的标准误
- 2.2 Stata 范例 2：采用 Bootstrap 获取 VIF 的标准误和置信区间
- 3. 主要参考文献
- Books
- Bootstrap 简介
- 聚类标准误-Bootstrap
- 关于我们
1. Bootstrap 简介
bootstrap 是一种崭新的增广样本统计方法，为解决小样本问题提供了很好的思路。它是非参数统计中一种重要的估计统计量方差进而进行区间估计的统计方法。对于回归模型：对于线性回归模型：
y t = X t β + u t , E ( u t ∣ X t ) = 0 , E ( u s u t = 0 ) ∀ s ≠ t y_t = X_t β+u_t, \\ E(u_t|X_t)=0,\ E(u_s u_t=0) \ ∀\ s≠t yt=Xtβ+ut,E(ut∣Xt)=0, E(usut=0) ∀ s=t
可以通过多种方法来建立 bootstrap 的数据生成过程 (DGP) 。所谓的 bootstrap DGP 是对未知的 「真实 DGP」的一种估计。如果 bootstrap DGP 在某种意义上接近真实的 DGP，那么由 bootstrap DGP 生成的数据将与真实 DGP 生成的数据相似(如果已知的话)。如果是这样，则进行模拟使用 bootstrap DGP 获得的 P 值与真实 P 值足够接近，可以进行准确的推理。
Bootstrap 的基本思想是：如果 观测样本 是从母体中随机抽取的，那么它将包含母体的全部的信息，那么我们不妨就把这个观测样本视为 “总体”。可以简单地概括为：既然样本是抽出来的，那我何不从样本中再抽样。
具体而言，Bootstrap 的第一步是生成一系列 bootstrap 经验样本 (Empirical Sample) (有时也被形象地称为 「伪样本」)，每个样本都是初始数据的一次有放回抽样。通过对 经验样本 的计算，获得统计量的分布。例如，要进行 1000 次 bootstrap，求平均值的置信区间，可以对每个经验样本 计算平均值。这样就获得了 1000 个平均值。对这 1000 个平均值的分位数进行计算， 即可获得置信区间。已经证明，在初始样本足够大且初始样本是从母体中随机抽取的情况下，bootstrap 抽样能够无偏接近总体的分布。
Bootstrap 的基本步骤如下：
- Step 1: 采用有放回抽样方法从原始样本中抽取一定数量的子样本。
- Step 2: 根据抽出的样本计算想要的统计量。
- Step 3: 重复前两步 K 次，得到 K 个统计量的估计值。
- Step 4: 根据 K 个估计值获得统计量的分布，并计算置信区间。
1.1 有放回抽样
所谓 「有放回抽样」 (Samping with replacement) 指的是在逐个抽取个体时，每次被抽到的个体放回总体中后，再进行下次抽取的抽样方法。
举个例子，对于由 0.1 和 0.3 这两个数字构成的观测样本而言， 记为 S 0 = ( 0.1 , 0.3 ) S_0 = (0.1, 0.3) S0=(0.1,0.3)。则采用有放回抽样 (Bootstrapping)，可以得到如下三种不同的经验样本： S 1 B S = ( 0.1 , 0.1 ) S_1^{BS} = (0.1,0.1) S1BS=(0.1,0.1), S 2 B S = ( 0.1 , 0.3 ) S_2^{BS} = (0.1,0.3) S2BS=(0.1,0.3), S 3 B S = ( 0.3 , 0.3 ) S_3^{BS} = (0.3,0.3) S3BS=(0.3,0.3)。
实际应用过程中，对于包含 N 个观察值的 观测样本 而言，Bootstrap 抽样得到的经验样本也会包含 N 个观察值。这意味着，在某个 经验样本 中，有些观察值可能被多次抽中，而有些观察值则可能一次都没有被抽中。例如，在上例中，经验样本 S 1 B S S_1^{BS} S1BS 中的观察值 0.1 被抽中了两次，而 0.3 则一次都没有被抽中。
1.2 标准差与标准误
简言之，统计量的标准差就称为 标准误。详情参见 「标准差与标准误 - 简书」，以及「标准差，标准与置信区间」。
2. 编写 bootstrap程序
Stata 的 bootstrap 命令非常方便，它不仅可以与估计命令（例如 OLS 回归和 logistic 回归）还与非估计命令（比如 summarize ）无缝衔接。bootstrap 可以自动执行自抽样过程，得到想要的统计量，并计算相关的统计指标（例如偏差和置信区间）。然而尽管这个命令非常方便，但在某些情况下想要获得的统计量却不能通过 bootstrap 实现。对于这些情况，您需要自行编写 bootstrap 程序。
本篇 Stata FAQ 将演示如何自行编写 bootstrap 程序。在第一个例子中，我们将 bootstrap 的结果与自行编写 bootstrap 程序的结果进行对比。通过比较, 可以发现自行编写 bootstrap 程序非常容易。接下来的一个示例将展示当 bootstrap 无法获得想要的统计量时,如何自行编写 bootstrap 程序。
为了使后续结果呈现更加美观，在执行 Stata 范例之前，我们可以先执行如下格式设定命令：
set cformat %4.3f //回归结果中系数的显示格式
set pformat %4.3f //回归结果中 p 值的显示格式
set sformat %4.2f //回归结果中 se值的显示格式
set showbaselevels off, permanently
set showemptycells off, permanently
set showomitted off, permanently
set fvlabel on, permanently
2.1 Stata 范例 1：OLS 回归的 RMSE 的标准误
在例一中，将通过使用 bootstrap 和自行编写的 bootstrap 程序来比较结果。我们使用 Stata 自带的 nlsw88.dta 数据中的年龄 (age)、种族 (race)、婚姻状况 (married) 和工作经验 (tenure) 对妇女工资 (wage) 进行回归，并通过 bootstrap 获得均方根误差 (rmse) 的标准误。对于 bootstrap 我们要求重复 100 次并且指定随机种子数，以确保结果开重现。
首先，进行初始回归。
. sysuse "nlsw88.dta", clear
. regress wage age race married tenure
结果如下：
Source | SS df MS Number of obs = 2,231
-------------+---------------------------------- F(4, 2226) = 26.36
Model | 3351.46097 4 837.865242 Prob > F = 0.0000
Residual | 70750.3667 2,226 31.7836328 R-squared = 0.0452
-------------+---------------------------------- Adj R-squared = 0.0435
Total | 74101.8276 2,230 33.2295191 Root MSE = 5.6377
------------------------------------------------------------------------------
wage | Coef. Std. Err. t P>|t| [95% Conf. Interval]
-------------+----------------------------------------------------------------
age | -0.107 0.039 -2.74 0.006 -0.184 -0.030
…… (Output omitted)
_cons | 12.842 1.608 7.99 0.000 9.689 15.996
------------------------------------------------------------------------------
此时的 RMSE = 5.6377。我们如何得到 RMSE 的标准误 (即 s.e.(RMSE)) 呢？显然，如果手头有足够的经费，我们可以从母体中另外抽取 100 份 样本 (Sample)，记为 S j ( j = 1 , 2 , ⋯ , 100 ) S_j (j = 1, 2, \cdots, 100) Sj(j=1,2,⋯,100)，经由每一个样本可以通过 OLS 获取对应的 RMSE，记为 R M S E j ( j = 1 , 2 , ⋯ , 100 ) RMSE_j (j = 1, 2, \cdots, 100) RMSEj(j=1,2,⋯,100)，进而得到 s . d . ( R M S E j ) s.d. (RMSE_j) s.d.(RMSEj)，即 R M S E j RMSE_j RMSEj 的标准差，而它就是我们所要求取的 R M S E RMSE RMSE 的标准误，即 s . d . ( R M S E j ) = s . e . ( R M S E ) s.d. (RMSE_j) = s.e.(RMSE) s.d.(RMSEj)=s.e.(RMSE)。
当然，现实中，我们通常无法获取这么多经费，而且也没有必要通过这种方式来估计 RMSE 的标准误。因为，如果我们手头的样本是从母体中随机抽取的，那么理论上可以证明 (Efron, 1993)，基于 Bootstrap 得到的 经验样本 (Empirical Sample) 可以很好地代替上述抽样样本 ( S j S_j Sj)。
或许各位读者已经明白我们要做的事情了：Bootstrap 其实就是一种抽样方法而已！
在本例中，nlsw88.dta 数据中涉及的 N = 2,231 个观察值记为原始样本 S 0 S_0 S0。执行 Bootstrap 的步骤如下：
- Step 1: 获取 Bootstrap 经验样本。从 S 0 S_0 S0 中有放回地抽取 N = 2,231 个观察值，得到经验样本 S 1 B S S_1^{BS} S1BS；
- Step 2: 使用经验样本进行估计。使用第 1 步中得到的经验样本 S 1 B S S_1^{BS} S1BS 进行 OLS 估计，得到 R M S E 1 B S RMSE_1^{BS} RMSE1BS；
- Step 3: 将第 1 步和第 2 步重复进行 K = 100 K = 100 K=100 次，得到 R M S E j B S ( j = 1 , 2 , ⋯ , 100 ) RMSE_j^{BS} (j = 1, 2, \cdots, 100) RMSEjBS(j=1,2,⋯,100)；
- Step 4: 计算 R M S E j B S RMSE_j^{BS} RMSEjBS 的标准差 s . d . ( R M S E j B S ) s.d. (RMSE_j^{BS}) s.d.(RMSEjBS)， 它就是 RMSE 这个统计量的抽样标准误。
下面，我们使用 bootstrap 命令来实现上述过程。这里，将种子值设定为 12345 (种子值可以是任何介于 0 和 0 and 2^31-1 (即 2,147,483,647) 之间的整数，详情参阅 help seed)，重复 100 次自抽样，得到 rmse 的标准误。
bootstrap rmse = e(rmse), reps(100) seed(12345): ///
regress wage age race married tenure
结果如下：
Bootstrap replications (100)
----+--- 1 ---+--- 2 ---+--- 3 ---+--- 4 ---+--- 5
.................................................. 50
.................................................. 100
Linear regression Number of obs = 2,231
Replications = 100
command: regress wage age race married tenure
rmse: e(rmse)
------------------------------------------------------------------------------
| Observed Bootstrap Normal-based
| Coef. Std. Err. z P>|z| [95% Conf. Interval]
-------------+----------------------------------------------------------------
rmse | 5.638 0.209 26.95 0.000 5.228 6.048
------------------------------------------------------------------------------
通过 estat bootstrap 得到上述 bootstrap 过程产生的所有信息。
estat bootstrap, all
Linear regression Number of obs = 2,231
Replications = 100
command: regress wage age race married tenure
rmse: e(rmse)
------------------------------------------------------------------------------
| Observed Bootstrap
| Coef. Bias Std. Err. [95% Conf. Interval]
-------------+----------------------------------------------------------------
rmse | 5.6376975 -.0200716 .2092082 5.227657 6.047738 (N)
| 5.248575 6.006687 (P)
| 5.251228 6.048207 (BC)
------------------------------------------------------------------------------
(N) normal confidence interval
(P) percentile confidence interval
(BC) bias-corrected confidence interval
编写 bootstrap 程序需要四步：
- 第一步，获得初始估计并将结果存储在 observe 矩阵中。此外，我们还必须记录分析中所使用的观测值个数。当我们计算 bootstrap 结果时，将需要使用这些数据。
- 第二步，我们编写了一个程序，我们将其称为 myboot，该程序通过重复抽样的方法对数据进行采样并返回需要的统计量。在此步骤中，我们首先利用
preserve命令保存数据，然后用bsample命令进行自抽样，其中bsample命令对原始数据进行重复抽样，这是 bootstrap 过程中必不可少的部分。对自抽样子样本进行回归，并使用return scalar输出需要的统计量。请注意，当使用program define myboot定义程序时，我们需要特别指定 rclass ，如果没有指定该项，将不会输出 bootstrap 统计量。 编写的 mybot 程序以restore结尾，该命令将使数据返回到 bootstrap 之前的原始状态。 - 第三步，使用前缀命令
simulate与 mybot 程序结合使用。此步骤指定与上一步中的相同的种子数和重复次数。 - 最后，我们使用
bstat命令来报告结果，其中存储在 observe 矩阵中的初始分析估计量和样本数分别放在 stat() 和 n() 选项中。
具体的 Stata 操作步骤如下：
sysuse "nlsw88.dta", clear
*Step 1
quietly regress wage age race married tenure
matrix observe = e(rmse)
*Step 2
*--------------------------------------begin----
capture program drop myboot
program define myboot, rclass
preserve
bsample
regress wage age race married tenure
return scalar rmse = e(rmse)
restore
end
*--------------------------------------over-----
*-
*-Note: 请选中上述 -begin- 和 -over- 之间的语句，
* 并按快捷键 `Ctrl + R`，以便把我们定义的 `myboot` 程序读入内存，
*Step 3
simulate rmse=r(rmse), reps(100) seed(12345): myboot
结果如下：
command: myboot
rmse: r(rmse)
Simulations (100)
----+--- 1 ---+--- 2 ---+--- 3 ---+--- 4 ---+--- 5
.................................................. 50
.................................................. 100
*Step 4
bstat, stat(observe) n(200)
Bootstrap results Number of obs = 200
Replications = 100
------------------------------------------------------------------------------
| Observed Bootstrap Normal-based
| Coef. Std. Err. z P>|z| [95% Conf. Interval]
-------------+----------------------------------------------------------------
rmse | 5.638 0.233 24.21 0.000 5.181 6.094
------------------------------------------------------------------------------
estat bootstrap, all
Bootstrap results Number of obs = 200
Replications = 100
------------------------------------------------------------------------------
| Observed Bootstrap
| Coef. Bias Std. Err. [95% Conf. Interval]
-------------+----------------------------------------------------------------
rmse | 5.6376975 -.0192764 .23284552 5.181329 6.094066 (N)
| 4.934809 5.973844 (P)
| 4.934809 5.973844 (BC)
------------------------------------------------------------------------------
(N) normal confidence interval
(P) percentile confidence interval
(BC) bias-corrected confidence interval
第四步的结果与上文使用 bootstrap 命令得到的结果一样。
2.2 Stata 范例 2：采用 Bootstrap 获取 VIF 的标准误和置信区间
在本例中，由于bootstrap 得到的统计量必须是可以直接通过 “analysis”得到的，否则bootstrap 将不能得到我们想要的统计量，因此就需要自行编写一个 bootstrap 手动操作程序来得到我们想要的统计量。如果想要查看内存中包含哪些统计量，则只需在 “analysis” 之后使用 ereturn list 或 return list 即可。ereturn list 和 return list 的区别在于 “analysis” 命令是否是估计命令。
假设我们想要通过 bootstrap 得到的统计量是方差膨胀因子（** VIF**），获得这个统计量首先需要运行regress,然后使用 estat vif 。然而在此情况下， 我们想要自抽样得到的统计量是通过后估计命令才能得到，而并不能直接从 regress 回归后获得，因此直接使用 bootstrap 并不能得到 VIF。因此，我们自行编写 bootstrap 程序来获得 VIF 的 bootstrap 估计。
sysuse "nlsw88.dta", clear
*-Step 1 进行初始回归计算 VIF
quietly regress wage age race married tenure
estat vif
初始回归后计算得到的 VIF 值结果如下：
Variable | VIF 1/VIF
-------------+----------------------
race | 1.04 0.958802
married | 1.04 0.962586
age | 1.01 0.990255
tenure | 1.01 0.991591
-------------+----------------------
Mean VIF | 1.03
通过 return list 查看 VIF 存储情况。
. return list
scalars:
r(vif_4) = 1.00848015716151
r(vif_3) = 1.009841055284585
r(vif_2) = 1.038868683195696
r(vif_1) = 1.042968550955925
macros:
r(name_4) : "tenure"
r(name_3) : "age"
r(name_2) : "married"
r(name_1) : "race"
新建一个 vif 矩阵存储计算得到的 VIF 值
matrix vif = ( r(vif_4), r(vif_3), r(vif_2), r(vif_1))
matrix list vif
结果如下：
vif[1,4]
c1 c2 c3 c4
r1 1.0084802 1.0098411 1.0388687 1.0429686
接下来，通过自行编写的 bootstrap 程序执行 bootstrap 。
*Step 2
capture program drop myboot2
program define myboot2, rclass
preserve
bsample
regress read female math write ses
estat vif
return scalar vif_4 = r(vif_4)
return scalar vif_3 = r(vif_3)
return scalar vif_2 = r(vif_2)
return scalar vif_1 = r(vif_1)
restore
end
*Step 3
simulate vif_4=r(vif_4) vif_3=r(vif_3) vif_2=r(vif_2) vif_1=r(vif_1), ///
reps(100) seed(12345): myboot2
command: myboot2
vif_4: r(vif_4)
vif_3: r(vif_3)
vif_2: r(vif_2)
vif_1: r(vif_1)
Simulations (100)
----+--- 1 ---+--- 2 ---+--- 3 ---+--- 4 ---+--- 5
.................................................. 50
.................................................. 100
bstat, stat(vif) n(200)
Bootstrap results Number of obs = 200
Replications = 100
------------------------------------------------------------------------------
| Observed Bootstrap Normal-based
| Coef. Std. Err. z P>|z| [95% Conf. Interval]
-------------+----------------------------------------------------------------
vif_4 | 1.00848 .0034994 288.19 0.000 1.001622 1.015339
vif_3 | 1.009841 .0038959 259.20 0.000 1.002205 1.017477
vif_2 | 1.038869 .0087988 118.07 0.000 1.021623 1.056114
vif_1 | 1.042969 .0093828 111.16 0.000 1.024579 1.061359
------------------------------------------------------------------------------
estat bootstrap, all
Bootstrap results Number of obs = 200
Replications = 100
------------------------------------------------------------------------------
| Observed Bootstrap
| Coef. Bias Std. Err. [95% Conf. Interval]
-------------+----------------------------------------------------------------
vif_4 | 1.0084802 .000774 .00349938 1.001622 1.015339 (N)
| 1.004008 1.016871 (P)
| 1.003899 1.015859 (BC)
vif_3 | 1.0098411 .0026392 .00389595 1.002205 1.017477 (N)
| 1.005817 1.02045 (P)
| 1.003658 1.01552 (BC)
vif_2 | 1.0388687 .0005919 .00879875 1.021623 1.056114 (N)
| 1.020785 1.058192 (P)
| 1.020785 1.058192 (BC)
vif_1 | 1.0429686 .0015362 .00938283 1.024579 1.061359 (N)
| 1.025008 1.064053 (P)
| 1.024641 1.063241 (BC)
------------------------------------------------------------------------------
(N) normal confidence interval
(P) percentile confidence interval
(BC) bias-corrected confidence interval
最终，我们通过自编 bootstrap 程序得到了 VIF 的 bootstrap 估计。
3. 主要参考文献
Books
- Efron, B., R. Tibshirani. An introduction to the bootstrap[M]. Chapman & Hall, 1997.
Bootstrap 简介
- Wehrens, R., H. Putter, L. Buydens (2000) Wehrens, R., H. Putter, L. Buydens, 2000, The bootstrap: A tutorial, Chemometrics and Intelligent Laboratory Systems, 54 (1): 35-52. 通俗易懂，比较浅显
- Davidson, R., E. Flachaire, 2008, The wild bootstrap, tamed at last, Journal of Econometrics, 146 (1): 162-169. [PDF]
- MacKinnon, J. G., 2006, Bootstrap methods in econometrics, Economic Record, 82: S2-S18.
- MacKinnon, J. G. 2013, Thirty years of heteroskedasticity-robust inference[C], Recent advances and future directions in causality, prediction, and specification analysis (Springer, 437-461.
- MacKinnon, J. G., M. D. Webb, 2017, Wild bootstrap inference for wildly different cluster sizes, Journal of Applied Econometrics, 32 (2): 233-254. [PDF]
- Davidson, R. (2012) Davidson, R., 2012, The bootstrap in econometrics, working paper, http://www.ncer.edu.au/events/documents/bootstrap_course.pdf. 这一篇实操性比较强
聚类标准误-Bootstrap
- Davidson, R., 2012, The bootstrap in econometrics, working paper, http://www.ncer.edu.au/events/documents/bootstrap_course.pdf.
关于我们
- Stata 连享会 由中山大学连玉君老师团队创办，定期分享实证分析经验。
- 欢迎赐稿： 欢迎赐稿至 StataChina@163。录用稿件达 三篇 以上，即可 免费 获得一期 Stata 现场培训资格。
- Stata连享会 (StataChina) 公众号关键词搜索/回复 功能已经上线。大家可以在公众号左下角点击键盘图标，输入简要关键词，以便快速呈现历史推文。常见关键词：
DID，RDD，PSM，面板，IV，合成控制法，plus，Profile, Bootstrap, 交乘项, 平方项, 工具, 软件, Sai2, gInk, Annotator, 直击面板数据, 连老师, 直播, 空间计量, 爬虫, 文本分析, 正则, Markdown幻灯片, marp, 盈余管理, MC……。
更多推荐
Stata：Bootstrap 简介
发布评论