add more r2 methods for glmmTMB

easystats · Nov 21, 2023 · 8329231 · 8329231
1 parent dd33955
commit 8329231
Show file tree

Hide file tree

Showing 5 changed files with 82 additions and 3 deletions.
diff --git a/NAMESPACE b/NAMESPACE
@@ -407,6 +407,7 @@ S3method(r2_coxsnell,coxph)
 S3method(r2_coxsnell,cpglm)
 S3method(r2_coxsnell,crch)
 S3method(r2_coxsnell,glm)
+S3method(r2_coxsnell,glmmTMB)
 S3method(r2_coxsnell,glmx)
 S3method(r2_coxsnell,logitmfx)
 S3method(r2_coxsnell,logitor)
@@ -463,6 +464,7 @@ S3method(r2_nagelkerke,coxph)
 S3method(r2_nagelkerke,cpglm)
 S3method(r2_nagelkerke,crch)
 S3method(r2_nagelkerke,glm)
+S3method(r2_nagelkerke,glmmTMB)
 S3method(r2_nagelkerke,glmx)
 S3method(r2_nagelkerke,logitmfx)
 S3method(r2_nagelkerke,logitor)

diff --git a/R/r2.R b/R/r2.R
@@ -489,17 +489,31 @@ r2.MixMod <- r2.merMod
 r2.rlmerMod <- r2.merMod
 
 #' @export
-r2.glmmTMB <- function(model, ci = NULL, tolerance = 1e-5, ...) {
+r2.glmmTMB <- function(model, ci = NULL, tolerance = 1e-5, verbose = TRUE, ...) {
   if (insight::is_mixed_model(model)) {
-    r2_nakagawa(model, ci = ci, tolerance = tolerance, ...)
+    return(r2_nakagawa(model, ci = ci, tolerance = tolerance, ...))
   } else {
+    if (!is.null(ci) && !is.na(ci)) {
+      return(.r2_ci(model, ci = ci, ...))
+    }
     mi <- insight::model_info(model, verbose = FALSE)
     if (mi$is_linear) {
-      .safe(.r2_lm_manual(model))
+      out <- .safe(.r2_lm_manual(model))
+    } else if (mi$is_logit && mi$is_bernoulli) {
+      out <- list(R2_Tjur = r2_tjur(model, model_info = mi, ...))
+      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) {
+      if (verbose) {
+        insight::format_warning("Can't calculate accurate R2 for binomial models that are not Bernoulli models.")
+      }
+      out <- NULL
     } else {
       insight::format_error("`r2()` does not support models of class `glmmTMB` without random effects and this link-function.") # nolint
     }
   }
+  out
 }
 
 #' @export

diff --git a/R/r2_coxsnell.R b/R/r2_coxsnell.R
@@ -89,6 +89,31 @@ r2_coxsnell.glm <- function(model, verbose = TRUE, ...) {
 r2_coxsnell.BBreg <- r2_coxsnell.glm
 
 
+#' @export
+r2_coxsnell.glmmTMB <- function(model, verbose = TRUE, ...) {
+  info <- list(...)$model_info
+  if (is.null(info)) {
+    info <- suppressWarnings(insight::model_info(model, verbose = FALSE))
+  }
+  if (info$is_binomial && !info$is_bernoulli) {
+    if (verbose) {
+      insight::format_alert("Can't calculate accurate R2 for binomial models that are not Bernoulli models.")
+    }
+    return(NULL)
+  } else {
+    dev <- stats::deviance(model)
+    # if no deviance, return NA
+    if (is.null(dev)) {
+      return(NULL)
+    }
+    null_dev <- stats::deviance(insight::null_model(model))
+    r2_coxsnell <- (1 - exp((dev - null_dev) / insight::n_obs(model, disaggregate = TRUE)))
+    names(r2_coxsnell) <- "Cox & Snell's R2"
+    r2_coxsnell
+  }
+}
+
+
 #' @export
 r2_coxsnell.nestedLogit <- function(model, ...) {
   n <- insight::n_obs(model, disaggregate = TRUE)

diff --git a/R/r2_nagelkerke.R b/R/r2_nagelkerke.R
@@ -78,6 +78,35 @@ r2_nagelkerke.glm <- function(model, verbose = TRUE, ...) {
 r2_nagelkerke.BBreg <- r2_nagelkerke.glm
 
 
+#' @export
+r2_nagelkerke.glmmTMB <- function(model, verbose = TRUE, ...) {
+  info <- list(...)$model_info
+  if (is.null(info)) {
+    info <- suppressWarnings(insight::model_info(model, verbose = FALSE))
+  }
+  if (info$is_binomial && !info$is_bernoulli) {
+    if (verbose) {
+      insight::format_warning("Can't calculate accurate R2 for binomial models that are not Bernoulli models.")
+    }
+    return(NULL)
+  } else {
+    dev <- stats::deviance(model)
+    # if no deviance, return NA
+    if (is.null(dev)) {
+      return(NULL)
+    }
+    null_dev <- stats::deviance(insight::null_model(model))
+    r2_cox <- (1 - exp((dev - null_dev) / insight::n_obs(model, disaggregate = TRUE)))
+    if (is.na(r2cox) || is.null(r2cox)) {
+      return(NULL)
+    }
+    r2_nagelkerke <- r2cox / (1 - exp(-null_dev / insight::n_obs(model, disaggregate = TRUE)))
+    names(r2_nagelkerke) <- "Nagelkerke's R2"
+    r2_nagelkerke
+  }
+}
+
+
 #' @export
 r2_nagelkerke.nestedLogit <- function(model, ...) {
   n <- insight::n_obs(model, disaggregate = TRUE)

diff --git a/tests/testthat/test-r2.R b/tests/testthat/test-r2.R
@@ -50,4 +50,13 @@ test_that("r2 glmmTMB, no ranef", {
   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)
 })