Detect custom MixMod families in model_info() (#901)

* Detect custom MixMod families in `model_info()`

* update readme

* remove badge

* fix

* fix?

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 = "[email protected]",
              comment = c(ORCID = "0000-0001-9560-6336", Twitter = "@bmwiernik")),
-      person(given = "Vincent", 
-             family = "Arel-Bundock", 
-             email = "[email protected]", 
+      person(given = "Vincent",
+             family = "Arel-Bundock",
+             email = "[email protected]",
              role = c("aut", "ctb"),
              comment = c(ORCID = "0000-0003-2042-7063")),
-      person(given = "Etienne", 
+      person(given = "Etienne",
              family = "Bacher",
-             email = "[email protected]", 
+             email = "[email protected]",
              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

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

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,

README.Rmd

@@ -17,7 +17,7 @@ options(knitr.kable.NA = "", digits = 2, width = 80)
 # insight <img src='man/figures/logo.png' align="right" height="139" />
 
 [![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!**
 

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)

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)
+})

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)
   })