Fix zero-inflation parameters in bootstrapped glmmTMB output

DESCRIPTION

@@ -1,7 +1,7 @@
 Type: Package
 Package: parameters
 Title: Processing of Model Parameters
-Version: 0.30.0.12
+Version: 0.30.0.13
 Authors@R:
     c(person(given = "Daniel",
              family = "Lüdecke",

NEWS.md

@@ -48,6 +48,11 @@
   per row. In the second kind, the response in the formula is not one column
   of the data, for example `log(y)` or `cbind()`.
 
+* Fixed issue where `bootstrap_parameters()` and `model_parameters()` with
+  `bootstrap = TRUE` mislabeled the zero-inflation parameters of *glmmTMB*
+  models. `bootstrap_model()` now adds the suffix `_zi` to these parameters,
+  and bootstrapped output no longer includes dispersion parameters (#1284).
+
 * The `p_adjust` argument in `model_parameters()` now accepts the methods from
   `stats::p.adjust.methods` in any letter case, for example `"bh"`. Before,
   this failed with an error. An invalid value for `p_adjust` no longer shows the

R/bootstrap_model.R

@@ -372,7 +372,7 @@
 # Find the grouping factors of a mixed model, the top factor that all of them
 # nest in, and the factors that need new labels for each drawn cluster.
 .bootstrap_cluster_groups <- function(model, data) {
   terms <- unique(insight::find_random(model, split_nested = FALSE, flatten = TRUE))
   columns <- unique(insight::find_random(model, split_nested = TRUE, flatten = TRUE))
 
   not_found <- setdiff(columns, colnames(data))
@@ -385,7 +385,7 @@
     ))
   }
 
   interaction_terms <- strsplit(terms[grepl(":", terms, fixed = TRUE)], ":", fixed = TRUE)
   fixed_predictors <- insight::find_predictors(model, effects = "fixed", flatten = TRUE)
   top <- if (length(columns) == 1L) {
     columns
@@ -429,7 +429,7 @@
 # nest in each of its factors. Other factors are nested if each of their levels
 # occurs within one level of the top factor.
 .bootstrap_top_factor <- function(columns, terms, data, fixed = NULL) {
   interaction_terms <- strsplit(terms[grepl(":", terms, fixed = TRUE)], ":", fixed = TRUE)
 
   nests_in <- function(inner, outer) {
     in_formula <- !inner %in% terms &&
@@ -463,7 +463,7 @@
   # with more than one candidate, a factor that is not a fixed effect, and then
   # the factor with the most clusters
   candidates <- columns[is_top]
   if (any(!candidates %in% fixed)) {
     candidates <- candidates[!candidates %in% fixed]
   }
   n_levels <- vapply(candidates, function(x) length(unique(data[[x]])), integer(1))
@@ -552,7 +552,12 @@
   out <- as.data.frame(results$t)
   out <- out[stats::complete.cases(out), ]
 
-  names(out) <- insight::find_parameters(model, effects = "fixed")$conditional
+  # zero-inflation parameters get the suffix "_zi", as in simulate_model()
+  fixed_names <- insight::find_parameters(model, effects = "fixed")
+  names(out) <- c(
+    fixed_names$conditional,
+    if (length(fixed_names$zero_inflated)) paste0(fixed_names$zero_inflated, "_zi")
+  )
   class(out) <- unique(c("bootstrap_model", "see_bootstrap_model", class(out)))
   attr(out, "original_model") <- model
   out

R/bootstrap_parameters.R

@@ -98,6 +98,7 @@ bootstrap_parameters.bootstrap_model <- function(
     ci_method = ci_method,
     ...
   )
+  out <- .add_component_bootstrap_glmmTMB(out, attr(model, "original_model"))
 
   class(out) <- c("bootstrap_parameters", "parameters_model", class(out))
   attr(out, "boot_samples") <- model
@@ -111,6 +112,29 @@ model_parameters.bootstrap_model <- bootstrap_parameters.bootstrap_model
 
 # utilities --------------------------------------------------------------------
 
+# restore names of zero-inflation parameters ("_zi" suffix) and add component
+.add_component_bootstrap_glmmTMB <- function(params, model) {
+  if (!inherits(model, "glmmTMB")) {
+    return(params)
+  }
+  fixed_names <- insight::find_parameters(model, effects = "fixed")
+  if (!length(fixed_names$zero_inflated)) {
+    return(params)
+  }
+  draw_names <- c(fixed_names$conditional, paste0(fixed_names$zero_inflated, "_zi"))
+  # user-added columns (e.g. linear combinations) keep their names
+  idx <- match(params$Parameter, draw_names)
+  found <- !is.na(idx)
+  fixed_parameters <- c(fixed_names$conditional, fixed_names$zero_inflated)
+  params$Parameter[found] <- fixed_parameters[idx[found]]
+  params$Component <- "conditional"
+  params$Component[found] <- rep(
+    c("conditional", "zero_inflated"),
+    c(length(fixed_names$conditional), length(fixed_names$zero_inflated))
+  )[idx[found]]
+  params
+}
+
 #' @keywords internal
 .summary_bootstrap <- function(data, test, centrality, ci, ci_method, ...) {
   # Is the p-value requested?

R/methods_glmmTMB.R

@@ -132,6 +132,9 @@
 #' For models where the dispersion parameter and the residual variance are
 #' the same, only the residual variance is shown in the output.
 #'
+#' With `bootstrap = TRUE`, dispersion parameters are not included in the
+#' output.
+#'
 #' @seealso [insight::standardize_names()] to
 #'   rename columns into a consistent, standardized naming scheme.
 #'
@@ -262,8 +265,21 @@
 
   if (effects %in% c("fixed", "all")) {
     if (bootstrap) {
+      if (component == "dispersion") {
+        insight::format_error(
+          "Bootstrapping does not return dispersion parameters. Use `component = \"all\"`, `\"conditional\"` or `\"zi\"`." # nolint
+        )
+      }
       # bootstrapped parameters
       params <- bootstrap_parameters(model, iterations = iterations, ci = ci, ...)
+      # keep requested component, and drop the component column for
+      # the conditional component, as for the regular parameters
+      if (!is.null(params$Component)) {
+        params <- .filter_component(params, component)
+        if (component == "conditional") {
+          params$Component <- NULL
+        }
+      }
       if (effects != "fixed") {
         effects <- "fixed"
         if (verbose) {
@@ -297,10 +313,12 @@
       params <- do.call(".extract_parameters_generic", fun_args)
     }
 
-    # add dispersion parameter
-    out <- .add_dispersion_param_glmmTMB(model, params, effects, component, ci, verbose)
-    params <- out$params
-    dispersion_param <- out$dispersion_param
+    # add dispersion parameter (not for bootstrapped draws, which have none)
+    if (!bootstrap) {
+      out <- .add_dispersion_param_glmmTMB(model, params, effects, component, ci, verbose)
+      params <- out$params
+      dispersion_param <- out$dispersion_param
+    }
 
     # exponentiate coefficients and SE/CI, if requested
     params <- .exponentiate_parameters(params, model, exponentiate)
@@ -344,11 +362,18 @@
     attributes(params) <- utils::modifyList(att, attributes(params))
   }
 
+  # bootstrapped intervals are quantiles of the draws
+  if (bootstrap) {
+    ci_method <- "quantile"
+  }
+
   params <- .add_model_parameters_attributes(
     params,
     model,
     ci = ci,
     exponentiate,
+    bootstrap = bootstrap,
+    iterations = iterations,
     ci_method = ci_method,
     p_adjust = p_adjust,
     verbose = verbose,
@@ -597,7 +622,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)
@@ -666,7 +691,7 @@
 .se_robust_glmmTMB <- function(
   model,
   component = "all",
   vcov,
   vcov_args = NULL,
   verbose = TRUE,
   ...

tests/testthat/test-bootstrap_parameters.R

@@ -33,7 +33,7 @@
   set.seed(123)
   out <- bootstrap_model(m, iterations = 20, type = "semiparametric")
   expect_named(out, c("(Intercept)", "Days"))
   expect_identical(nrow(out), 20L)
   set.seed(123)
   out <- bootstrap_parameters(m, iterations = 20, type = "semiparametric")
   expect_identical(out$Parameter, c("(Intercept)", "Days"))
@@ -96,7 +96,7 @@
     set.seed(7)
     draws <- boot::boot(seq_along(clusters), function(x, i) i, R = 20)$t
     expected <- apply(draws, 1, function(i) mean(cluster_means[clusters[i]]))
     expect_identical(nrow(out), 20L)
     expect_equal(out[["(Intercept)"]], unname(expected), tolerance = 1e-8)
   }
 })
@@ -171,7 +171,7 @@
   m <- nlme::lme(y ~ cond, random = ~ 1 | subject / cond, data = d)
   .expect_cluster_bootstrap(m, d, "subject")
   set.seed(42)
   expect_identical(nrow(bootstrap_model(m, iterations = 20)), 20L)
 })
 
 test_that("bootstrap_model errors for relabelled grouping factors that are fixed effects", {
@@ -403,3 +403,68 @@
     )
   )
 })
+
+test_that("bootstrap_parameters returns zero-inflation components of glmmTMB models", {
+  skip_if_not_installed("glmmTMB")
+  data(Salamanders, package = "glmmTMB")
+  m <- glmmTMB::glmmTMB(
+    count ~ mined + (1 | site),
+    ziformula = ~mined,
+    family = poisson,
+    data = Salamanders
+  )
+  set.seed(123)
+  b <- bootstrap_model(m, iterations = 5)
+  expect_named(b, c("(Intercept)", "minedno", "(Intercept)_zi", "minedno_zi"))
+  b$lin_comb <- b$minedno - b$minedno_zi
+  out <- bootstrap_parameters(b)
+  expect_identical(
+    out$Parameter,
+    c("(Intercept)", "minedno", "(Intercept)", "minedno", "lin_comb")
+  )
+  expect_identical(
+    out$Component,
+    c("conditional", "conditional", "zero_inflated", "zero_inflated", "conditional")
+  )
+  estimates <- insight::get_parameters(m)$Estimate
+  expect_lt(max(abs(out$Coefficient[1:4] - estimates)), 0.5)
+})
+
+test_that("model_parameters with bootstrap uses component for glmmTMB", {
+  skip_if_not_installed("glmmTMB")
+  data(Salamanders, package = "glmmTMB")
+  m <- glmmTMB::glmmTMB(
+    count ~ mined + (1 | site),
+    ziformula = ~mined,
+    dispformula = ~mined,
+    family = glmmTMB::nbinom2,
+    data = Salamanders
+  )
+  for (comp in c("all", "conditional", "zi")) {
+    expected <- model_parameters(m, component = comp, effects = "fixed")
+    if (comp == "all") {
+      expected <- expected[expected$Component != "dispersion", ]
+    }
+    set.seed(123)
+    out <- model_parameters(
+      m,
+      bootstrap = TRUE,
+      iterations = 5,
+      component = comp,
+      effects = "fixed"
+    )
+    expect_identical(out$Parameter, expected$Parameter)
+    expect_identical(out$Component, expected$Component)
+  }
+  printed <- evaluate_promise(print(out))
+  expect_match(paste(printed$messages, collapse = " "), "bootstrap", fixed = TRUE)
+  expect_error(
+    model_parameters(m, bootstrap = TRUE, iterations = 5, component = "dispersion"),
+    "Bootstrapping\\s+does\\s+not\\s+return\\s+dispersion\\s+parameters"
+  )
+
+  m <- glmmTMB::glmmTMB(count ~ mined + (1 | site), family = poisson, data = Salamanders)
+  set.seed(123)
+  out <- model_parameters(m, bootstrap = TRUE, iterations = 5, effects = "fixed")
+  expect_identical(out$Parameter, c("(Intercept)", "minedno"))
+})