add more r2 to non-mixed glmmTMB

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

@@ -43,20 +43,65 @@ test_that("r2 glm, ci", {
 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
+  # 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)
+  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)
 })