diff --git a/NEWS.md b/NEWS.md
index 960ce01e3..9c673a65e 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -1,5 +1,18 @@
 # insight (devel)
 
+## New supported models
+
+* Support for models of class `glmmTMB` fitted with the new `ordinal()`
+  family (*glmmTMB* >= 1.1.15). As for `ordinal::clm()`, the thresholds are
+  returned as conditional parameters ahead of the slopes, with standard errors
+  and statistics obtained via the delta method (the thresholds are estimated on
+  an internal, softmax-parameterized scale), while the intercept, which is
+  fixed to zero by *glmmTMB*, is dropped. `get_predicted()` returns
+  per-category probabilities for `predict = "expectation"` and the most likely
+  category for `predict = "classification"`, and `get_variance()` uses the
+  distribution-specific residual variance of the latent scale, as for
+  `ordinal::clmm()`.
+
 ## Bug fixes
 
 * `find_auxiliary()` no longer returns `"sigma"` for *brms* models that have no
diff --git a/R/compute_variances.R b/R/compute_variances.R
index 1e011d2d2..aa96bdd0c 100644
--- a/R/compute_variances.R
+++ b/R/compute_variances.R
@@ -432,17 +432,22 @@
       ),
       .badlink(faminfo$link_function, faminfo$family, verbose = verbose)
     )
-  } else if (faminfo$is_binomial) {
-    # binomial / bernoulli  ----
-    # --------------------------
+  } else if (faminfo$is_binomial || (faminfo$is_ordinal && inherits(model, "glmmTMB"))) {
+    # binomial / bernoulli / cumulative link  ----
+    # -------------------------------------------
 
     # we need this to adjust for "cbind()" outcomes
     y_factor <- .binomial_response_weight(model)
 
     # for observation level approximation, when we don't want the "fixed"
     # residual variance, pi^2/3, but the variance based on the distribution
-    # of the response
-    pmean <- .obs_level_variance(model_null, revar_null)
+    # of the response. Ordinal models have no meaningful pmean (glmmTMB fixes
+    # the null model's intercept to zero, the thresholds absorb it)
+    if (faminfo$is_ordinal) {
+      pmean <- NULL
+    } else {
+      pmean <- .obs_level_variance(model_null, revar_null)
+    }
 
     # sanity check - clmm-models are "binomial" but have no pmean
     if (is.null(pmean) && identical(approx_method, "observation_level")) {
diff --git a/R/find_parameters_mixed.R b/R/find_parameters_mixed.R
index 02090553f..0975b7b3e 100644
--- a/R/find_parameters_mixed.R
+++ b/R/find_parameters_mixed.R
@@ -50,7 +50,7 @@ find_parameters.glmmTMB <- function(
   )
 
   # installed
-  check_if_installed("lme4")
+  check_if_installed(c("lme4", "glmmTMB"))
 
   # we extract random effects only when really necessary, to save
   # computational time. In particular model with large sample and
@@ -73,6 +73,11 @@ find_parameters.glmmTMB <- function(
     ))
   }
 
+  # ordinal family: thresholds instead of the (fixed) intercept
+  if (.is_glmmtmb_ordinal(x)) {
+    l$conditional <- names(.glmmtmb_ordinal_conditional(x))
+  }
+
   .filter_parameters(l, effects = effects, component = component, flatten = flatten)
 }
 
diff --git a/R/get_parameters_mixed.R b/R/get_parameters_mixed.R
index fa426dc80..e5a2acca2 100644
--- a/R/get_parameters_mixed.R
+++ b/R/get_parameters_mixed.R
@@ -33,7 +33,7 @@
 #' get_parameters(m)
 #' @export
 get_parameters.glmmTMB <- function(x, effects = "fixed", component = "all", ...) {
-  check_if_installed("lme4")
+  check_if_installed(c("lme4", "glmmTMB"))
 
   effects <- validate_argument(effects, c("fixed", "random"))
   # fmt: skip
@@ -59,6 +59,11 @@ get_parameters.glmmTMB <- function(x, effects = "fixed", component = "all", ...)
     ))
   }
 
+  # ordinal family: thresholds instead of the (fixed) intercept
+  if (.is_glmmtmb_ordinal(x)) {
+    l$conditional <- .glmmtmb_ordinal_conditional(x)
+  }
+
   # ---- fixed effects (conditional model)
 
   fixed <- data.frame(
diff --git a/R/get_predicted_mixed.R b/R/get_predicted_mixed.R
index 994ec14fd..50973c3e7 100644
--- a/R/get_predicted_mixed.R
+++ b/R/get_predicted_mixed.R
@@ -103,6 +103,44 @@ get_predicted.glmmTMB <- function(
   verbose = TRUE,
   ...
 ) {
+  # ordinal family: "expectation" returns per-category probabilities (like
+  # `clm`), not `plogis()` of the linear predictor. As in `.get_predicted_args()`,
+  # a `type` argument takes precedence over `predict`
+  dots <- list(...)
+  if (.is_glmmtmb_ordinal(x)) {
+    if (is.null(dots$type)) {
+      requested <- predict[1]
+    } else {
+      requested <- dots$type[1]
+    }
+    ordinal_types <- c(
+      "expectation",
+      "expected",
+      "response",
+      "prediction",
+      "predicted",
+      "classification",
+      "probs"
+    )
+    if (isTRUE(requested %in% ordinal_types)) {
+      if (!is.null(iterations) && verbose) {
+        format_warning(
+          "Bootstrapped predictions are currently not supported for `glmmTMB` models with `ordinal()` family.",
+          "Ignoring the `iterations` argument."
+        )
+      }
+      return(.get_predicted_glmmtmb_ordinal(
+        x,
+        data = data,
+        predict = predict,
+        ci = ci,
+        include_random = include_random,
+        verbose = verbose,
+        dots = dots
+      ))
+    }
+  }
+
   # validation checks
   if (!is.null(predict) && predict %in% c("prediction", "predicted", "classification")) {
     predict <- "expectation"
diff --git a/R/get_statistic.R b/R/get_statistic.R
index 10b700036..d62304444 100644
--- a/R/get_statistic.R
+++ b/R/get_statistic.R
@@ -452,6 +452,25 @@ get_statistic.glmmTMB <- function(x, component = "all", ...) {
     choices = c("all", "conditional", "zi", "zero_inflated", "dispersion")
   )
 
+  # ordinal family: thresholds are not in the summary coefficient table,
+  # so compute Wald statistics from the (delta-method) covariance matrix
+  if (.is_glmmtmb_ordinal(x)) {
+    params <- get_parameters(x, effects = "fixed", component = "conditional")
+    se <- sqrt(diag(get_varcov(x, component = "conditional", verbose = FALSE)))
+    # match by name - non-estimated parameters (rank deficiency, `map`) are
+    # kept by get_parameters() but have no variance
+    stat <- data.frame(
+      Parameter = params$Parameter,
+      Statistic = params$Estimate / unname(se[params$Parameter]),
+      Component = "conditional",
+      stringsAsFactors = FALSE,
+      row.names = NULL
+    )
+    stat <- .filter_component(stat, component)
+    attr(stat, "statistic") <- find_statistic(x)
+    return(stat)
+  }
+
   cs <- compact_list(stats::coef(summary(x)))
   out <- lapply(names(cs), function(i) {
     data.frame(
diff --git a/R/get_varcov.R b/R/get_varcov.R
index 11a65bdc1..a7cd44ad2 100644
--- a/R/get_varcov.R
+++ b/R/get_varcov.R
@@ -593,15 +593,28 @@ get_varcov.glmmTMB <- function(
     c("conditional", "zero_inflated", "zi", "dispersion", "all", "full")
   )
 
-  if (is.null(vcov)) {
-    vc <- switch(
-      component,
-      conditional = .safe_vcov(x)[["cond"]],
-      zi = ,
-      zero_inflated = .safe_vcov(x)[["zi"]],
-      dispersion = .safe_vcov(x)[["disp"]],
-      stats::vcov(x, full = TRUE)
+  # ordinal family: supplied (robust) covariance matrices would be on the
+  # internal softmax scale for the thresholds and are not supported
+  if (.is_glmmtmb_ordinal(x) && !is.null(vcov)) {
+    format_error(
+      "The `vcov` argument is not supported for `glmmTMB` models with `ordinal()` family."
     )
+  }
+
+  if (is.null(vcov)) {
+    if (.is_glmmtmb_ordinal(x) && component %in% c("conditional", "all")) {
+      # thresholds (delta method) and estimated fixed effects
+      vc <- .glmmtmb_ordinal_varcov(x)
+    } else {
+      vc <- switch(
+        component,
+        conditional = .safe_vcov(x)[["cond"]],
+        zi = ,
+        zero_inflated = .safe_vcov(x)[["zi"]],
+        dispersion = .safe_vcov(x)[["disp"]],
+        stats::vcov(x, full = TRUE)
+      )
+    }
   } else {
     vc <- .get_varcov_sandwich(
       x,
@@ -1184,9 +1197,9 @@ get_varcov.LORgee <- get_varcov.gee
 
 # helper-functions -----------------------------------------------------
 
-.safe_vcov <- function(x) {
+.safe_vcov <- function(x, ...) {
   vc <- tryCatch(
-    suppressWarnings(stats::vcov(x)),
+    suppressWarnings(stats::vcov(x, ...)),
     error = function(e) e
   )
   if (inherits(vc, "error")) {
diff --git a/R/utils_glmmtmb_ordinal.R b/R/utils_glmmtmb_ordinal.R
new file mode 100644
index 000000000..ce0fdb81a
--- /dev/null
+++ b/R/utils_glmmtmb_ordinal.R
@@ -0,0 +1,187 @@
+# glmmTMB ordinal family ------------------------------------------------
+# =======================================================================
+
+# glmmTMB's `ordinal()` family stores the K-1 thresholds as family-specific
+# parameters ("psi", softmax-parameterized so that they stay ordered) rather
+# than as fixed effects, and fixes the intercept to zero via an internal map.
+# To match `ordinal::clm()`/`clmm()`, we report the thresholds (on the
+# threshold scale) as conditional parameters ahead of the slopes, drop the
+# non-estimated intercept, and transform the covariance matrix with the
+# delta method. `vcov(x, full = TRUE)` is on the internal psi scale.
+
+.is_glmmtmb_ordinal <- function(x) {
+  inherits(x, "glmmTMB") &&
+    identical(.safe(x$modelInfo$family$family), "ordinal")
+}
+
+
+# is the intercept fixed to zero by glmmTMB's internal map, i.e. not
+# estimated? This is the default unless the user supplied a beta map
+.glmmtmb_fixed_intercept <- function(x) {
+  cf <- .safe(lme4::fixef(x)$cond)
+  bmap <- .safe(x$obj$env$map$beta)
+  icpt <- which(names(cf) == "(Intercept)")
+  length(icpt) == 1L && !is.null(bmap) && is.na(bmap[icpt])
+}
+
+
+# named vector of thresholds followed by the estimated fixed effects
+.glmmtmb_ordinal_conditional <- function(x) {
+  cf <- lme4::fixef(x)$cond
+  if (.glmmtmb_fixed_intercept(x)) {
+    cf <- cf[names(cf) != "(Intercept)"]
+  }
+  c(glmmTMB::family_params(x), cf)
+}
+
+
+# covariance matrix of c(thresholds, fixed effects) on the threshold scale.
+# theta_j = qlogis(cumsum(softmax(c(psi, 0)))_j) is a joint function of all
+# psi elements, so the Jacobian is J[j, m] = s[m] * ((m <= j) - C_j) /
+# (C_j * (1 - C_j)), with s = softmax(c(psi, 0)) and C_j = cumsum(s)[j]
+.glmmtmb_ordinal_varcov <- function(x) {
+  check_if_installed(c("lme4", "glmmTMB"))
+  V <- .safe_vcov(x, full = TRUE)
+  thresholds <- glmmTMB::family_params(x)
+  cf <- lme4::fixef(x)$cond
+
+  # internal parameter ids ("cond1", "psi1", ...) are carried as names of the
+  # row labels; fall back to the labels if these are missing
+  ids <- names(rownames(V))
+  if (is.null(ids)) {
+    psi_i <- match(names(thresholds), rownames(V))
+    cond_i <- match(names(cf), rownames(V))
+  } else {
+    psi_i <- which(startsWith(ids, "psi"))
+    cond_i <- which(startsWith(ids, "cond"))
+  }
+  # non-estimated parameters (the mapped intercept) have NA variance
+  cond_i <- cond_i[!is.na(cond_i)]
+  cond_i <- cond_i[!is.na(diag(V)[cond_i])]
+
+  psi <- unname(x$fit$par[names(x$fit$par) == "psi"])
+  k <- length(psi)
+  w <- exp(c(psi, 0) - max(psi, 0))
+  s <- w / sum(w)
+  Cj <- cumsum(s)[seq_len(k)]
+  J_psi <- outer(
+    seq_len(k),
+    seq_len(k),
+    function(j, m) s[m] * ((m <= j) - Cj[j]) / (Cj[j] * (1 - Cj[j]))
+  )
+
+  J <- diag(k + length(cond_i))
+  J[seq_len(k), seq_len(k)] <- J_psi
+  idx <- c(psi_i, cond_i)
+  out <- J %*% V[idx, idx, drop = FALSE] %*% t(J)
+  nms <- c(names(thresholds), unname(rownames(V)[cond_i]))
+  dimnames(out) <- list(nms, nms)
+  out
+}
+
+
+# per-category probabilities (or the most likely category) from
+# `predict(type = "probs")`, reshaped to long format like `get_predicted.clm()`
+.get_predicted_glmmtmb_ordinal <- function(
+  x,
+  data = NULL,
+  predict = "expectation",
+  ci = NULL,
+  include_random = "default",
+  verbose = TRUE,
+  dots = list()
+) {
+  # `type` takes precedence over `predict`, as in `.get_predicted_args()`
+  if (is.null(dots$type)) {
+    requested <- predict[1]
+  } else {
+    requested <- dots$type[1]
+  }
+  classification <- identical(requested, "classification")
+
+  if (requested %in% c("prediction", "predicted") && verbose) {
+    format_warning(
+      "\"prediction\" is currently not supported by the `predict` argument for `glmmTMB` models.",
+      "Changing to `predict=\"expectation\"`."
+    )
+  }
+  if (classification && !is.null(ci)) {
+    if (verbose) {
+      format_warning("Confidence intervals are not available for classification.")
+    }
+    ci <- NULL
+  }
+
+  # sanitize input. This also handles the `predict`/`type` alert and the
+  # aliases; the type passed to `predict()` is always "probs"
+  my_args <- do.call(
+    .get_predicted_args,
+    c(
+      list(
+        x,
+        data = data,
+        predict = predict,
+        ci = ci,
+        include_random = include_random,
+        verbose = verbose
+      ),
+      dots
+    )
+  )
+  my_args$predict <- ifelse(classification, "classification", "expectation")
+
+  # remaining dot-arguments are forwarded to `predict()`, as for other
+  # families; those managed here take precedence
+  predict_args <- list(
+    x,
+    newdata = my_args$data,
+    type = "probs",
+    re.form = my_args$re.form,
+    allow.new.levels = my_args$allow_new_levels,
+    se.fit = !classification
+  )
+  dots[c("type", "newdata", "re.form", "allow.new.levels", "se.fit")] <- NULL
+  rez <- do.call(stats::predict, c(predict_args, dots))
+  if (classification) {
+    probs <- rez
+  } else {
+    probs <- rez$fit
+  }
+  resp_levels <- colnames(probs)
+
+  # the response is not a focal predictor of the long-format output
+  resp <- find_response(x)
+  if (!is.null(my_args$data) && !is.null(resp)) {
+    my_args$data <- my_args$data[, setdiff(colnames(my_args$data), resp), drop = FALSE]
+  }
+
+  if (classification) {
+    out <- factor(
+      resp_levels[max.col(probs, ties.method = "first")],
+      levels = resp_levels
+    )
+    return(.get_predicted_out(out, my_args = my_args))
+  }
+
+  # predictions matrix to long format (column-major, i.e. by response level)
+  out <- .get_predicted_out(probs, my_args = my_args)
+  ci_data <- data.frame(
+    Row = out$Row,
+    Response = out$Response,
+    SE = as.vector(rez$se.fit),
+    stringsAsFactors = FALSE
+  )
+
+  # intervals on the logit scale via the delta method, as in
+  # `ordinal:::predict.clm()`, so that they stay within [0, 1]
+  if (!is.null(ci) && !is.na(ci)) {
+    crit_val <- stats::qnorm((1 + ci) / 2)
+    p <- out$Predicted
+    se_logit <- ci_data$SE / (p * (1 - p))
+    ci_data$CI_low <- stats::plogis(stats::qlogis(p) - crit_val * se_logit)
+    ci_data$CI_high <- stats::plogis(stats::qlogis(p) + crit_val * se_logit)
+  }
+
+  attr(out, "ci_data") <- ci_data
+  out
+}
diff --git a/inst/WORDLIST b/inst/WORDLIST
index de9ee076a..811293458 100644
--- a/inst/WORDLIST
+++ b/inst/WORDLIST
@@ -191,6 +191,7 @@ sampleSelection
 serp
 sjPlot
 smicd
+softmax
 spaMM
 spatialreg
 specificities
diff --git a/tests/testthat/test-get_predicted.R b/tests/testthat/test-get_predicted.R
index 93599600d..23c560342 100644
--- a/tests/testthat/test-get_predicted.R
+++ b/tests/testthat/test-get_predicted.R
@@ -907,7 +907,7 @@ test_that("get_predicted works with brms-Wiener (cogmod-RT-choice)", {
     c(
       "Row",
       "Component",
-      "rt",
+      "RT",
       "Error",
       "Condition",
       "Predicted",
diff --git a/tests/testthat/test-glmmTMB-ordinal.R b/tests/testthat/test-glmmTMB-ordinal.R
new file mode 100644
index 000000000..7379d7715
--- /dev/null
+++ b/tests/testthat/test-glmmTMB-ordinal.R
@@ -0,0 +1,306 @@
+skip_if_not_installed("glmmTMB", minimum_version = "1.1.15")
+skip_if_not_installed("ordinal")
+skip_if_not_installed("lme4")
+
+data(wine, package = "ordinal")
+
+# fixed effects only, compared against ordinal::clm()
+m_tmb <- glmmTMB::glmmTMB(
+  rating ~ temp * contact,
+  data = wine,
+  family = glmmTMB::ordinal()
+)
+m_clm <- ordinal::clm(rating ~ temp * contact, data = wine)
+
+# mixed, compared against ordinal::clmm()
+m_tmb_mixed <- 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("model_info", {
+  expect_true(model_info(m_tmb)$is_ordinal)
+  expect_true(model_info(m_tmb)$is_logit)
+  expect_false(model_info(m_tmb)$is_linear)
+  expect_true(model_info(m_tmb_mixed)$is_mixed)
+})
+
+test_that("find_parameters: thresholds, no fixed intercept", {
+  expect_identical(find_parameters(m_tmb), find_parameters(m_clm))
+  expect_identical(
+    find_parameters(m_tmb_mixed, effects = "fixed"),
+    find_parameters(m_clmm)
+  )
+  expect_identical(
+    find_parameters(m_tmb_mixed),
+    list(
+      conditional = c("1|2", "2|3", "3|4", "4|5", "tempwarm", "contactyes"),
+      random = list(judge = "(Intercept)")
+    )
+  )
+})
+
+test_that("get_parameters: thresholds match ordinal", {
+  out <- get_parameters(m_tmb)
+  expect_identical(out$Parameter, get_parameters(m_clm)$Parameter)
+  expect_equal(out$Estimate, get_parameters(m_clm)$Estimate, tolerance = 1e-3)
+  expect_identical(unique(out$Component), "conditional")
+
+  out <- get_parameters(m_tmb_mixed)
+  expect_identical(out$Parameter, get_parameters(m_clmm)$Parameter)
+  expect_equal(out$Estimate, get_parameters(m_clmm)$Estimate, tolerance = 1e-3)
+
+  out <- get_parameters(m_tmb_mixed, effects = "random")
+  expect_named(out$random, "judge")
+})
+
+test_that("get_varcov: delta-method threshold SEs match ordinal", {
+  vc <- get_varcov(m_tmb)
+  expect_identical(dimnames(vc)[[1]], get_parameters(m_clm)$Parameter)
+  expect_equal(
+    sqrt(diag(vc)),
+    sqrt(diag(vcov(m_clm))),
+    tolerance = 1e-3,
+    ignore_attr = TRUE
+  )
+
+  vc <- get_varcov(m_tmb_mixed)
+  ref <- vcov(m_clmm)[1:6, 1:6]
+  expect_identical(dimnames(vc)[[1]], dimnames(ref)[[1]])
+  expect_equal(vc, ref, tolerance = 1e-2, ignore_attr = TRUE)
+
+  # "full" stays on glmmTMB's internal (softmax) scale for the thresholds
+  vc_full <- suppressWarnings(get_varcov(m_tmb, component = "full"))
+  expect_false(isTRUE(all.equal(
+    unname(sqrt(diag(vc_full))[c("1|2", "2|3")]),
+    unname(sqrt(diag(vc))[c("1|2", "2|3")])
+  )))
+})
+
+test_that("get_statistic matches ordinal", {
+  out <- get_statistic(m_tmb)
+  expect_identical(out$Parameter, get_statistic(m_clm)$Parameter)
+  expect_equal(out$Statistic, get_statistic(m_clm)$Statistic, tolerance = 1e-3)
+  expect_identical(attributes(out)$statistic, "z-statistic")
+
+  out <- get_statistic(m_tmb_mixed)
+  expect_equal(out$Statistic, get_statistic(m_clmm)$Statistic, tolerance = 1e-2)
+})
+
+test_that("get_predicted: per-category probabilities", {
+  out <- get_predicted(m_tmb, ci = 0.95, verbose = FALSE)
+  expect_s3_class(out, "data.frame")
+  expect_true(all(c("Row", "Response", "Predicted") %in% colnames(out)))
+  expect_shape(out, nrow = nrow(wine) * nlevels(wine$rating))
+
+  # matches glmmTMB's own predict(type = "probs")
+  pr <- predict(m_tmb, type = "probs", se.fit = TRUE)
+  expect_equal(out$Predicted, as.vector(pr$fit), tolerance = 1e-8)
+  ci_data <- attributes(out)$ci_data
+  expect_equal(ci_data$SE, as.vector(pr$se.fit), tolerance = 1e-8)
+
+  # intervals match ordinal::clm() and stay within [0, 1]
+  out_clm <- get_predicted(m_clm, ci = 0.95, verbose = FALSE)
+  ci_clm <- attributes(out_clm)$ci_data
+  expect_identical(colnames(out), setdiff(colnames(out_clm), "rating"))
+  expect_equal(out$Predicted, out_clm$Predicted, tolerance = 1e-3)
+  expect_equal(ci_data$SE, ci_clm$SE, tolerance = 1e-3)
+  expect_equal(ci_data$CI_low, ci_clm$CI_low, tolerance = 1e-3)
+  expect_equal(ci_data$CI_high, ci_clm$CI_high, tolerance = 1e-3)
+  expect_true(all(ci_data$CI_low >= 0 & ci_data$CI_high <= 1))
+  expect_false("rating" %in% colnames(out))
+
+  # data grid
+  dg <- get_datagrid(m_tmb, "temp", verbose = FALSE)
+  out <- get_predicted(m_tmb, data = dg, verbose = FALSE)
+  expect_shape(out, nrow = nrow(dg) * nlevels(wine$rating))
+  expect_true("temp" %in% colnames(out))
+  expect_equal(
+    out$Predicted,
+    as.vector(predict(m_tmb, newdata = dg, type = "probs")),
+    tolerance = 1e-8
+  )
+
+  # link scale is untouched
+  expect_equal(
+    as.vector(get_predicted(m_tmb, predict = "link", verbose = FALSE)),
+    unname(predict(m_tmb, type = "link")),
+    tolerance = 1e-8
+  )
+
+  # random effects
+  out <- get_predicted(m_tmb_mixed, verbose = FALSE)
+  expect_equal(
+    out$Predicted,
+    as.vector(predict(m_tmb_mixed, type = "probs")),
+    tolerance = 1e-8
+  )
+})
+
+test_that("get_predicted: type and include_random arguments", {
+  # explicit `type = "response"` means probabilities, as for clm
+  out <- get_predicted(m_tmb, predict = NULL, type = "response", verbose = FALSE)
+  expect_equal(
+    out$Predicted,
+    as.vector(predict(m_tmb, type = "probs")),
+    tolerance = 1e-8
+  )
+  out <- get_predicted(m_tmb, predict = NULL, type = "probs", verbose = FALSE)
+  expect_true("Response" %in% colnames(out))
+  out <- get_predicted(m_tmb, predict = "probs", verbose = FALSE)
+  expect_true("Response" %in% colnames(out))
+
+  # `type` takes precedence over the default `predict`
+  out <- suppressMessages(get_predicted(m_tmb, type = "link"))
+  expect_false("Response" %in% colnames(out))
+  expect_equal(
+    as.vector(out),
+    unname(predict(m_tmb, type = "link")),
+    tolerance = 1e-8
+  )
+
+  # glmmTMB's own other types pass through
+  expect_equal(
+    as.vector(get_predicted(
+      m_tmb,
+      predict = NULL,
+      type = "conditional",
+      verbose = FALSE
+    )),
+    unname(predict(m_tmb, type = "conditional")),
+    tolerance = 1e-8
+  )
+
+  # bootstrapping and prediction intervals are not implemented
+  expect_warning(get_predicted(m_tmb, iterations = 5), "Bootstrapped")
+  expect_warning(get_predicted(m_tmb, predict = "prediction"), "not supported")
+
+  # population-level predictions
+  out <- get_predicted(m_tmb_mixed, include_random = FALSE, verbose = FALSE)
+  expect_equal(
+    out$Predicted,
+    as.vector(predict(m_tmb_mixed, type = "probs", re.form = NA)),
+    tolerance = 1e-8
+  )
+  dg <- get_datagrid(m_tmb_mixed, "temp", verbose = FALSE)
+  out <- get_predicted(m_tmb_mixed, data = dg, verbose = FALSE)
+  expect_equal(
+    out$Predicted,
+    as.vector(predict(m_tmb_mixed, newdata = dg, type = "probs", re.form = NA)),
+    tolerance = 1e-8
+  )
+})
+
+test_that("probit link", {
+  m_probit <- glmmTMB::glmmTMB(
+    rating ~ temp + contact,
+    data = wine,
+    family = glmmTMB::ordinal(link = "probit")
+  )
+  m_clm_probit <- ordinal::clm(rating ~ temp + contact, data = wine, link = "probit")
+  expect_true(model_info(m_probit)$is_probit)
+  expect_equal(
+    get_parameters(m_probit)$Estimate,
+    get_parameters(m_clm_probit)$Estimate,
+    tolerance = 1e-3
+  )
+  expect_equal(
+    sqrt(diag(get_varcov(m_probit))),
+    sqrt(diag(vcov(m_clm_probit))),
+    tolerance = 1e-3,
+    ignore_attr = TRUE
+  )
+})
+
+test_that("get_variance matches clmm", {
+  out <- get_variance(m_tmb_mixed)
+  ref <- get_variance(m_clmm)
+  expect_equal(out$var.residual, pi^2 / 3, tolerance = 1e-6)
+  expect_equal(unlist(out), unlist(ref), tolerance = 1e-3)
+  # observation-level approximation falls back, as for clmm
+  expect_equal(
+    get_variance(m_tmb_mixed, approximation = "observation_level")$var.residual,
+    pi^2 / 3,
+    tolerance = 1e-6
+  )
+})
+
+test_that("non-estimated parameters", {
+  # rank deficient model: the dropped coefficient is kept (as NA) in
+  # get_parameters(), like for other glmmTMB models, but has no variance
+  wine2 <- wine
+  wine2$temp2 <- wine2$temp
+  m_rd <- suppressMessages(glmmTMB::glmmTMB(
+    rating ~ temp + temp2 + contact,
+    data = wine2,
+    family = glmmTMB::ordinal()
+  ))
+  params <- get_parameters(m_rd)
+  expect_identical(
+    params$Parameter,
+    c("1|2", "2|3", "3|4", "4|5", "tempwarm", "temp2warm", "contactyes")
+  )
+  expect_true(is.na(params$Estimate[params$Parameter == "temp2warm"]))
+  vc <- suppressWarnings(get_varcov(m_rd))
+  expect_false("temp2warm" %in% dimnames(vc)[[1]])
+  out <- suppressWarnings(get_statistic(m_rd))
+  expect_identical(out$Parameter, params$Parameter)
+  expect_true(is.na(out$Statistic[out$Parameter == "temp2warm"]))
+  ref <- coef(summary(m_rd))$cond[, "z value"]
+  ref <- ref[!is.na(ref)]
+  expect_equal(
+    out$Statistic[match(names(ref), out$Parameter)],
+    unname(ref),
+    tolerance = 1e-6
+  )
+
+  # supplied covariance matrices are rejected explicitly
+  expect_error(get_varcov(m_tmb, vcov = "HC0"), "not supported")
+
+  # errors from glmmTMB are wrapped in insight's message
+  m_disp <- glmmTMB::glmmTMB(
+    rating ~ temp + contact,
+    dispformula = ~contact,
+    data = wine,
+    family = glmmTMB::ordinal()
+  )
+  expect_error(get_varcov(m_disp), "Can't extract")
+})
+
+test_that("get_predicted: dot-arguments reach predict()", {
+  # `na.action` is not managed by insight and is forwarded to `predict()`:
+  # by default (`na.pass`), rows with missing predictors are kept as `NA`,
+  # with `na.omit` they are dropped
+  wine_na <- wine
+  wine_na$temp[1] <- NA
+  out <- get_predicted(
+    m_tmb,
+    data = wine_na,
+    predict = "classification",
+    verbose = FALSE
+  )
+  expect_length(out, nrow(wine_na))
+  expect_true(is.na(out[1]))
+  out <- get_predicted(
+    m_tmb,
+    data = wine_na,
+    predict = "classification",
+    na.action = stats::na.omit,
+    verbose = FALSE
+  )
+  expect_length(out, nrow(wine_na) - 1L)
+})
+
+test_that("get_predicted: classification", {
+  out <- get_predicted(m_tmb, predict = "classification", verbose = FALSE)
+  expect_s3_class(out, "factor")
+  expect_identical(levels(out), levels(wine$rating))
+  pr <- predict(m_tmb, type = "probs")
+  expect_identical(as.character(out), colnames(pr)[max.col(pr)])
+  expect_warning(
+    get_predicted(m_tmb, predict = "classification", ci = 0.95),
+    "not available"
+  )
+})