diff --git a/NAMESPACE b/NAMESPACE index e573d5ce0..0bb599f64 100644 --- 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 index 943ca2700..85fbd03ef 100644 --- 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 index cdb73dea7..d56803f5c 100644 --- 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 index 85bcc6e8f..e97e8fb6f 100644 --- 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 index 2d52f6b9a..069ee9ad2 100644 --- 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) })