Stata :Bootstrap 简介
作者: 吴雄(湘潭⼤学),童天天(中南财经政法⼤学)
Source:
原⽂链接:
⽂章⽬录
1. Bootstrap 简介
bootstrap 是⼀种崭新的增⼴样本统计⽅法,为解决⼩样本问题提供了很好的思路。它是⾮参数统计中⼀种重要的估计统计量⽅差进⽽进⾏区间估计的统计⽅法。对于回归模型:对于线性回归模型:
可以通过多种⽅法来建⽴ 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 个估计值获得统计量的分布,并计算置信区间。readplay
1.1 有放回抽样
所谓 「有放回抽样」 (Samping with replacement) 指的是在逐个抽取个体时,每次被抽到的个体放回总体中后,再进⾏下次抽取的抽样⽅法。
举个例⼦,对于由 0.1 和 0.3 这两个数字构成的观测样本⽽⾔, 记为 。则采⽤有放回抽样 (Bootstrapping),可以得到如下三种不同的经验样本:, , 。
y =t X β+t u ,
t E (u ∣X )=t t 0, E (u u =s t 0) ∀ s = t
S =0(0.1,0.3)S =1BS (0.1,0.1)S =2BS (0.1,0.3)S =3BS (0.3,0.3)
实际应⽤过程中,对于包含 N 个观察值的 观测样本 ⽽⾔,Bootstrap 抽样得到的经验样本也会包含 N 个观察值。这意味着,在某个 经验
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 范例之前,我们可以先执⾏如下格式设定命令:
t cformat %4.3f //回归结果中系数的显⽰格式
ordinal
t pformat %4.3f //回归结果中 p 值的显⽰格式
t sformat %4.2f //回归结果中 值的显⽰格式
t showbalevels off, permanently
t showemptycells off, permanently
t showomitted off, permanently
t fvlabel on, permanently
2.1 Stata 范例 1:OLS 回归的 RMSE 的标准误
在例⼀中,将通过使⽤ bootstrap 和⾃⾏编写的 bootstrap 程序来⽐较结果。我们使⽤ Stata ⾃带的 nlsw88.dta 数据中的年龄 (age)、种族 (race)、婚姻状况 (married) 和⼯作经验 (tenure) 对妇⼥⼯资 (wage) 进⾏回归,并通过 bootstrap 获得均⽅根误差 (rm) 的标准误。对于 bootstrap 我们要求重复 100 次并且指定随机种⼦数,以确保结果开重现。
⾸先,进⾏初始回归。
. sysu "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)
ccq
_cons | 12.842 1.608 7.99 0.000 9.689 15.996
------------------------------------------------------------------------------
此时的 RMSE = 5.6377。我们如何得到 RMSE 的标准误 (即 s.e.(RMSE)) 呢?显然,如果⼿头有⾜够的经费,我们可以从母体中另外抽取 100 份 样本 (Sample),记为 ,经由每⼀个样本可以通过 OLS 获取对应的 RMSE,记为 ,进⽽得到 ,即 的标准差,⽽它就是我们所要求取的 的标准误,即
。
当然,现实中,我们通常⽆法获取这么多经费,⽽且也没有必要通过这种⽅式来估计 RMSE 的标准误。因为,如果我们⼿头的样本是从母体中随机抽取的,那么理论上可以证明 (Efron, 1993),基于 Bootstrap 得到的 经验样本 (Empirical Sample) 可以很好地代替上述抽样样本 ()。
或许各位读者已经明⽩我们要做的事情了:Bootstrap 其实就是⼀种抽样⽅法⽽已!
在本例中,nlsw88.dta 数据中涉及的 N = 2,231 个观察值记为原始样本 。执⾏ Bootstrap 的步骤如下:
Step 1: 获取 Bootstrap 经验样本。从 中有放回地抽取 N = 2,231 个观察值,得到经验样本 ;
Step 2: 使⽤经验样本进⾏估计。使⽤第 1 步中得到的经验样本 进⾏ OLS 估计,得到 ;
Step 3: 将第 1 步和第 2 步重复进⾏ 次,得到 ;
Step 4: 计算 的标准差 , 它就是 RMSE 这个统计量的抽样标准误。下⾯,我们使⽤ bootstrap 命令来实现上述过程。这⾥,将种⼦值设定为 12345 (种⼦值可以是任何介于 0 和 0 and 2^31-1 (即2,147,483,647) 之间的整数,详情参阅 help ed ),重复 100 次⾃抽样,得到 rm 的标准误。
bootstrap rm = e(rm), reps(100) ed(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
rm: e(rm)
------------------------------------------------------------------------------stacker
| Obrved Bootstrap Normal-bad
| Coef. Std. Err. z P>|z| [95% Conf. Interval]
-------------+----------------------------------------------------------------
李阳rm | 5.638 0.209 26.95 0.000 5.228 6.048
------------------------------------------------------------------------------通过 estat bootstrap 得到上述 bootstrap 过程产⽣的所有信息。
estat bootstrap, all
S (j =j 1,2,⋯,100)RMSE (j =j 1,2,⋯,100)s .d .(RMSE )j RMSE j RMSE s .d .(RMSE )=j s .e .(RMSE )S j S 0S 0S 1BS
S 1BS RMSE 1BS K =100RMSE (j =j BS 1,2,⋯,100)RMSE j BS s .d .(RMSE )j BS
Linear regression Number of obs = 2,231
Replications = 100
command: regress wage age race married tenure
rm: e(rm)
------------------------------------------------------------------------------
| Obrved Bootstrap
| Coef. Bias Std. Err. [95% Conf. Interval]
-------------+----------------------------------------------------------------
rm | 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 程序需要四步:
loft什么意思第⼀步,获得初始估计并将结果存储在 obrve 矩阵中。此外,我们还必须记录分析中所使⽤的观测值个数。当我们计算 bootstrap 结果时,将需要使⽤这些数据。
第⼆步,我们编写了⼀个程序,我们将其称为 myboot,该程序通过重复抽样的⽅法对数据进⾏采样并返回需要的统计量。在此步骤中,我们⾸先利⽤ prerve 命令保存数据,然后⽤ bsample 命令进⾏⾃抽样,其中 bsample 命令对原始数据进⾏重复抽样,这是bootstrap 过程中必不可少的部分。对
⾃抽样⼦样本进⾏回归,并使⽤ return scalar 输出需要的统计量。请注意,当使⽤ program define myboot 定义程序时,我们需要特别指定 rclass ,如果没有指定该项,将不会输出 bootstrap 统计量。 编写的 mybot 程序以restore 结尾,该命令将使数据返回到 bootstrap 之前的原始状态。
第三步,使⽤前缀命令 simulate 与 mybot 程序结合使⽤。此步骤指定与上⼀步中的相同的种⼦数和重复次数。
最后,我们使⽤ bstat 命令来报告结果,其中存储在 obrve 矩阵中的初始分析估计量和样本数分别放在 stat() 和 n() 选项中。
具体的 Stata 操作步骤如下:
sysu "nlsw88.dta", clear
*Step 1
quietly regress wage age race married tenure
matrix obrve = e(rm)
*Step 2
*--------------------------------------begin----
capture program drop myboot
program define myboot, rclass
prerve
bsample
regress wage age race married tenure
return scalar rm = e(rm)
restore
end
*--------------------------------------over-----
fml*-
*-Note: 请选中上述 -begin- 和 -over- 之间的语句,
* 并按快捷键 `Ctrl + R`,以便把我们定义的 `myboot` 程序读⼊内存,
*Step 3
simulate rm=r(rm), reps(100) ed(12345): myboot
结果如下:
command: myboot
rm: r(rm)
Simulations (100)
----+--- 1 ---+--- 2 ---+--- 3 ---+--- 4 ---+--- 5 (50)
(100)
*Step 4
bstat, stat(obrve) n(200)
Bootstrap results Number of obs = 200
Replications = 100
------------------------------------------------------------------------------
| Obrved Bootstrap Normal-bad
| Coef. Std. Err. z P>|z| [95% Conf. Interval]
成人英语三级考试-------------+----------------------------------------------------------------
rm | 5.638 0.233 24.21 0.000 5.181 6.094
------------------------------------------------------------------------------
estat bootstrap, all
Bootstrap results Number of obs = 200
Replications = 100
------------------------------------------------------------------------------
| Obrved Bootstrap
| Coef. Bias Std. Err. [95% Conf. Interval]
-------------+----------------------------------------------------------------
rm | 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 的标准误和置信区间thirteen
在本例中,由于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 估计。
sysu "nlsw88.dta", clear
*-Step 1 进⾏初始回归计算 VIF
quietly regress wage age race married tenure
estat vif
初始回归后计算得到的 VIF 值结果如下: