Search code examples
rstatistics-bootstrap

Bootstrap Confidence Intervals for more than one statistics through boot.ci function


I want to get bootstrap confidence intervals for more than one statistics through boot.ci function. Here is my MWE.

I've two statistics in out and want to find the bootstrap confidence intervals for these two statistics. However, boot.ci function is providing the bootstrap confidence intervals for only first statistic (t1*) but not for the second statistic (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, conf = 0.95, type = c("norm", "basic", 
                                                     "perc", "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

Solution

  • The boot package is (IMO) a little clunky for regular use. The short answer is that you need to specify index (default value is 1) to boot.ci, e.g. boot.ci(boot.out,index=2). The long answer is that it would certainly be convenient to get the bootstrap CIs for all of the bootstrap statistics at once!

    Get all CI for a specified result slot:

    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))
    

    Results (maybe not what you want, but easy to reshape):

             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
    

    Alternatively, if you can live with getting one bootstrap method at a time, my version of the broom package on Github has this capability (I've submitted a pull request)

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