布袋 • 桃花影 | 白泉山书院读写会11月参评篇目

开始今天的内容之前,需要大吼一声放假啦!感恩节假期正式开始,并且为防止节假日人流带来的更严重的新冠传染,感恩节之后本学期所有的课程也都将转为线上,也就是说,直到春季学期开学,我都会蜗居在自己的公寓里。来美国那么久,除了学校所在的城市,都没去过其他地方。本来今年暑假打算出游,可惜被疫情打乱,真希望疫情能早点过去。

回到我们的计量话题,今天要讲的内容是bootstrap自抽样。我还记得我在统计课上第一次学bootstrap的时候,全程都在神游,以至于后来做作业时以为bootstrap是什么特别高级难懂的技术问题,很长一段时间内碰到bootstrap就头大。但是后来重新复习的时候,发现bootstrap不过是一个非常简单的重复抽样过程,无非是使用R语言执行相关命令符的时候稍微有那么一点点复杂。如果仅仅局限于线性回归,那么Stata的bootstrap命令符使用就会简单很多。

在之前关于渐近性的推送中,我们有提到过,如果样本容量足够大的话,我们认为渐近性分布接近正态分布,也就满足了关于线性回归的正态分布的重要假设。但是,如果我们的样本容量没有那么大,没办法满足渐近性的要求,那又该怎么办呢?这时候bootstrap就是一个很有用的替代方法,bootstrap本质就是对现有的样本进行重复抽样得到新样本,以这种方法来模拟扩大样本容量。bootstrap的适用范围不仅限于线性回归估计,任何统计学里的估计量都可以通过bootstrap进行量化并且计算standard error标准误。并且我们完全不用在意估计量是否呈正态分布,通过bootstrap进行重复抽样,当重复样本足够多的时候,我们就可以模拟出估计量的分布状态。

bootstrap

bootstrap的基本理论

在前文介绍中我们已经提到,bootstrap本质就是对样本自身的重复抽样。我以一个非常简单的例子来进行解释,我们都知道扔硬币任意一面朝上的概率为0.5,假设我们有一枚硬币,我们随机扔10次,扔100次,扔1000次,随着次数的增多,那么其中一面朝上的概率会越来越接近0.5。假如我们没办法扔那么多次,比如只能扔10次,得出的概率值很大可能不等于0.5,如果我们想知道根据这10次计算出来的概率的误差有多少,应该怎么办呢?渐近性肯定是没办法用了,只有10次,样本容量太小。这时就我们就可以利用bootstrap,我们将这10次的结果记录下来,正面朝上为1,反面朝上为0。之后,我们对样本进行重复抽样,容量仍为10,因为是重复抽样,这10个结果中有的结果可能被抽出来两次甚至多次。第一次重复抽样得到一个新的样本,计算正面朝上的概率值;第二次重复抽样得到一个新的样本,计算正面朝上的概率值;第三次重复抽样得到一个新的样本,计算正面朝上的概率值;一直到第次,计算正面朝上的概率值。将通过重复抽样得到的估计值按照大小从低到高排列,假如我们想要获得95%置信区间,那么我们只需掐头去尾,去掉序列中前2.5%以及后2.5%的估计值,剩下的就是我们的置信区间。比如我们重复了1000次,那么95%的置信区间,从左往右数第25个以及从右往左数第25个就是分位数。我们也可以绘制分布直方图以及密度图,查看估计值的分布形态。下图是我的统计学老师的笔记中解释bootstrap流程的图:

布袋 • 桃花影 | 白泉山书院读写会11月参评篇目

回到线性回归分析,假如我们想要获得回归系数的bootstrap standard error,那么我们只需对样本进行重复抽样,每次对重复抽样获得的样本进行一次回归,就可计算出一个回归系数估计值,而后不断进行重复抽样,假如进行1000次,我们就会获得1000个回归系数的估计值,使用standard error的计算公式,我们就可以计算出bootstrap standard error:

因为任何一个统计学的估计量都可以通过bootstrap进行量化,那么线性回归分析中的任何估计值我们都可以使用bootstrap,比如R方。

R语言手动实现bootstrap

在R语言中实现bootstrap稍稍有点难度,一般使用boot工具包执行相关的命令符,之所以有点难,是因为要先自己编写新的符合条件的命令符。我们先手动模拟下刚才我们所说的扔硬币过程,然后再介绍如何使用boot工具包的命令符:

### 加载工具包library(tidyverse)### 加载boot工具包library(boot)### 设定种子set.seed(123)### 生成10个硬币样本coin <- sample(c(0, 1), size = 10, replace = TRUE)### 计算平均值 也即正面朝上的概率 结果为0.4mean(coin)
布袋 • 桃花影 | 白泉山书院读写会11月参评篇目
### 生成空白向量 填充bootstrap样本的平均值bootstrap <- c()### 重复抽样1000次 得到1000个估计值for(i in 1:1000){  b <- coin[sample(c(1:10), size = 10, replace = TRUE)]  bootstrap[i] <- mean(b)}### 计算bootstrap估计值的平均值mean(bootstrap)
布袋 • 桃花影 | 白泉山书院读写会11月参评篇目

通过bootstrap样本得到的估计值的平均值为0.3916,和原始样本的估计值0.4非常接近。下面根据公式计算standard error:

sqrt(sum((bootstrap - mean(bootstrap))^2)/999)
布袋 • 桃花影 | 白泉山书院读写会11月参评篇目

bootstrap计算得出的standard error约为0.152,也就是说,如果从总体中随机抽取一个容量为10的样本,我们认为该样本的平均值和其真实值之间的偏差约为0.152。这里我们知道硬币朝上的概率为0.5,样本的估计值为0.4,而bootstrap估计的样本估计值和0.5的偏差约为0.152。我们也可查看bootstrap重复抽样估计值的分布:

hist(bootstrap)
布袋 • 桃花影 | 白泉山书院读写会11月参评篇目

boot工具包

boot是比较常用的计算bootstrap的工具包,核心命令符为boot()boot.ci()。第一个命令符就是用来执行bootstrap自抽样,第二个命令符用来计算bootstrap置信区间。

命令符boot()有三个最基本的参数。第一个参数data是告诉R对哪一个样本进行重复抽样;第三个参数R告诉R重复抽样多少次。注意这里第二个参数statistic输入的值必须是一个新的命令符,这个命令符能够生成一个估计量的值或者多个估计量的数字向量,并且这个命令符必须包括两个参数,第一个参数是样本数据,第二个参数是筛选样本。如果没有这两个参数,命令符boot()就没办法对原始样本进行重复抽样并进行估计。下面我们以刚才扔硬币的例子为例:

### 设置新的命令符 至少包括两个参数### 其中一个表示样本数据### 第二个参数表示筛选样本### 新定义的函数就是对原始样本重新抽样然后计算平均值### 结果为单个数字coin_mean <- function(data, index){  b <- data[index]  return(mean(b))}### 只要我们输入相应的参数值 上述命令符是可以单独运行的coin_mean(data = coin, index = sample(1:10, 10, TRUE))
布袋 • 桃花影 | 白泉山书院读写会11月参评篇目
### 使用boot()命令符执行bootstrapboot1 <- boot(data = coin, statistic = coin_mean, R = 1000)boot1
布袋 • 桃花影 | 白泉山书院读写会11月参评篇目

上图中的结果,其中original是指原始样本估计的平均值,我们最开始就已经知道是0.4,bias是指通过bootstrap得到的1000个估计值的平均值和原始样本平均值之间的差值,等下我们可以进行验证。std.error就是我们求出的标准误约为0.156。注意因为抽样是随机的,所以结果和我们刚才手动计算出来的不一样。这里命令符生成的结果是以list形式保存的,我们可以查看结果中都包括什么:

attributes(boot1)
布袋 • 桃花影 | 白泉山书院读写会11月参评篇目

这里面的t0是指原始样本估计的平均值,t是指通过bootstrap得到的那1000个估计值,以矩阵形式保存,这里为的矩阵。刚才我们说bias是这1000个估计值的平均值与t0的差值,我们可以手动验证:

### 查看矩阵head(boot1$t)
布袋 • 桃花影 | 白泉山书院读写会11月参评篇目
### 手动计算biasmean(boot1$t[,1]) - boot1$t0
布袋 • 桃花影 | 白泉山书院读写会11月参评篇目

R语言自带的plot()函数可以直接查看bootstrap估计值的分布以及检验其正态分布的QQ plot:

plot(boot1)
布袋 • 桃花影 | 白泉山书院读写会11月参评篇目

QQ plot是用来检验数据是否呈现正态分布的,横坐标是正态分布的分位数,纵坐标是样本数据分布的分位数,如果样本数据为正态分布,其分位数和正态分布的分位数一一对应,那么所有的点连在一起我们就可以看到一个45度斜向上的线。

使用boot.ci()可以查看bootstrap估计的置信区间:

### 这里参数type设定的值为percentile### perc计算的置信区间就是前文所讲述的寻找置信区间的方法### 除此之外命令符还提供其他类型的置信区间### 感兴趣的可以看命令符的帮助页面boot.ci(boot1, type = "perc")
布袋 • 桃花影 | 白泉山书院读写会11月参评篇目

下面我们用boot命令符来对线性回归进行检验。我们使用R自带的数据mtcars进行演示。回归的因变量为mpg每英里耗油量,自变量为wt车重和disp汽车排量。假设我们要计算回归系数的bootstrap standard error,首先写新的命令符,这个命令符能够返回估计的系数值:

### 这里的函数有三个参数值 同样至少有两个参数### 其中一个是样本 另外一个是筛选样本lmcoef <- function(data, index, formula){  ### 对样本数据进行重复抽样  b <- data[index, ]  ### 线性回归估计  lm1 <- lm(formula, data = b)  ### 返回系数估计值  return(coef(lm1))}### 验证新的函数可以使用 返回数字向量lmcoef(mtcars, index = sample(1:nrow(mtcars), size = nrow(mtcars), replace = T), mpg ~ wt + disp)
布袋 • 桃花影 | 白泉山书院读写会11月参评篇目
### 使用boot()执行bootstrap### 注意刚才新写的命令符里面的参数formula这里要输入对应的参数值boot2 <- boot(data = mtcars, statistic = lmcoef, R = 1000, formula = mpg ~ wt + disp)boot2
布袋 • 桃花影 | 白泉山书院读写会11月参评篇目
head(boot2$t)
布袋 • 桃花影 | 白泉山书院读写会11月参评篇目

因为我们估计的系数值有三个,所以bootstrap给出的结果也有三个,和刚才命令符lmcoef()返回的向量里面的结果排列顺序一致。t1*就是(Intercept)t2*就是wtt3*就是disp。此时结果里面保存的t应该是的矩阵,因为我们重复了1000次,估计了3个变量,矩阵中每一列均为该变量的1000个估计值。我们依旧可以使用plot()查看bootstrap的分布图,因为这时有三个估计值,所以我们添加一个新的参数index,1代表t1*,2代表t2*,3代表t3*。如果不给参数index设置值,会默认返还t1*的结果。假如我们想查看disp的分布图:

plot(boot2, index = 3)
布袋 • 桃花影 | 白泉山书院读写会11月参评篇目

计算回归系数的置信区间,同理,需要输入index参数值,假如我们想要查看wt的置信区间:

boot.ci(boot2, type = "perc", index = 2)
布袋 • 桃花影 | 白泉山书院读写会11月参评篇目

假如我们想要计算R方的bootstrap标准误,我们只需写一个能够返回R方的命令符然后再用boot()执行即可:

lmrsquared <- function(data, index, formula){  b <- data[index, ]  lm1 <- lm(formula, data = b)  return(summary(lm1)$r.squared)}boot3 <- boot(data = mtcars, statistic = lmrsquared, R = 1000, formula = mpg ~ wt + disp)boot3
布袋 • 桃花影 | 白泉山书院读写会11月参评篇目
boot.ci(boot3, type = "perc")
布袋 • 桃花影 | 白泉山书院读写会11月参评篇目

在上述所有的bootstrap估计中,我都将重复抽样次数设置为1000,这是因为对于使用percentile计算置信区间的方法,1000次重复抽样才能够给出足够多的数据模拟结果。Stata默认的重复抽样次数为50,这对于假设的正态分布估计已经足够。最后,简要讲下如何在Stata中对线性回归进行bootstrap重复抽样估计:

### 使用Stata自带的auto.dta数据sysuse auto### 假设同方差进行回归regress mpg weight displacement
布袋 • 桃花影 | 白泉山书院读写会11月参评篇目
### bootstrap### 皆为默认50次### 可以直接在回归命令符后面添加参数 设置bootstrap standard errorregress mpg weight displacement, vce(bootstrap)### 也可以通过bootstrap命令符bootstrap: regress mpg weight displacement
布袋 • 桃花影 | 白泉山书院读写会11月参评篇目
### 重复抽样1000次bootstrap, rep(1000): regress mpg weight displacement
布袋 • 桃花影 | 白泉山书院读写会11月参评篇目

以上,就是在R语言中执行bootstrap自抽样的相关简介。

来源:爱尔兰粉公爵,本文观点不代表自营销立场,网址:https://www.zyxiao.com/p/100053

发表评论

电子邮件地址不会被公开。 必填项已用*标注

侵权联系
分享本页
返回顶部