florianhartig / DHARMa

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

testOverdispersionParametric() not accepting glmer.nb #47

Closed DesmondCampbell closed 6 years ago

DesmondCampbell commented 6 years ago

Are you interested in the below problem with DHARMa testOverdispersionParametric() ?

The plot generated by DHARMa is

image

Regards Desmond


> testOverdispersionParametric( fitLmer ) # doesn't work

        Chisq test for overdispersion in GLMMs

data:  Error in cat("data:  ", x$data.name, "\n", sep = "") : 
  argument 2 (type 'language') cannot be handled by 'cat'
> 
> class(fitLmer)
[1] "glmerMod"
attr(,"package")
[1] "lme4"
> 
> fitLmer
Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
 Family: Negative Binomial(7.6604)  ( log )
Formula: value ~ 1 + tissue + PLATFORM + (1 | metabId) + (1 | metabId:tissue) +      (1 | specimenId) + (1 | specimenPlatform)
   Data: dfData
      AIC       BIC    logLik  deviance  df.resid 
 868136.3  868258.9 -434053.1  868106.3     26286 
Random effects:
 Groups           Name        Std.Dev.
 metabId:tissue   (Intercept) 1.0322  
 metabId          (Intercept) 2.1001  
 specimenPlatform (Intercept) 0.1092  
 specimenId       (Intercept) 0.1192  
Number of obs: 26301, groups:  metabId:tissue, 2201; metabId, 660; specimenPlatform, 240; specimenId, 60
Fixed Effects:
            (Intercept)              tissueBrain              tissueHeart             tissueKidney      tissueOuter Medulla        PLATFORMLC/MS Neg      PLATFORMLC/MS Polar        PLATFORMLC/MS Pos  PLATFORMLC/MS Pos Early  
               12.96533                 -0.01186                  0.43224                  1.16883                  1.09679                  1.34259                  0.62587                  2.95287                  2.80882  
 PLATFORMLC/MS Pos Late  
                3.13469  
convergence code 0; 1 optimizer warnings; 0 lme4 warnings 
>

> oSimRes <- NA; system.time( oSimRes <- simulateResiduals( fitLmer, n=1000) )
   user  system elapsed 
  14.06    1.33   15.40 
> 
> try(plotSimulatedResiduals( oSimRes ))
> 
> 
> testOverdispersion( oSimRes )

        DHARMa nonparametric overdispersion test via IQR of scaled residuals against IQR expected under uniform

data:  oSimRes
dispersion = 1.0662, p-value < 2.2e-16
alternative hypothesis: overdispersion

Warning message:
In testOverdispersion(oSimRes) :
  You have called the non-parametric test for overdispersion based on the scaled residuals. Simulations show that this test is less powerful for detecting overdispersion than the default uniform test on the sclaed residuals, and a lot less powerful than a parametric overdispersion test, or the non-parametric test on re-simulated residuals. The test you called is only implemented for testing / development purposes, there is no scenario where it would be preferred. See vignette for details.
>
florianhartig commented 6 years ago

Hi, yes, it seems the parametric overdispersion creates an error when used on glmer.nb. I have to say I didn't try this because I never intended this to be used with the neg binom.

Can I ask why you want to test overdispersion with neg binom? As the dispersion is adjusted, one wouldn't expect overdispersion here?

-> fix: just throw a warning?

-> I also noted some unexpected behavior with refit = T, will open a new ticket for this

Reproducible example

library(lme4)

testData = createData(sampleSize = 300, overdispersion = 5, randomEffectVariance = 0, family = poisson())
fittedModel <- glmer.nb(observedResponse ~ Environment1 + (1|group) , data = testData)
simulationOutput <- simulateResiduals(fittedModel = fittedModel)
plot(simulationOutput, rank = T)

# works
testOverdispersion(simulationOutput)

# works, but produces overdispersion! -> new ticket
simulationOutput2 <- simulateResiduals(fittedModel = fittedModel, refit = T)
testOverdispersion(simulationOutput2)

# doesn't work
testOverdispersionParametric(fittedModel)
florianhartig commented 6 years ago

the testOverdispersionParametric(fittedModel) was in the end just a string conversion problem, fixed with https://github.com/florianhartig/DHARMa/commit/9bb149d4bf2f56169ecc39014c1c5a06dcf8f565

see further notes on simulated testOverdispersion in https://github.com/florianhartig/DHARMa/issues/48