R中的Bootstrap数据分析 _R开发

本文概述

一点符号(对不起！)
百分位数CI

Bootstrap是一种使用样本数据推断总体的方法。布拉德利·埃夫隆(Bradley Efron)于1979年在论文中首次介绍了它。Bootstrap依赖于采样, 并从样本数据中进行替换。该技术可用于估计任何统计信息的标准误差, 并为其获取置信区间(CI)。如果CI没有封闭的形式, 或者具有非常复杂的形式, 则Bootstrap尤其有用。
【R中的Bootstrap数据分析】假设我们有一个n个元素的样本：X = {x1, x2, … , xn}, 我们对CI的一些统计量T = t(X)感兴趣。 Bootstrap框架很简单。我们只重复R次以下方案：对于第i次重复, 从可用样本中替换n个元素的样本(其中一些元素将被多次提取)。将此新样本称为第i个Bootstrap样本Xi, 并计算所需的统计量Ti = t(Xi)。
结果, 我们将获得统计数据的R值：T1, T2, …, TR。我们称它们为T的Bootstrap实现或T的Bootstrap分布。基于此, 我们可以计算T的CI。有多种实现方法。采取百分位数似乎是最简单的方法。
让我们(再次使用)著名的虹膜数据集。看一下前几行：

head(iris)##Sepal.Length Sepal.Width Petal.Length Petal.Width Species ## 15.13.51.40.2setosa ## 24.93.01.40.2setosa ## 34.73.21.30.2setosa ## 44.63.11.50.2setosa ## 55.03.61.40.2setosa ## 65.43.91.70.4setosa

假设我们要找到中位数Sepal.Length, 中位数Sepal.Width和Spearman的秩相关系数的CI。我们将使用R的启动程序包和一个名为… boot的函数。要使用其功能, 我们必须创建一个函数, 该函数可以根据重采样数据计算统计信息。它至少应具有两个参数：一个数据集和一个向量, 其中包含一个数据集中元素的索引, 这些元素被挑选来创建引导程序样本。
如果我们希望一次计算多个统计量的CI, 则我们的函数必须将它们作为单个向量返回。
对于我们的示例, 它可能看起来像这样：

library(boot)foo < - function(data, indices){ dt< -data[indices, ] c( cor(dt[, 1], dt[, 2], method='s'), median(dt[, 1]), median(dt[, 2]) ) }

foo从数据中选择所需的元素(将数字存储在索引中)并计算前两列的相关系数(方法=’ s’ 用于选择Spearman的秩系数, method =’ p’ 将产生Pearson系数)及其中位数。
我们还可以添加一些额外的参数, 例如。让用户选择od类型的相关系数：

foo < - function(data, indices, cor.type){ dt< -data[indices, ] c( cor(dt[, 1], dt[, 2], method=cor.type), median(dt[, 1]), median(dt[, 2]) ) }

为了使本教程更加通用, 我将使用foo的最新版本。
现在, 我们可以使用启动功能了。我们必须为其提供名称od数据集, 刚刚创建的函数, 重复次数(R)以及函数的其他任何参数(例如cor.type)。下面, 我使用set.seed来重现此示例。

set.seed(12345) myBootstrap < - boot(iris, foo, R=1000, cor.type='s')

引导功能返回一个名为… (是的, 你是对的！)引导的类的对象。它有两个有趣的元素。 $ t包含由引导程序(T的引导实现)生成的统计值的R值：

head(myBootstrap$t)##[, 1] [, 2] [, 3] ## [1, ] -0.26405188 5.703 ## [2, ] -0.12973299 5.803 ## [3, ] -0.07972066 5.753 ## [4, ] -0.16122705 6.003 ## [5, ] -0.20664808 6.003 ## [6, ] -0.12221170 5.803

$ t0包含原始, 完整数据集中的统计值：

myBootstrap$t0## [1] -0.16677775.80000003.0000000

将启动对象打印到控制台可提供一些进一步的信息：

myBootstrap## ## ORDINARY NONPARAMETRIC BOOTSTRAP ## ## ## Call: ## boot(data = http://www.srcmini.com/iris, statistic = foo, R = 1000, cor.type ="s") ## ## ## Bootstrap Statistics : ##originalbiasstd. error ## t1* -0.16677770.0025463910.07573983 ## t2*5.8000000 -0.0133500000.10295571 ## t3*3.00000000.0079000000.02726414

original与$ t0相同。偏差是引导实现的平均值(来自$ t的引导实现)之间的差, 称为T的引导估计和原始数据集中的值(来自$ t0的引导估计)。

colMeans(myBootstrap$t)-myBootstrap$t0## [1]0.002546391 -0.0133500000.007900000

标准误差是引导程序估计的标准误差, 等于引导程序实现的标准偏差。

apply(myBootstrap$t, 2, sd)## [1] 0.07573983 0.10295571 0.02726414

在开始使用CI之前, 始终值得一看引导实现的分布。我们可以使用plot函数, 通过索引告诉我们希望查看foo中计算的统计量。此处index = 1是萼片长度和宽度之间的Spearman相关系数, index = 2是中位数od sepal长度, index = 3是中位数od sepal宽度。

plot(myBootstrap, index=1)

文章图片
Bootstrap相关系数的分布似乎很正常。让我们为其找到CI。我们可以使用boot.ci。默认值为95％的配置项, 但可以使用conf参数进行更改。

boot.ci(myBootstrap, index=1)## Warning in boot.ci(myBootstrap, index = 1): bootstrap variances needed for ## studentized intervals## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS ## Based on 1000 bootstrap replicates ## ## CALL : ## boot.ci(boot.out = myBootstrap, index = 1) ## ## Intervals : ## LevelNormalBasic ## 95%(-0.3178, -0.0209 )(-0.3212, -0.0329 ) ## ## LevelPercentileBCa ## 95%(-0.3007, -0.0124 )(-0.3005, -0.0075 ) ## Calculations and Intervals on Original Scale

boot.ci提供5种类型的引导CI。其中之一, 学生间隔时间是唯一的。它需要估计Bootstrap方差。我们没有提供它, 因此R会显示警告：学生化间隔所需的引导程序差异。方差估计可以通过二级引导程序获得, 也可以通过折刀技术轻松获得。这在某种程度上超出了本教程的范围, 因此让我们集中讨论其余四种类型的引导CI。
如果我们不想全部看到它们, 可以在type参数中选择相关的参数。可能的值为norm, basic, stud, perc, bca或这些的向量。

boot.ci(myBootstrap, index=1, type=c('basic', 'perc'))## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS ## Based on 1000 bootstrap replicates ## ## CALL : ## boot.ci(boot.out = myBootstrap, type = c("basic", "perc"), index = 1) ## ## Intervals : ## LevelBasicPercentile ## 95%(-0.3212, -0.0329 )(-0.3007, -0.0124 ) ## Calculations and Intervals on Original Scale

boot.ci函数创建类的对象… (是的, 你猜对了！)bootci。就像在类型实参中使用的CI类型一样, 调用其元素。 $ norm是包含置信度水平和CI边界的3元素向量。

boot.ci(myBootstrap, index=1, type='norm')$norm##conf ## [1, ] 0.95 -0.3177714 -0.02087672

$ basic, $ stud, $ perc和$ bca是由5个元素组成的向量, 还包含用于计算CI的百分位数(我们稍后会再介绍)：

boot.ci(myBootstrap, index=1, type='basic')$basic##conf ## [1, ] 0.95 975.98 25.03 -0.3211981 -0.03285178

一点符号(对不起！) 要了解不同类型的CI是什么, 我们需要引入一些符号。所以让：

t?是引导估计(引导实现的平均值),
t0是原始数据集中的统计值,
这是引导估计的标准误差,
b是Bootstrap估计的偏差b =t?-t0
α为置信度, 通常为α= 0.95,
zα是标准正态分布的$ 1- \ frac \ alpha 2 $-分位数,
θα是引导实现的分布的α百分数。

百分位数CI 使用上述符号, 百分位数CI为：
(θ(1-α)/ 2, θ1-(1-α)/ 2)
所以这只需要相关的百分位数。而已。
典型的Wald CI为：
t0±zα?se?
但是在Bootstrap情况下, 我们应该对其进行校正以消除偏差。这样就变成：
$$ t_0-b \ pm z_ alpha \ cdot se ^ \ star \\ 2t_0-t ^ \ star \ pm alpha_ cdot se ^ \ star $$
通常不建议使用百分位数CI, 因为它在处理怪异的分布时表现不佳。基本配置项(也称为关键配置项或经验配置项)对此更为健壮。其背后的原理是计算每个引导复制和t0之间的差异, 并使用其分布的百分位数。可以找到完整的详细信息, 例如在L. Wasserman的所有统计资料中
基本CI的最终公式为：
(2t0-θ1-(1-α)/ 2, 2t0-θ(1-α)/ 2)
BCα来自偏差校正, 加速。它的公式不是很复杂, 但是有点不直观, 因此我将跳过它。如果你对详细信息感兴趣, 请参阅Thomas J. DiCiccio和Bradley Efron的文章。
方法名称中提到的加速要求使用引导实现的特定百分比。有时可能会出现一些极端的百分位数, 可能是异常值。在这种情况下, BCα可能不稳定。
让我们来看看BCαCI的花瓣宽度中位数。在原始数据集中, 该中值正好等于3。

boot.ci(myBootstrap, index=3)## Warning in boot.ci(myBootstrap, index = 3): bootstrap variances needed for ## studentized intervals## Warning in norm.inter(t, adj.alpha): extreme order statistics used as ## endpoints## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS ## Based on 1000 bootstrap replicates ## ## CALL : ## boot.ci(boot.out = myBootstrap, index = 3) ## ## Intervals : ## LevelNormalBasic ## 95%( 2.939, 3.046 )( 2.900, 3.000 ) ## ## LevelPercentileBCa ## 95%( 3.0, 3.1 )( 2.9, 2.9 ) ## Calculations and Intervals on Original Scale ## Warning : BCa Intervals used Extreme Quantiles ## Some BCa intervals may be unstable

我们得到BCαCI(2.9, 2.9)。这很奇怪, 但是幸运的是R告诫我们极端顺序统计用作端点。让我们看看这里到底发生了什么：

plot(myBootstrap, index=3)

文章图片
引导实现的分布并不常见。其中绝大多数(超过90％)为3秒。

table(myBootstrap$t[, 3])## ##2.9 2.953 3.053.1 3.153.2 ##11908246312

在这种情况下, BCα方法中的” 加速度” 很容易” 跳” 至极值。在此, 构成BCαCI基础的百分位数CI为(3, 3.1)。现在, 当我们将其” 向左” 扩展时, 必须使用0.002％。 “ 向右” 的扩展也达到极限：0.997％。
有时, 我们需要重新创建引导复制。如果我们可以使用R, 这不是问题。 set.seed函数解决了这个问题。在其他软件中重新创建引导程序复制将更加复杂。而且, 对于大的R, R中的重新计算也不是一种选择(例如, 由于时间不足)。
我们可以解决这个问题, 保存形成每个引导程序样本的原始数据集元素的索引。这就是boot.array函数(带有index = T参数)的作用。

tableOfIndices< -boot.array(myBootstrap, indices=T)

每行是一个引导程序样本。例如, 我们的第一个示例包含以下元素：

tableOfIndices[1, ]##[1] 10912 143654128 1056223 1023734555338311 ##[18] 14631 12285 141 118 116 11778 10013984959882456 ##[35] 13645967 126 118 104 1017013121921 14973 13367 ##[52]86688 131 105 121 145 12183706871 1117673 122 116 ##[69]85 144 114 13989 124642781785839175037 11776 ##[86]97 114641793 14165 12480 13754375770 1281055 ## [103]13535945 1161477 118 108 138507849 10449185 ## [120]28438274643355325962 11281196301430 ## [137]3885668597 107418763531 1332769

设置index = F(默认值), 我们将得到以下问题的答案：” 每个引导样本中原始数据集的每个元素出现了多少次？” 。例如, 在第一个样本中：数据集的第一个元素出现一次, 根本没有第二个元素出现, 第三个元素出现一次, 第四个元素出现两次, 依此类推。

tableOfAppearances< -boot.array(myBootstrap) tableOfAppearances[1, ]##[1] 1 0 1 2 0 1 0 1 0 1 2 2 3 2 0 0 2 1 1 0 1 0 1 1 0 0 2 2 0 2 2 1 1 1 1 ##[36] 0 3 2 1 0 1 0 1 0 1 0 0 0 3 2 0 0 2 1 3 1 1 1 4 0 0 2 0 3 2 1 2 1 1 3 ##[71] 1 0 2 1 0 3 1 3 0 1 1 1 1 0 5 1 0 2 1 0 0 0 1 0 0 1 2 1 0 1 1 1 0 2 2 ## [106] 0 1 1 1 0 1 1 0 2 0 3 2 3 0 0 2 2 0 2 0 1 0 1 0 0 1 0 2 0 0 1 1 1 1 0 ## [141] 2 0 1 1 1 1 0 0 1 0

这样的表使我们可以在R外部重新创建引导实现。或者, 当我们不想使用set.seed并再次执行所有计算时, 可以在R自身中重新创建引导实现。

onceAgain< -apply(tableOfIndices, 1, foo, data=http://www.srcmini.com/iris, cor.type='s')

让我们检查结果是否相同：

head(t(onceAgain))##[, 1] [, 2] [, 3] ## [1, ] -0.26405188 5.703 ## [2, ] -0.12973299 5.803 ## [3, ] -0.07972066 5.753 ## [4, ] -0.16122705 6.003 ## [5, ] -0.20664808 6.003 ## [6, ] -0.12221170 5.803head(myBootstrap$t)##[, 1] [, 2] [, 3] ## [1, ] -0.26405188 5.703 ## [2, ] -0.12973299 5.803 ## [3, ] -0.07972066 5.753 ## [4, ] -0.16122705 6.003 ## [5, ] -0.20664808 6.003 ## [6, ] -0.12221170 5.803all(t(onceAgain)==myBootstrap$t)## [1] TRUE

对, 他们是！
如果你想了解有关R中机器学习的更多信息, 请参加srcmini的机器学习工具箱课程。