diff --git a/NEWS.md b/NEWS.md index 3277fcda8..7112835b2 100644 --- a/NEWS.md +++ b/NEWS.md @@ -2,6 +2,14 @@ ## Bug fixes +* If the filter in `by` keeps one level, `estimate_contrasts()` with + `backend = "marginaleffects"` now returns contrasts of slopes for ordinal + models, such as `MASS::polr()` or `ordinal::clm()`. An example is + `estimate_contrasts(m, "x", by = "Infl = 'Low'")`. The output compares the + slopes of the response categories at that level. Before, calls with + `estimate = "typical"` stopped with the error "No contrasts to show", and + calls with `estimate = "average"` ignored the filter. + * Filters on the response categories of frequentist ordinal models with `backend = "emmeans"`, where `Sat` is the response: * `estimate_means()`, `estimate_contrasts()` and `estimate_slopes()` now diff --git a/R/format.R b/R/format.R index 87a9556ea..d857db73f 100644 --- a/R/format.R +++ b/R/format.R @@ -345,23 +345,43 @@ format.marginaleffects_contrasts <- function( # If we have, for example, `contrast = c("vs", "am"), by = "gear='5'"`, the # `by` column is the one with one unique value only, we thus have to update # `by`, and also `contrast` (the latter not(!) for numerics)... + all_by <- by by <- by[lengths(lapply(dgrid[by], unique)) > 1] # for slopes, we have our "levels" in the by-variable, because `by` indicates # the levels for which slopes are contrasted - thus, we replace contrast with by + one_level_slopes <- FALSE if (inherits(x, "estimate_slopes")) { contrast <- by + # models with response categories (e.g., ordinal models) have one slope + # per category. When the filter in `by` keeps one level, e.g. + # `by = "Infl='Low'"`, the slopes of the categories at that level are + # still contrasted. The labels from marginaleffects then hold only the + # response categories, e.g. "(Medium) - (Low)". Custom hypotheses, such + # as "b1=b2", and group formulas, such as `~pairwise | Infl`, are not + # labelled this way, so these still stop with an error below + one_level_slopes <- !length(contrast) && + length(all_by) > 0 && + all(all_by %in% colnames(dgrid)) && + all(lengths(lapply(dgrid[all_by], unique)) == 1) && + !.is_custom_comparison(comparison) && + !(inherits(comparison, "formula") && + grepl("|", insight::safe_deparse(comparison), fixed = TRUE)) && + .is_response_level_model(model, x) + if (one_level_slopes) { + contrast <- all_by + } by <- NULL + } else { + # for contrasts, we also filter variables with one unique value, but we + # keep numeric variables. When these are hold constant in the data grid, + # they are set to their mean value - meaning, they only have one unique + # value in the data grid, anyway. so we need to keep them. + keep_contrasts <- lengths(lapply(dgrid[contrast], unique)) > 1 | + vapply(dgrid[contrast], is.numeric, logical(1)) # nolint + contrast <- contrast[keep_contrasts] } - # for contrasts, we also filter variables with one unique value, but we - # keep numeric variables. When these are hold constant in the data grid, - # they are set to their mean value - meaning, they only have one unique - # value in the data grid, anyway. so we need to keep them. - keep_contrasts <- lengths(lapply(dgrid[contrast], unique)) > 1 | - vapply(dgrid[contrast], is.numeric, logical(1)) # nolint - contrast <- contrast[keep_contrasts] - # set to NULL, if all by-values have been removed here if (!length(by)) { by <- NULL @@ -373,8 +393,12 @@ format.marginaleffects_contrasts <- function( insight::format_error("No contrasts to show. Please adjust `contrast`.") } - # for more than one term, we have comma-separated levels. - if (length(contrast) > 1) { + # for slopes with a one-level filter in `by`, the two columns hold the + # contrasted response categories. For more than one term, we have + # comma-separated levels. + if (one_level_slopes) { + new_colnames <- c("Response1", "Response2") + } else if (length(contrast) > 1) { # levels may contain the separator char. to be 100% certain we extract # levels correctly, we now replace levels with a special "token", and later # replace those tokens with the original levels again @@ -487,8 +511,11 @@ format.marginaleffects_contrasts <- function( # all contrast terms that belong to one "contrast group", separated # by comma, and each the two "new contrast groups" go into separate # columns named "Level1" and "Level2". For one contrast term, we just - # need to rename the columns - if (length(contrast) == 1) { + # need to rename the columns. For slopes with a one-level filter in `by`, + # "Level1" and "Level2" get that level, next to the response categories + if (one_level_slopes) { + params <- .one_level_slope_labels(params, model, dgrid[all_by]) + } else if (length(contrast) == 1) { # we define this object to avoid scoping issues in data_rename(), # see https://github.com/easystats/modelbased/issues/401 rename_columns <- paste0(contrast, 1:2) @@ -522,13 +549,17 @@ format.marginaleffects_contrasts <- function( # level and the contrasted factor level into the same string (e.g., # "setosa 3"). We split those into separate "Response1" and "Response2" # columns, so contrasted factor levels and response levels are not - # confounded - params <- .split_contrast_response_levels(params, model, x) + # confounded. Slopes with a one-level filter in `by` already have these + # columns, because their labels hold only the response categories + if (!one_level_slopes) { + params <- .split_contrast_response_levels(params, model, x) + } # we need to update these variables, because these are the new column # names for contrasts and focal terms. We use `intersect()` here because # "Response1"/"Response2" are only present if the contrasted factor and - # response category levels were successfully split above + # response category levels were successfully split above, or for slopes + # with a one-level filter in `by` contrast <- intersect( c("Level1", "Response1", "Level2", "Response2"), colnames(params) @@ -594,6 +625,61 @@ format.marginaleffects_contrasts <- function( identical(comparison, "ratio") } +# TRUE for models where predictions are made separately for each level of the +# response (categorical, ordinal/cumulative, multinomial or multivariate models) + +.is_response_level_model <- function(model, x = NULL) { + if (is.null(model)) { + return(FALSE) + } + # prefer the model info already stored as attribute (consistent with how + # other format-methods retrieve model information), and only compute it + # ourselves as a fallback + model_info <- attr(x, "model_info", exact = TRUE) + if (is.null(model_info)) { + model_info <- .safe(insight::model_info(model, verbose = FALSE)) + } + isTRUE(model_info$is_multinomial) || + isTRUE(model_info$is_categorical) || + isTRUE(model_info$is_ordinal) || + isTRUE(model_info$is_cumulative) || + isTRUE(.safe(insight::is_multivariate(model))) +} + +# for contrasts of slopes of models with response categories, when the filter +# in `by` keeps one level: `params` has the response categories in "Response1" +# and "Response2", and "Level1" and "Level2" get the one level of `by`, from +# the data grid. Several `by` variables are combined with ", ", as in +# contrasts of more than one term + +.one_level_slope_labels <- function(params, model, by_grid) { + by_level <- paste( + vapply(by_grid, function(i) as.character(unique(i)), character(1)), + collapse = ", " + ) + # keep the order of the response levels, as in + # `.split_contrast_response_levels()`. Labels that are no response level, + # e.g. because marginaleffects removes parentheses from the labels, follow + # in the order in which they appear, so that no label becomes `NA` + labels <- unique(c(params$Response1, params$Response2)) + response <- .safe(insight::get_response(model, verbose = FALSE)) + if (is.factor(response)) { + known_levels <- levels(response) + } else if (is.null(dim(response))) { + known_levels <- unique(as.character(response)) + } else { + known_levels <- colnames(response) + } + response_levels <- c(intersect(known_levels, labels), setdiff(labels, known_levels)) + data.frame( + Level1 = factor(by_level), + Response1 = factor(params$Response1, levels = response_levels), + Level2 = factor(by_level), + Response2 = factor(params$Response2, levels = response_levels), + stringsAsFactors = FALSE + ) +} + # since we combine levels from different factors, we have to make # sure levels are unique across different terms. If not, paste # variable names to levels. We first find the intersection of all @@ -613,23 +699,10 @@ format.marginaleffects_contrasts <- function( return(params) } - # prefer the model info already stored as attribute (consistent with how - # other format-methods retrieve model information), and only compute it - # ourselves as a fallback - model_info <- attr(x, "model_info", exact = TRUE) - if (is.null(model_info)) { - model_info <- .safe(insight::model_info(model, verbose = FALSE)) - } - is_multivariate <- isTRUE(.safe(insight::is_multivariate(model))) - is_response_model <- isTRUE(model_info$is_multinomial) || - isTRUE(model_info$is_categorical) || - isTRUE(model_info$is_ordinal) || - isTRUE(model_info$is_cumulative) || - is_multivariate - - if (!is_response_model) { + if (!.is_response_level_model(model, x)) { return(params) } + is_multivariate <- isTRUE(.safe(insight::is_multivariate(model))) # response levels, used to detect and split off the response-category # label from the contrasted factor level. We keep the original order of diff --git a/tests/testthat/test-ordinal.R b/tests/testthat/test-ordinal.R index 8ee10011e..a28649be0 100644 --- a/tests/testthat/test-ordinal.R +++ b/tests/testthat/test-ordinal.R @@ -1338,6 +1338,76 @@ test_that("estimate_means, backend emmeans, ordinal, several by terms on the res }) +test_that("estimate_contrasts, ordinal slopes, by filter that keeps one level", { + skip_if_not_installed("ordinal") + d <- MASS::housing + set.seed(1) + d$x <- rnorm(nrow(d)) + models <- list( + polr = MASS::polr(Sat ~ x * Infl + Type, data = d, Hess = TRUE), + clm = ordinal::clm(Sat ~ x * Infl + Type, data = d) + ) + keys <- c("Level1", "Response1", "Level2", "Response2") + + for (m in models) { + for (estimate in c("typical", "average")) { + out_all <- estimate_contrasts(m, "x", by = "Infl", estimate = estimate) + + # one level: the contrasts between the response categories at that level + out <- estimate_contrasts(m, "x", by = "Infl='Low'", estimate = estimate) + expected <- out_all[out_all$Level1 == "Low" & out_all$Level2 == "Low", ] + expect_identical(nrow(out), 3L) + expect_identical( + lapply(as.data.frame(out)[keys], as.character), + lapply(as.data.frame(expected)[keys], as.character) + ) + expect_equal(out$Difference, expected$Difference, tolerance = 1e-8) + expect_equal(out$SE, expected$SE, tolerance = 1e-8) + + # two levels: the matching rows of the unfiltered output + out <- estimate_contrasts(m, "x", by = "Infl=c('Low','High')", estimate = estimate) + expect_identical(nrow(out), 15L) + out_keys <- do.call(paste, lapply(as.data.frame(out)[keys], as.character)) + all_keys <- do.call(paste, lapply(as.data.frame(out_all)[keys], as.character)) + matched <- out_all[match(out_keys, all_keys), ] + expect_false(anyNA(match(out_keys, all_keys))) + expect_equal(out$Difference, matched$Difference, tolerance = 1e-8) + expect_equal(out$SE, matched$SE, tolerance = 1e-8) + } + } + + # response levels with parentheses keep a label + d2 <- d + levels(d2$Sat) <- c("Low (bad)", "Mid", "High (good)") + m <- MASS::polr(Sat ~ x * Infl + Type, data = d2, Hess = TRUE) + out <- estimate_contrasts(m, "x", by = "Infl='Low'") + expect_identical(nrow(out), 3L) + expect_false(anyNA(out$Response1)) + expect_false(anyNA(out$Response2)) + expect_length(unique(c(as.character(out$Response1), as.character(out$Response2))), 3L) + + # custom hypotheses and group formulas are not labelled for one level + expect_error( + estimate_contrasts(models$polr, "x", by = "Infl='Low'", comparison = "b1=b2"), + "No contrasts to show", + fixed = TRUE + ) + expect_error( + estimate_contrasts(models$polr, "x", by = "Infl='Low'", comparison = ~ pairwise | Infl), + "No contrasts to show", + fixed = TRUE + ) + + # models without response categories still need two slopes + mt <- transform(mtcars, g = factor(am)) + expect_error( + estimate_contrasts(lm(mpg ~ wt * g, data = mt), "wt", by = "g='0'"), + "The filter in `by` leaves fewer than two slopes to compare", + fixed = TRUE + ) +}) + + # Wald statistic of probability contrasts of a polr model, by the delta method. # For `polr(Sat ~ Infl + Type)`, P(Sat <= k) = plogis(zeta_k - eta), with eta # the linear predictor. The category probabilities are differences of these