《Discovering Statistics Using R》笔记9-用bootstrap法计算Kendall’s tau相关系数的置信区间

笔记说明

读《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相关系数有统计学意义。

最后编辑于
©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 203,271评论 5 476
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 85,275评论 2 380
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 150,151评论 0 336
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 54,550评论 1 273
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 63,553评论 5 365
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 48,559评论 1 281
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 37,924评论 3 395
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 36,580评论 0 257
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 40,826评论 1 297
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 35,578评论 2 320
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 37,661评论 1 329
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 33,363评论 4 318
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 38,940评论 3 307
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 29,926评论 0 19
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 31,156评论 1 259
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 42,872评论 2 349
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 42,391评论 2 342