equivalence_test() bug when applied to glmmTMB models using the beta family (#1159)

* equivalence_test() bug when applied to glmmTMB models using the beta family
Fixes #1137

* add test

* version

* news

* fixes

* fix

* fix

* fix

* fix

* fix test

* trigger CI

Author: strengejacke

DESCRIPTION

@@ -1,7 +1,7 @@
 Type: Package
 Package: parameters
 Title: Processing of Model Parameters
-Version: 0.28.0.10
+Version: 0.28.0.13
 Authors@R:
     c(person(given = "Daniel",
              family = "Lüdecke",
@@ -231,4 +231,4 @@ Config/testthat/parallel: true
 Config/Needs/website: easystats/easystatstemplate
 Config/Needs/check: stan-dev/cmdstanr
 Config/rcmdcheck/ignore-inconsequential-notes: true
-Remotes: easystats/insight, easystats/bayestestR, easystats/modelbased, vincentarelbundock/marginaleffects
+Remotes: easystats/insight, easystats/bayestestR, easystats/modelbased#557, vincentarelbundock/marginaleffects

NEWS.md

@@ -8,6 +8,11 @@
 * `model_parameters()` for *marginaleffects* objects is now more robust in
   detecting Bayesian models.
 
+## Bug fixes
+
+* Fixed issue with `equivalence_test()` for models of class `glmmTMB` with
+  `beta_family()`.
+
 # parameters 0.28.0
 
 ## Breaking Changes

R/equivalence_test.R

@@ -736,6 +736,11 @@ equivalence_test.parameters_model <- function(x,
     {
       params <- insight::get_parameters(model)
 
+      # remove dispersion components
+      if (!is.null(params$Component)) {
+        params <- params[params$Component != "dispersion", ]
+      }
+
       # degrees of freedom
       dof <- insight::get_df(x = model, type = "wald")
 
@@ -745,6 +750,11 @@ equivalence_test.parameters_model <- function(x,
       # se
       se <- standard_error(model, vcov = vcov, vcov_args = vcov_args, ...)
 
+      # remove dispersion components
+      if (!is.null(se$Component)) {
+        se <- se[se$Component != "dispersion", ]
+      }
+
       stats::pt((range[1] - params$mu) / se$SE, df = dof, lower.tail = TRUE) +
         stats::pt((range[2] - params$mu) / se$SE, df = dof, lower.tail = FALSE)
     },

R/methods_marginaleffects.R

@@ -9,23 +9,22 @@ model_parameters.marginaleffects <- function(model,
                                              exponentiate = FALSE,
                                              verbose = TRUE,
                                              ...) {
-  insight::check_if_installed("marginaleffects", minimum_version = "0.28.0.21")
+  insight::check_if_installed("marginaleffects", minimum_version = "0.28.0.22")
 
   # Bayesian models have posterior draws as attribute
   is_bayesian <- !is.null(suppressWarnings(marginaleffects::get_draws(model, "PxD")))
 
   if (is_bayesian) {
     # Bayesian
-    tidy_model <- suppressWarnings(bayestestR::describe_posterior(
+    out <- suppressWarnings(bayestestR::describe_posterior(
       model,
       ci = ci,
       verbose = verbose,
       ...
     ))
   } else {
-    # handle non-Bayesian models
-    tidy_model <- marginaleffects::tidy(model, conf_level = ci, ...)
-
+    # non-Bayesian
+    out <- as.data.frame(model)
     # all columns in data grid and model data, we only want to keep "by" variables
     all_data_cols <- union(
       colnames(marginaleffects::components(model, "newdata")),
@@ -36,16 +35,12 @@ model_parameters.marginaleffects <- function(model,
       marginaleffects::components(model, "variable_names_by"),
       marginaleffects::components(model, "variable_names_by_hypothesis")
     )
-
-    ## FIXME: hack to workaround https://github.com/vincentarelbundock/marginaleffects/issues/1573
-    tidy_model <- .fix_duplicated_by_columns(tidy_model, by_cols)
-
     # remove redundant columns
     to_remove <- setdiff(all_data_cols, by_cols)
-    tidy_model <- tidy_model[, !colnames(tidy_model) %in% to_remove, drop = FALSE]
+    out <- out[, !colnames(out) %in% to_remove, drop = FALSE]
   }
 
-  out <- .rename_reserved_marginaleffects(tidy_model)
+  out <- .rename_reserved_marginaleffects(out)
 
   # need to standardize names for non-Bayesian models. Bayesian models have
   # been processed through describe_posterior() already
@@ -118,7 +113,7 @@ model_parameters.predictions <- function(model,
                                          exponentiate = FALSE,
                                          verbose = TRUE,
                                          ...) {
-  insight::check_if_installed("marginaleffects", minimum_version = "0.28.0.21")
+  insight::check_if_installed("marginaleffects", minimum_version = "0.28.0.22")
 
   # Bayesian models have posterior draws as attribute
   is_bayesian <- !is.null(suppressWarnings(marginaleffects::get_draws(model, "PxD")))
@@ -137,10 +132,6 @@ model_parameters.predictions <- function(model,
       marginaleffects::components(model, "variable_names_by"),
       marginaleffects::components(model, "variable_names_by_hypothesis")
     )
-
-    ## FIXME: hack to workaround https://github.com/vincentarelbundock/marginaleffects/issues/1573
-    model <- .fix_duplicated_by_columns(model, by_cols)
-
     # handle non-Bayesian models
     out <- .rename_reserved_marginaleffects(model)
     out <- datawizard::data_rename(out, "estimate", "predicted")

tests/testthat/test-equivalence_test.R

@@ -18,6 +18,36 @@ test_that("equivalence_test, robust", {
   expect_snapshot(print(x))
 })
 
+test_that("equivalence_test, robust", {
+  skip_on_cran()
+  skip_if_not_installed("glmmTMB")
+
+  set.seed(123)
+  sim_data <- data.frame(
+    x = rbeta(100, 2, 5),
+    y = sample(0:1, 100, replace = TRUE),
+    z = rnorm(100),
+    dataset = rep(1:10, each = 10)
+  )
+
+  mod <- glmmTMB::glmmTMB(
+    x ~ y * z + (1 | dataset),
+    data = sim_data,
+    family = glmmTMB::beta_family(link = "logit")
+  )
+  out <- equivalence_test(mod)
+  expect_identical(dim(out), c(4L, 10L))
+  expect_identical(
+    out$ROPE_Equivalence,
+    c("Rejected", "Undecided", "Accepted", "Undecided")
+  )
+  expect_equal(
+    out$SGPV,
+    c(0, 0.8726, 0.9739, 0.6741),
+    tolerance = 1e-3
+  )
+})
+
 test_that("equivalence_test, unequal rope-range", {
   data(iris)
   m <- lm(Sepal.Length ~ Species, data = iris)

tests/testthat/test-marginaleffects.R

@@ -1,4 +1,4 @@
-skip_if_not_installed("marginaleffects", minimum_version = "0.25.0")
+skip_if_not_installed("marginaleffects", minimum_version = "0.28.0.22")
 skip_if_not_installed("rstanarm")
 
 test_that("marginaleffects()", {
@@ -48,12 +48,12 @@ test_that("marginaleffects()", {
   model <- mgcv::gam(Sepal.Width ~ s(Petal.Length, by = Species), data = iris)
   mfx <- marginaleffects::avg_slopes(model, variables = "Petal.Length")
   out <- model_parameters(mfx)
-  expect_identical(dim(out), c(1L, 11L))
+  expect_identical(dim(out), c(1L, 10L))
   expect_named(
     out,
     c(
       "Parameter", "Comparison", "Coefficient", "SE", "Statistic",
-      "p", "S", "CI", "CI_low", "CI_high", "Predicted"
+      "p", "S", "CI", "CI_low", "CI_high"
     )
   )
   mfx <- marginaleffects::avg_slopes(model, variables = "Petal.Length", by = "Species")
@@ -181,7 +181,7 @@ test_that("predictions, bmrs with special response formula", {
   skip_if_offline()
   skip_if_not_installed("httr2")
   skip_if_not_installed("brms")
-  skip_if_not_installed("marginaleffects", minimum_version = "0.28.0.21")
+  skip_if_not_installed("marginaleffects", minimum_version = "0.28.0.22")
 
   m <- insight::download_model("brms_ipw_1")
   skip_if(is.null(m))
@@ -193,14 +193,13 @@ test_that("predictions, bmrs with special response formula", {
 
 
 test_that("modelbased, tidiers work", {
-  skip_if_not_installed("marginaleffects", minimum_version = "0.28.0.21")
+  skip_if_not_installed("marginaleffects", minimum_version = "0.28.0.22")
   skip_if_not_installed("modelbased", minimum_version = "0.12.0.17")
   skip_if(getRversion() < "4.5.0")
 
   data(penguins)
   m <- lm(bill_len ~ island * sex + bill_dep + species, data = penguins)
 
-  ## FIXME: Need to wait for https://github.com/vincentarelbundock/marginaleffects/issues/1573
   out <- modelbased::estimate_contrasts(m, "island", by = "sex", comparison = ratio ~ pairwise)
   expect_named(
     out,
@@ -266,7 +265,7 @@ test_that("predictions, using bayestestR #1063", {
   skip_if_offline()
   skip_if_not_installed("httr2")
   skip_if_not_installed("brms")
-  skip_if_not_installed("marginaleffects", minimum_version = "0.28.0.21")
+  skip_if_not_installed("marginaleffects", minimum_version = "0.28.0.22")
 
   m <- insight::download_model("brms_mixed_3")
   skip_if(is.null(m))