add more r2 to non-mixed glmmTMB (#657)

DESCRIPTION

@@ -1,7 +1,7 @@
 Type: Package
 Package: performance
 Title: Assessment of Regression Models Performance
-Version: 0.10.8.3
+Version: 0.10.8.4
 Authors@R: 
     c(person(given = "Daniel",
              family = "Lüdecke",

R/r2.R

@@ -490,25 +490,40 @@ r2.rlmerMod <- r2.merMod
 
 #' @export
 r2.glmmTMB <- function(model, ci = NULL, tolerance = 1e-5, verbose = TRUE, ...) {
+  # most models are mixed models
   if (insight::is_mixed_model(model)) {
     return(r2_nakagawa(model, ci = ci, tolerance = tolerance, ...))
   } else {
     if (!is.null(ci) && !is.na(ci)) {
       return(.r2_ci(model, ci = ci, ...))
     }
+    # calculate r2 for non-mixed glmmTMB models here -------------------------
     info <- insight::model_info(model, verbose = FALSE)
+
     if (info$is_linear) {
+      # for linear models, use the manual calculation
       out <- .safe(.r2_lm_manual(model))
     } else if (info$is_logit && info$is_bernoulli) {
+      # logistic regression with binary outcome
       out <- list(R2_Tjur = r2_tjur(model, model_info = info, ...))
       attr(out, "model_type") <- "Logistic"
       names(out$R2_Tjur) <- "Tjur's R2"
       class(out) <- c("r2_pseudo", class(out))
     } else if (info$is_binomial && !info$is_bernoulli) {
+      # currently, non-bernoulli binomial models are not supported
       if (verbose) {
         insight::format_warning("Can't calculate accurate R2 for binomial models that are not Bernoulli models.")
       }
       out <- NULL
+    } else if ((info$is_poisson && !info$is_zero_inflated) || info$is_exponential) {
+      # Poisson-regression or Gamma uses Nagelkerke's R2
+      out <- list(R2_Nagelkerke = r2_nagelkerke(model, ...))
+      names(out$R2_Nagelkerke) <- "Nagelkerke's R2"
+      attr(out, "model_type") <- "Generalized Linear"
+      class(out) <- c("r2_pseudo", class(out))
+    } else if (info$is_zero_inflated) {
+      # zero-inflated models use the default method
+      out <- r2_zeroinflated(model)
     } else {
       insight::format_error("`r2()` does not support models of class `glmmTMB` without random effects and this link-function.") # nolint
     }

R/r2_nagelkerke.R

@@ -102,7 +102,8 @@ r2_nagelkerke.glmmTMB <- function(model, verbose = TRUE, ...) {
     return(NULL)
   }
 
-  null_dev <- stats::deviance(insight::null_model(model))
+  null_mod <- suppressWarnings(insight::null_model(model))
+  null_dev <- stats::deviance(null_mod)
   r2cox <- (1 - exp((dev - null_dev) / insight::n_obs(model, disaggregate = TRUE)))
 
   if (is.na(r2cox) || is.null(r2cox)) {

R/r2_zeroinflated.R

@@ -63,7 +63,6 @@ r2_zeroinflated <- function(model, method = c("default", "correlation")) {
   k <- length(insight::find_parameters(model)[["conditional"]])
 
   y <- insight::get_response(model, verbose = FALSE)
-  # pred <- stats::predict(model, type = "response")
 
   var_fixed <- sum((stats::fitted(model) - mean(y))^2)
   var_resid <- sum(stats::residuals(model, type = "pearson")^2)

tests/testthat/test-r2.R

@@ -40,23 +40,74 @@ test_that("r2 glm, ci", {
   )
 })
 
-test_that("r2 glmmTMB, no ranef", {
-  skip_if_not_installed("glmmTMB")
-  data(Owls, package = "glmmTMB")
-  m <- glmmTMB::glmmTMB(NegPerChick ~ BroodSize + ArrivalTime, data = Owls)
-  out <- r2(m)
-  expect_equal(out$R2, 0.05597288, tolerance = 1e-3, ignore_attr = TRUE)
-  # validate against lm
-  m2 <- lm(NegPerChick ~ BroodSize + ArrivalTime, data = Owls)
-  out2 <- r2(m2)
-  expect_equal(out$R2, out2$R2, tolerance = 1e-3, ignore_attr = TRUE)
-  # binomial
-  data(mtcars)
-  m <- glmmTMB::glmmTMB(am ~ mpg, data = mtcars, family = binomial())
-  out <- r2(m)
-  expect_equal(out[[1]], 0.3677326, tolerance = 1e-3, ignore_attr = TRUE)
-  # validate against glm
-  m2 <- glm(am ~ mpg, data = mtcars, family = binomial())
-  out2 <- r2(m2)
-  expect_equal(out[[1]], out[[1]], tolerance = 1e-3, ignore_attr = TRUE)
-})
+# glmmTMB, non-mixed --------------------------------------------------------
+
+skip_if_not_installed("withr")
+withr::with_environment(
+  new.env(),
+  test_that("r2 glmmTMB, no ranef", {
+    skip_if_not_installed("glmmTMB")
+    data(Owls, package = "glmmTMB")
+    # linear ---------------------------------------------------------------
+    m <- glmmTMB::glmmTMB(NegPerChick ~ BroodSize + ArrivalTime, data = Owls)
+    out <- r2(m)
+    expect_equal(out$R2, 0.05597288, tolerance = 1e-3, ignore_attr = TRUE)
+    # validate against lm
+    m2 <- lm(NegPerChick ~ BroodSize + ArrivalTime, data = Owls)
+    out2 <- r2(m2)
+    expect_equal(out$R2, out2$R2, tolerance = 1e-3, ignore_attr = TRUE)
+    # binomial -------------------------------------------------------------
+    data(mtcars)
+    m <- glmmTMB::glmmTMB(am ~ mpg, data = mtcars, family = binomial())
+    out <- r2(m)
+    expect_equal(out[[1]], 0.3677326, tolerance = 1e-3, ignore_attr = TRUE)
+    # validate against glm
+    m2 <- glm(am ~ mpg, data = mtcars, family = binomial())
+    out2 <- r2(m2)
+    expect_equal(out[[1]], out2[[1]], tolerance = 1e-3, ignore_attr = TRUE)
+    # poisson --------------------------------------------------------------
+    d <<- data.frame(
+      counts = c(18, 17, 15, 20, 10, 20, 25, 13, 12),
+      outcome = gl(3, 1, 9),
+      treatment = gl(3, 3)
+    )
+    m <- glmmTMB::glmmTMB(counts ~ outcome + treatment, family = poisson(), data = d)
+    out <- r2(m)
+    expect_equal(out[[1]], 0.6571698, tolerance = 1e-3, ignore_attr = TRUE)
+    # validate against glm
+    m2 <- glm(counts ~ outcome + treatment, family = poisson(), data = d)
+    out2 <- r2(m2)
+    expect_equal(out[[1]], out2[[1]], tolerance = 1e-3, ignore_attr = TRUE)
+    # zero-inflated --------------------------------------------------------------
+    skip_if_not_installed("pscl")
+    data(bioChemists, package = "pscl")
+    m <- glmmTMB::glmmTMB(
+      art ~ fem + mar + kid5 + ment,
+      ziformula = ~ kid5 + phd,
+      family = poisson(),
+      data = bioChemists
+    )
+    out <- r2(m)
+    expect_equal(out[[1]], 0.1797549, tolerance = 1e-3, ignore_attr = TRUE)
+    # validate against pscl::zeroinfl
+    m2 <- pscl::zeroinfl(
+      art ~ fem + mar + kid5 + ment | kid5 + phd,
+      data = bioChemists
+    )
+    out2 <- r2(m2)
+    expect_equal(out[[1]], out2[[1]], tolerance = 1e-3, ignore_attr = TRUE)
+    # Gamma --------------------------------------------------------------
+    clotting <<- data.frame(
+      u = c(5, 10, 15, 20, 30, 40, 60, 80, 100),
+      lot1 = c(118, 58, 42, 35, 27, 25, 21, 19, 18),
+      lot2 = c(69, 35, 26, 21, 18, 16, 13, 12, 12)
+    )
+    m <- suppressWarnings(glmmTMB::glmmTMB(lot1 ~ log(u), data = clotting, family = Gamma()))
+    out <- r2(m)
+    expect_equal(out[[1]], 0.996103, tolerance = 1e-3, ignore_attr = TRUE)
+    # validate against glm
+    m2 <- glm(lot1 ~ log(u), data = clotting, family = Gamma())
+    out2 <- r2(m2)
+    expect_equal(out[[1]], out2[[1]], tolerance = 1e-3, ignore_attr = TRUE)
+  })
+)