checking diagonal elements of a covariance matrix is insufficient to …

…diagnose singularity (#882)

* warnings/approximation of observation-level variance in .get_variance_distributional
Fixes #877

* fix test

DESCRIPTION

@@ -1,7 +1,7 @@
 Type: Package
 Package: insight
 Title: Easy Access to Model Information for Various Model Objects
-Version: 0.20.1
+Version: 0.20.1.1
 Authors@R: 
     c(person(given = "Daniel",
              family = "Lüdecke",

NEWS.md

@@ -1,3 +1,9 @@
+# insight 0.20.2
+
+## General
+
+* Improved accuracy of singularity-checks in `get_variance()`.
+
 # insight 0.20.1
 
 ## Bug fixes
@@ -21,9 +27,6 @@
 
 * Fixed errors in CRAN checks.
 
-* Fixed issues in `compact_list()` for objects that contained variables of
-  class `vctrs`.
-
 # insight 0.19.11
 
 ## General

R/compute_variances.R

@@ -12,6 +12,9 @@
   ## Major revisions and adaption to more complex models and other packages
   ## by Daniel Lüdecke
 
+  # needed for singularity check
+  check_if_installed("performance", reason = "to check for singularity")
+
   faminfo <- model_info(x, verbose = FALSE)
 
   if (any(faminfo$family == "truncated_nbinom1")) {
@@ -36,7 +39,7 @@
 
   # Test for non-zero random effects ((near) singularity)
   no_random_variance <- FALSE
-  if (.is_singular(x, vals, tolerance = tolerance) && !(component %in% c("slope", "intercept"))) {
+  if (performance::check_singularity(x, tolerance = tolerance) && !(component %in% c("slope", "intercept"))) {
     if (verbose) {
       format_warning(
         sprintf("Can't compute %s. Some variance components equal zero. Your model may suffer from singularity (see `?lme4::isSingular` and `?performance::check_singularity`).", name_full),

R/helper_functions.R

@@ -300,47 +300,6 @@
 
 
 
-
-# checks if a mixed model fit is singular or not. Need own function,
-# because lme4::isSingular() does not work with glmmTMB
-.is_singular <- function(x, vals, tolerance = 1e-5) {
-  check_if_installed("lme4", reason = "to compute variances for mixed models")
-
-  tryCatch(
-    {
-      if (inherits(x, "glmmTMB")) {
-        eigen_values <- list()
-        for (component in c("cond", "zi")) {
-          for (i in seq_along(vals$vc[[component]])) {
-            eigen_values <- c(eigen_values, list(eigen(vals$vc[[component]][[i]], only.values = TRUE)$values))
-          }
-        }
-        is_si <- any(vapply(eigen_values, min, numeric(1), na.rm = TRUE) < tolerance)
-      } else if (inherits(x, c("clmm", "cpglmm"))) {
-        is_si <- any(sapply(vals$vc, function(.x) any(abs(diag(.x)) < tolerance)))
-      } else if (inherits(x, "merMod")) {
-        theta <- lme4::getME(x, "theta")
-        diag.element <- lme4::getME(x, "lower") == 0
-        is_si <- any(abs(theta[diag.element]) < tolerance)
-      } else if (inherits(x, "MixMod")) {
-        vc <- diag(x$D)
-        is_si <- any(sapply(vc, function(.x) any(abs(.x) < tolerance)))
-      } else if (inherits(x, "lme")) {
-        is_si <- any(abs(stats::na.omit(as.numeric(diag(vals$vc))) < tolerance))
-      } else {
-        is_si <- FALSE
-      }
-
-      is_si
-    },
-    error = function(x) {
-      FALSE
-    }
-  )
-}
-
-
-
 # Filter parameters from Stan-model fits
 .filter_pars <- function(l, parameters = NULL, is_mv = NULL) {
   if (!is.null(parameters)) {

tests/testthat/test-rlmer.R

@@ -5,12 +5,12 @@ skip_if_not(getRversion() >= "4.1.0")
 data(sleepstudy, package = "lme4")
 
 set.seed(123)
-sleepstudy$mygrp <- sample(1:5, size = 180, replace = TRUE)
+sleepstudy$mygrp <- sample.int(5, size = 180, replace = TRUE)
 sleepstudy$mysubgrp <- NA
 for (i in 1:5) {
   filter_group <- sleepstudy$mygrp == i
   sleepstudy$mysubgrp[filter_group] <-
-    sample(1:30,
+    sample.int(30,
       size = sum(filter_group),
       replace = TRUE
     )
@@ -121,18 +121,12 @@ test_that("link_inverse", {
 })
 
 test_that("get_data", {
-  expect_identical(colnames(get_data(m1)), c("Reaction", "Days", "Subject"))
-  expect_identical(colnames(get_data(m1, effects = "all")), c("Reaction", "Days", "Subject"))
-  expect_identical(colnames(get_data(m1, effects = "random")), "Subject")
-  expect_identical(
-    colnames(get_data(m2)),
-    c("Reaction", "Days", "mysubgrp", "mygrp", "Subject")
-  )
-  expect_identical(
-    colnames(get_data(m2, effects = "all")),
-    c("Reaction", "Days", "mysubgrp", "mygrp", "Subject")
-  )
-  expect_identical(colnames(get_data(m2, effects = "random")), c("mysubgrp", "mygrp", "Subject"))
+  expect_named(get_data(m1), c("Reaction", "Days", "Subject"))
+  expect_named(get_data(m1, effects = "all"), c("Reaction", "Days", "Subject"))
+  expect_named(get_data(m1, effects = "random"), "Subject")
+  expect_named(get_data(m2), c("Reaction", "Days", "mysubgrp", "mygrp", "Subject"))
+  expect_named(get_data(m2, effects = "all"), c("Reaction", "Days", "mysubgrp", "mygrp", "Subject"))
+  expect_named(get_data(m2, effects = "random"), c("mysubgrp", "mygrp", "Subject"))
 })
 
 test_that("find_formula", {
@@ -203,16 +197,13 @@ test_that("get_response", {
 })
 
 test_that("get_predictors", {
-  expect_identical(colnames(get_predictors(m1)), "Days")
-  expect_identical(colnames(get_predictors(m2)), "Days")
+  expect_named(get_predictors(m1), "Days")
+  expect_named(get_predictors(m2), "Days")
 })
 
 test_that("get_random", {
-  expect_identical(colnames(get_random(m1)), "Subject")
-  expect_identical(
-    colnames(get_random(m2)),
-    c("mysubgrp", "mygrp", "Subject")
-  )
+  expect_named(get_random(m1), "Subject")
+  expect_named(get_random(m2), c("mysubgrp", "mygrp", "Subject"))
 })
 
 test_that("clean_names", {
@@ -279,13 +270,22 @@ test_that("get_variance", {
     tolerance = 1e-3
   )
 
+  expect_warning(
+    {
+      out <- get_variance(m2)
+    },
+    regex = "Can't compute random effect"
+  )
   expect_equal(
-    get_variance(m2),
+    out,
     list(
-      var.fixed = 914.841369705921, var.random = 1406.78220075798,
-      var.residual = 809.318117542254, var.distribution = 809.318117542254,
-      var.dispersion = 0,
-      var.intercept = c(`mysubgrp:mygrp` = 0, Subject = 1390.66848960835, mygrp = 16.1137111496379)
+      var.fixed = 914.841369705924, var.residual = 809.318117542254,
+      var.distribution = 809.318117542254, var.dispersion = 0,
+      var.intercept = c(
+        `mysubgrp:mygrp` = 0,
+        Subject = 1390.66848960834,
+        mygrp = 16.1137111496375
+      )
     ),
     tolerance = 1e-3
   )