From 7f5c28a60a9521824c8a7e133b1b777305892544 Mon Sep 17 00:00:00 2001 From: Jeffrey Girard Date: Mon, 12 Oct 2026 23:30:44 -0500 Subject: [PATCH 1/4] Show response level names for clmm models with the emmeans backend --- DESCRIPTION | 2 +- NEWS.md | 4 ++ R/get_emmeans.R | 14 ++++ tests/testthat/test-ordinal.R | 121 ++++++++++++++++++++++++++++++++++ 4 files changed, 140 insertions(+), 1 deletion(-) diff --git a/DESCRIPTION b/DESCRIPTION index e14b5c72a..fddf7b278 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,7 +1,7 @@ Type: Package Package: modelbased Title: Estimation of Model-Based Predictions, Contrasts and Means -Version: 0.17.0.9 +Version: 0.17.0.10 Authors@R: c(person(given = "Dominique", family = "Makowski", diff --git a/NEWS.md b/NEWS.md index e6e975aa5..c86f5ddcd 100644 --- a/NEWS.md +++ b/NEWS.md @@ -33,6 +33,10 @@ ## Bug fixes +* `estimate_means()` with `backend = "emmeans"` now shows the response level + names for models from `ordinal::clmm()`. Before, the `Response` column showed + `1`, `2`, `3`, and the `Threshold` column showed `1|2`, `2|3`. + * Fixed issue in `estimate_contrasts()` with wrong assignment of estimates in custom comparisons (e.g., when `comparison = "(b1 - b12) = (b4 - b18)"`), when some of the indexed coefficients contained the same digit (in this case, diff --git a/R/get_emmeans.R b/R/get_emmeans.R index ce0c5f07c..0ec93db65 100644 --- a/R/get_emmeans.R +++ b/R/get_emmeans.R @@ -303,6 +303,20 @@ get_emmeans <- function( pseudo_factor <- .emmeans_ordinal_pseudo_factor(model, ordinal_mode) if (!is.null(pseudo_factor) && pseudo_factor %in% colnames(means)) { new_name <- if (ordinal_mode == "prob") "Response" else "Threshold" + # for clmm models, emmeans codes the response categories and thresholds + # as 1, 2, 3, ... - we replace these codes by the names from the model + if (inherits(model, "clmm")) { + level_names <- if (ordinal_mode == "prob") model$y.levels else names(model$alpha) + codes <- as.character(means[[pseudo_factor]]) + if (ordinal_mode != "prob") { + # threshold codes look like "1|2" + codes <- sub("\\|.*", "", codes) + } + means[[pseudo_factor]] <- factor( + level_names[as.integer(codes)], + levels = level_names + ) + } colnames(means)[colnames(means) == pseudo_factor] <- new_name } means diff --git a/tests/testthat/test-ordinal.R b/tests/testthat/test-ordinal.R index af8ea05af..ce09d090f 100644 --- a/tests/testthat/test-ordinal.R +++ b/tests/testthat/test-ordinal.R @@ -238,6 +238,127 @@ test_that("estimate_means, backend emmeans, ordinal, transformed response", { }) +test_that("estimate_means, backend emmeans, clmm shows response level names", { + skip_if_not_installed("emmeans") + skip_if_not_installed("ordinal") + data(wine, package = "ordinal") + # emmeans codes the response categories and thresholds of clmm models as + # 1, 2, 3, ..., so we check that each code is replaced by its level name + text_levels <- c("vbad", "bad", "ok", "good", "vgood") + number_levels <- c("2", "4", "6", "8", "10") + wine$rating_text <- factor(wine$rating, labels = text_levels, ordered = TRUE) + wine$rating_number <- factor(as.integer(wine$rating) * 2, ordered = TRUE) + + for (response in c("rating_text", "rating_number")) { + response_levels <- if (response == "rating_text") text_levels else number_levels + threshold_names <- paste(response_levels[-5], response_levels[-1], sep = "|") + f <- stats::as.formula(paste(response, "~ temp + contact + (1 | judge)")) + m <- ordinal::clmm(f, data = wine) + expect_identical(m$y.levels, response_levels) + expect_identical(names(m$alpha), threshold_names) + + for (by in list("temp", c("temp", "contact"))) { + # response categories + out <- estimate_means(m, by, backend = "emmeans") + expected <- as.data.frame(suppressMessages(emmeans::emmeans( + m, + stats::as.formula(paste("~", paste(by, collapse = "+"), "|", response)), + mode = "prob" + ))) + expect_identical(nrow(out), nrow(expected)) + expect_identical( + as.character(out$Response), + response_levels[as.integer(as.character(expected[[response]]))] + ) + expect_identical(levels(out$Response), response_levels) + expect_equal(out$Probability, expected$prob, tolerance = 1e-6) + + # thresholds + modes <- c( + cum.prob = "Probability", + exc.prob = "Probability", + linear.predictor = "Linear_predictor" + ) + for (mode in names(modes)) { + out <- estimate_means(m, by, predict = mode, backend = "emmeans") + expected <- as.data.frame(suppressMessages(emmeans::emmeans( + m, + stats::as.formula(paste("~", paste(by, collapse = "+"), "| cut")), + mode = mode + ))) + code <- as.integer(sub("\\|.*", "", as.character(expected$cut))) + expect_identical(nrow(out), nrow(expected)) + expect_identical(as.character(out$Threshold), threshold_names[code]) + expect_identical(levels(out$Threshold), threshold_names) + estimate_column <- setdiff( + colnames(expected), + c(by, "cut", "SE", "df", "asymp.LCL", "asymp.UCL") + ) + expect_equal(out[[modes[[mode]]]], expected[[estimate_column]], tolerance = 1e-6) + } + } + } +}) + + +test_that("estimate_means, backend emmeans, ordinal level names unchanged", { + skip_if_not_installed("emmeans") + skip_if_not_installed("ordinal") + skip_if_not_installed("glmmTMB", minimum_version = "1.1.15.2") + data(housing, package = "MASS") + housing$SatNum <- as.integer(housing$Sat) + data(wine, package = "ordinal") + + sat_levels <- c("Low", "Medium", "High") + models <- list( + polr = list( + MASS::polr(Sat ~ Infl + Type + Cont, weights = Freq, data = housing, Hess = TRUE), + sat_levels + ), + polr_transformed = list( + MASS::polr(factor(SatNum) ~ Infl + Type + Cont, weights = Freq, data = housing, Hess = TRUE), + c("1", "2", "3") + ), + clm = list( + ordinal::clm(Sat ~ Infl + Type + Cont, weights = Freq, data = housing), + sat_levels + ), + glmmTMB = list( + glmmTMB::glmmTMB( + Sat ~ Infl + Type + Cont, + weights = Freq, + data = housing, + family = glmmTMB::ordinal() + ), + sat_levels + ), + clmm = list( + ordinal::clmm(rating ~ temp + contact + (1 | judge), data = wine), + c("1", "2", "3", "4", "5") + ) + ) + + for (i in names(models)) { + m <- models[[i]][[1]] + response_levels <- models[[i]][[2]] + by <- if (i == "clmm") "temp" else "Type" + k <- length(response_levels) + threshold_names <- paste(response_levels[-k], response_levels[-1], sep = "|") + + out <- estimate_means(m, by, predict = "prob", backend = "emmeans") + expect_identical(unique(as.character(out$Response)), response_levels, info = i) + for (mode in c("cum.prob", "exc.prob", "linear.predictor")) { + out <- estimate_means(m, by, predict = mode, backend = "emmeans") + expect_identical( + unique(as.character(out$Threshold)), + threshold_names, + info = paste(i, mode) + ) + } + } +}) + + test_that("estimate_means, print bracl", { skip_if_not_installed("brglm2") # required for the penguins dataset, which was added in R 4.5.0 From 508b57eed3f02d016a74e41e7f691eb48a2eb84e Mon Sep 17 00:00:00 2001 From: Jeffrey Girard Date: Mon, 12 Oct 2026 23:38:58 -0500 Subject: [PATCH 2/4] Build clmm threshold names from the response levels --- R/get_emmeans.R | 28 +++++++++++++++++--------- tests/testthat/test-ordinal.R | 37 +++++++++++++++++++++++++++++++++++ 2 files changed, 56 insertions(+), 9 deletions(-) diff --git a/R/get_emmeans.R b/R/get_emmeans.R index 0ec93db65..17069a906 100644 --- a/R/get_emmeans.R +++ b/R/get_emmeans.R @@ -306,16 +306,26 @@ get_emmeans <- function( # for clmm models, emmeans codes the response categories and thresholds # as 1, 2, 3, ... - we replace these codes by the names from the model if (inherits(model, "clmm")) { - level_names <- if (ordinal_mode == "prob") model$y.levels else names(model$alpha) - codes <- as.character(means[[pseudo_factor]]) - if (ordinal_mode != "prob") { - # threshold codes look like "1|2" - codes <- sub("\\|.*", "", codes) + response_levels <- model$y.levels + n_levels <- length(response_levels) + values <- as.character(means[[pseudo_factor]]) + if (ordinal_mode == "prob") { + level_names <- response_levels + codes <- as.character(seq_len(n_levels)) + } else { + # threshold codes look like "1|2". We don't use `names(model$alpha)`, + # because these are not the thresholds for `threshold = "symmetric"` + # or `"equidistant"` + level_names <- paste(response_levels[-n_levels], response_levels[-1], sep = "|") + codes <- paste(seq_len(n_levels - 1), seq_len(n_levels)[-1], sep = "|") + } + # only replace codes, in case emmeans returns level names in the future + if (setequal(values, codes)) { + means[[pseudo_factor]] <- factor( + level_names[match(values, codes)], + levels = level_names + ) } - means[[pseudo_factor]] <- factor( - level_names[as.integer(codes)], - levels = level_names - ) } colnames(means)[colnames(means) == pseudo_factor] <- new_name } diff --git a/tests/testthat/test-ordinal.R b/tests/testthat/test-ordinal.R index ce09d090f..8d306e441 100644 --- a/tests/testthat/test-ordinal.R +++ b/tests/testthat/test-ordinal.R @@ -301,6 +301,43 @@ test_that("estimate_means, backend emmeans, clmm shows response level names", { }) +test_that("estimate_means, backend emmeans, clmm threshold names, all threshold types", { + skip_if_not_installed("emmeans") + skip_if_not_installed("ordinal") + data(wine, package = "ordinal") + text_levels <- c("vbad", "bad", "ok", "good", "vgood") + threshold_names <- c("vbad|bad", "bad|ok", "ok|good", "good|vgood") + wine$rating_text <- factor(wine$rating, labels = text_levels, ordered = TRUE) + + # for "symmetric" and "equidistant", the model parameters are not one per + # threshold, but emmeans still returns one row per threshold + for (threshold in c("flexible", "symmetric", "equidistant")) { + m <- ordinal::clmm( + rating_text ~ temp + (1 | judge), + data = wine, + threshold = threshold + ) + for (mode in c("cum.prob", "exc.prob", "linear.predictor")) { + out <- estimate_means(m, "temp", predict = mode, backend = "emmeans") + expected <- as.data.frame(suppressMessages(emmeans::emmeans( + m, + ~ temp | cut, + mode = mode + ))) + code <- as.integer(sub("\\|.*", "", as.character(expected$cut))) + expect_identical( + as.character(out$Threshold), + threshold_names[code], + info = paste(threshold, mode) + ) + expect_identical(levels(out$Threshold), threshold_names, info = paste(threshold, mode)) + } + out <- estimate_means(m, "temp", backend = "emmeans") + expect_identical(levels(out$Response), text_levels, info = threshold) + } +}) + + test_that("estimate_means, backend emmeans, ordinal level names unchanged", { skip_if_not_installed("emmeans") skip_if_not_installed("ordinal") From 3e851e229f50c4a89c0e59f84b42b9e7a8a206e1 Mon Sep 17 00:00:00 2001 From: Jeffrey Girard Date: Mon, 12 Oct 2026 23:40:03 -0500 Subject: [PATCH 3/4] Test the symmetric2 threshold type for clmm --- R/get_emmeans.R | 4 ++-- tests/testthat/test-ordinal.R | 6 +++--- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/R/get_emmeans.R b/R/get_emmeans.R index 17069a906..1fb8eb10b 100644 --- a/R/get_emmeans.R +++ b/R/get_emmeans.R @@ -314,8 +314,8 @@ get_emmeans <- function( codes <- as.character(seq_len(n_levels)) } else { # threshold codes look like "1|2". We don't use `names(model$alpha)`, - # because these are not the thresholds for `threshold = "symmetric"` - # or `"equidistant"` + # because these are not the thresholds for `threshold = "symmetric"`, + # `"symmetric2"` or `"equidistant"` level_names <- paste(response_levels[-n_levels], response_levels[-1], sep = "|") codes <- paste(seq_len(n_levels - 1), seq_len(n_levels)[-1], sep = "|") } diff --git a/tests/testthat/test-ordinal.R b/tests/testthat/test-ordinal.R index 8d306e441..88a4687f0 100644 --- a/tests/testthat/test-ordinal.R +++ b/tests/testthat/test-ordinal.R @@ -309,9 +309,9 @@ test_that("estimate_means, backend emmeans, clmm threshold names, all threshold threshold_names <- c("vbad|bad", "bad|ok", "ok|good", "good|vgood") wine$rating_text <- factor(wine$rating, labels = text_levels, ordered = TRUE) - # for "symmetric" and "equidistant", the model parameters are not one per - # threshold, but emmeans still returns one row per threshold - for (threshold in c("flexible", "symmetric", "equidistant")) { + # for "symmetric", "symmetric2" and "equidistant", the model parameters are + # not one per threshold, but emmeans still returns one row per threshold + for (threshold in c("flexible", "symmetric", "symmetric2", "equidistant")) { m <- ordinal::clmm( rating_text ~ temp + (1 | judge), data = wine, From 7bf806461d39ce42826d6228eb0fbfde13e62a79 Mon Sep 17 00:00:00 2001 From: Jeffrey Girard Date: Mon, 12 Oct 2026 23:43:53 -0500 Subject: [PATCH 4/4] Guard clmm level mapping and test reversed numeric levels --- R/get_emmeans.R | 5 +++-- tests/testthat/test-ordinal.R | 34 +++++++++++++++++++++++++--------- 2 files changed, 28 insertions(+), 11 deletions(-) diff --git a/R/get_emmeans.R b/R/get_emmeans.R index 1fb8eb10b..ba4b47ced 100644 --- a/R/get_emmeans.R +++ b/R/get_emmeans.R @@ -319,8 +319,9 @@ get_emmeans <- function( level_names <- paste(response_levels[-n_levels], response_levels[-1], sep = "|") codes <- paste(seq_len(n_levels - 1), seq_len(n_levels)[-1], sep = "|") } - # only replace codes, in case emmeans returns level names in the future - if (setequal(values, codes)) { + # we only replace the values when they are exactly the codes, and when + # each code gets its own name (level names can contain "|") + if (setequal(values, codes) && !anyDuplicated(level_names)) { means[[pseudo_factor]] <- factor( level_names[match(values, codes)], levels = level_names diff --git a/tests/testthat/test-ordinal.R b/tests/testthat/test-ordinal.R index 88a4687f0..a2a93000e 100644 --- a/tests/testthat/test-ordinal.R +++ b/tests/testthat/test-ordinal.R @@ -246,11 +246,22 @@ test_that("estimate_means, backend emmeans, clmm shows response level names", { # 1, 2, 3, ..., so we check that each code is replaced by its level name text_levels <- c("vbad", "bad", "ok", "good", "vgood") number_levels <- c("2", "4", "6", "8", "10") + reversed_levels <- c("5", "4", "3", "2", "1") wine$rating_text <- factor(wine$rating, labels = text_levels, ordered = TRUE) wine$rating_number <- factor(as.integer(wine$rating) * 2, ordered = TRUE) + wine$rating_reversed <- factor( + as.character(wine$rating), + levels = reversed_levels, + ordered = TRUE + ) + all_levels <- list( + rating_text = text_levels, + rating_number = number_levels, + rating_reversed = reversed_levels + ) - for (response in c("rating_text", "rating_number")) { - response_levels <- if (response == "rating_text") text_levels else number_levels + for (response in names(all_levels)) { + response_levels <- all_levels[[response]] threshold_names <- paste(response_levels[-5], response_levels[-1], sep = "|") f <- stats::as.formula(paste(response, "~ temp + contact + (1 | judge)")) m <- ordinal::clmm(f, data = wine) @@ -341,7 +352,6 @@ test_that("estimate_means, backend emmeans, clmm threshold names, all threshold test_that("estimate_means, backend emmeans, ordinal level names unchanged", { skip_if_not_installed("emmeans") skip_if_not_installed("ordinal") - skip_if_not_installed("glmmTMB", minimum_version = "1.1.15.2") data(housing, package = "MASS") housing$SatNum <- as.integer(housing$Sat) data(wine, package = "ordinal") @@ -360,7 +370,17 @@ test_that("estimate_means, backend emmeans, ordinal level names unchanged", { ordinal::clm(Sat ~ Infl + Type + Cont, weights = Freq, data = housing), sat_levels ), - glmmTMB = list( + clmm = list( + ordinal::clmm(rating ~ temp + contact + (1 | judge), data = wine), + c("1", "2", "3", "4", "5") + ) + ) + # ordinal family in glmmTMB requires glmmTMB >= 1.1.15.2 + if ( + requireNamespace("glmmTMB", quietly = TRUE) && + utils::packageVersion("glmmTMB") >= "1.1.15.2" + ) { + models$glmmTMB <- list( glmmTMB::glmmTMB( Sat ~ Infl + Type + Cont, weights = Freq, @@ -368,12 +388,8 @@ test_that("estimate_means, backend emmeans, ordinal level names unchanged", { family = glmmTMB::ordinal() ), sat_levels - ), - clmm = list( - ordinal::clmm(rating ~ temp + contact + (1 | judge), data = wine), - c("1", "2", "3", "4", "5") ) - ) + } for (i in names(models)) { m <- models[[i]][[1]]