From fcc5e3b7579b6e3e44b07e3712eb3fbc7e48a516 Mon Sep 17 00:00:00 2001 From: Daniel Date: Sat, 29 Aug 2026 12:04:10 +0200 Subject: [PATCH 1/7] glmmTMB ordinal support Fixes #1220 --- DESCRIPTION | 2 +- NEWS.md | 4 ++++ R/find_parameters_mixed.R | 7 +++++++ R/get_parameters_mixed.R | 4 ++++ 4 files changed, 16 insertions(+), 1 deletion(-) diff --git a/DESCRIPTION b/DESCRIPTION index b130a53e6..9b2236068 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,7 +1,7 @@ Type: Package Package: insight Title: Easy Access to Model Information for Various Model Objects -Version: 1.5.3.1 +Version: 1.5.3.2 Authors@R: c(person(given = "Daniel", family = "Lüdecke", diff --git a/NEWS.md b/NEWS.md index 4fc070cb9..f8944ba31 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,5 +1,9 @@ # insight (devel) +## Changes + +* *insight* now supports the new `ordinal` family from *glmmTMB* models. + ## Bug fixes * Fixed a bug where `get_statistic()` incorrectly extracted degrees of freedom diff --git a/R/find_parameters_mixed.R b/R/find_parameters_mixed.R index 02090553f..2ccadcee2 100644 --- a/R/find_parameters_mixed.R +++ b/R/find_parameters_mixed.R @@ -73,6 +73,13 @@ find_parameters.glmmTMB <- function( )) } + # handle ordinal models - no intercept + if (identical(stats::family(x)$family, "ordinal")) { + l <- lapply(l, function(i) { + setdiff(i, "(Intercept)") + }) + } + .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..4f5f933f2 100644 --- a/R/get_parameters_mixed.R +++ b/R/get_parameters_mixed.R @@ -105,6 +105,10 @@ get_parameters.glmmTMB <- function(x, effects = "fixed", component = "all", ...) zero_inflated = fixedzi, dispersion = fixeddisp ) + # handle ordinal models - no intercept + if (identical(stats::family(x)$family, "ordinal")) { + out <- out[-1, ] + } text_remove_backticks(out) } else if (effects == "random") { switch( From 35888b522e828e34684dc9c32a99f65b03a3fb38 Mon Sep 17 00:00:00 2001 From: Daniel Date: Mon, 31 Aug 2026 11:00:48 +0200 Subject: [PATCH 2/7] include thresholds --- R/find_parameters_mixed.R | 2 ++ R/get_parameters_mixed.R | 5 +++++ R/get_varcov.R | 18 +++++++++++++++++- 3 files changed, 24 insertions(+), 1 deletion(-) diff --git a/R/find_parameters_mixed.R b/R/find_parameters_mixed.R index 2ccadcee2..a13ca294a 100644 --- a/R/find_parameters_mixed.R +++ b/R/find_parameters_mixed.R @@ -78,6 +78,8 @@ find_parameters.glmmTMB <- function( l <- lapply(l, function(i) { setdiff(i, "(Intercept)") }) + # add threshold names + l$conditional <- c(names(glmmTMB::family_params(x)), l$conditional) } .filter_parameters(l, effects = effects, component = component, flatten = flatten) diff --git a/R/get_parameters_mixed.R b/R/get_parameters_mixed.R index 4f5f933f2..c7b1ef0b5 100644 --- a/R/get_parameters_mixed.R +++ b/R/get_parameters_mixed.R @@ -59,6 +59,11 @@ get_parameters.glmmTMB <- function(x, effects = "fixed", component = "all", ...) )) } + # handle ordinal models - add thresholds + if (identical(stats::family(x)$family, "ordinal")) { + l$conditional <- c(glmmTMB::family_params(x), l$conditional) + } + # ---- fixed effects (conditional model) fixed <- data.frame( diff --git a/R/get_varcov.R b/R/get_varcov.R index 11a65bdc1..431e8418e 100644 --- a/R/get_varcov.R +++ b/R/get_varcov.R @@ -593,6 +593,12 @@ get_varcov.glmmTMB <- function( c("conditional", "zero_inflated", "zi", "dispersion", "all", "full") ) + # handle ordinal models - we need full varcov here + is_ordinal <- identical(stats::family(x)$family, "ordinal") + if (is_ordinal) { + component <- "full" + } + if (is.null(vcov)) { vc <- switch( component, @@ -631,10 +637,20 @@ get_varcov.glmmTMB <- function( # drop theta parameters theta_parms <- startsWith(colnames(vc), "theta_") - if (any(theta_parms) && component != "full") { + if (any(theta_parms) && (is_ordinal || component != "full")) { vc <- vc[!theta_parms, !theta_parms, drop = FALSE] } + # reorder for ordinal, to be in line with order of parameters + if (is_ordinal) { + # find parameters to re-order the vcov-matrix + params <- intersect(find_parameters(x)$conditional, colnames(vc)) + # check if dimensions still match + if (length(params) == ncol(vc)) { + vc <- vc[params, params] + } + } + .process_vcov(vc, verbose, ...) } From 556aa2d3471b4084efe83f6af2b1bb9d2faf0832 Mon Sep 17 00:00:00 2001 From: Jeffrey Girard Date: Tue, 1 Sep 2026 17:46:04 -0500 Subject: [PATCH 3/7] support glmmTMB ordinal() family (#1220) Thresholds are stored by glmmTMB as softmax-parameterized family parameters (psi) and the intercept is fixed to zero via an internal map. To match ordinal::clm()/clmm(), find_parameters() and get_parameters() now return the thresholds ahead of the slopes and drop the intercept, get_varcov() transforms the covariance matrix to the threshold scale with the delta method, get_statistic() derives Wald statistics from it (it previously errored on the row-count mismatch), and get_predicted() returns per-category probabilities via predict(type = "probs") for "expectation" and the most likely category for "classification". --- NEWS.md | 13 +- R/compute_variances.R | 15 +- R/find_parameters_mixed.R | 10 +- R/get_parameters_mixed.R | 10 +- R/get_predicted_mixed.R | 38 ++++ R/get_statistic.R | 19 ++ R/get_varcov.R | 43 ++-- R/utils_glmmtmb_ordinal.R | 182 +++++++++++++++++ tests/testthat/test-glmmTMB-ordinal.R | 279 ++++++++++++++++++++++++++ 9 files changed, 561 insertions(+), 48 deletions(-) create mode 100644 R/utils_glmmtmb_ordinal.R create mode 100644 tests/testthat/test-glmmTMB-ordinal.R diff --git a/NEWS.md b/NEWS.md index f8944ba31..5fc35b8ad 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,8 +1,17 @@ # insight (devel) -## Changes +## New supported models -* *insight* now supports the new `ordinal` family from *glmmTMB* 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 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 a13ca294a..b8bc52eb5 100644 --- a/R/find_parameters_mixed.R +++ b/R/find_parameters_mixed.R @@ -73,13 +73,9 @@ find_parameters.glmmTMB <- function( )) } - # handle ordinal models - no intercept - if (identical(stats::family(x)$family, "ordinal")) { - l <- lapply(l, function(i) { - setdiff(i, "(Intercept)") - }) - # add threshold names - l$conditional <- c(names(glmmTMB::family_params(x)), l$conditional) + # 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 c7b1ef0b5..48c2a34ba 100644 --- a/R/get_parameters_mixed.R +++ b/R/get_parameters_mixed.R @@ -59,9 +59,9 @@ get_parameters.glmmTMB <- function(x, effects = "fixed", component = "all", ...) )) } - # handle ordinal models - add thresholds - if (identical(stats::family(x)$family, "ordinal")) { - l$conditional <- c(glmmTMB::family_params(x), l$conditional) + # ordinal family: thresholds instead of the (fixed) intercept + if (.is_glmmtmb_ordinal(x)) { + l$conditional <- .glmmtmb_ordinal_conditional(x) } # ---- fixed effects (conditional model) @@ -110,10 +110,6 @@ get_parameters.glmmTMB <- function(x, effects = "fixed", component = "all", ...) zero_inflated = fixedzi, dispersion = fixeddisp ) - # handle ordinal models - no intercept - if (identical(stats::family(x)$family, "ordinal")) { - out <- out[-1, ] - } text_remove_backticks(out) } else if (effects == "random") { switch( 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 431e8418e..e3d2d8ed9 100644 --- a/R/get_varcov.R +++ b/R/get_varcov.R @@ -593,21 +593,20 @@ get_varcov.glmmTMB <- function( c("conditional", "zero_inflated", "zi", "dispersion", "all", "full") ) - # handle ordinal models - we need full varcov here - is_ordinal <- identical(stats::family(x)$family, "ordinal") - if (is_ordinal) { - component <- "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) - ) + 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, @@ -637,20 +636,10 @@ get_varcov.glmmTMB <- function( # drop theta parameters theta_parms <- startsWith(colnames(vc), "theta_") - if (any(theta_parms) && (is_ordinal || component != "full")) { + if (any(theta_parms) && component != "full") { vc <- vc[!theta_parms, !theta_parms, drop = FALSE] } - # reorder for ordinal, to be in line with order of parameters - if (is_ordinal) { - # find parameters to re-order the vcov-matrix - params <- intersect(find_parameters(x)$conditional, colnames(vc)) - # check if dimensions still match - if (length(params) == ncol(vc)) { - vc <- vc[params, params] - } - } - .process_vcov(vc, verbose, ...) } @@ -1200,9 +1189,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..a6776c617 --- /dev/null +++ b/R/utils_glmmtmb_ordinal.R @@ -0,0 +1,182 @@ +# 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) { + 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") + + rez <- stats::predict( + x, + newdata = my_args$data, + type = "probs", + re.form = my_args$re.form, + allow.new.levels = my_args$allow_new_levels, + se.fit = !classification + ) + 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/tests/testthat/test-glmmTMB-ordinal.R b/tests/testthat/test-glmmTMB-ordinal.R new file mode 100644 index 000000000..cf199df5e --- /dev/null +++ b/tests/testthat/test-glmmTMB-ordinal.R @@ -0,0 +1,279 @@ +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_identical(nrow(out), 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_identical(nrow(out), 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 + ) + + # 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: 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" + ) +}) From f331d21c445f242f067f411b318222eecab70908 Mon Sep 17 00:00:00 2001 From: Jeffrey Girard Date: Tue, 8 Sep 2026 11:00:14 -0500 Subject: [PATCH 4/7] add softmax to WORDLIST, use expect_shape() in ordinal tests Clears the two lint-changed-files diagnostics and the one spelling error introduced by the glmmTMB ordinal() support. The remaining lints and spelling words in these checks predate this branch. --- inst/WORDLIST | 1 + tests/testthat/test-glmmTMB-ordinal.R | 4 ++-- 2 files changed, 3 insertions(+), 2 deletions(-) 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-glmmTMB-ordinal.R b/tests/testthat/test-glmmTMB-ordinal.R index cf199df5e..5b52b63cd 100644 --- a/tests/testthat/test-glmmTMB-ordinal.R +++ b/tests/testthat/test-glmmTMB-ordinal.R @@ -93,7 +93,7 @@ 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_identical(nrow(out), nrow(wine) * nlevels(wine$rating)) + 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) @@ -115,7 +115,7 @@ test_that("get_predicted: per-category probabilities", { # data grid dg <- get_datagrid(m_tmb, "temp", verbose = FALSE) out <- get_predicted(m_tmb, data = dg, verbose = FALSE) - expect_identical(nrow(out), nrow(dg) * nlevels(wine$rating)) + expect_shape(out, nrow = nrow(dg) * nlevels(wine$rating)) expect_true("temp" %in% colnames(out)) expect_equal( out$Predicted, From ac174b2143e79f83f17106a8939f8847c96b4929 Mon Sep 17 00:00:00 2001 From: Jeffrey Girard Date: Fri, 11 Sep 2026 01:17:50 -0500 Subject: [PATCH 5/7] address review: move NEWS entry to devel, reject vcov for ordinal, forward dots insight 1.5.4 was released without the ordinal support, so the NEWS entry moves from the released section to (devel). get_varcov() now rejects supplied covariance matrices for the ordinal family explicitly (they would be on the internal softmax scale), and the ordinal prediction path forwards remaining dot-arguments to predict(), as the other families do. --- NEWS.md | 26 +++++++++++++------------- R/get_varcov.R | 8 ++++++++ R/utils_glmmtmb_ordinal.R | 6 +++++- tests/testthat/test-glmmTMB-ordinal.R | 27 +++++++++++++++++++++++++++ 4 files changed, 53 insertions(+), 14 deletions(-) diff --git a/NEWS.md b/NEWS.md index 55b4bcf0c..9c673a65e 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,18 +1,5 @@ # insight (devel) -## Bug fixes - -* `find_auxiliary()` no longer returns `"sigma"` for *brms* models that have no - residual standard deviation, but auxiliary parameters whose names merely - contain `"sigma"` (e.g. `"sigmadrift"` or `"sigmabias"` from custom families). - -* `clean_parameters()` no longer assigns auxiliary parameters whose names - contain `"sigma"` (like `"sigmabias"`) to the `"sigma"` component. These are - now returned as their own component, which also fixes the related grouping in - `parameters::model_parameters()`. - -# insight 1.5.4 - ## New supported models * Support for models of class `glmmTMB` fitted with the new `ordinal()` @@ -28,6 +15,19 @@ ## Bug fixes +* `find_auxiliary()` no longer returns `"sigma"` for *brms* models that have no + residual standard deviation, but auxiliary parameters whose names merely + contain `"sigma"` (e.g. `"sigmadrift"` or `"sigmabias"` from custom families). + +* `clean_parameters()` no longer assigns auxiliary parameters whose names + contain `"sigma"` (like `"sigmabias"`) to the `"sigma"` component. These are + now returned as their own component, which also fixes the related grouping in + `parameters::model_parameters()`. + +# insight 1.5.4 + +## Bug fixes + * Fixed a bug where `get_statistic()` incorrectly extracted degrees of freedom (df) instead of the t value for `rlmerMod` models (*robustlmm*) after Satterthwaite degrees of freedom were cached (e.g., following `emmeans` diff --git a/R/get_varcov.R b/R/get_varcov.R index e3d2d8ed9..a7cd44ad2 100644 --- a/R/get_varcov.R +++ b/R/get_varcov.R @@ -593,6 +593,14 @@ get_varcov.glmmTMB <- function( c("conditional", "zero_inflated", "zi", "dispersion", "all", "full") ) + # 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 diff --git a/R/utils_glmmtmb_ordinal.R b/R/utils_glmmtmb_ordinal.R index a6776c617..2d45f5f20 100644 --- a/R/utils_glmmtmb_ordinal.R +++ b/R/utils_glmmtmb_ordinal.R @@ -129,7 +129,9 @@ ) my_args$predict <- ifelse(classification, "classification", "expectation") - rez <- stats::predict( + # 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", @@ -137,6 +139,8 @@ 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 { diff --git a/tests/testthat/test-glmmTMB-ordinal.R b/tests/testthat/test-glmmTMB-ordinal.R index 5b52b63cd..7379d7715 100644 --- a/tests/testthat/test-glmmTMB-ordinal.R +++ b/tests/testthat/test-glmmTMB-ordinal.R @@ -256,6 +256,9 @@ test_that("non-estimated parameters", { 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, @@ -266,6 +269,30 @@ test_that("non-estimated parameters", { 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") From be2e80d3287d183bec73a1e040af82cffa928f3b Mon Sep 17 00:00:00 2001 From: Daniel Date: Wed, 16 Sep 2026 11:55:29 +0200 Subject: [PATCH 6/7] minor edits --- R/find_parameters_mixed.R | 2 +- R/get_parameters_mixed.R | 2 +- R/utils_glmmtmb_ordinal.R | 1 + 3 files changed, 3 insertions(+), 2 deletions(-) diff --git a/R/find_parameters_mixed.R b/R/find_parameters_mixed.R index b8bc52eb5..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 diff --git a/R/get_parameters_mixed.R b/R/get_parameters_mixed.R index 48c2a34ba..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 diff --git a/R/utils_glmmtmb_ordinal.R b/R/utils_glmmtmb_ordinal.R index 2d45f5f20..ce0fdb81a 100644 --- a/R/utils_glmmtmb_ordinal.R +++ b/R/utils_glmmtmb_ordinal.R @@ -40,6 +40,7 @@ # 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 From b04d48f14013123451abaaef2ee4bd7115bb1323 Mon Sep 17 00:00:00 2001 From: Daniel Date: Wed, 16 Sep 2026 12:18:22 +0200 Subject: [PATCH 7/7] fix test --- tests/testthat/test-get_predicted.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) 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",