florianhartig / DHARMa

Diagnostics for HierArchical Regession Models
http://florianhartig.github.io/DHARMa/
211 stars 22 forks source link

simulateResiduals(), with betabinomial error, generates warning message #68

Closed jhmaindonald closed 6 years ago

jhmaindonald commented 6 years ago

Here is an example:

codling1988 <- subset(DAAG::codling, year==1988)
codling1988$gp <- factor(codling1988$gp)
library(glmmTMB)
library(bbmle)
library(DHARMa)
Rho0.TMB <- glmmTMB(cbind(dead,tot-dead)~0+Cultivar/dose+(1|gp),
                   family=betabinomial(link='logit'),
                   data=subset(codling1988,dose>0))
Rho2.TMB <- update(Rho0.TMB, dispformula=~Cultivar+poly(dose,2),
                   formula=cbind(dead,tot-dead)~0+Cultivar/dose+(1|gp))
simRef <- simulateResiduals(Rho2.TMB)
## Warning message:
## In securityAssertion("nKcase") :
##  Message from DHARMa package: a security assertion was not met.
##  . . .
plot(simRef)  ## Plot suggests "Can do much better!"

Note: The model Rho2.TMB is not a good fit. One gets a better fit, thus (one has to fiddle the form of basis function to get the fit to converge):

Rho0.TMBclog <- glmmTMB(cbind(dead,tot-dead)~0+Cultivar/dose+(1|gp),
                        family=list(family='betabinomial',link='cloglog'),
                        data=subset(codling1988,dose>0))
  ## In spite of warning messages, fit seems OK
X3 <- mgcv::smoothCon(s(dose,k=4), data=subset(codling1988,dose>0), knots=NULL)[[1]]$X[,-3]
  ## Omit column for constant term,, so that formula does not duplicate it
Rho3.TMBclog <- update(Rho0.TMBclog, dispformula=~Cultivar+splines::ns(dose,3),
                          formula=cbind(dead,tot-dead)~0+Cultivar+X3+(1|gp))
  ## Updates without warning messages
simRef <- simulateResiduals(Rho2.TMB)
  ## Same warning message as before
plot(simRef)  ## Plot looks much better
bbmle::AICtab(Rho3.TMBclog, Rho2.TMB)

I have named output from simulateResiduals() simRef, because the intent is to create reference distributions. This distinguishes the object from the output from simulate().

The relevant part of the output from sessionInfo() is:

> sessionInfo()
R version 3.5.0 (2018-04-23)
Platform: x86_64-apple-darwin15.6.0 (64-bit)
Running under: macOS High Sierra 10.13.2

Matrix products: default
BLAS: /System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/Versions/A/libBLAS.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRlapack.dylib

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

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

other attached packages:
 [1] DHARMa_0.1.6.5   lme4_1.1-17      Matrix_1.2-14    mgcv_1.8-23     
 [5] nlme_3.1-137     lattice_0.20-35  bbmle_1.0.20     glmmTMB_0.2.1.0 
 [9] magrittr_1.5     kableExtra_0.8.0 knitr_1.20  

loaded via a namespace (and not attached):
. . .
[31] TMB_1.7.13          munsell_0.4.3       compiler_3.5.0 
. . .
florianhartig commented 6 years ago

Hi, thanks for reporting this, it's super helpful that people try out the new glmmTMB link before this is on CRAN, so that I can fix stuff.

I'll look into this and report back.

florianhartig commented 6 years ago

Note: I think results were fine, the reason for the warning was only that I hadn't include the betabinomial in the list of families that support k/n coding.