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