From c0cde7b4c0e2b4f250331001464df0a3709bb159 Mon Sep 17 00:00:00 2001 From: Jeffrey Girard Date: Mon, 19 Oct 2026 19:02:39 -0500 Subject: [PATCH 1/8] Add tests for latent-scale effect sizes of ordinal models and for effect sizes of reordered contrasts --- ...test-estimate_contrasts_effectsize_order.R | 14 ++ tests/testthat/test-ordinal.R | 151 ++++++++++++++++++ 2 files changed, 165 insertions(+) create mode 100644 tests/testthat/test-estimate_contrasts_effectsize_order.R diff --git a/tests/testthat/test-estimate_contrasts_effectsize_order.R b/tests/testthat/test-estimate_contrasts_effectsize_order.R new file mode 100644 index 000000000..c1343109b --- /dev/null +++ b/tests/testthat/test-estimate_contrasts_effectsize_order.R @@ -0,0 +1,14 @@ +skip_on_cran() +skip_if_not_installed("emmeans") + +test_that("estimate_contrasts, effectsize emmeans, effect sizes follow their contrasts", { + # the levels are not in alphabetical order, so the contrast rows are sorted + # in another order than the rows of emmeans + d <- iris + d$Species <- factor(d$Species, levels = c("virginica", "setosa", "versicolor")) + g <- lm(Sepal.Length ~ Species, data = d) + out <- estimate_contrasts(g, "Species", effectsize = "emmeans", backend = "emmeans") + expect_identical(nrow(out), 3L) + # e.g. versicolor - setosa: 0.930 / sigma(g) = 0.930 / 0.5148 = 1.807 + expect_equal(out$partial_d, out$Difference / sigma(g), tolerance = 1e-6) +}) diff --git a/tests/testthat/test-ordinal.R b/tests/testthat/test-ordinal.R index a28649be0..9863151a1 100644 --- a/tests/testthat/test-ordinal.R +++ b/tests/testthat/test-ordinal.R @@ -1986,3 +1986,154 @@ test_that("estimate_contrasts, backend emmeans, ordinal, slopes with a by filter c("Difference", "SE") ) }) + + +test_that("estimate_contrasts, backend emmeans, ordinal, latent-scale effect sizes", { + skip_if_not_installed("emmeans") + skip_if_not_installed("ordinal") + data(housing, package = "MASS") + + # the latent errors of the logit link follow the standard logistic + # distribution, with variance pi^2 / 3 (?stats::Logistic), and those of the + # probit link the standard normal distribution, with SD 1 (?stats::Normal) + sd_logit <- pi / sqrt(3) + + m <- MASS::polr(Sat ~ Infl + Type, weights = Freq, data = housing, Hess = TRUE) + out <- estimate_contrasts( + m, + "Infl", + predict = "latent", + effectsize = "emmeans", + backend = "emmeans" + ) + expect_identical(nrow(out), 3L) + # sigma is taken as known (`edf = Inf`), so the effect sizes and their + # confidence limits are the differences and their limits divided by sigma + expect_equal(out$partial_d, out$Difference / sd_logit, tolerance = 1e-6) + expect_equal(out$es_CI_low, out$CI_low / sd_logit, tolerance = 1e-6) + expect_equal(out$es_CI_high, out$CI_high / sd_logit, tolerance = 1e-6) + # with no interaction, the latent differences to `Low` are the coefficients, + # e.g. coef(m)["InflHigh"] / (pi / sqrt(3)) = 1.2373 / 1.8138 = 0.6822 + row_ml <- out$Level1 == "Medium" & out$Level2 == "Low" + row_hl <- out$Level1 == "High" & out$Level2 == "Low" + expect_equal( + out$partial_d[row_ml], + unname(coef(m)["InflMedium"]) / sd_logit, + tolerance = 1e-6 + ) + expect_equal( + out$partial_d[row_hl], + unname(coef(m)["InflHigh"]) / sd_logit, + tolerance = 1e-6 + ) + + mp <- MASS::polr( + Sat ~ Infl + Type, + weights = Freq, + data = housing, + Hess = TRUE, + method = "probit" + ) + out_p <- estimate_contrasts( + mp, + "Infl", + predict = "latent", + effectsize = "emmeans", + backend = "emmeans" + ) + expect_identical(nrow(out_p), 3L) + expect_equal(out_p$partial_d, out_p$Difference, tolerance = 1e-6) + expect_equal(out_p$es_CI_low, out_p$CI_low, tolerance = 1e-6) + expect_equal(out_p$es_CI_high, out_p$CI_high, tolerance = 1e-6) + row_ml <- out_p$Level1 == "Medium" & out_p$Level2 == "Low" + row_hl <- out_p$Level1 == "High" & out_p$Level2 == "Low" + expect_equal(out_p$partial_d[row_ml], unname(coef(mp)["InflMedium"]), tolerance = 1e-6) + expect_equal(out_p$partial_d[row_hl], unname(coef(mp)["InflHigh"]), tolerance = 1e-6) + + # clm fits of the same models give the same effect sizes + mc <- ordinal::clm(Sat ~ Infl + Type, weights = Freq, data = housing) + out_c <- estimate_contrasts( + mc, + "Infl", + predict = "latent", + effectsize = "emmeans", + backend = "emmeans" + ) + expect_identical(nrow(out_c), 3L) + expect_identical(as.character(out_c$Level1), as.character(out$Level1)) + expect_identical(as.character(out_c$Level2), as.character(out$Level2)) + expect_equal(out_c$partial_d, out$partial_d, tolerance = 1e-3) + + mcp <- ordinal::clm(Sat ~ Infl + Type, weights = Freq, data = housing, link = "probit") + out_cp <- estimate_contrasts( + mcp, + "Infl", + predict = "latent", + effectsize = "emmeans", + backend = "emmeans" + ) + expect_identical(as.character(out_cp$Level1), as.character(out_p$Level1)) + expect_identical(as.character(out_cp$Level2), as.character(out_p$Level2)) + expect_equal(out_cp$partial_d, out_p$partial_d, tolerance = 1e-3) +}) + + +test_that("estimate_contrasts, backend emmeans, ordinal, effect size errors", { + skip_if_not_installed("emmeans") + skip_if_not_installed("ordinal") + data(housing, package = "MASS") + m <- MASS::polr(Sat ~ Infl + Type, weights = Freq, data = housing, Hess = TRUE) + mc <- ordinal::clm(Sat ~ Infl + Type, weights = Freq, data = housing) + + # effect sizes need the latent scale + for (model in list(m, mc)) { + err <- expect_error( + estimate_contrasts(model, "Infl", effectsize = "emmeans", backend = "emmeans") + ) + expect_match(conditionMessage(err), "predict = \"latent\"", fixed = TRUE) + } + + # only the logit and probit links have a known latent SD here + mcl <- MASS::polr( + Sat ~ Infl + Type, + weights = Freq, + data = housing, + Hess = TRUE, + method = "cloglog" + ) + mcc <- ordinal::clm(Sat ~ Infl + Type, weights = Freq, data = housing, link = "cloglog") + for (model in list(mcl, mcc)) { + err <- expect_error( + estimate_contrasts( + model, + "Infl", + predict = "latent", + effectsize = "emmeans", + backend = "emmeans" + ) + ) + expect_match(conditionMessage(err), "logit", fixed = TRUE) + expect_match(conditionMessage(err), "probit", fixed = TRUE) + } + + # a scale formula makes the latent SD differ between observations + ms <- ordinal::clm(Sat ~ Infl + Type, scale = ~Cont, weights = Freq, data = housing) + err <- expect_error( + estimate_contrasts( + ms, + "Infl", + predict = "latent", + effectsize = "emmeans", + backend = "emmeans" + ) + ) + expect_match(conditionMessage(err), "scale", fixed = TRUE) + + # the other effect sizes still need a Gaussian model + for (es in c("marginal", "boot")) { + err <- expect_error( + estimate_contrasts(m, "Infl", predict = "latent", effectsize = es, backend = "emmeans") + ) + expect_match(conditionMessage(err), "non-Gaussian", fixed = TRUE) + } +}) From 3991d8041f77af3bb4cd5452ed7f58c4bcb337fb Mon Sep 17 00:00:00 2001 From: Jeffrey Girard Date: Mon, 19 Oct 2026 19:02:39 -0500 Subject: [PATCH 2/8] Latent-scale effect sizes for polr and clm models, and effect sizes matched to their contrasts --- R/estimate_contrasts.R | 1 + R/estimate_contrasts_effectsize.R | 54 +++++++++++++++++++++++++++++-- R/get_emcontrasts.R | 12 ++++++- 3 files changed, 63 insertions(+), 4 deletions(-) diff --git a/R/estimate_contrasts.R b/R/estimate_contrasts.R index 50427fecb..d16a59796 100644 --- a/R/estimate_contrasts.R +++ b/R/estimate_contrasts.R @@ -564,6 +564,7 @@ estimate_contrasts.default <- function( backend = backend ) } + attr(out, "emmeans_rows") <- NULL # sanity check - did method return standard errors? .check_standard_errors( diff --git a/R/estimate_contrasts_effectsize.R b/R/estimate_contrasts_effectsize.R index 242cb9d7a..c8268615a 100644 --- a/R/estimate_contrasts_effectsize.R +++ b/R/estimate_contrasts_effectsize.R @@ -33,11 +33,15 @@ )) } + # For ordinal models, `effectsize = "emmeans"` standardizes the contrasts of + # the latent variable by the SD of its error distribution + latent_ordinal <- effectsize == "emmeans" && inherits(model, c("polr", "clm")) + # Check if the model's response variable follows a Gaussian (normal) # distribution. Effect size calculations implemented in this function are # designed for Gaussian models and may not be appropriate or meaningful for # other types of distributions. - if (insight::get_family(model)$family != "gaussian") { + if (!latent_ordinal && insight::get_family(model)$family != "gaussian") { insight::format_error(paste0( "We strongly recommend not using the `effectsize` ", "argument with non-Gaussian models." @@ -47,15 +51,28 @@ switch( effectsize, emmeans = { + if (latent_ordinal) { + sigma <- .ordinal_latent_sigma(model, attributes(estimated)$predict) + edf <- Inf + } else { + sigma <- insight::get_sigma(model, verbose = FALSE) + edf <- insight::get_df(model) + } eff <- emmeans::eff_size( estimated, - sigma = insight::get_sigma(model, verbose = FALSE), - edf = insight::get_df(model), + sigma = sigma, + edf = edf, method = "identity" ) eff <- as.data.frame(eff) eff <- eff[c(2, 5:6)] names(eff) <- c("partial_d", "es_CI_low", "es_CI_high") + # the contrasts are sorted in another order than the emmeans grid + emmeans_rows <- attributes(contrasts_results)$emmeans_rows + if (!is.null(emmeans_rows)) { + eff <- eff[emmeans_rows, , drop = FALSE] + rownames(eff) <- NULL + } contrasts_results <- cbind(contrasts_results, eff) }, marginal = { @@ -112,6 +129,37 @@ } +# SD of the error distribution of the latent variable of ordinal models, which +# is the standard logistic (SD pi / sqrt(3)) for the logit link, and the +# standard normal (SD 1) for the probit link +.ordinal_latent_sigma <- function(model, predict) { + if (!identical(predict, "latent")) { + insight::format_error(paste0( + "Effect sizes for ordinal models are only possible on the latent scale. ", + "Please use `predict = \"latent\"`." + )) + } + if (inherits(model, "clm") && !is.null(model$formulas$scale)) { + insight::format_error( + "Effect sizes are not possible for ordinal models with a `scale` formula." + ) + } + link <- if (inherits(model, "polr")) model$method else model$link + switch( + link, + logistic = , + logit = pi / sqrt(3), + probit = 1, + insight::format_error(paste0( + "Effect sizes for ordinal models are only possible for the `logit` and ", + "`probit` links, but the model uses the `", + link, + "` link." + )) + ) +} + + .get_group_variable <- function(estimated, backend) { if (backend == "emmeans") { names(estimated@model.info$xlev) diff --git a/R/get_emcontrasts.R b/R/get_emcontrasts.R index 417a5466a..b3303ed15 100644 --- a/R/get_emcontrasts.R +++ b/R/get_emcontrasts.R @@ -267,6 +267,7 @@ get_emcontrasts <- function( .format_emmeans_contrasts <- function(model, estimated, ci, p_adjust, ...) { predict <- attributes(estimated)$predict m_info <- insight::model_info(model, response = 1) + emmeans_rows <- NULL # Summarize and clean if (m_info$is_bayesian) { @@ -287,6 +288,9 @@ get_emcontrasts <- function( out <- as.data.frame(estimated) } else { estimated_data <- as.data.frame(estimated) + # `merge()` sorts the rows, so we keep the row of each contrast in the + # emmeans grid, to match other emmeans results (e.g. effect sizes) + estimated_data$.emmeans_row <- seq_len(nrow(estimated_data)) out <- as.data.frame(merge( estimated_data, stats::confint(estimated, level = ci, adjust = p_adjust) @@ -305,6 +309,8 @@ get_emcontrasts <- function( out <- out[match(row_keys(estimated_data), row_keys(out)), , drop = FALSE] rownames(out) <- NULL } + emmeans_rows <- out$.emmeans_row + out$.emmeans_row <- NULL } # the response categories or thresholds are renamed first, because the # response can have the name of an output column (e.g. "response") @@ -367,5 +373,9 @@ get_emcontrasts <- function( } # add posterior draws? - .add_posterior_draws_emmeans(attributes(estimated), out) + out <- .add_posterior_draws_emmeans(attributes(estimated), out) + if (!is.null(emmeans_rows)) { + attr(out, "emmeans_rows") <- emmeans_rows + } + out } From 45cd6ac7b9f07dd4c098fc9b11b4342273f12de4 Mon Sep 17 00:00:00 2001 From: Jeffrey Girard Date: Mon, 19 Oct 2026 19:06:22 -0500 Subject: [PATCH 3/8] Document latent-scale effect sizes for ordinal models, add NEWS, bump version --- DESCRIPTION | 2 +- NEWS.md | 16 ++++++++++++++++ R/estimate_contrasts.R | 12 ++++++++++-- man/estimate_contrasts.Rd | 12 ++++++++++-- 4 files changed, 37 insertions(+), 5 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index 72c0fb47f..6c94bc4dd 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,7 +1,7 @@ Type: Package Package: modelbased Title: Estimation of Model-Based Predictions, Contrasts and Means -Version: 0.18.0.3 +Version: 0.18.0.4 Authors@R: c(person(given = "Dominique", family = "Makowski", diff --git a/NEWS.md b/NEWS.md index 7112835b2..bd6c729b0 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,5 +1,16 @@ # modelbased (devel) +## Changes + +* `estimate_contrasts()` with `backend = "emmeans"` and + `effectsize = "emmeans"` now returns effect sizes for `polr` and `clm` + models with `predict = "latent"` and a logit or probit link. The contrasts + of the latent variable are divided by the SD of its error distribution: + `pi / sqrt(3)` for the logit link and `1` for the probit link. Other values + of `predict`, other links and `clm` models with a `scale` formula give an + error. `effectsize = "marginal"` and `effectsize = "boot"` still give an + error for ordinal models. + ## Bug fixes * If the filter in `by` keeps one level, `estimate_contrasts()` with @@ -10,6 +21,11 @@ `estimate = "typical"` stopped with the error "No contrasts to show", and calls with `estimate = "average"` ignored the filter. +* `estimate_contrasts()` with `effectsize = "emmeans"` now puts each effect + size on the row of its own contrast. Before, if the factor levels were not + in alphabetical order, the contrasts were sorted in another order than in + *emmeans*, and some effect sizes were on the wrong rows. + * 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/estimate_contrasts.R b/R/estimate_contrasts.R index d16a59796..e037cb0fb 100644 --- a/R/estimate_contrasts.R +++ b/R/estimate_contrasts.R @@ -103,7 +103,9 @@ #' specified post-process tests. See 'Examples'. #' @param effectsize Desired measure of standardized effect size, one of #' `"emmeans"`, `"marginal"`, or `"boot"`. Default is `NULL`, i.e. no effect -#' size will be computed. +#' size will be computed. Effect sizes need a Gaussian model without random +#' effects. Only `"emmeans"` also works for `polr` and `clm` models with +#' `predict = "latent"` and a logit or probit link. #' @param es_type Specifies the type of effect-size measure to estimate when #' using `effectsize = "boot"`. One of `"unstandardized"`, `"cohens.d"`, #' `"hedges.g"`, `"cohens.d.sigma"`, `"r"`, or `"akp.robust.d"`. See `effect.type` @@ -247,7 +249,13 @@ #' this works when the contrasts are the only predictors in the model, but not #' when there are covariates. The response variance accounted for by the #' covariates should not be removed from the SD used to standardize. Otherwise, -#' _d_ will be overestimated. +#' _d_ will be overestimated. For `polr` and `clm` models with +#' `predict = "latent"`, sigma is the SD of the error distribution of the +#' latent variable: `pi / sqrt(3)` (the standard logistic distribution) for +#' the logit link, and `1` (the standard normal distribution) for the probit +#' link. Sigma is taken as known, with `edf = Inf`, so the effect sizes and +#' their confidence limits are the contrasts and their limits divided by +#' sigma. #' #' `effectsize = "marginal"` uses the following formula to compute effect #' size: `d_adj <- difference * (1- R2)/ sigma`. This standardizes diff --git a/man/estimate_contrasts.Rd b/man/estimate_contrasts.Rd index 645aad7da..898c6c91a 100644 --- a/man/estimate_contrasts.Rd +++ b/man/estimate_contrasts.Rd @@ -319,7 +319,9 @@ to the output. You can reshape them to a long format by running \item{effectsize}{Desired measure of standardized effect size, one of \code{"emmeans"}, \code{"marginal"}, or \code{"boot"}. Default is \code{NULL}, i.e. no effect -size will be computed.} +size will be computed. Effect sizes need a Gaussian model without random +effects. Only \code{"emmeans"} also works for \code{polr} and \code{clm} models with +\code{predict = "latent"} and a logit or probit link.} \item{iterations}{This argument has two distinct effects, depending on the context. When \code{effectsize = "boot"}, it defines the number of bootstrap @@ -543,7 +545,13 @@ aren't part of the contrast. this works when the contrasts are the only predictors in the model, but not when there are covariates. The response variance accounted for by the covariates should not be removed from the SD used to standardize. Otherwise, -\emph{d} will be overestimated. +\emph{d} will be overestimated. For \code{polr} and \code{clm} models with +\code{predict = "latent"}, sigma is the SD of the error distribution of the +latent variable: \code{pi / sqrt(3)} (the standard logistic distribution) for +the logit link, and \code{1} (the standard normal distribution) for the probit +link. Sigma is taken as known, with \code{edf = Inf}, so the effect sizes and +their confidence limits are the contrasts and their limits divided by +sigma. \code{effectsize = "marginal"} uses the following formula to compute effect size: \code{d_adj <- difference * (1- R2)/ sigma}. This standardizes From 6190710be674fe8a971d02ee6d7b384ddac27381 Mon Sep 17 00:00:00 2001 From: Jeffrey Girard Date: Mon, 19 Oct 2026 19:11:47 -0500 Subject: [PATCH 4/8] Select emmeans effect-size columns by name, so a by variable no longer shifts them; correct effect-size docs --- NEWS.md | 4 +++- R/estimate_contrasts.R | 10 ++++++---- R/estimate_contrasts_effectsize.R | 8 ++++++-- man/estimate_contrasts.Rd | 10 ++++++---- .../test-estimate_contrasts_effectsize_order.R | 13 +++++++++++++ tests/testthat/test-ordinal.R | 18 +++++++++++++++++- 6 files changed, 51 insertions(+), 12 deletions(-) diff --git a/NEWS.md b/NEWS.md index bd6c729b0..d88d165fb 100644 --- a/NEWS.md +++ b/NEWS.md @@ -24,7 +24,9 @@ * `estimate_contrasts()` with `effectsize = "emmeans"` now puts each effect size on the row of its own contrast. Before, if the factor levels were not in alphabetical order, the contrasts were sorted in another order than in - *emmeans*, and some effect sizes were on the wrong rows. + *emmeans*, and some effect sizes were on the wrong rows. With a `by` + variable, the effect-size columns held the values of the `by` variable and + of other *emmeans* columns instead of the effect sizes. * Filters on the response categories of frequentist ordinal models with `backend = "emmeans"`, where `Sat` is the response: diff --git a/R/estimate_contrasts.R b/R/estimate_contrasts.R index e037cb0fb..a28bdd933 100644 --- a/R/estimate_contrasts.R +++ b/R/estimate_contrasts.R @@ -105,7 +105,8 @@ #' `"emmeans"`, `"marginal"`, or `"boot"`. Default is `NULL`, i.e. no effect #' size will be computed. Effect sizes need a Gaussian model without random #' effects. Only `"emmeans"` also works for `polr` and `clm` models with -#' `predict = "latent"` and a logit or probit link. +#' `predict = "latent"` and a logit or probit link, but not for `clm` models +#' with a `scale` formula. #' @param es_type Specifies the type of effect-size measure to estimate when #' using `effectsize = "boot"`. One of `"unstandardized"`, `"cohens.d"`, #' `"hedges.g"`, `"cohens.d.sigma"`, `"r"`, or `"akp.robust.d"`. See `effect.type` @@ -253,9 +254,10 @@ #' `predict = "latent"`, sigma is the SD of the error distribution of the #' latent variable: `pi / sqrt(3)` (the standard logistic distribution) for #' the logit link, and `1` (the standard normal distribution) for the probit -#' link. Sigma is taken as known, with `edf = Inf`, so the effect sizes and -#' their confidence limits are the contrasts and their limits divided by -#' sigma. +#' link. Sigma is taken as known, with `edf = Inf`, so the effect sizes are +#' the contrasts divided by sigma. Their confidence limits are the unadjusted +#' confidence limits of the contrasts divided by sigma, also when `p_adjust` +#' is not `"none"`. #' #' `effectsize = "marginal"` uses the following formula to compute effect #' size: `d_adj <- difference * (1- R2)/ sigma`. This standardizes diff --git a/R/estimate_contrasts_effectsize.R b/R/estimate_contrasts_effectsize.R index c8268615a..0478801a1 100644 --- a/R/estimate_contrasts_effectsize.R +++ b/R/estimate_contrasts_effectsize.R @@ -65,9 +65,13 @@ method = "identity" ) eff <- as.data.frame(eff) - eff <- eff[c(2, 5:6)] + # select columns by name, because `by` variables add columns + es_col <- intersect(c("effect.size", "estimate"), colnames(eff))[1] + ci_low <- intersect(c("lower.CL", "asymp.LCL"), colnames(eff))[1] + ci_high <- intersect(c("upper.CL", "asymp.UCL"), colnames(eff))[1] + eff <- eff[c(es_col, ci_low, ci_high)] names(eff) <- c("partial_d", "es_CI_low", "es_CI_high") - # the contrasts are sorted in another order than the emmeans grid + # the contrasts can be sorted in another order than the emmeans grid emmeans_rows <- attributes(contrasts_results)$emmeans_rows if (!is.null(emmeans_rows)) { eff <- eff[emmeans_rows, , drop = FALSE] diff --git a/man/estimate_contrasts.Rd b/man/estimate_contrasts.Rd index 898c6c91a..48fefccec 100644 --- a/man/estimate_contrasts.Rd +++ b/man/estimate_contrasts.Rd @@ -321,7 +321,8 @@ to the output. You can reshape them to a long format by running \code{"emmeans"}, \code{"marginal"}, or \code{"boot"}. Default is \code{NULL}, i.e. no effect size will be computed. Effect sizes need a Gaussian model without random effects. Only \code{"emmeans"} also works for \code{polr} and \code{clm} models with -\code{predict = "latent"} and a logit or probit link.} +\code{predict = "latent"} and a logit or probit link, but not for \code{clm} models +with a \code{scale} formula.} \item{iterations}{This argument has two distinct effects, depending on the context. When \code{effectsize = "boot"}, it defines the number of bootstrap @@ -549,9 +550,10 @@ covariates should not be removed from the SD used to standardize. Otherwise, \code{predict = "latent"}, sigma is the SD of the error distribution of the latent variable: \code{pi / sqrt(3)} (the standard logistic distribution) for the logit link, and \code{1} (the standard normal distribution) for the probit -link. Sigma is taken as known, with \code{edf = Inf}, so the effect sizes and -their confidence limits are the contrasts and their limits divided by -sigma. +link. Sigma is taken as known, with \code{edf = Inf}, so the effect sizes are +the contrasts divided by sigma. Their confidence limits are the unadjusted +confidence limits of the contrasts divided by sigma, also when \code{p_adjust} +is not \code{"none"}. \code{effectsize = "marginal"} uses the following formula to compute effect size: \code{d_adj <- difference * (1- R2)/ sigma}. This standardizes diff --git a/tests/testthat/test-estimate_contrasts_effectsize_order.R b/tests/testthat/test-estimate_contrasts_effectsize_order.R index c1343109b..e3d8e48f1 100644 --- a/tests/testthat/test-estimate_contrasts_effectsize_order.R +++ b/tests/testthat/test-estimate_contrasts_effectsize_order.R @@ -12,3 +12,16 @@ test_that("estimate_contrasts, effectsize emmeans, effect sizes follow their con # e.g. versicolor - setosa: 0.930 / sigma(g) = 0.930 / 0.5148 = 1.807 expect_equal(out$partial_d, out$Difference / sigma(g), tolerance = 1e-6) }) + + +test_that("estimate_contrasts, effectsize emmeans, effect sizes with a by variable", { + d <- mtcars + d$cyl <- factor(d$cyl) + d$am <- factor(d$am) + m <- lm(mpg ~ cyl * am, data = d) + out <- estimate_contrasts(m, "cyl", by = "am", effectsize = "emmeans", backend = "emmeans") + expect_identical(nrow(out), 6L) + # e.g. cyl 6 - cyl 4 at am 0: -3.775 / sigma(m) = -3.775 / 3.032 = -1.245 + expect_type(out$partial_d, "double") + expect_equal(out$partial_d, out$Difference / sigma(m), tolerance = 1e-6) +}) diff --git a/tests/testthat/test-ordinal.R b/tests/testthat/test-ordinal.R index 9863151a1..5012276ae 100644 --- a/tests/testthat/test-ordinal.R +++ b/tests/testthat/test-ordinal.R @@ -2075,6 +2075,22 @@ test_that("estimate_contrasts, backend emmeans, ordinal, latent-scale effect siz expect_identical(as.character(out_cp$Level1), as.character(out_p$Level1)) expect_identical(as.character(out_cp$Level2), as.character(out_p$Level2)) expect_equal(out_cp$partial_d, out_p$partial_d, tolerance = 1e-3) + + # a `by` variable adds a column to the emmeans effect sizes + mi <- MASS::polr(Sat ~ Infl * Type, weights = Freq, data = housing, Hess = TRUE) + out_by <- estimate_contrasts( + mi, + "Infl", + by = "Type", + predict = "latent", + effectsize = "emmeans", + backend = "emmeans" + ) + expect_identical(nrow(out_by), 12L) + expect_type(out_by$partial_d, "double") + expect_equal(out_by$partial_d, out_by$Difference / sd_logit, tolerance = 1e-6) + expect_equal(out_by$es_CI_low, out_by$CI_low / sd_logit, tolerance = 1e-6) + expect_equal(out_by$es_CI_high, out_by$CI_high / sd_logit, tolerance = 1e-6) }) @@ -2093,7 +2109,7 @@ test_that("estimate_contrasts, backend emmeans, ordinal, effect size errors", { expect_match(conditionMessage(err), "predict = \"latent\"", fixed = TRUE) } - # only the logit and probit links have a known latent SD here + # only the logit and probit links are supported here mcl <- MASS::polr( Sat ~ Infl + Type, weights = Freq, From e17716852e6ddb3a66e3f085162e205b3b11df44 Mon Sep 17 00:00:00 2001 From: Jeffrey Girard Date: Mon, 19 Oct 2026 19:28:09 -0500 Subject: [PATCH 5/8] Effect-size limits at the ci level, rescale-aware latent sigma, and clear errors for nominal clm models and joint tests --- NEWS.md | 22 +++--- R/estimate_contrasts.R | 23 +++--- R/estimate_contrasts_effectsize.R | 24 ++++++- man/estimate_contrasts.Rd | 19 +++-- ...test-estimate_contrasts_effectsize_order.R | 35 +++++++++ tests/testthat/test-ordinal.R | 71 +++++++++++++++++++ 6 files changed, 168 insertions(+), 26 deletions(-) diff --git a/NEWS.md b/NEWS.md index d88d165fb..f63e0e95b 100644 --- a/NEWS.md +++ b/NEWS.md @@ -7,9 +7,9 @@ models with `predict = "latent"` and a logit or probit link. The contrasts of the latent variable are divided by the SD of its error distribution: `pi / sqrt(3)` for the logit link and `1` for the probit link. Other values - of `predict`, other links and `clm` models with a `scale` formula give an - error. `effectsize = "marginal"` and `effectsize = "boot"` still give an - error for ordinal models. + of `predict`, other links and `clm` models with a `scale` or `nominal` + formula give an error. `effectsize = "marginal"` and `effectsize = "boot"` + still give an error for ordinal models. ## Bug fixes @@ -22,11 +22,17 @@ calls with `estimate = "average"` ignored the filter. * `estimate_contrasts()` with `effectsize = "emmeans"` now puts each effect - size on the row of its own contrast. Before, if the factor levels were not - in alphabetical order, the contrasts were sorted in another order than in - *emmeans*, and some effect sizes were on the wrong rows. With a `by` - variable, the effect-size columns held the values of the `by` variable and - of other *emmeans* columns instead of the effect sizes. + size on the row of its own contrast. Before, the contrasts were sometimes sorted + in another order than in *emmeans*, for example with factor levels that + were not in alphabetical order or with `comparison = "poly"`, and then some + effect sizes were on the wrong rows. With a `by` variable, the effect-size + columns held the values of the `by` variable and of other *emmeans* columns + instead of the effect sizes. + +* `estimate_contrasts()` with `effectsize = "emmeans"` now computes the + confidence limits of the effect sizes at the level of `ci`. Before, they + were always 95% limits. With `comparison = "joint"`, effect sizes now give + a clear error. * Filters on the response categories of frequentist ordinal models with `backend = "emmeans"`, where `Sat` is the response: diff --git a/R/estimate_contrasts.R b/R/estimate_contrasts.R index a28bdd933..0736aabc0 100644 --- a/R/estimate_contrasts.R +++ b/R/estimate_contrasts.R @@ -106,7 +106,8 @@ #' size will be computed. Effect sizes need a Gaussian model without random #' effects. Only `"emmeans"` also works for `polr` and `clm` models with #' `predict = "latent"` and a logit or probit link, but not for `clm` models -#' with a `scale` formula. +#' with a `scale` or `nominal` formula. Effect sizes are not possible for joint +#' tests. #' @param es_type Specifies the type of effect-size measure to estimate when #' using `effectsize = "boot"`. One of `"unstandardized"`, `"cohens.d"`, #' `"hedges.g"`, `"cohens.d.sigma"`, `"r"`, or `"akp.robust.d"`. See `effect.type` @@ -245,8 +246,9 @@ #' aren't part of the contrast. #' #' `effectsize = "emmeans"` uses [emmeans::eff_size] with -#' `sigma = stats::sigma(model)`, `edf = stats::df.residual(model)` and -#' `method = "identity"`. This standardizes using the MSE (sigma). Some believe +#' `method = "identity"`. For Gaussian models, it uses +#' `sigma = stats::sigma(model)` and `edf = stats::df.residual(model)`. This +#' standardizes using the MSE (sigma). Some believe #' this works when the contrasts are the only predictors in the model, but not #' when there are covariates. The response variance accounted for by the #' covariates should not be removed from the SD used to standardize. Otherwise, @@ -254,10 +256,13 @@ #' `predict = "latent"`, sigma is the SD of the error distribution of the #' latent variable: `pi / sqrt(3)` (the standard logistic distribution) for #' the logit link, and `1` (the standard normal distribution) for the probit -#' link. Sigma is taken as known, with `edf = Inf`, so the effect sizes are -#' the contrasts divided by sigma. Their confidence limits are the unadjusted -#' confidence limits of the contrasts divided by sigma, also when `p_adjust` -#' is not `"none"`. +#' link. With a `rescale` argument for *emmeans*, sigma is multiplied by the +#' absolute value of `rescale[2]`, so the effect sizes do not change. Sigma is +#' taken as known, with `edf = Inf`, so the effect sizes are the contrasts +#' divided by sigma, and their confidence limits are the unadjusted confidence +#' limits of the contrasts divided by sigma. For all models, the confidence +#' limits of the effect sizes are not adjusted for multiple comparisons, also +#' when `p_adjust` is not `"none"`. #' #' `effectsize = "marginal"` uses the following formula to compute effect #' size: `d_adj <- difference * (1- R2)/ sigma`. This standardizes @@ -571,7 +576,9 @@ estimate_contrasts.default <- function( effectsize = effectsize, bootstraps = bootstraps, bootES_type = es_type, - backend = backend + backend = backend, + ci = ci, + rescale = list(...)[["rescale"]] ) } attr(out, "emmeans_rows") <- NULL diff --git a/R/estimate_contrasts_effectsize.R b/R/estimate_contrasts_effectsize.R index 0478801a1..7601e44e9 100644 --- a/R/estimate_contrasts_effectsize.R +++ b/R/estimate_contrasts_effectsize.R @@ -5,11 +5,17 @@ effectsize, bootstraps, bootES_type, - backend + backend, + ci = 0.95, + rescale = NULL ) { # Add standardized effect size insight::validate_argument(effectsize, c("none", "emmeans", "marginal", "boot")) + if (effectsize != "none" && isTRUE(attributes(estimated)$joint_test)) { + insight::format_error("Effect sizes are not possible for joint tests.") + } + if (effectsize == "emmeans" && backend != "emmeans") { insight::format_error( "`effectsize = \"emmeans\"` only possible with `backend = \"emmeans\"`" @@ -40,7 +46,8 @@ # Check if the model's response variable follows a Gaussian (normal) # distribution. Effect size calculations implemented in this function are # designed for Gaussian models and may not be appropriate or meaningful for - # other types of distributions. + # other types of distributions. Ordinal models with `effectsize = "emmeans"` + # are the exception, because the latent variable has a known error SD. if (!latent_ordinal && insight::get_family(model)$family != "gaussian") { insight::format_error(paste0( "We strongly recommend not using the `effectsize` ", @@ -53,6 +60,11 @@ emmeans = { if (latent_ordinal) { sigma <- .ordinal_latent_sigma(model, attributes(estimated)$predict) + # `rescale` (passed to emmeans) scales the latent variable, and so its + # error SD + if (!is.null(rescale)) { + sigma <- sigma * abs(rescale[2]) + } edf <- Inf } else { sigma <- insight::get_sigma(model, verbose = FALSE) @@ -64,7 +76,7 @@ edf = edf, method = "identity" ) - eff <- as.data.frame(eff) + eff <- as.data.frame(stats::confint(eff, level = ci, adjust = "none")) # select columns by name, because `by` variables add columns es_col <- intersect(c("effect.size", "estimate"), colnames(eff))[1] ci_low <- intersect(c("lower.CL", "asymp.LCL"), colnames(eff))[1] @@ -148,6 +160,12 @@ "Effect sizes are not possible for ordinal models with a `scale` formula." ) } + # nominal effects make the thresholds differ between groups + if (inherits(model, "clm") && !is.null(model$formulas$nominal)) { + insight::format_error( + "Effect sizes are not possible for ordinal models with a `nominal` formula." + ) + } link <- if (inherits(model, "polr")) model$method else model$link switch( link, diff --git a/man/estimate_contrasts.Rd b/man/estimate_contrasts.Rd index 48fefccec..330b82ad1 100644 --- a/man/estimate_contrasts.Rd +++ b/man/estimate_contrasts.Rd @@ -322,7 +322,8 @@ to the output. You can reshape them to a long format by running size will be computed. Effect sizes need a Gaussian model without random effects. Only \code{"emmeans"} also works for \code{polr} and \code{clm} models with \code{predict = "latent"} and a logit or probit link, but not for \code{clm} models -with a \code{scale} formula.} +with a \code{scale} or \code{nominal} formula. Effect sizes are not possible for joint +tests.} \item{iterations}{This argument has two distinct effects, depending on the context. When \code{effectsize = "boot"}, it defines the number of bootstrap @@ -541,8 +542,9 @@ sigma adding back in the response variance associated with the variables that aren't part of the contrast. \code{effectsize = "emmeans"} uses \link[emmeans:eff_size]{emmeans::eff_size} with -\code{sigma = stats::sigma(model)}, \code{edf = stats::df.residual(model)} and -\code{method = "identity"}. This standardizes using the MSE (sigma). Some believe +\code{method = "identity"}. For Gaussian models, it uses +\code{sigma = stats::sigma(model)} and \code{edf = stats::df.residual(model)}. This +standardizes using the MSE (sigma). Some believe this works when the contrasts are the only predictors in the model, but not when there are covariates. The response variance accounted for by the covariates should not be removed from the SD used to standardize. Otherwise, @@ -550,10 +552,13 @@ covariates should not be removed from the SD used to standardize. Otherwise, \code{predict = "latent"}, sigma is the SD of the error distribution of the latent variable: \code{pi / sqrt(3)} (the standard logistic distribution) for the logit link, and \code{1} (the standard normal distribution) for the probit -link. Sigma is taken as known, with \code{edf = Inf}, so the effect sizes are -the contrasts divided by sigma. Their confidence limits are the unadjusted -confidence limits of the contrasts divided by sigma, also when \code{p_adjust} -is not \code{"none"}. +link. With a \code{rescale} argument for \emph{emmeans}, sigma is multiplied by the +absolute value of \code{rescale[2]}, so the effect sizes do not change. Sigma is +taken as known, with \code{edf = Inf}, so the effect sizes are the contrasts +divided by sigma, and their confidence limits are the unadjusted confidence +limits of the contrasts divided by sigma. For all models, the confidence +limits of the effect sizes are not adjusted for multiple comparisons, also +when \code{p_adjust} is not \code{"none"}. \code{effectsize = "marginal"} uses the following formula to compute effect size: \code{d_adj <- difference * (1- R2)/ sigma}. This standardizes diff --git a/tests/testthat/test-estimate_contrasts_effectsize_order.R b/tests/testthat/test-estimate_contrasts_effectsize_order.R index e3d8e48f1..c30084717 100644 --- a/tests/testthat/test-estimate_contrasts_effectsize_order.R +++ b/tests/testthat/test-estimate_contrasts_effectsize_order.R @@ -25,3 +25,38 @@ test_that("estimate_contrasts, effectsize emmeans, effect sizes with a by variab expect_type(out$partial_d, "double") expect_equal(out$partial_d, out$Difference / sigma(m), tolerance = 1e-6) }) + + +test_that("estimate_contrasts, effectsize emmeans, polynomial contrasts and ci level", { + g <- lm(Sepal.Length ~ Species, data = iris) + # the polynomial contrasts are sorted as "linear", "quadratic" + out <- estimate_contrasts( + g, + "Species", + comparison = "poly", + effectsize = "emmeans", + backend = "emmeans" + ) + expect_identical(nrow(out), 2L) + expect_equal(out$partial_d, out$Difference / sigma(g), tolerance = 1e-6) + + # the effect-size limits are at the level of `ci` + e <- emmeans::eff_size( + emmeans::emmeans(g, "Species"), + sigma = sigma(g), + edf = df.residual(g) + ) + e80 <- as.data.frame(confint(e, level = 0.8, adjust = "none")) + out80 <- estimate_contrasts(g, "Species", effectsize = "emmeans", backend = "emmeans", ci = 0.8) + # emmeans contrasts the levels in the other direction, so the rows are + # matched by the absolute effect size, and the half-widths are compared + idx <- match(round(abs(out80$partial_d), 6), round(abs(e80$effect.size), 6)) + expect_false(anyNA(idx)) + expect_equal( + (out80$es_CI_high - out80$es_CI_low) / 2, + ((e80$upper.CL - e80$lower.CL) / 2)[idx], + tolerance = 1e-6 + ) + out95 <- estimate_contrasts(g, "Species", effectsize = "emmeans", backend = "emmeans") + expect_true(all(out80$es_CI_high - out80$es_CI_low < out95$es_CI_high - out95$es_CI_low)) +}) diff --git a/tests/testthat/test-ordinal.R b/tests/testthat/test-ordinal.R index 5012276ae..706e23333 100644 --- a/tests/testthat/test-ordinal.R +++ b/tests/testthat/test-ordinal.R @@ -2153,3 +2153,74 @@ test_that("estimate_contrasts, backend emmeans, ordinal, effect size errors", { expect_match(conditionMessage(err), "non-Gaussian", fixed = TRUE) } }) + + +test_that("estimate_contrasts, backend emmeans, ordinal, effect sizes with ci, p_adjust and rescale", { + skip_if_not_installed("emmeans") + skip_if_not_installed("ordinal") + data(housing, package = "MASS") + sd_logit <- pi / sqrt(3) + m <- MASS::polr(Sat ~ Infl + Type, weights = Freq, data = housing, Hess = TRUE) + out <- estimate_contrasts(m, "Infl", predict = "latent", effectsize = "emmeans", backend = "emmeans") + + # the effect-size limits are at the level of `ci` + out80 <- estimate_contrasts( + m, + "Infl", + predict = "latent", + effectsize = "emmeans", + backend = "emmeans", + ci = 0.8 + ) + expect_equal(out80$es_CI_low, out80$CI_low / sd_logit, tolerance = 1e-6) + expect_equal(out80$es_CI_high, out80$CI_high / sd_logit, tolerance = 1e-6) + expect_true(all(out80$es_CI_low > out$es_CI_low)) + + # the effect-size limits are not adjusted for multiple comparisons + out_tukey <- estimate_contrasts( + m, + "Infl", + predict = "latent", + effectsize = "emmeans", + backend = "emmeans", + p_adjust = "tukey" + ) + expect_equal(out_tukey$partial_d, out$partial_d, tolerance = 1e-6) + expect_equal(out_tukey$es_CI_low, out$CI_low / sd_logit, tolerance = 1e-6) + expect_equal(out_tukey$es_CI_high, out$CI_high / sd_logit, tolerance = 1e-6) + expect_true(all(out_tukey$CI_low < out$CI_low)) + + # `rescale` scales the latent variable and its error SD alike + out_r <- estimate_contrasts( + m, + "Infl", + predict = "latent", + effectsize = "emmeans", + backend = "emmeans", + rescale = c(0, 2) + ) + expect_equal(out_r$Difference, 2 * out$Difference, tolerance = 1e-6) + expect_equal(out_r$partial_d, out$partial_d, tolerance = 1e-6) + expect_equal(out_r$es_CI_low, out$es_CI_low, tolerance = 1e-6) + + # nominal effects make the thresholds differ between groups + mn <- ordinal::clm(Sat ~ Infl, nominal = ~Type, weights = Freq, data = housing) + err <- expect_error( + estimate_contrasts(mn, "Type", predict = "latent", effectsize = "emmeans", backend = "emmeans") + ) + expect_match(gsub("\\s+", " ", conditionMessage(err)), "with a `nominal` formula", fixed = TRUE) + + # joint tests have no contrasts to standardize + err <- expect_error( + estimate_contrasts( + m, + "Infl", + by = "Type", + comparison = "joint", + predict = "latent", + effectsize = "emmeans", + backend = "emmeans" + ) + ) + expect_match(conditionMessage(err), "not possible for joint tests", fixed = TRUE) +}) From 0a1186f1b37ec9a615e3b35313cd94c9a95c3082 Mon Sep 17 00:00:00 2001 From: Jeffrey Girard Date: Mon, 19 Oct 2026 19:36:18 -0500 Subject: [PATCH 6/8] Correct effect-size docs for Gaussian sigma and negative rescale; test polynomial contrasts with 6 levels --- R/estimate_contrasts.R | 9 ++++++--- man/estimate_contrasts.Rd | 9 ++++++--- .../test-estimate_contrasts_effectsize_order.R | 15 +++++++++------ 3 files changed, 21 insertions(+), 12 deletions(-) diff --git a/R/estimate_contrasts.R b/R/estimate_contrasts.R index 0736aabc0..434a36aef 100644 --- a/R/estimate_contrasts.R +++ b/R/estimate_contrasts.R @@ -246,8 +246,9 @@ #' aren't part of the contrast. #' #' `effectsize = "emmeans"` uses [emmeans::eff_size] with -#' `method = "identity"`. For Gaussian models, it uses -#' `sigma = stats::sigma(model)` and `edf = stats::df.residual(model)`. This +#' `method = "identity"`. For Gaussian models, it uses the residual SD as +#' `sigma` and the residual degrees of freedom as `edf` (from +#' `insight::get_sigma()` and `insight::get_df()`). This #' standardizes using the MSE (sigma). Some believe #' this works when the contrasts are the only predictors in the model, but not #' when there are covariates. The response variance accounted for by the @@ -257,7 +258,9 @@ #' latent variable: `pi / sqrt(3)` (the standard logistic distribution) for #' the logit link, and `1` (the standard normal distribution) for the probit #' link. With a `rescale` argument for *emmeans*, sigma is multiplied by the -#' absolute value of `rescale[2]`, so the effect sizes do not change. Sigma is +#' absolute value of `rescale[2]`, so the size of the effect sizes does not +#' change. A negative `rescale[2]` reverses the sign of the contrasts and of +#' the effect sizes. Sigma is #' taken as known, with `edf = Inf`, so the effect sizes are the contrasts #' divided by sigma, and their confidence limits are the unadjusted confidence #' limits of the contrasts divided by sigma. For all models, the confidence diff --git a/man/estimate_contrasts.Rd b/man/estimate_contrasts.Rd index 330b82ad1..add336fe3 100644 --- a/man/estimate_contrasts.Rd +++ b/man/estimate_contrasts.Rd @@ -542,8 +542,9 @@ sigma adding back in the response variance associated with the variables that aren't part of the contrast. \code{effectsize = "emmeans"} uses \link[emmeans:eff_size]{emmeans::eff_size} with -\code{method = "identity"}. For Gaussian models, it uses -\code{sigma = stats::sigma(model)} and \code{edf = stats::df.residual(model)}. This +\code{method = "identity"}. For Gaussian models, it uses the residual SD as +\code{sigma} and the residual degrees of freedom as \code{edf} (from +\code{insight::get_sigma()} and \code{insight::get_df()}). This standardizes using the MSE (sigma). Some believe this works when the contrasts are the only predictors in the model, but not when there are covariates. The response variance accounted for by the @@ -553,7 +554,9 @@ covariates should not be removed from the SD used to standardize. Otherwise, latent variable: \code{pi / sqrt(3)} (the standard logistic distribution) for the logit link, and \code{1} (the standard normal distribution) for the probit link. With a \code{rescale} argument for \emph{emmeans}, sigma is multiplied by the -absolute value of \code{rescale[2]}, so the effect sizes do not change. Sigma is +absolute value of \code{rescale[2]}, so the size of the effect sizes does not +change. A negative \code{rescale[2]} reverses the sign of the contrasts and of +the effect sizes. Sigma is taken as known, with \code{edf = Inf}, so the effect sizes are the contrasts divided by sigma, and their confidence limits are the unadjusted confidence limits of the contrasts divided by sigma. For all models, the confidence diff --git a/tests/testthat/test-estimate_contrasts_effectsize_order.R b/tests/testthat/test-estimate_contrasts_effectsize_order.R index c30084717..b9c67ea79 100644 --- a/tests/testthat/test-estimate_contrasts_effectsize_order.R +++ b/tests/testthat/test-estimate_contrasts_effectsize_order.R @@ -28,17 +28,20 @@ test_that("estimate_contrasts, effectsize emmeans, effect sizes with a by variab test_that("estimate_contrasts, effectsize emmeans, polynomial contrasts and ci level", { - g <- lm(Sepal.Length ~ Species, data = iris) - # the polynomial contrasts are sorted as "linear", "quadratic" + # with 6 levels, the polynomial contrasts ("linear" to "degree 5") are + # sorted in another order than in emmeans + w <- lm(weight ~ feed, data = chickwts) out <- estimate_contrasts( - g, - "Species", + w, + "feed", comparison = "poly", effectsize = "emmeans", backend = "emmeans" ) - expect_identical(nrow(out), 2L) - expect_equal(out$partial_d, out$Difference / sigma(g), tolerance = 1e-6) + expect_identical(nrow(out), 5L) + expect_equal(out$partial_d, out$Difference / sigma(w), tolerance = 1e-6) + + g <- lm(Sepal.Length ~ Species, data = iris) # the effect-size limits are at the level of `ci` e <- emmeans::eff_size( From 0a50fc078b64cb15129ed6861f1bd0007a841ae7 Mon Sep 17 00:00:00 2001 From: Jeffrey Girard Date: Mon, 19 Oct 2026 19:49:05 -0500 Subject: [PATCH 7/8] Name omnibus tests in the effect-size error, accept a partly matched rescale, and clarify effect-size docs and NEWS --- NEWS.md | 8 ++++++-- R/estimate_contrasts.R | 6 +++--- R/estimate_contrasts_effectsize.R | 18 +++++++++++++++++- man/estimate_contrasts.Rd | 4 ++-- tests/testthat/test-ordinal.R | 25 ++++++++++++++++++++++++- 5 files changed, 52 insertions(+), 9 deletions(-) diff --git a/NEWS.md b/NEWS.md index f63e0e95b..99522ee36 100644 --- a/NEWS.md +++ b/NEWS.md @@ -31,8 +31,12 @@ * `estimate_contrasts()` with `effectsize = "emmeans"` now computes the confidence limits of the effect sizes at the level of `ci`. Before, they - were always 95% limits. With `comparison = "joint"`, effect sizes now give - a clear error. + were always 95% limits. + +* `estimate_contrasts()` with an `effectsize` and `comparison = "joint"`, or + `comparison = "omnibus"` with `backend = "marginaleffects"`, now gives a + clear error. Before, the emmeans effect sizes gave an internal error, and + `effectsize = "marginal"` returned `NA` effect sizes. * Filters on the response categories of frequentist ordinal models with `backend = "emmeans"`, where `Sat` is the response: diff --git a/R/estimate_contrasts.R b/R/estimate_contrasts.R index 434a36aef..5d63d1fad 100644 --- a/R/estimate_contrasts.R +++ b/R/estimate_contrasts.R @@ -107,7 +107,7 @@ #' effects. Only `"emmeans"` also works for `polr` and `clm` models with #' `predict = "latent"` and a logit or probit link, but not for `clm` models #' with a `scale` or `nominal` formula. Effect sizes are not possible for joint -#' tests. +#' or omnibus tests. #' @param es_type Specifies the type of effect-size measure to estimate when #' using `effectsize = "boot"`. One of `"unstandardized"`, `"cohens.d"`, #' `"hedges.g"`, `"cohens.d.sigma"`, `"r"`, or `"akp.robust.d"`. See `effect.type` @@ -260,7 +260,7 @@ #' link. With a `rescale` argument for *emmeans*, sigma is multiplied by the #' absolute value of `rescale[2]`, so the size of the effect sizes does not #' change. A negative `rescale[2]` reverses the sign of the contrasts and of -#' the effect sizes. Sigma is +#' the effect sizes. For these models, sigma is #' taken as known, with `edf = Inf`, so the effect sizes are the contrasts #' divided by sigma, and their confidence limits are the unadjusted confidence #' limits of the contrasts divided by sigma. For all models, the confidence @@ -581,7 +581,7 @@ estimate_contrasts.default <- function( bootES_type = es_type, backend = backend, ci = ci, - rescale = list(...)[["rescale"]] + rescale = .rescale_argument(...) ) } attr(out, "emmeans_rows") <- NULL diff --git a/R/estimate_contrasts_effectsize.R b/R/estimate_contrasts_effectsize.R index 7601e44e9..ccafaaccf 100644 --- a/R/estimate_contrasts_effectsize.R +++ b/R/estimate_contrasts_effectsize.R @@ -13,7 +13,7 @@ insight::validate_argument(effectsize, c("none", "emmeans", "marginal", "boot")) if (effectsize != "none" && isTRUE(attributes(estimated)$joint_test)) { - insight::format_error("Effect sizes are not possible for joint tests.") + insight::format_error("Effect sizes are not possible for joint or omnibus tests.") } if (effectsize == "emmeans" && backend != "emmeans") { @@ -182,6 +182,22 @@ } +# the `rescale` argument for emmeans, which emmeans also accepts partly +# matched, e.g. `resc = c(0, 2)` +.rescale_argument <- function(...) { + dots <- list(...) + dot_names <- names(dots) + if (is.null(dot_names)) { + return(NULL) + } + i <- which(nzchar(dot_names) & startsWith("rescale", dot_names)) + if (length(i) == 0) { + return(NULL) + } + dots[[i[1]]] +} + + .get_group_variable <- function(estimated, backend) { if (backend == "emmeans") { names(estimated@model.info$xlev) diff --git a/man/estimate_contrasts.Rd b/man/estimate_contrasts.Rd index add336fe3..6050cc9cf 100644 --- a/man/estimate_contrasts.Rd +++ b/man/estimate_contrasts.Rd @@ -323,7 +323,7 @@ size will be computed. Effect sizes need a Gaussian model without random effects. Only \code{"emmeans"} also works for \code{polr} and \code{clm} models with \code{predict = "latent"} and a logit or probit link, but not for \code{clm} models with a \code{scale} or \code{nominal} formula. Effect sizes are not possible for joint -tests.} +or omnibus tests.} \item{iterations}{This argument has two distinct effects, depending on the context. When \code{effectsize = "boot"}, it defines the number of bootstrap @@ -556,7 +556,7 @@ the logit link, and \code{1} (the standard normal distribution) for the probit link. With a \code{rescale} argument for \emph{emmeans}, sigma is multiplied by the absolute value of \code{rescale[2]}, so the size of the effect sizes does not change. A negative \code{rescale[2]} reverses the sign of the contrasts and of -the effect sizes. Sigma is +the effect sizes. For these models, sigma is taken as known, with \code{edf = Inf}, so the effect sizes are the contrasts divided by sigma, and their confidence limits are the unadjusted confidence limits of the contrasts divided by sigma. For all models, the confidence diff --git a/tests/testthat/test-ordinal.R b/tests/testthat/test-ordinal.R index 706e23333..21abe3572 100644 --- a/tests/testthat/test-ordinal.R +++ b/tests/testthat/test-ordinal.R @@ -2202,6 +2202,29 @@ test_that("estimate_contrasts, backend emmeans, ordinal, effect sizes with ci, p expect_equal(out_r$Difference, 2 * out$Difference, tolerance = 1e-6) expect_equal(out_r$partial_d, out$partial_d, tolerance = 1e-6) expect_equal(out_r$es_CI_low, out$es_CI_low, tolerance = 1e-6) + # a partly matched name reaches emmeans too + out_resc <- estimate_contrasts( + m, + "Infl", + predict = "latent", + effectsize = "emmeans", + backend = "emmeans", + resc = c(0, 2) + ) + expect_equal(out_resc$Difference, 2 * out$Difference, tolerance = 1e-6) + expect_equal(out_resc$partial_d, out$partial_d, tolerance = 1e-6) + # a negative scale reverses the sign of the contrasts and effect sizes + out_neg <- estimate_contrasts( + m, + "Infl", + predict = "latent", + effectsize = "emmeans", + backend = "emmeans", + rescale = c(1, -3) + ) + expect_equal(out_neg$Difference, -3 * out$Difference, tolerance = 1e-6) + expect_equal(out_neg$partial_d, -out$partial_d, tolerance = 1e-6) + expect_equal(out_neg$es_CI_low, -out$es_CI_high, tolerance = 1e-6) # nominal effects make the thresholds differ between groups mn <- ordinal::clm(Sat ~ Infl, nominal = ~Type, weights = Freq, data = housing) @@ -2222,5 +2245,5 @@ test_that("estimate_contrasts, backend emmeans, ordinal, effect sizes with ci, p backend = "emmeans" ) ) - expect_match(conditionMessage(err), "not possible for joint tests", fixed = TRUE) + expect_match(gsub("\\s+", " ", conditionMessage(err)), "not possible for joint or omnibus tests", fixed = TRUE) }) From 371d45ab73ef69d5a382ec24b2c4e7b79ba63dec Mon Sep 17 00:00:00 2001 From: Jeffrey Girard Date: Wed, 21 Oct 2026 13:38:47 -0500 Subject: [PATCH 8/8] Shorten effect-size docs, NEWS and comments --- NEWS.md | 37 +++++++------------ R/estimate_contrasts.R | 35 +++++++----------- R/estimate_contrasts_effectsize.R | 19 ++++------ man/estimate_contrasts.Rd | 35 +++++++----------- ...test-estimate_contrasts_effectsize_order.R | 9 ++--- tests/testthat/test-ordinal.R | 7 +--- 6 files changed, 51 insertions(+), 91 deletions(-) diff --git a/NEWS.md b/NEWS.md index 99522ee36..122f29d2e 100644 --- a/NEWS.md +++ b/NEWS.md @@ -2,14 +2,10 @@ ## Changes -* `estimate_contrasts()` with `backend = "emmeans"` and - `effectsize = "emmeans"` now returns effect sizes for `polr` and `clm` - models with `predict = "latent"` and a logit or probit link. The contrasts - of the latent variable are divided by the SD of its error distribution: - `pi / sqrt(3)` for the logit link and `1` for the probit link. Other values - of `predict`, other links and `clm` models with a `scale` or `nominal` - formula give an error. `effectsize = "marginal"` and `effectsize = "boot"` - still give an error for ordinal models. +* `effectsize = "emmeans"` in `estimate_contrasts()` now also works for `polr` + and `clm` models with `predict = "latent"` and a logit or probit link. The + contrasts of the latent variable are divided by the SD of the latent error, + i.e. `pi / sqrt(3)` (logit) or `1` (probit). ## Bug fixes @@ -21,22 +17,15 @@ `estimate = "typical"` stopped with the error "No contrasts to show", and calls with `estimate = "average"` ignored the filter. -* `estimate_contrasts()` with `effectsize = "emmeans"` now puts each effect - size on the row of its own contrast. Before, the contrasts were sometimes sorted - in another order than in *emmeans*, for example with factor levels that - were not in alphabetical order or with `comparison = "poly"`, and then some - effect sizes were on the wrong rows. With a `by` variable, the effect-size - columns held the values of the `by` variable and of other *emmeans* columns - instead of the effect sizes. - -* `estimate_contrasts()` with `effectsize = "emmeans"` now computes the - confidence limits of the effect sizes at the level of `ci`. Before, they - were always 95% limits. - -* `estimate_contrasts()` with an `effectsize` and `comparison = "joint"`, or - `comparison = "omnibus"` with `backend = "marginaleffects"`, now gives a - clear error. Before, the emmeans effect sizes gave an internal error, and - `effectsize = "marginal"` returned `NA` effect sizes. +* Fixed issue where `effectsize = "emmeans"` in `estimate_contrasts()` returned + effect sizes in the wrong rows (e.g., when factor levels were not sorted + alphabetically) or in the wrong columns (when `by` was specified). + +* The confidence intervals of effect sizes from `effectsize = "emmeans"` now + use the level of `ci`. Before, they were always 95% intervals. + +* `estimate_contrasts()` now gives an informative error when effect sizes are + requested for joint or omnibus tests. * Filters on the response categories of frequentist ordinal models with `backend = "emmeans"`, where `Sat` is the response: diff --git a/R/estimate_contrasts.R b/R/estimate_contrasts.R index 5d63d1fad..f49b783d8 100644 --- a/R/estimate_contrasts.R +++ b/R/estimate_contrasts.R @@ -103,11 +103,9 @@ #' specified post-process tests. See 'Examples'. #' @param effectsize Desired measure of standardized effect size, one of #' `"emmeans"`, `"marginal"`, or `"boot"`. Default is `NULL`, i.e. no effect -#' size will be computed. Effect sizes need a Gaussian model without random -#' effects. Only `"emmeans"` also works for `polr` and `clm` models with -#' `predict = "latent"` and a logit or probit link, but not for `clm` models -#' with a `scale` or `nominal` formula. Effect sizes are not possible for joint -#' or omnibus tests. +#' size will be computed. Effect sizes are only available for Gaussian models +#' without random effects, and (only for `"emmeans"`) for ordinal models on the +#' latent scale (see 'Details'). #' @param es_type Specifies the type of effect-size measure to estimate when #' using `effectsize = "boot"`. One of `"unstandardized"`, `"cohens.d"`, #' `"hedges.g"`, `"cohens.d.sigma"`, `"r"`, or `"akp.robust.d"`. See `effect.type` @@ -246,26 +244,19 @@ #' aren't part of the contrast. #' #' `effectsize = "emmeans"` uses [emmeans::eff_size] with -#' `method = "identity"`. For Gaussian models, it uses the residual SD as -#' `sigma` and the residual degrees of freedom as `edf` (from -#' `insight::get_sigma()` and `insight::get_df()`). This -#' standardizes using the MSE (sigma). Some believe +#' `sigma = stats::sigma(model)`, `edf = stats::df.residual(model)` and +#' `method = "identity"`. This standardizes using the MSE (sigma). Some believe #' this works when the contrasts are the only predictors in the model, but not #' when there are covariates. The response variance accounted for by the #' covariates should not be removed from the SD used to standardize. Otherwise, -#' _d_ will be overestimated. For `polr` and `clm` models with -#' `predict = "latent"`, sigma is the SD of the error distribution of the -#' latent variable: `pi / sqrt(3)` (the standard logistic distribution) for -#' the logit link, and `1` (the standard normal distribution) for the probit -#' link. With a `rescale` argument for *emmeans*, sigma is multiplied by the -#' absolute value of `rescale[2]`, so the size of the effect sizes does not -#' change. A negative `rescale[2]` reverses the sign of the contrasts and of -#' the effect sizes. For these models, sigma is -#' taken as known, with `edf = Inf`, so the effect sizes are the contrasts -#' divided by sigma, and their confidence limits are the unadjusted confidence -#' limits of the contrasts divided by sigma. For all models, the confidence -#' limits of the effect sizes are not adjusted for multiple comparisons, also -#' when `p_adjust` is not `"none"`. +#' _d_ will be overestimated. +#' +#' For `polr` and `clm` models with `predict = "latent"` and a logit or probit +#' link, `effectsize = "emmeans"` uses the SD of the latent error as sigma, i.e. +#' `pi / sqrt(3)` (logit) or `1` (probit), with `edf = Inf`. A `rescale` +#' argument passed to *emmeans* also rescales sigma. `clm` models with a `scale` +#' or `nominal` formula are not supported. The confidence intervals of the +#' effect sizes are not adjusted for multiple comparisons. #' #' `effectsize = "marginal"` uses the following formula to compute effect #' size: `d_adj <- difference * (1- R2)/ sigma`. This standardizes diff --git a/R/estimate_contrasts_effectsize.R b/R/estimate_contrasts_effectsize.R index ccafaaccf..2272626aa 100644 --- a/R/estimate_contrasts_effectsize.R +++ b/R/estimate_contrasts_effectsize.R @@ -39,15 +39,14 @@ )) } - # For ordinal models, `effectsize = "emmeans"` standardizes the contrasts of - # the latent variable by the SD of its error distribution + # ordinal models are standardized by the SD of the latent error latent_ordinal <- effectsize == "emmeans" && inherits(model, c("polr", "clm")) # Check if the model's response variable follows a Gaussian (normal) # distribution. Effect size calculations implemented in this function are # designed for Gaussian models and may not be appropriate or meaningful for - # other types of distributions. Ordinal models with `effectsize = "emmeans"` - # are the exception, because the latent variable has a known error SD. + # other types of distributions. Ordinal models on the latent scale are an + # exception for `effectsize = "emmeans"`. if (!latent_ordinal && insight::get_family(model)$family != "gaussian") { insight::format_error(paste0( "We strongly recommend not using the `effectsize` ", @@ -60,8 +59,7 @@ emmeans = { if (latent_ordinal) { sigma <- .ordinal_latent_sigma(model, attributes(estimated)$predict) - # `rescale` (passed to emmeans) scales the latent variable, and so its - # error SD + # `rescale` also rescales the latent error SD if (!is.null(rescale)) { sigma <- sigma * abs(rescale[2]) } @@ -83,7 +81,7 @@ ci_high <- intersect(c("upper.CL", "asymp.UCL"), colnames(eff))[1] eff <- eff[c(es_col, ci_low, ci_high)] names(eff) <- c("partial_d", "es_CI_low", "es_CI_high") - # the contrasts can be sorted in another order than the emmeans grid + # match effect sizes to the contrasts, which may be sorted differently emmeans_rows <- attributes(contrasts_results)$emmeans_rows if (!is.null(emmeans_rows)) { eff <- eff[emmeans_rows, , drop = FALSE] @@ -145,9 +143,7 @@ } -# SD of the error distribution of the latent variable of ordinal models, which -# is the standard logistic (SD pi / sqrt(3)) for the logit link, and the -# standard normal (SD 1) for the probit link +# SD of the latent error of ordinal models: pi / sqrt(3) for logit, 1 for probit .ordinal_latent_sigma <- function(model, predict) { if (!identical(predict, "latent")) { insight::format_error(paste0( @@ -182,8 +178,7 @@ } -# the `rescale` argument for emmeans, which emmeans also accepts partly -# matched, e.g. `resc = c(0, 2)` +# `rescale` argument for emmeans from `...`, also when partly matched .rescale_argument <- function(...) { dots <- list(...) dot_names <- names(dots) diff --git a/man/estimate_contrasts.Rd b/man/estimate_contrasts.Rd index 6050cc9cf..b507015a2 100644 --- a/man/estimate_contrasts.Rd +++ b/man/estimate_contrasts.Rd @@ -319,11 +319,9 @@ to the output. You can reshape them to a long format by running \item{effectsize}{Desired measure of standardized effect size, one of \code{"emmeans"}, \code{"marginal"}, or \code{"boot"}. Default is \code{NULL}, i.e. no effect -size will be computed. Effect sizes need a Gaussian model without random -effects. Only \code{"emmeans"} also works for \code{polr} and \code{clm} models with -\code{predict = "latent"} and a logit or probit link, but not for \code{clm} models -with a \code{scale} or \code{nominal} formula. Effect sizes are not possible for joint -or omnibus tests.} +size will be computed. Effect sizes are only available for Gaussian models +without random effects, and (only for \code{"emmeans"}) for ordinal models on the +latent scale (see 'Details').} \item{iterations}{This argument has two distinct effects, depending on the context. When \code{effectsize = "boot"}, it defines the number of bootstrap @@ -542,26 +540,19 @@ sigma adding back in the response variance associated with the variables that aren't part of the contrast. \code{effectsize = "emmeans"} uses \link[emmeans:eff_size]{emmeans::eff_size} with -\code{method = "identity"}. For Gaussian models, it uses the residual SD as -\code{sigma} and the residual degrees of freedom as \code{edf} (from -\code{insight::get_sigma()} and \code{insight::get_df()}). This -standardizes using the MSE (sigma). Some believe +\code{sigma = stats::sigma(model)}, \code{edf = stats::df.residual(model)} and +\code{method = "identity"}. This standardizes using the MSE (sigma). Some believe this works when the contrasts are the only predictors in the model, but not when there are covariates. The response variance accounted for by the covariates should not be removed from the SD used to standardize. Otherwise, -\emph{d} will be overestimated. For \code{polr} and \code{clm} models with -\code{predict = "latent"}, sigma is the SD of the error distribution of the -latent variable: \code{pi / sqrt(3)} (the standard logistic distribution) for -the logit link, and \code{1} (the standard normal distribution) for the probit -link. With a \code{rescale} argument for \emph{emmeans}, sigma is multiplied by the -absolute value of \code{rescale[2]}, so the size of the effect sizes does not -change. A negative \code{rescale[2]} reverses the sign of the contrasts and of -the effect sizes. For these models, sigma is -taken as known, with \code{edf = Inf}, so the effect sizes are the contrasts -divided by sigma, and their confidence limits are the unadjusted confidence -limits of the contrasts divided by sigma. For all models, the confidence -limits of the effect sizes are not adjusted for multiple comparisons, also -when \code{p_adjust} is not \code{"none"}. +\emph{d} will be overestimated. + +For \code{polr} and \code{clm} models with \code{predict = "latent"} and a logit or probit +link, \code{effectsize = "emmeans"} uses the SD of the latent error as sigma, i.e. +\code{pi / sqrt(3)} (logit) or \code{1} (probit), with \code{edf = Inf}. A \code{rescale} +argument passed to \emph{emmeans} also rescales sigma. \code{clm} models with a \code{scale} +or \code{nominal} formula are not supported. The confidence intervals of the +effect sizes are not adjusted for multiple comparisons. \code{effectsize = "marginal"} uses the following formula to compute effect size: \code{d_adj <- difference * (1- R2)/ sigma}. This standardizes diff --git a/tests/testthat/test-estimate_contrasts_effectsize_order.R b/tests/testthat/test-estimate_contrasts_effectsize_order.R index b9c67ea79..3b9a27906 100644 --- a/tests/testthat/test-estimate_contrasts_effectsize_order.R +++ b/tests/testthat/test-estimate_contrasts_effectsize_order.R @@ -2,8 +2,7 @@ skip_on_cran() skip_if_not_installed("emmeans") test_that("estimate_contrasts, effectsize emmeans, effect sizes follow their contrasts", { - # the levels are not in alphabetical order, so the contrast rows are sorted - # in another order than the rows of emmeans + # levels not in alphabetical order d <- iris d$Species <- factor(d$Species, levels = c("virginica", "setosa", "versicolor")) g <- lm(Sepal.Length ~ Species, data = d) @@ -28,8 +27,7 @@ test_that("estimate_contrasts, effectsize emmeans, effect sizes with a by variab test_that("estimate_contrasts, effectsize emmeans, polynomial contrasts and ci level", { - # with 6 levels, the polynomial contrasts ("linear" to "degree 5") are - # sorted in another order than in emmeans + # with 6 levels, polynomial contrasts are sorted differently than in emmeans w <- lm(weight ~ feed, data = chickwts) out <- estimate_contrasts( w, @@ -51,8 +49,7 @@ test_that("estimate_contrasts, effectsize emmeans, polynomial contrasts and ci l ) e80 <- as.data.frame(confint(e, level = 0.8, adjust = "none")) out80 <- estimate_contrasts(g, "Species", effectsize = "emmeans", backend = "emmeans", ci = 0.8) - # emmeans contrasts the levels in the other direction, so the rows are - # matched by the absolute effect size, and the half-widths are compared + # emmeans contrasts in the other direction, so compare absolute values idx <- match(round(abs(out80$partial_d), 6), round(abs(e80$effect.size), 6)) expect_false(anyNA(idx)) expect_equal( diff --git a/tests/testthat/test-ordinal.R b/tests/testthat/test-ordinal.R index 21abe3572..1d59dd7c1 100644 --- a/tests/testthat/test-ordinal.R +++ b/tests/testthat/test-ordinal.R @@ -1993,9 +1993,7 @@ test_that("estimate_contrasts, backend emmeans, ordinal, latent-scale effect siz skip_if_not_installed("ordinal") data(housing, package = "MASS") - # the latent errors of the logit link follow the standard logistic - # distribution, with variance pi^2 / 3 (?stats::Logistic), and those of the - # probit link the standard normal distribution, with SD 1 (?stats::Normal) + # SD of the latent error: pi / sqrt(3) for logit, 1 for probit sd_logit <- pi / sqrt(3) m <- MASS::polr(Sat ~ Infl + Type, weights = Freq, data = housing, Hess = TRUE) @@ -2007,8 +2005,7 @@ test_that("estimate_contrasts, backend emmeans, ordinal, latent-scale effect siz backend = "emmeans" ) expect_identical(nrow(out), 3L) - # sigma is taken as known (`edf = Inf`), so the effect sizes and their - # confidence limits are the differences and their limits divided by sigma + # with `edf = Inf`, effect sizes and limits are the contrasts divided by sigma expect_equal(out$partial_d, out$Difference / sd_logit, tolerance = 1e-6) expect_equal(out$es_CI_low, out$CI_low / sd_logit, tolerance = 1e-6) expect_equal(out$es_CI_high, out$CI_high / sd_logit, tolerance = 1e-6)