improved check_singularity (#801)

Author: strengejacke

DESCRIPTION

@@ -1,7 +1,7 @@
 Type: Package
 Package: performance
 Title: Assessment of Regression Models Performance
-Version: 0.13.0
+Version: 0.13.0.1
 Authors@R:
     c(person(given = "Daniel",
              family = "Lüdecke",
@@ -137,6 +137,7 @@ Suggests:
     qqplotr (>= 0.0.6),
     randomForest,
     RcppEigen,
+    reformulas,
     rempsyc,
     rmarkdown,
     rstanarm,

NEWS.md

@@ -1,3 +1,13 @@
+# performance (devel)
+
+## Changes
+
+* Singularity checks with `check_singularity()` are now more efficient and also
+  include the random effects for the dispersion component (from package *glmmTMB*).
+  Furthermore, a `check` argument allows to check for general singularity (for
+  the full model), or can return singularity checks for each random effects term
+  separately.
+
 # performance 0.13.0
 
 ## Breaking changes

R/check_singularity.R

@@ -5,8 +5,10 @@
 #'
 #' @param x A mixed model.
 #' @param tolerance Indicates up to which value the convergence result is
-#'    accepted. The larger `tolerance` is, the stricter the test
-#'    will be.
+#'   accepted. The larger `tolerance` is, the stricter the test
+#'   will be.
+#' @param check Indicates whether singularity check should be carried out for
+#'   the full model (`"model"`, the default), or per random effects term (`"terms"`).
 #' @param ... Currently not used.
 #'
 #' @return `TRUE` if the model fit is singular.
@@ -95,7 +97,10 @@
 #'   Reaction ~ Days + (1 | mygrp / mysubgrp) + (1 | Subject),
 #'   data = sleepstudy
 #' )
+#' # any singular fits?
 #' check_singularity(model)
+#' # singular fit for which particular random effects terms?
+#' check_singularity(model, check = "terms")
 #'
 #' \dontrun{
 #' # Fixing singularity issues using priors in glmmTMB
@@ -129,35 +134,61 @@ check_singularity <- function(x, tolerance = 1e-5, ...) {
 
 
 #' @export
-check_singularity.merMod <- function(x, tolerance = 1e-5, ...) {
-  insight::check_if_installed("lme4")
-
-  theta <- lme4::getME(x, "theta")
-  # diagonal elements are identifiable because they are fitted
-  #  with a lower bound of zero ...
-  diag.element <- lme4::getME(x, "lower") == 0
-  any(abs(theta[diag.element]) < tolerance)
+check_singularity.merMod <- function(x, tolerance = 1e-5, check = "model", ...) {
+  insight::check_if_installed(c("lme4", "reformulas"))
+
+  check <- insight::validate_argument(check, c("model", "terms"))
+  result <- list()
+  vv <- lme4::VarCorr(x)
+
+  re_names <- vapply(
+    reformulas::findbars(stats::formula(x)),
+    insight::safe_deparse,
+    FUN.VALUE = character(1)
+  )
+  result <- vapply(
+    vv,
+    function(x) det(x) < tolerance,
+    FUN.VALUE = logical(1)
+  )
+
+  switch(check,
+    model = any(unlist(result, use.names = FALSE)),
+    insight::compact_list(result)
+  )
 }
 
 #' @export
 check_singularity.rlmerMod <- check_singularity.merMod
 
 
+#' @rdname check_singularity
 #' @export
-check_singularity.glmmTMB <- function(x, tolerance = 1e-5, ...) {
-  insight::check_if_installed("lme4")
+check_singularity.glmmTMB <- function(x, tolerance = 1e-5, check = "model", ...) {
+  insight::check_if_installed(c("lme4", "reformulas"))
 
-  eigen_values <- list()
+  check <- insight::validate_argument(check, c("model", "terms"))
+  result <- list()
   vv <- lme4::VarCorr(x)
-  for (component in c("cond", "zi")) {
-    for (i in seq_along(vv[[component]])) {
-      eigen_values <- c(
-        eigen_values,
-        list(eigen(vv[[component]][[i]], only.values = TRUE)$values)
-      )
-    }
+
+  for (component in c("cond", "zi", "disp")) {
+    re_names <- vapply(
+      reformulas::findbars(stats::formula(x, component = component)),
+      insight::safe_deparse,
+      FUN.VALUE = character(1)
+    )
+    result[[component]] <- vapply(
+      vv[[component]],
+      function(x) det(x) < tolerance,
+      FUN.VALUE = logical(1)
+    )
+    names(result[[component]]) <- re_names
   }
-  any(vapply(eigen_values, min, numeric(1), na.rm = TRUE) < tolerance)
+
+  switch(check,
+    model = any(unlist(result, use.names = FALSE)),
+    insight::compact_list(result)
+  )
 }
 
 

inst/WORDLIST

@@ -261,6 +261,7 @@ brms
 brmsfit
 cauchy
 clusterable
+codecov
 confounder
 confounders
 concurvity

man/check_singularity.Rd

@@ -16,6 +16,13 @@ test_that("check_singularity, lme4", {
     data = sleepstudy
   ))
   expect_true(check_singularity(model))
+
+  # term-wise
+  out <- check_singularity(model, check = "terms")
+  expect_identical(
+    out,
+    c(`mysubgrp:mygrp` = TRUE, Subject = FALSE, mygrp = FALSE)
+  )
 })
 
 
@@ -36,4 +43,22 @@ test_that("check_singularity", {
     m2 <- glmmTMB::glmmTMB(y ~ 1 + (x | f), data = dd, REML = FALSE)
   }))
   expect_true(check_singularity(m2))
+
+  data(Salamanders, package = "glmmTMB")
+  m <- glmmTMB::glmmTMB(
+    count ~ spp + mined + (1 | site),
+    data = Salamanders[Salamanders$count > 0, , drop = FALSE],
+    family = glmmTMB::truncated_nbinom2(),
+    ziformula = ~ spp + (1 | site),
+    dispformula = ~ spp + (1 | site)
+  )
+  out <- check_singularity(m, check = "terms")
+  expect_identical(
+    out,
+    list(
+      cond = c(`1 | site` = TRUE),
+      zi = c(`1 | site` = TRUE),
+      disp = c(`1 | site` = TRUE)
+    )
+  )
 })