通过boot.ci函数获取多个统计信息的Bootstrap Confidence Interv

发布时间：2020-12-17 21:21:25 所属栏目：安全来源：网络整理

导读：我希望通过boot.ci函数获得多个统计信息的bootstrap置信区间.这是我的MWE. 我有两个统计数据,想要找到这两个统计信息的bootstrap置信区间.但是,boot.ci函数仅为第一个统计量(t1 *)提供引导置信区间,但不为第二个统计量(t2 *)提供自举置信区间. set.seed(1234

我希望通过boot.ci函数获得多个统计信息的bootstrap置信区间.这是我的MWE.

我有两个统计数据,想要找到这两个统计信息的bootstrap置信区间.但是,boot.ci函数仅为第一个统计量(t1 *)提供引导置信区间,但不为第二个统计量(t2 *)提供自举置信区间.

set.seed(12345)
df <- rnorm(n=10,mean = 0,sd = 1)


Boot.fun <- 
  function(data,idx) {
    data1 <- sample(data[idx],replace=TRUE)
    m1 <- mean(data1)
    sd1 <- sd(data1)
    out <- cbind(m1,sd1)
    return(out)
  }

Boot.fun(data = df)

library(boot)
boot.out <- boot(df,Boot.fun,R = 20)
boot.out

RDINARY NONPARAMETRIC BOOTSTRAP


Call:
  boot(data = df,statistic = Boot.fun,R = 20)


Bootstrap Statistics :
  original     bias    std. error
t1* -0.4815861  0.3190424   0.2309631
t2*  0.9189246 -0.1998455   0.2499412

boot.ci(boot.out=boot.out,conf = 0.95,type = c("norm","basic","perc","bca"))

BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
Based on 20 bootstrap replicates

CALL : 
  boot.ci(boot.out = boot.out,"bca"))

Intervals : 
  Level      Normal              Basic         
95%   (-1.2533,-0.3479 )   (-1.1547,-0.4790 )  

Level     Percentile            BCa          
95%   (-0.4842,0.1916 )   (-0.4842,-0.4629 )  
Calculations and Intervals on Original Scale
Warning : Basic Intervals used Extreme Quantiles
Some basic intervals may be unstable
Warning : Percentile Intervals used Extreme Quantiles
Some percentile intervals may be unstable
Warning : BCa Intervals used Extreme Quantiles
Some BCa intervals may be unstable
Warning messages:
  1: In norm.inter(t,(1 + c(conf,-conf))/2) :
  extreme order statistics used as endpoints
2: In norm.inter(t,alpha) : extreme order statistics used as endpoints
3: In norm.inter(t,adj.alpha) :
  extreme order statistics used as endpoints

解决方法

引导包(IMO)对于常规使用来说有点笨拙.简短的回答是你需要为boot.ci指定索引(默认值为1),例如boot.ci(boot.out,索引= 2).很长的答案是,一次获得所有引导程序统计信息的引导CI肯定会很方便！

获取指定结果槽的所有CI：

getCI <- function(x,w) {
   b1 <- boot.ci(x,index=w)
   ## extract info for all CI types
   tab <- t(sapply(b1[-(1:3)],function(x) tail(c(x),2)))
   ## combine with metadata: CI method,index
   tab <- cbind(w,rownames(tab),as.data.frame(tab))
   colnames(tab) <- c("index","method","lwr","upr")
   tab
}
## do it for both parameters
do.call(rbind,lapply(1:2,getCI,x=boot.out))

结果(可能不是你想要的,但很容易重塑)：

index  method        lwr        upr
normal       1  normal -1.2533079 -0.3479490
basic        1   basic -1.1547310 -0.4789996
percent      1 percent -0.4841726  0.1915588
bca          1     bca -0.4841726 -0.4628899
normal1      2  normal  0.6288945  1.6086459
basic1       2   basic  0.5727462  1.4789105
percent1     2 percent  0.3589388  1.2651031
bca1         2     bca  0.6819394  1.2651031

或者,如果您可以一次获得一个引导方法,我在Github上的扫帚包版本具有此功能(我已经提交了拉取请求)

## devtools::install_github("bbolker/broom")
library(broom)
tidy(boot.out,conf.int=TRUE,conf.method="perc")

（编辑：李大同）

【声明】本站内容均来自网络，其相关言论仅代表作者个人观点，不代表本站立场。若无意侵犯到您的权利，请及时与联系站长删除相关内容!