Fix standard errors for glmmTMB ordinal() models (needs easystats/insight#1221)

DESCRIPTION

@@ -1,7 +1,7 @@
 Type: Package
 Package: parameters
 Title: Processing of Model Parameters
-Version: 0.29.3
+Version: 0.29.3.1
 Authors@R:
     c(person(given = "Daniel",
              family = "Lüdecke",
@@ -85,7 +85,7 @@ Depends:
 Imports:
     bayestestR (>= 0.18.1),
     datawizard (>= 1.3.1),
-    insight (>= 1.5.3),
+    insight (>= 1.5.4),
     graphics,
     methods,
     stats,
@@ -240,3 +240,4 @@ Config/Needs/website: easystats/easystatstemplate
 Config/Needs/check: stan-dev/cmdstanr
 Config/rcmdcheck/ignore-inconsequential-notes: true
 Config/roxygen2/version: 8.0.0
+Remotes: easystats/insight, glmmTMB/glmmTMB/glmmTMB

NEWS.md

@@ -1,3 +1,14 @@
+# parameters (devel)
+
+## Bug fixes
+
+* `standard_error()` and hence `model_parameters()` returned wrong (recycled)
+  standard errors for models of class `glmmTMB` fitted with the new `ordinal()`
+  family, because the thresholds are not part of the summary coefficient table.
+  Standard errors are now taken from `insight::get_varcov()`, which returns
+  them on the threshold scale via the delta method (requires *insight* > 1.5.4
+  and *glmmTMB* >= 1.1.15).
+
 # parameters 0.29.3
 
 ## Changes

R/methods_glmmTMB.R

@@ -596,7 +596,7 @@
 
   if (effects == "random") {
     .se_random_effects_glmmTMB(model)
   } else if (!is.null(vcov)) {
     .se_robust_glmmTMB(model, component, vcov, vcov_args, verbose, ...)
   } else {
     .se_fixed_effects_glmmTMB(model, component, method, verbose)
@@ -617,6 +617,27 @@
     return(se_kenward(model, component = "conditional"))
   }
 
+  # ordinal family: the thresholds are not part of the summary coefficient
+  # table, and their standard errors require the delta method (glmmTMB
+  # estimates them on an internal softmax scale) - use insight's covariance
+  # matrix, which is on the threshold scale and aligned with the parameters
+  if (identical(insight::get_family(model)$family, "ordinal")) {
+    params <- insight::find_parameters(
+      model,
+      effects = "fixed",
+      component = "conditional",
+      flatten = TRUE
+    )
+    vc <- insight::get_varcov(model, component = "conditional", verbose = FALSE)
+    se_vec <- sqrt(diag(vc))
+    se <- .data_frame(
+      Parameter = params,
+      SE = unname(se_vec[params]),
+      Component = "conditional"
+    )
+    return(.filter_component(se, component))
+  }
+
   cs <- suppressWarnings(insight::compact_list(stats::coef(summary(model))))
   x <- lapply(names(cs), function(i) {
     .data_frame(
@@ -644,7 +665,7 @@
 .se_robust_glmmTMB <- function(
   model,
   component = "all",
   vcov,
   vcov_args = NULL,
   verbose = TRUE,
   ...

tests/testthat/test-glmmTMB-ordinal.R

@@ -0,0 +1,46 @@
+skip_on_cran()
+skip_if_not_installed("glmmTMB", minimum_version = "1.1.15")
+skip_if_not_installed("insight", minimum_version = "1.5.4.1")
+skip_if_not_installed("ordinal")
+skip_if_not_installed("withr")
+
+data(wine, package = "ordinal")
+
+withr::with_options(list(parameters_exponentiate = FALSE), {
+  m_tmb <- glmmTMB::glmmTMB(
+    rating ~ temp + contact + (1 | judge),
+    data = wine,
+    family = glmmTMB::ordinal()
+  )
+  m_clmm <- ordinal::clmm(rating ~ temp + contact + (1 | judge), data = wine)
+
+  test_that("standard_error, glmmTMB ordinal", {
+    out <- standard_error(m_tmb)
+    expect_identical(
+      out$Parameter,
+      c("1|2", "2|3", "3|4", "4|5", "tempwarm", "contactyes")
+    )
+    expect_equal(
+      out$SE,
+      unname(sqrt(diag(vcov(m_clmm)))[out$Parameter]),
+      tolerance = 1e-3
+    )
+    expect_identical(unique(out$Component), "conditional")
+    expect_message(
+      expect_null(standard_error(m_tmb, component = "zi")),
+      "no zero-inflation"
+    )
+  })
+
+  test_that("model_parameters, glmmTMB ordinal matches clmm", {
+    out <- model_parameters(m_tmb)
+    ref <- model_parameters(m_clmm)
+    expect_identical(out$Parameter, ref$Parameter)
+    for (col in c("Coefficient", "SE", "CI_low", "CI_high", "z", "p")) {
+      expect_equal(out[[col]], ref[[col]], tolerance = 1e-2, ignore_attr = TRUE)
+    }
+    # random effects are still reported
+    out <- model_parameters(m_tmb, effects = "all")
+    expect_true("SD (Intercept)" %in% out$Parameter)
+  })
+})