An update on this question: it's now an issue on github, at https://github.com/lme4/lme4/issues/231 Thanks to Ben Bolker for following this up - see the comments for updates on progress. Paul Johnson
On 14 Aug 2014, at 14:15, Paul Johnson <Paul.Johnson at glasgow.ac.uk> wrote:
Hi, I?ve been using confint.merMod to get 95% CIs by parametric bootstrapping, and have been getting the warnings of the type "In bootMer(object, bootFun, nsim = nsim, ...) : some bootstrap runs failed (30/100)?. (I realise than 100 samples are too few - this is just an example). The function still returns CIs using the bootstrap samples that didn?t fail, but I don?t think these CIs are valid as the failures are very unlikely to be a random sample. Here?s an example using the grouseticks data set that comes with lme4 (see ?grouseticks). I chose this data set because it?s a classic example data set so presumably isn?t pathological, but I?ve been encountering the same problem with other apparently healthy glmer fits.
library(lme4) full_mod1 <- glmer(TICKS ~ YEAR+scale(HEIGHT)+(1|BROOD)+(1|INDEX)+(1|LOCATION), family="poisson", data=grouseticks) set.seed(12345) confint(full_mod1, method = "boot", boot.type = "basic", nsim = 100, parallel = "multicore", ncpus = 8)
Computing bootstrap confidence intervals ...
2.5 % 97.5 %
sd_(Intercept)|INDEX 0.4624700492 0.6550336
sd_(Intercept)|BROOD 0.5740796527 1.0482547
sd_(Intercept)|LOCATION 0.2497343230 1.0574176
(Intercept) -0.0002611865 0.7208569
YEAR96 0.7770837606 1.5857202
YEAR97 -1.4648325648 -0.4484355
scale(HEIGHT) -1.1335522182 -0.5879858
Warning message:
In bootMer(object, bootFun, nsim = nsim, ...) :
some bootstrap runs failed (30/100)
(NB I scaled height to prevent a warning but it makes no difference to the likelihood.)
I tried to replicate the problem using bootMer directly to produce fixed effect estimates from 10 samples:
bootMer(full_mod1, fixef, nsim = 10)$t
(Intercept) YEAR96 YEAR97 scale(HEIGHT) [1,] 0.6387534 0.7467317 -1.1291369 -0.8741711 [2,] NA NA NA NA [3,] 0.1920928 1.3666922 -0.5821328 -0.6742498 [4,] NA NA NA NA [5,] NA NA NA NA [6,] 0.5312426 0.9639612 -1.1041643 -0.7550492 [7,] NA NA NA NA [8,] 0.2328530 1.1016917 -0.8010541 -0.8480263 [9,] 0.2685594 1.2806690 -0.7452719 -0.8302675 [10,] 0.1726928 1.1423894 -0.8989542 -1.0344128 Warning message: In bootMer(full_mod1, fixef, nsim = 10) : some bootstrap runs failed (4/10) ?and get the same problem, not surprisingly. Here?s a manual version of what (I believe) bootMer is doing:
simTICKS.tab <- simulate(full_mod1, nsim=10) t(apply(simTICKS.tab, 2, function(simTICKS) fixef(glmer(simTICKS ~ YEAR+scale(HEIGHT)+(1|BROOD)+(1|INDEX)+(1|LOCATION), family="poisson", data=grouseticks))))
(Intercept) YEAR96 YEAR97 scale(HEIGHT) sim_1 0.1793276 1.2265976 -0.5753168 -0.6435871 sim_2 0.3918824 1.4583529 -0.8958396 -0.8513233 sim_3 0.2565916 1.1230342 -0.9369728 -1.0817019 sim_4 0.3507285 1.1865904 -1.0654895 -0.9193403 sim_5 0.5688213 0.7291811 -0.9734654 -0.8575988 sim_6 0.6160919 1.0552260 -1.5139770 -1.0050926 sim_7 0.1594954 1.1342809 -1.0624134 -0.8107930 sim_8 0.7956874 0.8290606 -1.2491060 -1.0024299 sim_9 0.3081545 1.2677562 -1.0398100 -0.9458849 sim_10 0.4325345 0.9390844 -0.9241206 -0.7341589 ?no errors, not even a warning. The discrepancy between bootMer and the DIY version above appears to arise from the bootMer fitting function, refit, being less tolerant of warning signs (non-convergence?) than glmer, which can be seen by putting refit into the DIY bootstrap:
t(apply(simTICKS.tab, 2, function(simTICKS) fixef(refit(full_mod1, simTICKS))))
Warning messages: 1: In pwrssUpdate(pp, resp, tolPwrss, GQmat, compDev, fac, verbose) : Cholmod warning 'not positive definite' at file:../Cholesky/t_cholmod_rowfac.c, line 431 2: In pwrssUpdate(pp, resp, tolPwrss, GQmat, compDev, fac, verbose) : Cholmod warning 'not positive definite' at file:../Cholesky/t_cholmod_rowfac.c, line 431 Error in t(apply(simTICKS.tab, 2, function(simTICKS) fixef(refit(full_mod1, : error in evaluating the argument 'x' in selecting a method for function 't': Error in t(apply(simTICKS.tab, 2, function(simTICKS) fixef(refit(full_mod1, : (maxstephalfit) PIRLS step-halvings failed to reduce deviance in pwrssUpdate My questions are: Should I trust the error- and warning-free DIY bootstrap results? I would say ?yes?, on the (slightly vague) basis that glmer produced no warnings, and that the original model fits well to what is a pretty information-rich data set. So is bootMer being oversensitive? I?m aware that there have been some teething problems with the optimisers and convergence checks in lme4 1.0+ ? are these ongoing? I?m writing a tutorial document in which I?d like to demonstrate parametric bootstrapping for CIs. I could use the DIY approach, but I?d much rather use confint, so it would be great if there was a simple fix to this problem. Thanks in advance, Paul Johnson (the Glasgow one, not the Kansas or Oxford ones)
sessionInfo()
R version 3.1.1 (2014-07-10) Platform: x86_64-apple-darwin13.1.0 (64-bit) locale: [1] en_GB.UTF-8/en_GB.UTF-8/en_GB.UTF-8/C/en_GB.UTF-8/en_GB.UTF-8 attached base packages: [1] parallel stats graphics grDevices utils datasets methods base other attached packages: [1] lme4_1.1-7 Rcpp_0.11.2 Matrix_1.1-4 loaded via a namespace (and not attached): [1] boot_1.3-11 grid_3.1.1 lattice_0.20-29 MASS_7.3-33 minqa_1.2.3 nlme_3.1-117 nloptr_1.0.0 splines_3.1.1 tools_3.1.1
_______________________________________________ R-sig-mixed-models at r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models