`boot` tidier doesn't work for some `conf.method` types

[IndrajeetPatil]

Tidier for boot object works inconsistently.

Here is a custom function to get confidence interval for Spearman's rho using boot package.

# setup
set.seed(123)
library(tidyverse)

# custom function
cor_ci <- function(data,
                   x,
                   y,
                   method = "spearman",
                   exact = FALSE,
                   continuity = TRUE,
                   alternative = "two.sided",
                   nboot = 100,
                   conf.level = 0.95,
                   conf.type = "norm",
                   ...) {
  # creating a dataframe from entered data
  data <- dplyr::select(
    .data = data,
    x = !!rlang::enquo(x),
    y = !!rlang::enquo(y)
  ) %>%
    dplyr::filter(.data = ., !is.na(x), !is.na(y)) %>%
    tibble::as_tibble(x = .)

  # function to obtain 95% CI for xi
  corci <- function(formula,
                      x,
                      y,
                      method = method,
                      exact = exact,
                      continuity = continuity,
                      alternative = alternative,
                      indices) {
    # allows boot to select sample
    d <- data[indices, ]
    # allows boot to select sample
    boot_df <-
      broom::tidy(
        stats::cor.test(
          formula = stats::as.formula(~ x + y),
          data = d,
          method = method,
          exact = exact,
          continuity = continuity,
          alternative = alternative,
          na.action = na.omit
        )
      )

    # return the value of interest: effect size
    return(boot_df$estimate)
  }

  # save the bootstrapped results to an object
  bootobj <- boot::boot(
    data = data,
    statistic = corci,
    R = nboot,
    x = x,
    y = y,
    method = method,
    exact = exact,
    continuity = continuity,
    alternative = alternative,
    parallel = "multicore",
    ...
  )

  # tidy dataframe with confidence intervals
  df <- broom::tidy(
    bootobj,
    conf.int = TRUE,
    conf.level = conf.level,
    conf.method = conf.type
  )

  # get 95% CI from the bootstrapped object
  bootci <- boot::boot.ci(
    boot.out = bootobj,
    conf = conf.level,
    type = conf.type
  )

  # return both of these objects
  list(df, bootci)
}

• basic: works

cor_ci(mtcars, wt, mpg, conf.type = "basic")

#> [[1]]
#> # A tibble: 1 x 5
#>   statistic    bias std.error conf.low conf.high
#>       <dbl>   <dbl>     <dbl>    <dbl>     <dbl>
#> 1    -0.886 0.00665    0.0542    -1.03    -0.818
#> 
#> [[2]]
#> BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
#> Based on 100 bootstrap replicates
#> 
#> CALL : 
#> boot::boot.ci(boot.out = bootobj, conf = conf.level, type = conf.type)
#> 
#> Intervals : 
#> Level      Basic         
#> 95%   (-1.0316, -0.8181 )  
#> Calculations and Intervals on Original Scale
#> Some basic intervals may be unstable

• bca: works

cor_ci(mtcars, wt, mpg, conf.type = "bca")

#> [[1]]
#> # A tibble: 1 x 5
#>   statistic    bias std.error conf.low conf.high
#>       <dbl>   <dbl>     <dbl>    <dbl>     <dbl>
#> 1    -0.886 0.00909    0.0420   -0.954    -0.783
#> 
#> [[2]]
#> BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
#> Based on 100 bootstrap replicates
#> 
#> CALL : 
#> boot::boot.ci(boot.out = bootobj, conf = conf.level, type = conf.type)
#> 
#> Intervals : 
#> Level       BCa          
#> 95%   (-0.9543, -0.7826 )  
#> Calculations and Intervals on Original Scale
#> Some BCa intervals may be unstable

• percentile: works

cor_ci(mtcars, wt, mpg, conf.type = "perc")

#> [[1]]
#> # A tibble: 1 x 5
#>   statistic   bias std.error conf.low conf.high
#>       <dbl>  <dbl>     <dbl>    <dbl>     <dbl>
#> 1    -0.886 0.0130    0.0655   -0.955    -0.689
#> 
#> [[2]]
#> BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
#> Based on 100 bootstrap replicates
#> 
#> CALL : 
#> boot::boot.ci(boot.out = bootobj, conf = conf.level, type = conf.type)
#> 
#> Intervals : 
#> Level     Percentile     
#> 95%   (-0.9550, -0.6893 )  
#> Calculations and Intervals on Original Scale
#> Some percentile intervals may be unstable

• normal : doesn't work properly; confidence intervals are available but tidier outputs NAs instead

cor_ci(mtcars, wt, mpg, conf.type = "norm")

#> [[1]]
#> # A tibble: 1 x 5
#>   statistic   bias std.error conf.low conf.high
#>       <dbl>  <dbl>     <dbl>    <dbl>     <dbl>
#> 1    -0.886 0.0104    0.0531       NA        NA
#> 
#> [[2]]
#> BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
#> Based on 100 bootstrap replicates
#> 
#> CALL : 
#> boot::boot.ci(boot.out = bootobj, conf = conf.level, type = conf.type)
#> 
#> Intervals : 
#> Level      Normal        
#> 95%   (-1.0009, -0.7928 )  
#> Calculations and Intervals on Original Scale

• Student's t: doesn't work

cor_ci(mtcars, wt, mpg, conf.type = "stud")

#> Warning in FUN(X[[i]], ...): bootstrap variances needed for studentized
#> intervals
#> Error in dimnames(x) <- dn: length of 'dimnames' [2] not equal to array extent

^{Created on 2019-01-12 by the reprex package (v0.2.1)}

Session info

devtools::session_info()
#> - Session info ----------------------------------------------------------
#>  setting  value                                             
#>  version  R Under development (unstable) (2018-11-30 r75724)
#>  os       Windows 10 x64                                    
#>  system   x86_64, mingw32                                   
#>  ui       RTerm                                             
#>  language (EN)                                              
#>  collate  English_United States.1252                        
#>  ctype    English_United States.1252                        
#>  tz       Asia/Calcutta                                     
#>  date     2019-01-12                                        
#> 
#> - Packages --------------------------------------------------------------
#>  package     * version     date       lib
#>  assertthat    0.2.0       2017-04-11 [1]
#>  backports     1.1.3       2018-12-14 [1]
#>  boot          1.3-20      2017-08-06 [2]
#>  broom         0.5.1.9000  2019-01-11 [1]
#>  callr         3.1.1       2018-12-21 [1]
#>  cellranger    1.1.0       2016-07-27 [1]
#>  cli           1.0.1.9000  2018-10-30 [1]
#>  colorspace    1.3-2       2016-12-14 [1]
#>  crayon        1.3.4       2017-09-16 [1]
#>  desc          1.2.0       2018-10-30 [1]
#>  devtools      2.0.1       2018-10-26 [1]
#>  digest        0.6.18      2018-10-10 [1]
#>  dplyr       * 0.8.0       2019-01-06 [1]
#>  evaluate      0.12        2018-10-09 [1]
#>  fansi         0.4.0       2018-11-05 [1]
#>  forcats     * 0.3.0       2018-02-19 [1]
#>  fs            1.2.6       2018-08-23 [1]
#>  generics      0.0.2       2018-11-29 [1]
#>  ggplot2     * 3.1.0.9000  2018-12-15 [1]
#>  glue          1.3.0       2018-07-17 [1]
#>  gtable        0.2.0       2016-02-26 [1]
#>  haven         2.0.0       2018-11-22 [1]
#>  highr         0.7         2018-06-09 [1]
#>  hms           0.4.2       2018-03-10 [1]
#>  htmltools     0.3.6       2017-04-28 [1]
#>  httr          1.4.0       2018-12-11 [1]
#>  jsonlite      1.6         2018-12-07 [1]
#>  knitr         1.21        2018-12-10 [1]
#>  lazyeval      0.2.1       2017-10-29 [1]
#>  lubridate     1.7.4       2018-04-11 [1]
#>  magrittr      1.5         2014-11-22 [1]
#>  memoise       1.1.0       2017-04-21 [1]
#>  modelr        0.1.2       2018-05-11 [1]
#>  munsell       0.5.0       2018-06-12 [1]
#>  pillar        1.3.1       2018-12-15 [1]
#>  pkgbuild      1.0.2       2018-10-16 [1]
#>  pkgconfig     2.0.2       2018-08-16 [1]
#>  pkgload       1.0.2       2018-10-29 [1]
#>  prettyunits   1.0.2       2015-07-13 [1]
#>  processx      3.2.1       2018-12-05 [1]
#>  ps            1.3.0       2018-12-21 [1]
#>  purrr       * 0.2.99.9000 2019-01-04 [1]
#>  R6            2.3.0       2018-10-04 [1]
#>  Rcpp          1.0.0       2018-11-07 [1]
#>  readr       * 1.3.1       2018-12-21 [1]
#>  readxl        1.2.0       2018-12-19 [1]
#>  remotes       2.0.2       2018-10-30 [1]
#>  rlang         0.3.1       2019-01-08 [1]
#>  rmarkdown     1.11        2018-12-08 [1]
#>  rprojroot     1.3-2       2018-01-03 [1]
#>  rvest         0.3.2       2016-06-17 [1]
#>  scales        1.0.0       2018-08-09 [1]
#>  sessioninfo   1.1.1       2018-11-05 [1]
#>  stringi       1.2.4       2018-07-20 [1]
#>  stringr     * 1.3.1       2018-05-10 [1]
#>  testthat      2.0.1       2018-10-13 [1]
#>  tibble      * 2.0.0       2019-01-04 [1]
#>  tidyr       * 0.8.2       2018-10-28 [1]
#>  tidyselect    0.2.5       2018-10-11 [1]
#>  tidyverse   * 1.2.1       2017-11-14 [1]
#>  usethis       1.4.0.9000  2018-12-12 [1]
#>  utf8          1.1.4       2018-05-24 [1]
#>  withr         2.1.2       2018-03-15 [1]
#>  xfun          0.4         2018-10-23 [1]
#>  xml2          1.2.0       2018-01-24 [1]
#>  yaml          2.2.0       2018-07-25 [1]
#>  source                            
#>  CRAN (R 3.5.1)                    
#>  CRAN (R 3.6.0)                    
#>  CRAN (R 3.6.0)                    
#>  Github (tidymodels/broom@e1465ea) 
#>  CRAN (R 3.6.0)                    
#>  CRAN (R 3.5.1)                    
#>  Github (r-lib/cli@56538e3)        
#>  CRAN (R 3.5.1)                    
#>  CRAN (R 3.5.1)                    
#>  Github (r-lib/desc@7c12d36)       
#>  CRAN (R 3.5.1)                    
#>  CRAN (R 3.5.1)                    
#>  Github (tidyverse/dplyr@59d4600)  
#>  CRAN (R 3.5.1)                    
#>  Github (brodieG/fansi@ab11e9c)    
#>  CRAN (R 3.5.1)                    
#>  CRAN (R 3.5.1)                    
#>  CRAN (R 3.6.0)                    
#>  Github (tidyverse/ggplot2@9292832)
#>  CRAN (R 3.5.1)                    
#>  CRAN (R 3.5.1)                    
#>  CRAN (R 3.6.0)                    
#>  CRAN (R 3.5.1)                    
#>  CRAN (R 3.5.1)                    
#>  CRAN (R 3.5.1)                    
#>  CRAN (R 3.6.0)                    
#>  CRAN (R 3.6.0)                    
#>  CRAN (R 3.6.0)                    
#>  CRAN (R 3.5.1)                    
#>  CRAN (R 3.5.1)                    
#>  CRAN (R 3.5.1)                    
#>  CRAN (R 3.5.1)                    
#>  CRAN (R 3.5.1)                    
#>  CRAN (R 3.5.1)                    
#>  CRAN (R 3.6.0)                    
#>  CRAN (R 3.5.1)                    
#>  CRAN (R 3.5.1)                    
#>  CRAN (R 3.6.0)                    
#>  CRAN (R 3.5.1)                    
#>  CRAN (R 3.6.0)                    
#>  CRAN (R 3.6.0)                    
#>  Github (tidyverse/purrr@3b4a013)  
#>  CRAN (R 3.5.1)                    
#>  CRAN (R 3.6.0)                    
#>  CRAN (R 3.6.0)                    
#>  CRAN (R 3.6.0)                    
#>  CRAN (R 3.6.0)                    
#>  CRAN (R 3.6.0)                    
#>  CRAN (R 3.6.0)                    
#>  CRAN (R 3.5.1)                    
#>  CRAN (R 3.5.1)                    
#>  CRAN (R 3.5.1)                    
#>  CRAN (R 3.6.0)                    
#>  CRAN (R 3.6.0)                    
#>  CRAN (R 3.5.1)                    
#>  CRAN (R 3.5.1)                    
#>  Github (tidyverse/tibble@d06dd5c) 
#>  CRAN (R 3.6.0)                    
#>  CRAN (R 3.5.1)                    
#>  CRAN (R 3.5.1)                    
#>  Github (r-lib/usethis@923dd75)    
#>  CRAN (R 3.5.1)                    
#>  CRAN (R 3.5.1)                    
#>  CRAN (R 3.6.0)                    
#>  CRAN (R 3.5.1)                    
#>  CRAN (R 3.5.1)                    
#> 
#> [1] C:/Users/inp099/Documents/R/win-library/3.6
#> [2] C:/Program Files/R/R-devel/library

[IndrajeetPatil]

boot.ci also allows type = "all", which doesn't seem to work as expected with the tidier.

# setup
set.seed(123)
library(boot)

# data
clotting <- data.frame(
  u = c(5, 10, 15, 20, 30, 40, 60, 80, 100),
  lot1 = c(118, 58, 42, 35, 27, 25, 21, 19, 18),
  lot2 = c(69, 35, 26, 21, 18, 16, 13, 12, 12)
)

# model
g1 <- glm(lot2 ~ log(u), data = clotting, family = Gamma)

# bootstrapping function
bootfun <- function(d, i) {
  coef(update(g1, data = d[i, ]))
}

# using the custom bootstrapping function
bootres <- boot::boot(
  data = clotting,
  statistic = bootfun,
  R = 999
)

# getting all confidence intervals
broom::tidy(
  x = bootres,
  conf.int = TRUE,
  conf.method = "all"
)
#> Warning in FUN(X[[i]], ...): bootstrap variances needed for studentized
#> intervals

#> Warning in FUN(X[[i]], ...): bootstrap variances needed for studentized
#> intervals
#> # A tibble: 2 x 6
#>   term        statistic      bias std.error conf.low conf.high
#>   <chr>           <dbl>     <dbl>     <dbl> <list>   <list>   
#> 1 (Intercept)   -0.0239 -0.00184    0.00378 <NULL>   <NULL>   
#> 2 log(u)         0.0236  0.000559   0.00114 <NULL>   <NULL>

^{Created on 2019-01-12 by the reprex package (v0.2.1)}

[asbates]

For the studentized interval, boot.ci needs a way to compute the variance of the statistic of interest. Typically, this is done by having the function statistic that gets passed to boot return a vector where the first element is the statistic and the second is the variance of the statistic. In this case, corci would need to end with something like return(boot_df$estimate, boot_df$variance) instead of return(boot_df$estimate)

[alexpghayes]

What would tidy() do when conf.method = "all"? The standard result with an additional method column indexing method?

[IndrajeetPatil]

Yeah, but I am not sure how helpful that would be. Just wanted to bring up the possibility that this particular argument can get this particular input and the method should have some way to handle it. If not, it can print a helpful error message.

[alexpghayes]

Good point. Let's error on conf.method = "all" and use arg_match() on the conf.method options if we aren't already.

[ankrauska]

I'm going to work on this today at ChiRUnconf19

[ankrauska]

I fixed the issue for "norm", but not the issue with “stud”. Instead, I changed it so that “stud” is no longer supported, and so “all” is also no longer supported.

The lapply() does not work in lines 107-110 for conf.method = “stud”, if anyone wants to work on expanding support to “stud” and “all”.

[alexpghayes]

Quote: "Supporting conf.method = "stud" would not qualify as a beginner friendly issue"

[github-actions]

This issue has been automatically locked. If you believe you have found a related problem, please file a new issue (with a reprex: https://reprex.tidyverse.org) and link to this issue.