Error for simulateResiduals() when refit=T

[RegalPlatypus]

I'm afraid it might not be easy to give reproducible code for this, but I'll give as much information as I can. I receive no warnings in the model itself, but simulateResiduals(binary.glmer, n=500, refit=T) gives the following error:

Error in pwrssUpdate(pp, resp, tol = tolPwrss, GQmat = GQmat, compDev = compDev,  : 
  Downdated VtV is not positive definite
In addition: There were 12 warnings (use warnings() to see them)

Any suggestions would be appreciated. I always get this warning with this model, but I've run similar models on similar data structures that worked alright.

Model:
binary.glmer <- glmer(Status ~ Line + (1|Day), data=ParBin.data, family="binomial", nAGQ=25)

Data:

> print(head(ParBin.data))
     Line TreeID Day Status
1    Zoar    3.1   1      1
2    Zoar    3.3   1      0
3 Cropper    4.1   1      1
4 Cropper    4.1   1      1
5 Cropper    4.2   1      1
6 Cropper    4.2   1      1

> str(ParBin.data)
'data.frame':	306 obs. of  4 variables:
 $ Line  : Factor w/ 7 levels "BC3F3","Cropper",..: 7 7 2 2 2 2 2 2 2 2 ...
 $ TreeID: Factor w/ 37 levels "1.1","3.1","3.3",..: 2 3 5 5 6 6 6 6 8 9 ...
 $ Day   : Factor w/ 3 levels "1","2","3": 1 1 1 1 1 1 1 1 1 1 ...
 $ Status: int  1 0 1 1 1 1 1 0 0 1 ...

sessionInfo():

R version 3.3.2 (2016-10-31)
Platform: x86_64-apple-darwin13.4.0 (64-bit)
Running under: macOS Sierra 10.12.1

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

attached base packages:
[1] grid      stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] DHARMa_0.1.3       lsmeans_2.25       estimability_1.2   lmtest_0.9-34      zoo_1.7-14        
 [6] lme4_1.1-12        Matrix_1.2-8       EnvStats_2.2.1     boot_1.3-18        aod_1.3           
[11] multcompView_0.1-7 multcomp_1.4-6     TH.data_1.0-8      MASS_7.3-45        survival_2.40-1   
[16] mvtnorm_1.0-5      ggedit_0.0.2       miniUI_0.1.1       dplyr_0.5.0        colourpicker_0.3  
[21] plyr_1.8.4         plotly_4.5.6       ggplot2_2.2.1      scales_0.4.1       gridExtra_2.2.1   
[26] reshape2_1.4.2     shinyBS_0.61       shinyAce_0.2.1     shiny_1.0.0       

loaded via a namespace (and not attached):
 [1] purrr_0.2.2       splines_3.3.2     lattice_0.20-34   colorspace_1.3-2  htmltools_0.3.5   viridisLite_0.1.3
 [7] base64enc_0.1-3   nloptr_1.0.4      DBI_0.5-1         foreach_1.4.3     stringr_1.1.0     munsell_0.4.3    
[13] gtable_0.2.0      htmlwidgets_0.8   coda_0.19-1       codetools_0.2-15  labeling_0.3      httpuv_1.3.3     
[19] Rcpp_0.12.9       xtable_1.8-2      jsonlite_1.2      mime_0.5          digest_0.6.12     stringi_1.1.2    
[25] gap_1.1-16        tools_3.3.2       sandwich_2.3-4    magrittr_1.5      lazyeval_0.2.0    tibble_1.2       
[31] tidyr_0.6.1       qrnn_1.1.3        iterators_1.0.8   assertthat_0.1    minqa_1.2.4       httr_1.2.1       
[37] rstudioapi_0.6    R6_2.2.0          nlme_3.1-130     

[florianhartig]

Thanks for reporting this - I assume if you could dput your data, you would have done it, right?

But anyway, also without data, the problem is obviously that, for at least one of the 500 simulated datasets, the refit with lme4 does not converge and produces an error. This is fundamentally an lme4 problem, but I should changed the code to catch these errors, so that the failing model gets removed from the refit, and DHARMa doesn't stop. I have noted this under #18, but I don't have the time to do this immediately.

In the meantime, 3 remarks:

a) I don't really see a statistical reason to use refit = T in your case, so l would suggest to simply use refit = F

b) The problem may disappear when changing the optimizer, use control=glmerControl(optimizer="bobyqa")

c) The fact that the error occurs probably has something to do with the structure of your data. Either your dataset is very small, or the distribution of your predictors are such that the optimizer gets problems (do you have enough observations per RE?). You will find lots of case studies = error reports if you search, e.g. lme4/lme4#133 .

[RegalPlatypus]

a) Thanks for the feedback! If refit=T could be written with the changes you mentioned (such that the failed simulations are removed), I would like to use the testOverdispersion() function.

b) Using the bobyqa optimizer, refit=T worked, albeit with >50 warnings. Are the results still valid, or would you recommend trying a different method of estimating over/underdispersion?

c) This is very likely the case. There are a lot of missing observations both within the "Day" random effect and for individual trees in the fixed effects. It's far from ideal; I'm just trying to make the best of analysis possible out of the data I have.

Thank you for developing this package! It definitely makes diagnostics of GLMMs much much easier, especially for thick-skulled ecologists like me who simply can't take every statistics course offered in the course of an MS.

[florianhartig]

a) it looks like you are running a binomial with 1/0 data. With 1/0 data, you generally don't look for dispersion, because your data are either too high or too low, they cannot be "more spread out than expected" (definition of overdispersion) - the testOverdispersion() function will still react to problems with 1/0 data (it tests fitted likelihoods against expected), but if it shows a problem, you should call it a signal of misfit, not overdispersion. If you want to test for overdispersion in a 1/0 binomial, you need to group your data.

But again, asfaiks, in your case with one predictor, the standard DHARMa plot is totally sufficient for you - I'm going on the record here saying: if you see no clear problems there + no spatial autocorrelation, you're good!

b) I guess the warnings are convergence warnings - results are probably still fine, but one cannot be sure.

Glad you find the package useful!

[florianhartig]

Hi Regal, I have implemented something, but it is difficult to test for me because I haven't managed to find / create a dataset that makes lme4 crash in that way. Having your dataset via email (confidentially) would be a great help.

[RegalPlatypus]

Will do!

[florianhartig]

Implemented in d5cfe8f

Regal, if you install the development version (see main page), does this solve your problem?

[florianhartig]

Hi Regal, this should be fixed in DHARMa 0.1.5 on CRAN - a feedback if this solves your problem would be appreciated. Anyway, I will close this now.