Include changes from latest glmmTMB (#933)

* Include changes from latest glmmTMB * rename betadisp * nbinom * fix * run on glmmtmb dev * Bell family * fix for tweedie * be safe * fix for parameters * add test * news, desc * use family_params * docs
easystats · Sep 27, 2024 · fca4b12 · fca4b12
1 parent 68564af
commit fca4b12
Show file tree

Hide file tree

Showing 12 changed files with 89 additions and 30 deletions.
diff --git a/DESCRIPTION b/DESCRIPTION
@@ -1,7 +1,7 @@
 Type: Package
 Package: insight
 Title: Easy Access to Model Information for Various Model Objects
-Version: 0.20.4.4
+Version: 0.20.4.5
 Authors@R:
     c(person(given = "Daniel",
              family = "Lüdecke",
@@ -128,7 +128,7 @@ Suggests:
     geoR,
     ggeffects,
     GLMMadaptive,
-    glmmTMB,
+    glmmTMB (>= 1.1.10),
     glmtoolbox,
     gmnl,
     grDevices,
@@ -218,4 +218,3 @@ Config/testthat/edition: 3
 Config/testthat/parallel: true
 Config/Needs/website: easystats/easystatstemplate
 Config/Needs/check: stan-dev/cmdstanr
-Remotes: easystats/bayestestR, easystats/parameters
diff --git a/NEWS.md b/NEWS.md
@@ -12,6 +12,8 @@
 
 * Fixed rendering issue of the example in `?insight::display`.
 
+* Fixed issues due to recent changes in the *glmmTMB* package.
+
 # insight 0.20.4
 
 ## New supported models

diff --git a/R/compute_variances.R b/R/compute_variances.R
@@ -848,6 +848,7 @@
       nbinom = ,
       nbinom1 = ,
       nbinom2 = ,
+      nbinom12 = ,
       negbinomial = ,
       tweedie = ,
       `negative binomial` = exp(mu + 0.5 * as.vector(revar_null)),
@@ -888,6 +889,7 @@
           nbinom = ,
           nbinom1 = ,
           nbinom2 = ,
+          nbinom12 = ,
           negbinomial = ,
           `negative binomial` = ,
           `zero-inflated negative binomial` = .variance_family_nbinom(model, mu, sig, faminfo),
@@ -924,6 +926,7 @@
           nbinom = ,
           nbinom1 = ,
           nbinom2 = ,
+          nbinom12 = ,
           quasipoisson = ,
           negbinomial = ,
           `negative binomial` = ,
@@ -967,6 +970,7 @@
           nbinom = ,
           nbinom1 = ,
           nbinom2 = ,
+          nbinom12 = ,
           negbinomial = ,
           genpois = ,
           `negative binomial` = ((1 / mu) + (1 / dispersion_param))^-1,
@@ -979,6 +983,7 @@
           nbinom = ,
           nbinom1 = ,
           nbinom2 = ,
+          nbinom12 = ,
           negbinomial = ,
           genpois = ,
           `negative binomial` = (1 / mu) + (1 / dispersion_param),
@@ -1069,12 +1074,8 @@
     phi <- model@phi
     p <- model@p - 2
   } else {
-    if ("psi" %in% names(model$fit$par)) {
-      psi <- model$fit$par["psi"] # glmmmTMB >= 1.1.5
-    } else {
-      psi <- model$fit$par["thetaf"]
-    }
-    p <- unname(stats::plogis(psi) + 1)
+    check_if_installed("glmmTMB")
+    p <- unname(unlist(glmmTMB::family_params(model)))
   }
   phi * mu^p
 }
@@ -1092,6 +1093,14 @@
       return(rep(1e-16, length(mu)))
     }
     mu * (1 + sig)
+  } else if (identical(faminfo$family, "nbinom12")) {
+    # nbinom12-family from glmmTMB requires psi-parameter
+    if ("psi" %in% names(model$fit$par)) {
+      psi <- model$fit$par["psi"]
+    } else {
+      format_error("Could not extract psi-parameter for the distributional variance for nbinom12-family.")
+    }
+    stats::family(model)$variance(mu, sig, psi)
   } else {
     stats::family(model)$variance(mu, sig)
   }
@@ -1104,19 +1113,31 @@
 # ----------------------------------------------
 .variance_zinb <- function(model, sig, faminfo, family_var) {
   if (inherits(model, "glmmTMB")) {
+    if (identical(faminfo$family, "nbinom12")) {
+      # nbinom12-family from glmmTMB requires psi-parameter
+      if ("psi" %in% names(model$fit$par)) {
+        psi <- model$fit$par["psi"]
+      } else {
+        format_error("Could not extract psi-parameter for the distributional variance for nbinom12-family.")
+      }
+    }
     v <- stats::family(model)$variance
     # zi probability
     p <- stats::predict(model, type = "zprob")
     # mean of conditional distribution
     mu <- stats::predict(model, type = "conditional")
     # sigma
-    betad <- model$fit$par["betad"]
+    betadisp <- model$fit$par["betadisp"]
     k <- switch(faminfo$family,
-      gaussian = exp(0.5 * betad),
-      Gamma = exp(-0.5 * betad),
-      exp(betad)
+      gaussian = exp(0.5 * betadisp),
+      Gamma = exp(-0.5 * betadisp),
+      exp(betadisp)
     )
-    pvar <- (1 - p) * v(mu, k) + mu^2 * (p^2 + p)
+    if (identical(faminfo$family, "nbinom12")) {
+      pvar <- (1 - p) * v(mu, k, psi) + mu^2 * (p^2 + p)
+    } else {
+      pvar <- (1 - p) * v(mu, k) + mu^2 * (p^2 + p)
+    }
   } else if (inherits(model, "MixMod")) {
     v <- family_var
     p <- stats::plogis(stats::predict(model, type_pred = "link", type = "zero_part"))

diff --git a/R/find_formula.R b/R/find_formula.R
@@ -1117,6 +1117,7 @@ find_formula.glmmTMB <- function(x, verbose = TRUE, ...) {
   f.zi <- stats::formula(x, component = "zi")
   f.disp <- stats::formula(x, component = "disp")
 
+  # check for "empty" formulas
   if (identical(safe_deparse(f.zi), "~0") || identical(safe_deparse(f.zi), "~1")) {
     f.zi <- NULL
   }
@@ -1125,7 +1126,7 @@ find_formula.glmmTMB <- function(x, verbose = TRUE, ...) {
     f.disp <- NULL
   }
 
-
+  # extract random parts of formula
   f.random <- lapply(.findbars(f.cond), function(.x) {
     f <- safe_deparse(.x)
     stats::as.formula(paste0("~", f))
@@ -1147,16 +1148,30 @@ find_formula.glmmTMB <- function(x, verbose = TRUE, ...) {
     f.zirandom <- f.zirandom[[1]]
   }
 
+  f.disprandom <- lapply(.findbars(f.disp), function(.x) {
+    f <- safe_deparse(.x)
+    if (f == "NULL") {
+      return(NULL)
+    }
+    stats::as.formula(paste0("~", f))
+  })
+
+  if (length(f.disprandom) == 1L) {
+    f.disprandom <- f.disprandom[[1]]
+  }
 
+  # extract fixed effects parts
   f.cond <- stats::as.formula(.get_fixed_effects(f.cond))
   if (!is.null(f.zi)) f.zi <- stats::as.formula(.get_fixed_effects(f.zi))
+  if (!is.null(f.disp)) f.disp <- stats::as.formula(.get_fixed_effects(f.disp))
 
   f <- compact_list(list(
     conditional = f.cond,
     random = f.random,
     zero_inflated = f.zi,
     zero_inflated_random = f.zirandom,
-    dispersion = f.disp
+    dispersion = f.disp,
+    dispersion_random = f.disprandom
   ))
   .find_formula_return(f, verbose = verbose)
 }

diff --git a/R/format_table.R b/R/format_table.R
@@ -55,7 +55,7 @@
 #'
 #' @note `options(insight_use_symbols = TRUE)` overrides the `use_symbols` argument
 #'   and always displays symbols, if possible.
-#' @examplesIf require("rstanarm", warn.conflicts = FALSE) && require("parameters", , warn.conflicts = FALSE)
+#' @examplesIf require("rstanarm", warn.conflicts = FALSE) && require("parameters", warn.conflicts = FALSE) && packageVersion("parameters") > "0.22.2"
 #' format_table(head(iris), digits = 1)
 #'
 #' m <- lm(Sepal.Length ~ Species * Sepal.Width, data = iris)

diff --git a/R/get_residuals.R b/R/get_residuals.R
@@ -285,11 +285,11 @@ print.insight_residuals <- function(x, ...) {
     # mean of conditional distribution
     mu <- stats::predict(model, type = "conditional")
     # sigma
-    betad <- model$fit$par["betad"]
+    betadisp <- model$fit$par["betadisp"]
     k <- switch(faminfo$family,
-      gaussian = exp(0.5 * betad),
-      Gamma = exp(-0.5 * betad),
-      exp(betad)
+      gaussian = exp(0.5 * betadisp),
+      Gamma = exp(-0.5 * betadisp),
+      exp(betadisp)
     )
     pvar <- (1 - p) * v(mu, k) + mu^2 * (p^2 + p)
     pred <- stats::predict(model, type = "response") ## (1 - p) * mu

diff --git a/R/get_variances.R b/R/get_variances.R
@@ -135,7 +135,7 @@
 #' - Bernoulli (logistic) regression
 #' - Binomial regression (with other than binary outcomes)
 #' - Poisson and Quasi-Poisson regression
-#' - Negative binomial regression (including nbinom1 and nbinom2 families)
+#' - Negative binomial regression (including nbinom1, nbinom2 and nbinom12 families)
 #' - Gaussian regression (linear models)
 #' - Gamma regression
 #' - Tweedie regression
@@ -149,6 +149,7 @@
 #' - Compound Poisson regression
 #' - Generalized Poisson regression
 #' - Log-normal regression
+#' - Skew-normal regression
 #'
 #' Extracting variance components for models with zero-inflation part is not
 #' straightforward, because it is not definitely clear how the distribution-specific

diff --git a/R/utils_model_info.R b/R/utils_model_info.R
@@ -28,7 +28,7 @@
   # poisson family --------
 
   poisson_fam <-
-    fitfam %in% c("poisson", "quasipoisson", "genpois", "ziplss") |
+    fitfam %in% c("poisson", "quasipoisson", "genpois", "ziplss", "bell") |
       grepl("poisson", fitfam_lower, fixed = TRUE)
 
 
@@ -345,7 +345,7 @@
     dirichlet_fam || is.ordinal || zero.inf || is.censored || is.survival || is_binomtest ||
     is.categorical || hurdle || is.multinomial || is_chi2test || is_proptest || is_xtab) {
     linear_model <- FALSE
-  } else if (!(fitfam %in% c("Student's-t", "t Family", "gaussian", "Gaussian", "lognormal")) && !grepl("(\\st)$", fitfam)) {
+  } else if (!(fitfam %in% c("Student's-t", "t Family", "gaussian", "Gaussian", "lognormal", "skewnormal")) && !grepl("(\\st)$", fitfam)) {
     linear_model <- FALSE
   }
   if (!linear_model && is.survival && fitfam == "gaussian") {

diff --git a/man/format_table.Rd b/man/format_table.Rd
diff --git a/man/get_variance.Rd b/man/get_variance.Rd
diff --git a/tests/testthat/test-r2_nakagawa_negbin_zi.R b/tests/testthat/test-r2_nakagawa_negbin_zi.R
@@ -127,3 +127,19 @@ test_that("glmmTMB, Nbinom1 zero-inflated", {
   expect_equal(out$R2_conditional, 0.6051817, tolerance = 1e-4, ignore_attr = TRUE)
   expect_equal(out$R2_marginal, 0.5173316, tolerance = 1e-4, ignore_attr = TRUE)
 })
+
+
+test_that("glmmTMB, Nbinom12, zero-inflated", {
+  # results are in line with other nbinom families
+  skip_if_not_installed("glmmTMB", minimum_version = "1.1.10")
+  data(Salamanders, package = "glmmTMB")
+  m <- glmmTMB::glmmTMB(
+    count ~ mined + spp + (1 | site),
+    ziformula = ~mined,
+    family = glmmTMB::nbinom12(),
+    data = Salamanders, REML = TRUE
+  )
+  out <- performance::r2_nakagawa(m)
+  expect_equal(out$R2_conditional, 0.7857598, tolerance = 1e-4, ignore_attr = TRUE)
+  expect_equal(out$R2_marginal, 0.6742057, tolerance = 1e-4, ignore_attr = TRUE)
+})
diff --git a/tests/testthat/test-r2_nakagawa_tweedie.R b/tests/testthat/test-r2_nakagawa_tweedie.R
@@ -84,15 +84,19 @@ test_that("cplm, tweedie", {
   # results for glmmTMB are very close to cplm
   m0 <- glmmTMB::glmmTMB(
     Parasite ~ (1 | Population) + (1 | Container),
-    family = glmmTMB::tweedie(1),
+    family = glmmTMB::tweedie(),
+    start = list(psi = -12),
+    map = list(psi = factor(NA)),
     data = DataAll
   )
   m <- glmmTMB::glmmTMB(
     Parasite ~ Sex + Treatment + Habitat + (1 | Population) + (1 | Container),
-    family = glmmTMB::tweedie(1),
+    family = glmmTMB::tweedie(),
+    start = list(psi = -12),
+    map = list(psi = factor(NA)),
     data = DataAll
   )
   out3 <- performance::r2_nakagawa(m, null_model = m0)
-  expect_equal(out$R2_conditional, out3$R2_conditional, tolerance = 1e-2, ignore_attr = TRUE)
-  expect_equal(out$R2_marginal, out3$R2_marginal, tolerance = 1e-2, ignore_attr = TRUE)
+  expect_equal(out$R2_conditional, out3$R2_conditional, tolerance = 1e-1, ignore_attr = TRUE)
+  expect_equal(out$R2_marginal, out3$R2_marginal, tolerance = 1e-1, ignore_attr = TRUE)
 })