笔记说明
读《Discovering Statistics Using R》第六章 Correlation中的6.5.7节Bootstrapping correlations做的笔记。在书中bootstrap最初出现于5.8.4节。本笔记主要介绍bootstrap基本概念和使用bootstrap计算Kendall’s tau相关系数的置信区间。
示例数据
设我们想要验证一个理论:创造力更强的人可以讲出更厉害的故事。有这么一个比赛“the World's Biggest Liar competition”每年举办一次。作者收集了68个参赛者的比赛名次数据并让他们做了一份考察创造力的量表,满分60分。数据在这里:The Biggest Liar.dat
library(rio)
liarData <- import("data/The Biggest Liar.dat")
str(liarData)
## 'data.frame': 68 obs. of 3 variables:
## $ Creativity: int 53 36 31 43 30 41 32 54 47 50 ...
## $ Position : int 1 3 4 2 4 1 4 1 2 2 ...
## $ Novice : int 0 1 0 0 1 0 0 1 1 0 ...
Position即为比赛名次,Creativity即为创造力评分。
由于position变量为定序变量,而Pearson相关系数要求数据为定距变量,不适合使用Pearson相关系数。名次数据排序相同的情况有点多,比较适合使用Kendall’s tau相关系数。本系列笔记8中可以看到使用cor.test()
等方法无法计算Kendall’s tau相关系数的置信区间,我们使用bootstrap方法来计算Kendall’s tau相关系数的置信区间。
仍然是先做一个散点图看一下数据情况:
#散点图
library(ggplot2)
scatter <- ggplot(liarData, aes(Creativity, Position)) + geom_point()
Bootstrap
bootstrap是一种非参数的参数估计方法。也叫做自助法。
参数估计中的一个重要问题是我们不知道抽样分布(即统计量的分布)的形状。传统参数估计的解决方法:
- 如果样本数据符合正态分布,那么我们可以合理推断抽样分布符合正态分布。
- 中心极限定理告诉我们,在大样本下,抽样分布符合正态分布。
bootstrap法将样本数据作为一个总体,从这个总体中有放回地重复抽取很多更小的样本,样本量和原样本数据相同,新抽到的样本称为bootstrap样本,根据bootstrap样本的信息估计抽样分布的特征:在每个bootstrap样本中计算要估计的统计量,进而可以模拟一个抽样分布,用模拟出的抽样分布的标准差来作为待估计统计量的标准误的估计值。有了标准误的估计值就可以计算置信区间和进行假设检验了。
我们主要借助boot
包中的boot()
完成bootstrap. boot()
的一般用法为:
object <- boot(data = , statistic =, R = , ...)
生成一个bootstrap对象
-
data
指定使用的数据集 -
statistic
指定 计算统计量的函数。该函数的第一个参数应为原始数据,第二个参数应为用于选择数据行、确定bootstrap样本的indices。 -
R
即重复抽取bootstrap样本的次数。为使bootstrap估计的结果稳健,重复抽样数应该足够大,《Discovering Statistics Using R》作者一般将重复数设为2000。
用boot()
进行bootstrap法时的关键在于statistic=
参数对应的计算统计量的函数。这里我们要对Kendall’s tau相关系数这个统计量进行bootstrap,使用的函数为:
bootTau <-function(liarData,i) {
cor(liarData$Position[i], liarData$Creativity[i],
use = "complete.obs", method = "kendall")
}
这个函数的作用就是输入bootstrap样本(用i来确定),返回对应的样本Kendall’s tau相关系数。
定义好bootTau
后我们就可以使用boot()
来生成bootstrap对象:
library(boot)
set.seed(1234)
boot_kendall <- boot(liarData, bootTau, 2000)
boot_kendall
由于bootstrap过程涉及随机抽样,每次的计算结果都会不一样(重复次数足够大时结果会比较稳定),用set.seed()
确定随机种子数,以保证可重复性(即每次运行代码结果都是一样的)。
## ORDINARY NONPARAMETRIC BOOTSTRAP
##
##
## Call:
## boot(data = liarData, statistic = bootTau, R = 2000)
##
##
## Bootstrap Statistics :
## original bias std. error
## t1* -0.3002413 -0.00220919 0.09730457
从结果可以看到bootstrap法对Kendall’s tau相关系数的估计值为-0.300,标准误为0.097。
对于生成的bootstrap对象,可以用plot()
来查看bootstrap得到的抽样分布
plot(boot_kendall)
可以看出来bootstrap得到的Kendall’s tau相关系数的抽样分布比较符合正态分布。
对于生成的bootstrap对象,可以用
boot.ci()
来生成统计量的置信区间。boot.ci(bootobject, conf, type)
-
bootobject
即为boot()
生成的bootstrap对象 -
conf
指定置信区间长度,默认为0.95 -
type
指定置信区间计算方法,可能值为 "norm","basic", "stud", "perc", "bca","all",默认为"all",即输出全部5种方法的结果。具体这五种方法不在此展开,其中"bca"法对应adjusted bootstrap percentile (BCa) method,会根据偏差(bootstrap对象中显示的bias)对区间做简单调整。在大部分情况中是更可取的(此结论来自《R in action》第二版)
boot.ci(boot_kendall)
## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 2000 bootstrap replicates
##
## CALL :
## boot.ci(boot.out = boot_kendall)
##
## Intervals :
## Level Normal Basic
## 95% (-0.4887, -0.1073 ) (-0.4934, -0.1060 )
##
## Level Percentile BCa
## 95% (-0.4945, -0.1071 ) (-0.4887, -0.1005 )
## Calculations and Intervals on Original Scale
## Warning message:
## In boot.ci(boot_kendall) :
## bootstrap variances needed for studentized intervals
结果中给出了四种方法计算得到的置信区间( "stud"对应结果未能正常输出)。四个置信区间都没有跨0,表示Kendall’s tau相关系数有统计学意义。