From b2f8fd5586fff583e1615032eafdbd73060f36ee Mon Sep 17 00:00:00 2001
From: Daniel <mail@danielluedecke.de>
Date: Sun, 22 Dec 2024 13:26:17 +0100
Subject: [PATCH] fix, add tests

---
 R/r2.R                            | 13 ++++++++++--
 R/r2_coxsnell.R                   | 11 +++++++++-
 R/r2_mcfadden.R                   | 10 ++++++++-
 tests/testthat/test-r2_mcfadden.R | 35 +++++++++++++++++++++++++++++++
 4 files changed, 65 insertions(+), 4 deletions(-)

diff --git a/R/r2.R b/R/r2.R
index 946e385b3..50f62177e 100644
--- a/R/r2.R
+++ b/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
@@ -527,8 +528,16 @@ r2.glmmTMB <- function(model, ci = NULL, tolerance = 1e-5, verbose = TRUE, ...)
       names(out$R2_Tjur) <- "Tjur's R2"
       class(out) <- c("r2_pseudo", class(out))
     } else if (info$is_betabinomial) {
-      # betabinomial default to mcfadden, see pscl:::pR2Work
-      out <- r2_mcfadden(model)
+      # 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) {
diff --git a/R/r2_coxsnell.R b/R/r2_coxsnell.R
index 826faca4e..5215c283b 100644
--- a/R/r2_coxsnell.R
+++ b/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_betabinomial && !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)
diff --git a/R/r2_mcfadden.R b/R/r2_mcfadden.R
index f3ab0a38c..2da19787c 100644
--- a/R/r2_mcfadden.R
+++ b/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_betabinomial && !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)
diff --git a/tests/testthat/test-r2_mcfadden.R b/tests/testthat/test-r2_mcfadden.R
index ad64e92dd..37f433641 100644
--- a/tests/testthat/test-r2_mcfadden.R
+++ b/tests/testthat/test-r2_mcfadden.R
@@ -26,3 +26,38 @@ test_that("r2_mcfadden", {
     }
   )
 })
+
+
+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 = "betabinomial"
+  )[[1]]
+  dd$success <- round(runif(nrow(dd), 0, dd$y))
+
+  m <- glmmTMB::glmmTMB(
+    y/10 ~ 1 + x,
+    data = dd,
+    weights = rep(10, nrow(dd)),
+    family = "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 = dd,
+    weights = rep(10, nrow(dd)),
+    family = "betabinomial"
+  )
+  expect_warning(r2(m), regex = "calculate accurate")
+  expect_warning(r2_mcfadden(m), regex = "calculate accurate")
+})