Stata：Bootstrap 简介

作者: stata连享会 | 来源:发表于2019-09-15 00:02 被阅读0次

Stata：Bootstrap 简介
stata -- bootstrap
Stata语言编程 | 介绍篇
2019-02-17
BootStrap简介
Bootstrap简介
Bootstrap 简介
bootstrap简介
Bootstrap简介
一、bootstrap简介

作者：吴雄(湘潭大学)，童天天(中南财经政法大学) ，连玉君 (中山大学)

Stata 连享会：知乎 | 简书 | 码云

Source: The Bootstrap in Stata

连享会：内生性问题及估计方法专题

连享会-内生性专题现场班-2019.11.14-17

Stata连享会计量专题 || 精品课程 || 简书推文 || 公众号合集

点击查看完整推文列表

1. Bootstrap 简介

bootstrap 是一种崭新的增广样本统计方法，为解决小样本问题提供了很好的思路。它是非参数统计中一种重要的估计统计量方差进而进行区间估计的统计方法。对于回归模型：对于线性回归模型：

$y_t = X_t β+u_t, \\ E(u_t|X_t)=0,\ E(u_s u_t=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 个估计值获得统计量的分布，并计算置信区间。

有放回抽样

所谓 「有放回抽样」 (Samping with replacement) 指的是在逐个抽取个体时，每次被抽到的个体放回总体中后，再进行下次抽取的抽样方法。

举个例子，对于由 0.1 和 0.3 这两个数字构成的观测样本而言，记为 $S_0 = (0.1, 0.3)$ 。则采用有放回抽样 (Bootstrapping)，可以得到如下三种不同的经验样本： $S_1^{BS} = (0.1,0.1)$ , $S_2^{BS} = (0.1,0.3)$ , $S_3^{BS} = (0.3,0.3)$ 。

实际应用过程中，对于包含 N 个观察值的 观测样本 而言，Bootstrap 抽样得到的经验样本也会包含 N 个观察值。这意味着，在某个 经验样本 中，有些观察值可能被多次抽中，而有些观察值则可能一次都没有被抽中。例如，在上例中，经验样本 $S_1^{BS}$ 中的观察值 0.1 被抽中了两次，而 0.3 则一次都没有被抽中。

标准差与标准误

简言之，统计量的标准差就称为 标准误。详情参见「标准差与标准误 - 简书」，以及「标准差，标准与置信区间」。

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, \cdots, 100)$ ，经由每一个样本可以通过 OLS 获取对应的 RMSE，记为 $RMSE_j (j = 1, 2, \cdots, 100)$ ，进而得到 $s.d. (RMSE_j)$ ，即 $RMSE_j$ 的标准差，而它就是我们所要求取的 $RMSE$ 的标准误，即 $s.d. (RMSE_j) = s.e.(RMSE)$ 。

当然，现实中，我们通常无法获取这么多经费，而且也没有必要通过这种方式来估计 RMSE 的标准误。因为，如果我们手头的样本是从母体中随机抽取的，那么理论上可以证明 (Efron, 1993)，基于 Bootstrap 得到的 经验样本 (Empirical Sample) 可以很好地代替上述抽样样本 ( $S_j$ )。

或许各位读者已经明白我们要做的事情了：Bootstrap 其实就是一种抽样方法而已！

在本例中，nlsw88.dta 数据中涉及的 N = 2,231 个观察值记为原始样本 $S_0$ 。执行 Bootstrap 的步骤如下：

Step 1: 获取 Bootstrap 经验样本。从 $S_0$ 中有放回地抽取 N = 2,231 个观察值，得到经验样本 $S_1^{BS}$ ；
Step 2: 使用经验样本进行估计。使用第 1 步中得到的经验样本 $S_1^{BS}$ 进行 OLS 估计，得到 $RMSE_1^{BS}$ ；
Step 3: 将第 1 步和第 2 步重复进行 $K = 100$ 次，得到 $RMSE_j^{BS} (j = 1, 2, \cdots, 100)$ ；
Step 4: 计算 $RMSE_j^{BS}$ 的标准差 $s.d. (RMSE_j^{BS})$ ，它就是 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 估计。

主要参考文献

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. 这一篇实操性比较强