Search code examples
rtidyversehierarchical-clusteringstatistics-bootstrap

How to use tryCatch() during bootstrapping with modelr/tidyverse?


This question is similar to Using tryCatch() to catch a bootstrap loop, but I'm having trouble applying the suggestions to my case using the tidyverse method of bootstrapping. I am trying to obtain confidence interval estimates for coefficients from a complicated mixed effects model, but some models fail during the bootstrap, halting the entire bootstrapping process. I want to ignore these failed runs (but also count and store them to know which they are) and continue the bootstrap. I am also open to suggestions using the boot package with tryCatch. Example using the diamonds dataset:

diamonds <- diamonds 
diamonds$clarity <- factor(diamonds$clarity, ordered=FALSE)
diamonds$cut <- factor(diamonds$cut, ordered=FALSE)
diamonds$color <- factor(diamonds$color, ordered=FALSE)
diamonds <- diamonds[diamonds$price <= 500,] # truncate the data set for faster processing 

Random non-sensical model, but that runs with no errors:

my_mod <- glmmTMB(carat ~
     cut*x*poly(depth,3) + table + price + 
     (1|color) + (1|clarity),
     REML=TRUE,
     data = diamonds)
summary(my_mod)

I want to sample with replacement at the level of 'clarity' (i.e., by a cluster, not the observation).

set.seed(30)
my_boot <- diamonds %>% 
modelr::bootstrap(n = 20, id = 'clarity') %>%
group_by(clarity) %>%
mutate(fit = map(strap,
               ~glmmTMB(carat ~
     cut*x*poly(depth,3) +
     table + price + 
     (1|color) + (1|clarity),
     REML=T,
     data = data.frame(.))))

 Warning message:
 In fitTMB(TMBStruc) :
 Model convergence problem; false convergence (8). See vignette('troubleshooting')

Solution

  • Just wrap tryCatch around the call to glmmTMB. And write error and warning functions. The error messages of the type

    Timing stopped at: 0.484 0.071 0.556
    Timing stopped at: 0.919 0.052 0.972
    Timing stopped at: 0.542 0.04 0.583

    will still be displayed but my_boot$fit will have the errors or warnings.

    set.seed(30)
    my_boot <- diamonds %>% 
      modelr::bootstrap(n = 20, id = 'clarity') %>%
      group_by(clarity) %>%
      mutate(fit = map(strap,
                       ~tryCatch(glmmTMB(carat ~
                                  cut*x*poly(depth,3) +
                                  table + price + 
                                  (1|color) + (1|clarity),
                                REML=T,
                                data = data.frame(.)),
                                error = function(e) e,
                                warning = function(w) w)))
    

    You can then check the results with

    err <- sapply(my_boot$fit, inherits, "error")
    warn <- sapply(my_boot$fit, inherits, "warning")
    ok <- !(err | warn)
    

    And use these vectors to subset the list my_boot. For instance:

    my_boot[ok]
    

    And the same for err and warn.