224 addition of glmmtmb to improve convergency for poisson in mixed effect models

DESCRIPTION

@@ -42,6 +42,7 @@ Depends:
 Imports: 
     doParallel,
     foreach,
+    glmmTMB,
     graph,
     lme4,
     mclogit,

R/build_score_cache.R

@@ -38,6 +38,7 @@
 #' @param epsilon Defaults to 1e-8. Positive convergence tolerance \eqn{\epsilon} that is directly passed to the \code{control} argument of \code{mclogit::mblogit} as \code{mclogit.control}. Only for \code{method='mle', group.var=...}.
 #' @param ncores The number of cores to parallelize to, see \sQuote{Details}. If >0, the number of CPU cores to be used. -1 for all available -1 core. Only for \code{method="mle"}.
 #' @param cluster.type The type of cluster to be used, see \code{?parallel::makeCluster}. \code{abn} then defaults to \code{"PSOCK"} on Windows and \code{"FORK"} on Unix-like systems. With "FORK" the child process are started with \code{rscript_args = "--no-environ"} to avoid loading the whole workspace into each child.
+#' @param only_glmmTMB_poisson logical, if TRUE only use \code{glmmTMB} to fit Poisson nodes with random effects. This is useful if \code{glmer} fails due to convergence issues. Default is FALSE.
 #' @param seed a non-negative integer which sets the seed in \code{set.seed(seed)}.
 #'
 #' @details
@@ -64,6 +65,7 @@
 #'                         trace.mblogit = FALSE,
 #'                         catcov.mblogit = "free",
 #'                         epsilon = 1e-6,
+#'                         only_glmmTMB_poisson=FALSE,
 #'                         seed = 9062019L)
 #' ctrlbayes <- abn::build.control(method = "bayes",
 #'                            max.mode.error = 10,
@@ -121,6 +123,7 @@ build.control <-
             trace.mblogit = FALSE,
             catcov.mblogit = "free",
             epsilon = 1e-6,
+            only_glmmTMB_poisson=FALSE,
             seed = 9062019L) {
     if (method == "bayes") {
       ctrl <-
@@ -167,6 +170,7 @@ build.control <-
           trace.mblogit = trace.mblogit,
           catcov.mblogit = catcov.mblogit,
           epsilon = epsilon,
+          only_glmmTMB_poisson=only_glmmTMB_poisson,
           seed = seed
         )
     } else {

R/build_score_cache_mle.R

@@ -100,7 +100,10 @@ forLoopContent <-
                 # if fit is still NULL, do not modify model further as this would change the structure of the dag but return very low score (further down)
             },
             poisson = {
-                if (verbose) {message(paste("using glmer with model:", deparse1(model)))} else NA
+              if (verbose) {message(paste("using glmer with model:", deparse1(model)))} else NA
+
+              fit <- NULL
+              if(control[["only_glmmTMB_poisson"]] == FALSE) {
                 tryCatch({
                   fit <- lme4::glmer(model, data = data.df.grouping, family = "poisson")
                 }, error=function(e)NULL)
@@ -118,8 +121,19 @@ forLoopContent <-
                                                                                            ftol = control[["ftol_abs"]])))
                   }, error=function(e)NULL)
                 }
+              } else if(control[["only_glmmTMB_poisson"]] == TRUE) {
+                if (is.null(fit)){
+                  # try glmmTMB as alternative
+                  if (verbose) {message(paste("trying glmmTMB with model:", deparse1(model)))} else NA
+                  tryCatch({
+                    fit <- glmmTMB::glmmTMB(model, data = data.df.grouping, family = "poisson")
+                  }, error=function(e)NULL)
+                }
+              } else {
+                stop("Invalid 'only_glmmTMB_poisson' argument. Must be one of TRUE or FALSE.")
+              }
 
-                # if fit is still NULL, do not modify model further as this would change the structure of the dag but return very low score (further down)
+              # if fit is still NULL, do not modify model further as this would change the structure of the dag but return very low score (further down)
             },
             multinomial = {
               if (length(parents.names) == 0){
@@ -216,7 +230,7 @@ forLoopContent <-
     # collect values to return
     if(!is.null(fit)){
       if (verbose) {
-        message("Sccessfully fitted local model.")
+        message("Successfully fitted local model.")
       }
       fit_loglik <- logLik(fit)
       fit_aic <- AIC(fit)
@@ -321,7 +335,7 @@ forLoopContent <-
     # Prepare return values
     if (!is.null(fit)) {
       if (verbose) {
-        message("Sccessfully fitted local model.")
+        message("Successfully fitted local model.")
       }
       c(fit$loglik,
         fit$aic,

tests/testthat/test-build_score_cache_mle.R

@@ -416,3 +416,36 @@ test_that("forLoopContent() prints local model when verbose.", {
     skip("`forLoopContent()` is tested mainly on Unix-like systems.")
   }
 })
+
+test_that("Poisson nodes step into calling glmmTMB.", {
+  if(.Platform$OS.type == "unix") {
+    suppressMessages({
+      # Suppress messages that are not related to the test but are printed when verbose
+      capture.output({
+        # load(file = "tests/testthat/testdata/forLoopContent_data.Rdata")
+        load(file = "testdata/forLoopContent_data.Rdata")
+
+        verbose <- TRUE
+        # with group.var
+        expect_message(
+          forLoopContent(row.num = 9,
+                         mycache = mycache,
+                         data.dists = data.dists,
+                         data.df.multi = data.df.multi,
+                         adj.vars = adj.vars,
+                         data.df = data.df,
+                         data.df.lvl = data.df.lvl,
+                         group.var = group.var,
+                         group.ids = group.ids,
+                         control = build.control(method = "mle", only_glmmTMB_poisson=TRUE),
+                         n = nvars,
+                         verbose = verbose),
+          regexp = "trying glmmTMB with model"
+        )
+      },
+      file = "/dev/null")
+    })
+  } else {
+    skip("`forLoopContent()` is tested mainly on Unix-like systems.")
+  }
+})