CODEWITHSUNDEEP
rtidyversehierarchical-clusteringstatistics-bootstrap
S.Bird
S.Bird

Reputation: 112

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

Upvotes: 0

Views: 296

Answers (1)

Rui Barradas
Rui Barradas

Reputation: 76641

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.

Upvotes: 1

Related Questions