Error for R2 for beta binomial glmmTMB model (#788)

Author: strengejacke

DESCRIPTION

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

NAMESPACE

@@ -473,6 +473,7 @@ S3method(r2_mcfadden,clm)
 S3method(r2_mcfadden,clm2)
 S3method(r2_mcfadden,cpglm)
 S3method(r2_mcfadden,glm)
+S3method(r2_mcfadden,glmmTMB)
 S3method(r2_mcfadden,glmx)
 S3method(r2_mcfadden,logitmfx)
 S3method(r2_mcfadden,logitor)

NEWS.md

@@ -16,6 +16,9 @@
 
 * Increased accuracy for `check_convergence()` for *glmmTMB* models.
 
+* `r2()` and `r2_mcfadden()` now support beta-binomial (non-mixed) models from
+  package *glmmTMB*.
+
 ## Bug fixes
 
 * `check_outliers()` did not warn that no numeric variables were found when only

R/r2.R

@@ -516,6 +516,7 @@ r2.glmmTMB <- function(model, ci = NULL, tolerance = 1e-5, verbose = TRUE, ...)
     }
     # calculate r2 for non-mixed glmmTMB models here -------------------------
     info <- insight::model_info(model, verbose = FALSE)
+    matrix_response <- grepl("cbind", insight::find_response(model), fixed = TRUE)
 
     if (info$is_linear) {
       # for linear models, use the manual calculation
@@ -526,6 +527,17 @@ r2.glmmTMB <- function(model, ci = NULL, tolerance = 1e-5, verbose = TRUE, ...)
       attr(out, "model_type") <- "Logistic"
       names(out$R2_Tjur) <- "Tjur's R2"
       class(out) <- c("r2_pseudo", class(out))
+    } else if (info$is_betabinomial) {
+      # currently, beta-binomial models without proportion response are not supported
+      if (matrix_response) {
+        if (verbose) {
+          insight::format_warning("Can't calculate accurate R2 for beta-binomial models with matrix-response formulation.")
+        }
+        out <- NULL
+      } else {
+        # betabinomial default to mcfadden, see pscl:::pR2Work
+        out <- r2_mcfadden(model)
+      }
     } else if (info$is_binomial && !info$is_bernoulli) {
       # currently, non-bernoulli binomial models are not supported
       if (verbose) {

R/r2_coxsnell.R

@@ -69,13 +69,22 @@ r2_coxsnell.glm <- function(model, verbose = TRUE, ...) {
   if (is.null(info)) {
     info <- suppressWarnings(insight::model_info(model, verbose = FALSE))
   }
+  matrix_response <- grepl("cbind", insight::find_response(model), fixed = TRUE)
+
   # Cox & Snell's R2 is not defined for binomial models that are not Bernoulli models
-  if (info$is_binomial && !info$is_bernoulli && class(model)[1] == "glm") {
+  if (info$is_binomial && !info$is_betabinomial && !info$is_bernoulli && class(model)[1] %in% c("glm", "glmmTMB")) {
     if (verbose) {
       insight::format_alert("Can't calculate accurate R2 for binomial models that are not Bernoulli models.")
     }
     return(NULL)
   }
+  # currently, beta-binomial models without proportion response are not supported
+  if (info$is_betabinomial && matrix_response) {
+    if (verbose) {
+      insight::format_warning("Can't calculate accurate R2 for beta-binomial models with matrix-response formulation.")
+    }
+    return(NULL)
+  }
   # if no deviance, return NULL
   if (is.null(model$deviance)) {
     return(NULL)
@@ -96,7 +105,7 @@ r2_coxsnell.glmmTMB <- function(model, verbose = TRUE, ...) {
     info <- suppressWarnings(insight::model_info(model, verbose = FALSE))
   }
   # Cox & Snell's R2 is not defined for binomial models that are not Bernoulli models
-  if (info$is_binomial && !info$is_bernoulli) {
+  if (info$is_binomial && !info$is_bernoulli && !info$is_betabinomial) {
     if (verbose) {
       insight::format_alert("Can't calculate accurate R2 for binomial models that are not Bernoulli models.")
     }

R/r2_mcfadden.R

@@ -63,13 +63,21 @@ r2_mcfadden.glm <- function(model, verbose = TRUE, ...) {
   if (is.null(info)) {
     info <- suppressWarnings(insight::model_info(model, verbose = FALSE))
   }
+  matrix_response <- grepl("cbind", insight::find_response(model), fixed = TRUE)
 
-  if (info$is_binomial && !info$is_bernoulli && class(model)[1] == "glm") {
+  if (info$is_binomial && !info$is_betabinomial && !info$is_bernoulli && class(model)[1] %in% c("glm", "glmmTMB")) {
     if (verbose) {
       insight::format_warning("Can't calculate accurate R2 for binomial models that are not Bernoulli models.")
     }
     return(NULL)
   }
+  # currently, beta-binomial models without proportion response are not supported
+  if (info$is_betabinomial && matrix_response) {
+    if (verbose) {
+      insight::format_warning("Can't calculate accurate R2 for beta-binomial models with matrix-response formulation.")
+    }
+    return(NULL)
+  }
 
   l_null <- insight::get_loglikelihood(stats::update(model, ~1))
   .r2_mcfadden(model, l_null)
@@ -99,6 +107,9 @@ r2_mcfadden.brmultinom <- r2_mcfadden.glm
 #' @export
 r2_mcfadden.censReg <- r2_mcfadden.glm
 
+#' @export
+r2_mcfadden.glmmTMB <- r2_mcfadden.glm
+
 #' @export
 r2_mcfadden.truncreg <- r2_mcfadden.glm
 

performance.Rproj

@@ -1,4 +1,5 @@
 Version: 1.0
+ProjectId: af6facf3-033e-40d4-ac22-2830774814a9
 
 RestoreWorkspace: No
 SaveWorkspace: No

tests/testthat/test-r2_mcfadden.R

@@ -26,3 +26,45 @@ test_that("r2_mcfadden", {
     }
   )
 })
+
+skip_if_not_installed("withr")
+
+withr::with_environment(
+  new.env(),
+  {
+    test_that("r2_mcfadden, glmmTMB-beta-binomial", {
+      skip_if_not_installed("glmmTMB")
+      set.seed(101)
+      dd <- data.frame(x = rnorm(200))
+      dd$y <- glmmTMB::simulate_new(
+        ~ 1 + x,
+        newdata = dd,
+        newparams = list(beta = c(0,1), betadisp = -1),
+        weights = rep(10, nrow(dd)),
+        family = glmmTMB::betabinomial()
+      )[[1]]
+      dd$success <- round(runif(nrow(dd), 0, dd$y))
+      d <<- dd
+
+      m <- glmmTMB::glmmTMB(
+        y/10 ~ 1 + x,
+        data = d,
+        weights = rep(10, nrow(d)),
+        family = glmmTMB::betabinomial()
+      )
+      out1 <- r2(m)
+      out2 <- r2_mcfadden(m)
+      expect_equal(out1$R2, out2$R2, tolerance = 1e-4, ignore_attr = TRUE)
+      expect_equal(out1$R2, 0.06892733, tolerance = 1e-4, ignore_attr = TRUE)
+
+      m <- glmmTMB::glmmTMB(
+        cbind(y, success) ~ 1 + x,
+        data = d,
+        weights = rep(10, nrow(d)),
+        family = glmmTMB::betabinomial()
+      )
+      expect_warning(r2(m), regex = "calculate accurate")
+      expect_warning(r2_mcfadden(m), regex = "calculate accurate")
+    })
+  }
+)