From afe04c1a2409abb007b17706b846b29b24fd3751 Mon Sep 17 00:00:00 2001 From: Daniel Date: Fri, 12 Jul 2024 19:07:41 +0200 Subject: [PATCH] Detect custom MixMod families in `model_info()` (#901) * Detect custom MixMod families in `model_info()` * update readme * remove badge * fix * fix? --- DESCRIPTION | 22 ++-- NEWS.md | 3 + R/model_info.R | 1 + README.Rmd | 2 +- README.md | 8 -- tests/testthat/test-GLMMadaptive.R | 110 ++++++++++++++++++ .../testthat/test-r2_nakagawa_ordered_beta.R | 7 +- 7 files changed, 132 insertions(+), 21 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index a9e04bf8c2..70914217b9 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,8 +1,8 @@ Type: Package Package: insight Title: Easy Access to Model Information for Various Model Objects -Version: 0.20.1.13 -Authors@R: +Version: 0.20.1.14 +Authors@R: c(person(given = "Daniel", family = "Lüdecke", role = c("aut", "cre"), @@ -33,14 +33,14 @@ Authors@R: role = c("aut", "ctb"), email = "brenton@wiernik.org", comment = c(ORCID = "0000-0001-9560-6336", Twitter = "@bmwiernik")), - person(given = "Vincent", - family = "Arel-Bundock", - email = "vincent.arel-bundock@umontreal.ca", + person(given = "Vincent", + family = "Arel-Bundock", + email = "vincent.arel-bundock@umontreal.ca", role = c("aut", "ctb"), comment = c(ORCID = "0000-0003-2042-7063")), - person(given = "Etienne", + person(given = "Etienne", family = "Bacher", - email = "etienne.bacher@protonmail.com", + email = "etienne.bacher@protonmail.com", role = c("aut", "ctb"), comment = c(ORCID = "0000-0002-9271-5075")), person(given = "Alex", @@ -76,13 +76,13 @@ Description: A tool to provide an easy, intuitive and consistent License: GPL-3 URL: https://easystats.github.io/insight/ BugReports: https://github.com/easystats/insight/issues -Depends: +Depends: R (>= 3.6) -Imports: +Imports: methods, stats, utils -Suggests: +Suggests: AER, afex, aod, @@ -203,7 +203,7 @@ Suggests: tweedie, VGAM, withr -VignetteBuilder: +VignetteBuilder: knitr Encoding: UTF-8 Language: en-US diff --git a/NEWS.md b/NEWS.md index a6e5f89097..3d3580f809 100644 --- a/NEWS.md +++ b/NEWS.md @@ -37,6 +37,9 @@ * `get_transformation()` now returns more transformations for power-transformed response variables. +* `model_info()` for `MixMod` objects from package *GLMMadaptive* now recognize + zero-inflation and hurdle models for custom families. + ## Bug fixes * `null_model()` now correctly handles zero-inflated models from package diff --git a/R/model_info.R b/R/model_info.R index f6aea308d7..61769a4cd6 100644 --- a/R/model_info.R +++ b/R/model_info.R @@ -515,6 +515,7 @@ model_info.MixMod <- function(x, verbose = TRUE, ...) { .make_family( x = x, fitfam = faminfo$family, + zero.inf = !is.null(stats::formula(x, type = "zi_fixed")), logit.link = faminfo$link == "logit", link.fun = faminfo$link, verbose = verbose, diff --git a/README.Rmd b/README.Rmd index b6cd247d0a..3433c08231 100644 --- a/README.Rmd +++ b/README.Rmd @@ -17,7 +17,7 @@ options(knitr.kable.NA = "", digits = 2, width = 80) # insight [![DOI](https://joss.theoj.org/papers/10.21105/joss.01412/status.svg)](https://doi.org/10.21105/joss.01412) -[![downloads](https://cranlogs.r-pkg.org/badges/insight)](https://cranlogs.r-pkg.org/) [![total](https://cranlogs.r-pkg.org/badges/grand-total/insight)](https://cranlogs.r-pkg.org/) [![status](https://tinyverse.netlify.com/badge/insight)](https://CRAN.R-project.org/package=insight) [![lifecycle](https://img.shields.io/badge/lifecycle-maturing-blue.svg)](https://lifecycle.r-lib.org/articles/stages.html) +[![downloads](https://cranlogs.r-pkg.org/badges/insight)](https://cranlogs.r-pkg.org/) [![total](https://cranlogs.r-pkg.org/badges/grand-total/insight)](https://cranlogs.r-pkg.org/) **Gain insight into your models!** diff --git a/README.md b/README.md index 1714160f21..f520cb6835 100644 --- a/README.md +++ b/README.md @@ -4,8 +4,6 @@ [![DOI](https://joss.theoj.org/papers/10.21105/joss.01412/status.svg)](https://doi.org/10.21105/joss.01412) [![downloads](https://cranlogs.r-pkg.org/badges/insight)](https://cranlogs.r-pkg.org/) [![total](https://cranlogs.r-pkg.org/badges/grand-total/insight)](https://cranlogs.r-pkg.org/) -[![status](https://tinyverse.netlify.com/badge/insight)](https://CRAN.R-project.org/package=insight) -[![lifecycle](https://img.shields.io/badge/lifecycle-maturing-blue.svg)](https://lifecycle.r-lib.org/articles/stages.html) **Gain insight into your models!** @@ -242,9 +240,6 @@ print_params <- function(model) { m1 <- lm(Sepal.Length ~ Petal.Width, data = iris) print_params(m1) #> [1] "My parameters are (Intercept), Petal.Width, thank you for your attention!" -``` - -``` r # obviously, something is missing in the output m2 <- mgcv::gam(Sepal.Length ~ Petal.Width + s(Petal.Length), data = iris) @@ -270,9 +265,6 @@ print_params <- function(model) { m1 <- lm(Sepal.Length ~ Petal.Width, data = iris) print_params(m1) #> [1] "My parameters are (Intercept), Petal.Width, thank you for your attention!" -``` - -``` r m2 <- mgcv::gam(Sepal.Length ~ Petal.Width + s(Petal.Length), data = iris) print_params(m2) diff --git a/tests/testthat/test-GLMMadaptive.R b/tests/testthat/test-GLMMadaptive.R index 4e57ebb41e..cdf95fac4f 100644 --- a/tests/testthat/test-GLMMadaptive.R +++ b/tests/testthat/test-GLMMadaptive.R @@ -343,3 +343,113 @@ test_that("find_algorithm", { list(algorithm = "quasi-Newton", optimizer = "optim") ) }) + +test_that("detect custom families", { + skip_on_cran() + set.seed(1234) + n <- 100 # number of subjects + K <- 8 # number of measurements per subject + t_max <- 5 # maximum follow-up time + + # we construct a data frame with the design: + # everyone has a baseline measurement, and then measurements at random follow-up times + DF <- data.frame( + id = rep(seq_len(n), each = K), + time = c(replicate(n, c(0, sort(runif(K - 1, 0, t_max))))), + sex = rep(gl(2, n / 2, labels = c("male", "female")), each = K) + ) + + # design matrices for the fixed and random effects non-zero part + X <- model.matrix(~ sex * time, data = DF) + Z <- model.matrix(~time, data = DF) + # design matrices for the fixed and random effects zero part + X_zi <- model.matrix(~sex, data = DF) + Z_zi <- model.matrix(~1, data = DF) + + betas <- c(-2.13, -0.25, 0.24, -0.05) # fixed effects coefficients non-zero part + sigma <- 0.5 # standard deviation error terms non-zero part + gammas <- c(-1.5, 0.5) # fixed effects coefficients zero part + D11 <- 0.5 # variance of random intercepts non-zero part + D22 <- 0.1 # variance of random slopes non-zero part + D33 <- 0.4 # variance of random intercepts zero part + + # we simulate random effects + b <- cbind(rnorm(n, sd = sqrt(D11)), rnorm(n, sd = sqrt(D22)), rnorm(n, sd = sqrt(D33))) + # linear predictor non-zero part + eta_y <- as.vector(X %*% betas + rowSums(Z * b[DF$id, 1:2, drop = FALSE])) + # linear predictor zero part + eta_zi <- as.vector(X_zi %*% gammas + rowSums(Z_zi * b[DF$id, 3, drop = FALSE])) + # we simulate log-normal longitudinal data + DF$y <- exp(rnorm(n * K, mean = eta_y, sd = sigma)) + # we set the zeros from the logistic regression + DF$y[as.logical(rbinom(n * K, size = 1, prob = plogis(eta_zi)))] <- 0 + + hurdle.lognormal <- function() { + stats <- make.link("identity") + log_dens <- function(y, eta, mu_fun, phis, eta_zi) { + sigma <- exp(phis) + # binary indicator for y > 0 + ind <- y > 0 + # non-zero part + eta <- as.matrix(eta) + eta_zi <- as.matrix(eta_zi) + out <- eta + out[ind, ] <- plogis(eta_zi[ind, ], lower.tail = FALSE, log.p = TRUE) + + dnorm(x = log(y[ind]), mean = eta[ind, ], sd = sigma, log = TRUE) + # zero part + out[!ind, ] <- plogis(eta_zi[!ind, ], log.p = TRUE) + attr(out, "mu_y") <- eta + out + } + score_eta_fun <- function(y, mu, phis, eta_zi) { + sigma <- exp(phis) + # binary indicator for y > 0 + ind <- y > 0 + # non-zero part + eta <- as.matrix(mu) + out <- eta + out[!ind, ] <- 0 + out[ind, ] <- (log(y[ind]) - eta[ind, ]) / sigma^2 + out + } + score_eta_zi_fun <- function(y, mu, phis, eta_zi) { + ind <- y > 0 + probs <- plogis(as.matrix(eta_zi)) + out <- 1 - probs + out[ind, ] <- -probs[ind, ] + out + } + score_phis_fun <- function(y, mu, phis, eta_zi) { + sigma <- exp(phis) + # binary indicator for y > 0 + ind <- y > 0 + # non-zero part + eta <- as.matrix(mu) + out <- eta + out[!ind, ] <- 0 + out[ind, ] <- -1 + (log(y[ind]) - eta[ind, ])^2 / sigma^2 + out + } + simulate <- function(n, mu, phis, eta_zi) { + y <- rlnorm(n = n, meanlog = mu, sdlog = exp(phis)) + y[as.logical(rbinom(n, 1, plogis(eta_zi)))] <- 0 + y + } + structure( + list( + family = "two-part log-normal", link = stats$name, + linkfun = stats$linkfun, linkinv = stats$linkinv, log_dens = log_dens, + score_eta_fun = score_eta_fun, score_eta_zi_fun = score_eta_zi_fun, + score_phis_fun = score_phis_fun, simulate = simulate + ), + class = "family" + ) + } + km1 <- GLMMadaptive::mixed_model(y ~ sex * time, + random = ~ 1 | id, data = DF, + family = hurdle.lognormal(), n_phis = 1, + zi_fixed = ~sex + ) + out <- model_info(km1) + expect_true(out$is_zero_inflated) +}) diff --git a/tests/testthat/test-r2_nakagawa_ordered_beta.R b/tests/testthat/test-r2_nakagawa_ordered_beta.R index b5b4948f1c..2b3cac3ff9 100644 --- a/tests/testthat/test-r2_nakagawa_ordered_beta.R +++ b/tests/testthat/test-r2_nakagawa_ordered_beta.R @@ -22,7 +22,12 @@ withr::with_environment( data = sleepstudy, family = glmmTMB::ordbeta() ) - out <- suppressWarnings(performance::r2_nakagawa(m, verbose = FALSE)) + mnull <- glmmTMB::glmmTMB( + y ~ 1 + (Days | Subject), + data = sleepstudy, + family = glmmTMB::ordbeta() + ) + out <- suppressWarnings(performance::r2_nakagawa(m, null_model = mnull, verbose = FALSE)) expect_equal(out$R2_marginal, 0.2799715, ignore_attr = TRUE, tolerance = 1e-4) expect_equal(out$R2_conditional, 0.8504158, ignore_attr = TRUE, tolerance = 1e-4) })