From ded5c5d798450cbc38908ae8c6efaf6a6ad8e432 Mon Sep 17 00:00:00 2001 From: Daniel Date: Thu, 23 Nov 2023 11:25:52 +0100 Subject: [PATCH] add more r2 to non-mixed glmmTMB --- R/r2.R | 15 ++++++++++++ R/r2_nagelkerke.R | 3 ++- R/r2_zeroinflated.R | 1 - tests/testthat/test-r2.R | 49 ++++++++++++++++++++++++++++++++++++++-- 4 files changed, 64 insertions(+), 4 deletions(-) diff --git a/R/r2.R b/R/r2.R index f1399ce0d..dacdef614 100644 --- a/R/r2.R +++ b/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 } diff --git a/R/r2_nagelkerke.R b/R/r2_nagelkerke.R index 87b8cb3ab..00dfe36d2 100644 --- a/R/r2_nagelkerke.R +++ b/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)) { diff --git a/R/r2_zeroinflated.R b/R/r2_zeroinflated.R index b1e057bb6..23d859587 100644 --- a/R/r2_zeroinflated.R +++ b/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) diff --git a/tests/testthat/test-r2.R b/tests/testthat/test-r2.R index 069ee9ad2..2f450ccc9 100644 --- a/tests/testthat/test-r2.R +++ b/tests/testthat/test-r2.R @@ -43,6 +43,7 @@ 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) @@ -50,7 +51,7 @@ 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 + # binomial ------------------------------------------------------------- data(mtcars) m <- glmmTMB::glmmTMB(am ~ mpg, data = mtcars, family = binomial()) out <- r2(m) @@ -58,5 +59,49 @@ test_that("r2 glmmTMB, no ranef", { # 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) })