R 中的自举实际上是如何工作的?

机器算法验证 r 回归 引导程序 回归系数
2022-02-10 22:09:15

我一直在研究 R 中的引导包,虽然我找到了许多关于如何使用它的很好的入门书,但我还没有找到任何可以准确描述“幕后”发生的事情的东西。例如,在这个例子中,该指南展示了如何使用标准回归系数作为引导回归的起点,但没有解释引导过程实际上是在做什么来推导引导回归系数。似乎正在发生某种迭代过程,但我似乎无法弄清楚到底发生了什么。

2个回答

引导程序有几种“风格”或形式(例如,非参数、参数、残差重采样等等)。示例中的引导程序称为非参数引导程序案例重采样(有关回归应用程序,请参见此处此处此处此处)。基本思想是将样本视为总体,并通过替换重复从中抽取新样本。所有原始观测值被抽入新样本的概率相同。然后您计算并存储感兴趣的统计数据,这可能是使用新抽取的样本的平均值、中位数或回归系数. 这是重复的n次。在每次迭代中,您的原始样本中的一些观察被多次提取,而一些观察可能根本不被提取。n迭代,你有n感兴趣的统计数据的存储引导估计(例如,如果n=1000并且感兴趣的统计数据是平均值,您有 1000 个自举估计平均值)。最后,汇总统计数据,例如平均值、中位数和标准差n计算引导估计。

自举通常用于:

  1. 置信区间的计算(和标准误差的估计)
  2. 点估计偏差的估计

几种基于 bootstrap 样本计算置信区间的方法(本文提供了解释和指导)。计算 95% 置信区间的一种非常简单的方法是计算引导样本的经验 2.5th 和 97.5th 百分位数(这个区间称为引导百分位数区间;参见下面的代码)。简单的百分位区间法在实践中很少使用,因为有更好的方法,例如偏差校正和加速引导 (BCa)。BCa 间隔针对自举分布中的偏差和偏度进行调整。

偏差被简单地估计为平均值之间的差异n存储的引导样本和原始估计。

让我们复制网站上的示例,但使用我们自己的循环结合我上面概述的想法(重复绘制并替换):

#-----------------------------------------------------------------------------
# Load packages
#-----------------------------------------------------------------------------

require(ggplot2)
require(pscl)
require(MASS)
require(boot)

#-----------------------------------------------------------------------------
# Load data
#-----------------------------------------------------------------------------

zinb <- read.csv("http://www.ats.ucla.edu/stat/data/fish.csv")
zinb <- within(zinb, {
  nofish <- factor(nofish)
  livebait <- factor(livebait)
  camper <- factor(camper)
})

#-----------------------------------------------------------------------------
# Calculate zero-inflated regression
#-----------------------------------------------------------------------------

m1 <- zeroinfl(count ~ child + camper | persons, data = zinb,
               dist = "negbin", EM = TRUE)

#-----------------------------------------------------------------------------
# Store the original regression coefficients
#-----------------------------------------------------------------------------

original.estimates <- as.vector(t(do.call(rbind, coef(summary(m1)))[, 1:2]))

#-----------------------------------------------------------------------------
# Set the number of replications
#-----------------------------------------------------------------------------

n.sim <- 2000

#-----------------------------------------------------------------------------
# Set up a matrix to store the results
#-----------------------------------------------------------------------------

store.matrix <- matrix(NA, nrow=n.sim, ncol=12)

#-----------------------------------------------------------------------------
# The loop
#-----------------------------------------------------------------------------

set.seed(123)

for(i in 1:n.sim) {

  #-----------------------------------------------------------------------------
  # Draw the observations WITH replacement
  #-----------------------------------------------------------------------------

  data.new <- zinb[sample(1:dim(zinb)[1], dim(zinb)[1], replace=TRUE),]

  #-----------------------------------------------------------------------------
  # Calculate the model with this "new" data
  #-----------------------------------------------------------------------------

  m <- zeroinfl(count ~ child + camper | persons,
                data = data.new, dist = "negbin",
                start = list(count = c(1.3711, -1.5152, 0.879),
                             zero = c(1.6028, -1.6663)))

  #-----------------------------------------------------------------------------
  # Store the results
  #-----------------------------------------------------------------------------

  store.matrix[i, ] <- as.vector(t(do.call(rbind, coef(summary(m)))[, 1:2]))

}


#-----------------------------------------------------------------------------
# Save the means, medians and SDs of the bootstrapped statistics
#-----------------------------------------------------------------------------

boot.means <- colMeans(store.matrix, na.rm=T)

boot.medians <- apply(store.matrix,2,median, na.rm=T)

boot.sds <- apply(store.matrix,2,sd, na.rm=T)

#-----------------------------------------------------------------------------
# The bootstrap bias is the difference between the mean bootstrap estimates
# and the original estimates
#-----------------------------------------------------------------------------

boot.bias <- colMeans(store.matrix, na.rm=T) - original.estimates

#-----------------------------------------------------------------------------
# Basic bootstrap CIs based on the empirical quantiles
#-----------------------------------------------------------------------------

conf.mat <- matrix(apply(store.matrix, 2 ,quantile, c(0.025, 0.975), na.rm=T),
ncol=2, byrow=TRUE)
colnames(conf.mat) <- c("95%-CI Lower", "95%-CI Upper")

这是我们的汇总表:

#-----------------------------------------------------------------------------
# Set up summary data frame
#-----------------------------------------------------------------------------

summary.frame <- data.frame(mean=boot.means, median=boot.medians,
sd=boot.sds, bias=boot.bias, "CI_lower"=conf.mat[,1], "CI_upper"=conf.mat[,2])

summary.frame

      mean  median       sd       bias CI_lower CI_upper
1   1.2998  1.3013  0.39674 -0.0712912  0.51960   2.0605
2   0.2527  0.2486  0.03208 -0.0034461  0.19898   0.3229
3  -1.5662 -1.5572  0.26220 -0.0509239 -2.12900  -1.0920
4   0.2005  0.1986  0.01949  0.0049019  0.16744   0.2418
5   0.9544  0.9252  0.48915  0.0753405  0.03493   1.9025
6   0.2702  0.2688  0.02043  0.0009583  0.23272   0.3137
7  -0.8997 -0.9082  0.22174  0.0856793 -1.30664  -0.4380
8   0.1789  0.1781  0.01667  0.0029513  0.14494   0.2140
9   2.0683  1.7719  1.59102  0.4654898  0.44150   8.0471
10  4.0209  0.8270 13.23434  3.1845710  0.58114  57.6417
11 -2.0969 -1.6717  1.56311 -0.4306844 -8.43440  -1.1156
12  3.8660  0.6435 13.27525  3.1870642  0.33631  57.6062

一些解释

  • bootstrap 估计的均值与原始估计之间的差异就是所谓的“偏差”boot
  • boot所谓“std.error”的输出是自举估计的标准偏差

将其与以下输出进行比较boot

#-----------------------------------------------------------------------------
# Compare with boot output and confidence intervals
#-----------------------------------------------------------------------------

set.seed(10)
res <- boot(zinb, f, R = 2000, parallel = "snow", ncpus = 4)

res

Bootstrap Statistics :
       original       bias    std. error
t1*   1.3710504 -0.076735010  0.39842905
t2*   0.2561136 -0.003127401  0.03172301
t3*  -1.5152609 -0.064110745  0.26554358
t4*   0.1955916  0.005819378  0.01933571
t5*   0.8790522  0.083866901  0.49476780
t6*   0.2692734  0.001475496  0.01957823
t7*  -0.9853566  0.083186595  0.22384444
t8*   0.1759504  0.002507872  0.01648298
t9*   1.6031354  0.482973831  1.58603356
t10*  0.8365225  3.240981223 13.86307093
t11* -1.6665917 -0.453059768  1.55143344
t12*  0.6793077  3.247826469 13.90167954

perc.cis <- matrix(NA, nrow=dim(res$t)[2], ncol=2)
    for( i in 1:dim(res$t)[2] ) {
  perc.cis[i,] <- boot.ci(res, conf=0.95, type="perc", index=i)$percent[4:5] 
}
colnames(perc.cis) <- c("95%-CI Lower", "95%-CI Upper")

perc.cis 

      95%-CI Lower 95%-CI Upper
 [1,]      0.52240       2.1035
 [2,]      0.19984       0.3220
 [3,]     -2.12820      -1.1012
 [4,]      0.16754       0.2430
 [5,]      0.04817       1.9084
 [6,]      0.23401       0.3124
 [7,]     -1.29964      -0.4314
 [8,]      0.14517       0.2149
 [9,]      0.29993       8.0463
[10,]      0.57248      56.6710
[11,]     -8.64798      -1.1088
[12,]      0.33048      56.6702

#-----------------------------------------------------------------------------
# Our summary table
#-----------------------------------------------------------------------------

summary.frame

      mean  median       sd       bias CI_lower CI_upper
1   1.2998  1.3013  0.39674 -0.0712912  0.51960   2.0605
2   0.2527  0.2486  0.03208 -0.0034461  0.19898   0.3229
3  -1.5662 -1.5572  0.26220 -0.0509239 -2.12900  -1.0920
4   0.2005  0.1986  0.01949  0.0049019  0.16744   0.2418
5   0.9544  0.9252  0.48915  0.0753405  0.03493   1.9025
6   0.2702  0.2688  0.02043  0.0009583  0.23272   0.3137
7  -0.8997 -0.9082  0.22174  0.0856793 -1.30664  -0.4380
8   0.1789  0.1781  0.01667  0.0029513  0.14494   0.2140
9   2.0683  1.7719  1.59102  0.4654898  0.44150   8.0471
10  4.0209  0.8270 13.23434  3.1845710  0.58114  57.6417
11 -2.0969 -1.6717  1.56311 -0.4306844 -8.43440  -1.1156
12  3.8660  0.6435 13.27525  3.1870642  0.33631  57.6062

将“bias”列和“std.error”与我们自己的汇总表的“sd”列进行比较。我们的 95% 置信区间与使用百分位数方法计算的置信区间非常相似boot.ci(但并非全部:查看索引为 9 的参数下限)。

您应该关注boot作为“统计”参数传递给的函数,并注意它是如何构造的。

f <- function(data, i) {
  require(pscl)
  m <- zeroinfl(count ~ child + camper | persons,
    data = data[i, ], dist = "negbin",
    start = list(count = c(1.3711, -1.5152, 0.879), zero = c(1.6028, -1.6663)))
  as.vector(t(do.call(rbind, coef(summary(m)))[, 1:2]))
}

“data”参数将接收整个数据帧,但“i”参数将接收由“boot”生成并取自 1:NROW(data) 的行索引样本。正如您从该代码中看到的那样,“i”然后用于创建一个新样本,该样本被传递给zeroinl然后仅返回其结果的选定部分。

假设“i”是 {1,2,3,3,3,6,7,7,10}。“[”函数将仅返回具有 3 个第 3 行副本和 2 个第 7 行副本的行。这将是单个zeroinl()计算的基础,然后系数将boot作为该过程复制的结果返回。此类复制的数量由“R”参数控制。

由于statistic在这种情况下仅返回回归系数,因此该boot函数将这些累积的系数作为“t”的值返回。可以通过其他引导包功能进行进一步的比较。