Simulated LRT sometimes shows strange behaviour

[florianhartig]

When comparing Poisson with NegBin in the presence of many REs, type I error rates are not properly controlled.

Not 100% sure if this code actually runs through, this is from the last time I played around with this problem

library(DHARMa)

library(lme4)
library(glmmTMB)

getP <- function(X=NULL){
  dat <- DHARMa::createData(sampleSize = 500, overdispersion = 0, numGroups = 100)
  #m0 =lme4:: glmer(observedResponse ~ Environment1 + (1|group), family = "poisson", data = dat)
  m0 =glmmTMB::glmmTMB(observedResponse ~ Environment1 + (1|group), family = poisson, data = dat)
  #summary(m0)
  a = logLik(m0)
  #m1 = lme4::glmer.nb(observedResponse ~ Environment1 + (1|group), data = dat)  
  m1 =glmmTMB::glmmTMB(observedResponse ~ Environment1 + (1|group), family = nbinom1, data = dat)
  #summary(m1)
  #getME(m1, "glmer.nb.theta")
  b = logLik(m1)
  x = DHARMa::simulateLRT(m0,m1, plot = F, n = 250)
  out = c(a, b, x$p.value)
  return(out)
}

getP()

# out <- replicate(10, getP())

parReplicate<-function(n, expr, simplify = "array", ncores = 1, libraries = NULL){
  cl = parallel::makeCluster (ncores)
  on.exit(parallel::stopCluster (cl))
  parallel::clusterExport(cl, varlist = "libraries")
  parallel::clusterEvalQ(cl, lapply(libraries, require, character.only = TRUE) )
  res <- parallel::parSapply (cl, X = 1:n, expr) 
  return(t(res))
}

res <- parReplicate(300, getP, ncores = 30, libraries = c("DHARMa", "lme4"))

hist(res[,3])