Better support for models of class sdmTMB (#1075)

* Better support for models of class `sdmTMB`

* supported list

* supported model list

* readme

* more support

* fix

* fix

* find_predictors

* find_predictors

* model info

* more

* disable old tests

* get_predicted for sdmTMB

* docs

* remove unused code

* add to suggest

* styler

* fix

* fix

* new range option

* add getvarcov

* tests

* fix

* remove commented code

* skip on ubuntu

Author: strengejacke

DESCRIPTION

@@ -195,6 +195,7 @@ Suggests:
     rstudioapi,
     RWiener,
     sandwich,
+    sdmTMB,
     serp,
     speedglm,
     splines,

NAMESPACE

@@ -178,6 +178,7 @@ S3method(find_formula,poissonmfx)
 S3method(find_formula,probitmfx)
 S3method(find_formula,rlmerMod)
 S3method(find_formula,rma)
+S3method(find_formula,sdmTMB)
 S3method(find_formula,selection)
 S3method(find_formula,sem)
 S3method(find_formula,stanmvreg)
@@ -304,6 +305,7 @@ S3method(find_parameters,rms)
 S3method(find_parameters,rqs)
 S3method(find_parameters,rqss)
 S3method(find_parameters,scam)
+S3method(find_parameters,sdmTMB)
 S3method(find_parameters,selection)
 S3method(find_parameters,sem)
 S3method(find_parameters,sim)
@@ -329,6 +331,7 @@ S3method(find_predictors,default)
 S3method(find_predictors,fixest)
 S3method(find_predictors,insight_formula)
 S3method(find_predictors,logitr)
+S3method(find_predictors,sdmTMB)
 S3method(find_predictors,selection)
 S3method(find_random,afex_aov)
 S3method(find_random,default)
@@ -577,6 +580,7 @@ S3method(get_dispersion,model_fit)
 S3method(get_family,default)
 S3method(get_family,list)
 S3method(get_family,model_fit)
+S3method(get_family,sdmTMB)
 S3method(get_intercept,default)
 S3method(get_intercept,stanreg)
 S3method(get_loglikelihood,afex_aov)
@@ -743,6 +747,7 @@ S3method(get_parameters,rms)
 S3method(get_parameters,rqs)
 S3method(get_parameters,rqss)
 S3method(get_parameters,scam)
+S3method(get_parameters,sdmTMB)
 S3method(get_parameters,selection)
 S3method(get_parameters,sem)
 S3method(get_parameters,sim)
@@ -794,6 +799,7 @@ S3method(get_predicted,prcomp)
 S3method(get_predicted,principal)
 S3method(get_predicted,rlm)
 S3method(get_predicted,rma)
+S3method(get_predicted,sdmTMB)
 S3method(get_predicted,stanreg)
 S3method(get_predicted,survreg)
 S3method(get_predicted,zeroinfl)
@@ -963,6 +969,7 @@ S3method(get_statistic,rq)
 S3method(get_statistic,rqs)
 S3method(get_statistic,rqss)
 S3method(get_statistic,scam)
+S3method(get_statistic,sdmTMB)
 S3method(get_statistic,selection)
 S3method(get_statistic,sem)
 S3method(get_statistic,summary.lm)
@@ -1055,6 +1062,7 @@ S3method(get_varcov,probitmfx)
 S3method(get_varcov,robmixglm)
 S3method(get_varcov,rq)
 S3method(get_varcov,rqs)
+S3method(get_varcov,sdmTMB)
 S3method(get_varcov,selection)
 S3method(get_varcov,tobit)
 S3method(get_varcov,truncreg)
@@ -1205,6 +1213,7 @@ S3method(link_function,robmixglm)
 S3method(link_function,rq)
 S3method(link_function,rqs)
 S3method(link_function,rqss)
+S3method(link_function,sdmTMB)
 S3method(link_function,serp)
 S3method(link_function,speedglm)
 S3method(link_function,speedlm)
@@ -1336,6 +1345,7 @@ S3method(link_inverse,robmixglm)
 S3method(link_inverse,rq)
 S3method(link_inverse,rqs)
 S3method(link_inverse,rqss)
+S3method(link_inverse,sdmTMB)
 S3method(link_inverse,serp)
 S3method(link_inverse,speedglm)
 S3method(link_inverse,speedlm)
@@ -1492,6 +1502,7 @@ S3method(model_info,robmixglm)
 S3method(model_info,rq)
 S3method(model_info,rqs)
 S3method(model_info,rqss)
+S3method(model_info,sdmTMB)
 S3method(model_info,serp)
 S3method(model_info,speedglm)
 S3method(model_info,speedlm)

NEWS.md

@@ -1,5 +1,12 @@
 # insight (devel)
 
+## Changes
+
+* The `range` argument in `get_datagrid()` gets a new option, `"pretty"`, to
+  create a range of pretty values.
+
+* Better support for models of class `sdmTMB`.
+
 ## Bug fixes
 
 * Fixed issue in `clean_names()` for *brms* models with `mm()` in formula.

R/compute_variances.R

@@ -163,7 +163,6 @@
   if (component %in% c("residual", "distribution", "all")) {
     var.distribution <- .compute_variance_distribution(
       model,
-      var_cor = mixed_effects_info$vc,
       faminfo,
       model_null = model_null,
       revar_null = var.random_null,
@@ -586,7 +585,6 @@
 # different values for the log/delta/trigamma approximation.
 # -----------------------------------------------------------------------------
 .compute_variance_distribution <- function(model,
-                                           var_cor,
                                            faminfo,
                                            model_null = NULL,
                                            revar_null = NULL,

R/find_formula.R

@@ -1296,6 +1296,36 @@ find_formula.glmmTMB <- function(x, verbose = TRUE, ...) {
 }
 
 
+#' @export
+find_formula.sdmTMB <- function(x, verbose = TRUE, ...) {
+  f.cond <- stats::formula(x)
+
+  # sanity check, we might have a list of formulas, only need first element
+  if (is.list(f.cond)) {
+    f.cond <- f.cond[[1]]
+  }
+
+  # extract random parts of formula
+  f.random <- lapply(.findbars(f.cond), function(.x) {
+    f <- safe_deparse(.x)
+    stats::as.formula(paste0("~", f))
+  })
+
+  if (length(f.random) == 1L) {
+    f.random <- f.random[[1]]
+  }
+
+  # extract fixed effects parts
+  f.cond <- stats::as.formula(.get_fixed_effects(f.cond))
+
+  f <- compact_list(list(
+    conditional = f.cond,
+    random = f.random
+  ))
+  .find_formula_return(f, verbose = verbose)
+}
+
+
 #' @export
 find_formula.nlmerMod <- function(x, verbose = TRUE, ...) {
   f.random <- lapply(.findbars(stats::formula(x)), function(.x) {

R/find_parameters_other.R

@@ -63,6 +63,39 @@ find_parameters.glmgee <- function(x, component = "all", flatten = FALSE, ...) {
 }
 
 
+#' @export
+find_parameters.sdmTMB <- function(x, component = "all", flatten = FALSE, ...) {
+  delta_comp <- isTRUE(x$family$delta)
+  valid_comp <- compact_character(c("all", "conditional", ifelse(delta_comp, "delta", "")))
+  component <- validate_argument(component, valid_comp)
+
+  cf <- suppressMessages(stats::coef(x, model = 1))
+  conditional <- names(cf)
+
+  if (delta_comp) {
+    cf <- suppressMessages(stats::coef(x, model = 2))
+    delta <- names(cf)
+  }
+
+  if (delta_comp) {
+    out <- list(
+      conditional = text_remove_backticks(conditional),
+      delta = text_remove_backticks(delta)
+    )
+  } else {
+    out <- list(conditional = text_remove_backticks(conditional))
+  }
+
+  .filter_parameters(
+    out,
+    effects = "all",
+    component = component,
+    flatten = flatten,
+    recursive = FALSE
+  )
+}
+
+
 #' @export
 find_parameters.betareg <- function(x, component = "all", flatten = FALSE, ...) {
   component <- validate_argument(

R/find_predictors.R

@@ -187,7 +187,7 @@ find_predictors.default <- function(x,
   if (flatten) {
     unique(unlist(l, use.names = FALSE))
   } else {
-    l
+    compact_list(l)
   }
 }
 
@@ -295,7 +295,7 @@ find_predictors.afex_aov <- function(x,
   if (flatten) {
     unique(unlist(l, use.names = FALSE))
   } else {
-    l
+    compact_list(l)
   }
 }
 
@@ -364,6 +364,44 @@ find_predictors.brmsfit <- function(x,
 }
 
 
+#' @export
+find_predictors.sdmTMB <- function(x,
+                                   effects = "fixed",
+                                   flatten = FALSE,
+                                   verbose = TRUE,
+                                   ...) {
+  effects <- validate_argument(effects, c("fixed", "random", "all"))
+  elements <- .get_elements(effects, component = "conditional", model = x)
+
+  f <- find_formula(x, verbose = verbose)
+  f <- .prepare_predictors(x, f, elements)
+
+  # random effects are returned as list, so we need to unlist here
+  l <- .return_vars(f, x)
+
+  if (is_empty_object(l) || is_empty_object(compact_list(l))) {
+    return(NULL)
+  }
+
+  # add time variable
+  l$time <- x$call$time
+
+  # add random slope, if not yet present
+  if (object_has_names(l, "random") && effects == "all") {
+    random_slope <- unlist(find_random_slopes(x), use.names = FALSE)
+    all_predictors <- unlist(unique(l), use.names = FALSE)
+    rs_not_in_pred <- unique(setdiff(random_slope, all_predictors))
+    if (length(rs_not_in_pred)) l$random <- c(rs_not_in_pred, l$random)
+  }
+
+  if (flatten) {
+    unique(unlist(l, use.names = FALSE))
+  } else {
+    compact_list(l)
+  }
+}
+
+
 #' @export
 find_predictors.insight_formula <- function(x, flatten = FALSE, verbose = TRUE, ...) {
   is_mv <- is_multivariate(x)

R/find_statistic.R

@@ -140,6 +140,7 @@ find_statistic.default <- function(x, ...) {
     "qr", "QRNLMM", "QRLMM",
     "Rchoice", "riskRegression", "robmixglm", "rma", "rma.mv", "rma.uni", "rrvglm",
     "Sarlm", "sem", "SemiParBIV", "serp", "slm", "slopes", "survreg", "svy_vglm",
+    "sdmTMB",
     "test_mediation", "tobit",
     "vglm",
     "wbgee",

R/get_datagrid.R

@@ -124,6 +124,8 @@
 #'     length as indicated in `length` are generated. For numeric predictors not
 #'     specified at first in `by`, mean and -1/+1 SD around the mean are
 #'     returned. For factors, all levels are returned.
+#'   - `"pretty"` will create a range "pretty" values, using [`pretty()`], where
+#'     the value in `length` is used for the `n` argument in `pretty()`.
 #'
 #'   `range` can also be a vector of different values (see 'Examples'). In this
 #'   case, `range` must be of same length as numeric target variables. If
@@ -1181,7 +1183,7 @@ get_datagrid.comparisons <- get_datagrid.slopes
                            ...) {
   range <- validate_argument(
     tolower(range),
-    c("range", "iqr", "ci", "hdi", "eti", "sd", "mad", "grid")
+    c("range", "iqr", "ci", "hdi", "eti", "sd", "mad", "grid", "pretty")
   )
 
   # bayestestR only for some options
@@ -1239,6 +1241,11 @@ get_datagrid.comparisons <- get_datagrid.slopes
     return(out)
   }
 
+  # If Range is a range of pretty values
+  if (range == "pretty") {
+    return(pretty(x, n = length))
+  }
+
   # If Range is an interval
   if (range == "iqr") { # nolint
     mini <- stats::quantile(x, (1 - ci) / 2, ...)

R/get_family.R

@@ -47,6 +47,26 @@ get_family.model_fit <- function(x, ...) {
   get_family(x$fit, ...)
 }
 
+#' @export
+get_family.sdmTMB <- function(x, ...) {
+  check_if_installed("sdmTMB")
+  f <- x$family
+  if (length(f$family) > 1) {
+    f <- compact_list(list(
+      family = f$family[2],
+      link = f$link[2],
+      linkfun = f[[2]]$linkfun,
+      linkinv = f[[2]]$linkinv,
+      mu.eta = f[[2]]$mu.eta,
+      valideta = f[[2]]$valideta,
+      name = f[[2]]$name,
+      initialize = f[[2]]$initialize
+    ))
+    class(f) <- "family"
+  }
+  f
+}
+
 
 .get_family <- function(x, ...) {
   info <- model_info(x, response = 1, verbose = FALSE)