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,

R/fitabn.R

@@ -26,6 +26,7 @@
 #'                        epsilon = 1e-6,
 #'                        ncores = 2,
 #'                        cluster.type = "PSOCK",
+#'                        only_glmmTMB_poisson = FALSE,
 #'                        seed = 9062019L)
 #' ctrlbayes <- abn::fit.control(method = "bayes",
 #'                          mean = 0,
@@ -102,6 +103,7 @@ fit.control <-
            trace.mblogit = FALSE,
            catcov.mblogit = "free",
            epsilon = 1e-6,
+           only_glmmTMB_poisson = FALSE,
            seed = 9062019L) {
     if (method == "bayes") {
       ctrl <-
@@ -155,6 +157,7 @@ fit.control <-
           trace.mblogit = trace.mblogit,
           catcov.mblogit = catcov.mblogit,
           epsilon = epsilon,
+          only_glmmTMB_poisson = only_glmmTMB_poisson,
           seed = seed
         )
     } else {