From 884a4d755096aa0bcfc9b55619923136a35e2c9a Mon Sep 17 00:00:00 2001 From: Jeffrey Girard Date: Wed, 21 Oct 2026 22:20:47 -0500 Subject: [PATCH] Remove redundant tests Keep one regression test per fixed behavior from the recent ordinal, p_adjust, focal-name and footer PRs. Drop parity loops over sibling model classes and inputs, extra model fits, custom matching helpers and narrative comments. No R code changes. --- tests/testthat/test-brms.R | 36 +- .../test-estimate_contrasts_slopes_filter.R | 195 +- tests/testthat/test-estimate_slopes.R | 118 +- tests/testthat/test-focal-name-clashes.R | 615 +----- tests/testthat/test-ordinal.R | 1842 ++++------------- tests/testthat/test-p_adjust.R | 259 +-- tests/testthat/test-table_footer.R | 167 +- tests/testthat/test-visualisation_recipe.R | 230 +- 8 files changed, 756 insertions(+), 2706 deletions(-) diff --git a/tests/testthat/test-brms.R b/tests/testthat/test-brms.R index 246fb0085..93fb1fec1 100644 --- a/tests/testthat/test-brms.R +++ b/tests/testthat/test-brms.R @@ -87,17 +87,19 @@ test_that("estimate_means - brms, multivariate", { }) -test_that("estimate_contrasts - brms, multivariate, many observations", { +test_that("estimate_contrasts - brms, multivariate, many observations and response names with dots", { # the response values must not be used as response levels, else the - # regular expressions built from them grow too long and fail + # regular expressions built from them grow too long and fail. brms drops the + # dot from response names ("y.1" -> "y1"), and the contrast labels use the + # brms names set.seed(123) d <- data.frame( - y1 = rnorm(600), - y2 = rnorm(600), + y.1 = rnorm(600), + y.2 = rnorm(600), g = factor(sample(c("a", "b", "c"), 600, replace = TRUE)) ) m <- suppressWarnings(brms::brm( - brms::mvbind(y1, y2) ~ g, + brms::mvbind(y.1, y.2) ~ g, data = d, refresh = 0, iter = 500, @@ -106,35 +108,13 @@ test_that("estimate_contrasts - brms, multivariate, many observations", { )) out <- estimate_contrasts(m, "g", resp = "y1") - expect_identical(nrow(out), 3L) expect_identical(as.character(out$Level1), c("b", "c", "c")) expect_identical(as.character(out$Level2), c("a", "a", "b")) out <- estimate_contrasts(m, "g") expect_identical(nrow(out), 15L) - expect_true(all(c("Response1", "Response2") %in% colnames(out))) expect_identical(levels(out$Response1), c("y1", "y2")) -}) - - -test_that("estimate_contrasts - brms, multivariate, response names with dots", { - # brms drops the dot from response names ("Sepal.Length" -> "SepalLength"), - # and the contrast labels use the brms names - m <- suppressWarnings(brms::brm( - brms::mvbind(Sepal.Length, Petal.Length) ~ Species, - data = iris, - refresh = 0, - iter = 500, - chains = 1, - seed = 123 - )) - - out <- estimate_contrasts(m, "Species") - expect_identical(nrow(out), 15L) - expect_true(all(c("Response1", "Response2") %in% colnames(out))) - expect_identical(levels(out$Response1), c("SepalLength", "PetalLength")) - expect_true(all(as.character(out$Level1) %in% levels(iris$Species))) - expect_true(all(as.character(out$Level2) %in% levels(iris$Species))) + expect_true(all(as.character(out$Level1) %in% levels(d$g))) }) diff --git a/tests/testthat/test-estimate_contrasts_slopes_filter.R b/tests/testthat/test-estimate_contrasts_slopes_filter.R index d61c04998..b219a096e 100644 --- a/tests/testthat/test-estimate_contrasts_slopes_filter.R +++ b/tests/testthat/test-estimate_contrasts_slopes_filter.R @@ -6,154 +6,60 @@ mt$g <- factor(mt$cyl) mt$a <- factor(mt$am) m2 <- lm(mpg ~ wt * g, data = mt) m3 <- lm(mpg ~ wt * g * a, data = mt) -m6 <- lm(mpg ~ wt * g * a + hp, data = mt) -# filters for `by`, the matching `at` for emmeans, and the levels that each -# filter allows -filters <- list( - F1 = list( - by = c("g='4'", "a"), - at = list(g = "4"), - allowed = list(g = "4", a = c("0", "1")) - ), - F2 = list( - by = c("g=c('4','6')", "a='1'"), - at = list(g = c("4", "6"), a = "1"), - allowed = list(g = c("4", "6"), a = "1") - ), - F3 = list( - by = c("g", "a='1'"), - at = list(a = "1"), - allowed = list(g = c("4", "6", "8"), a = "1") - ), - F4 = list( - by = "g=c('8','4')", - at = list(g = c("8", "4")), - allowed = list(g = c("8", "4")) - ) -) - -# identify the slope in a label like "4, 1", "6" or "g4 a1"; variables that -# the label leaves out take the one value their filter allows -.slope_key <- function(label, allowed) { - vars <- names(allowed) - tokens <- strsplit(trimws(label), ",\\s*|\\s+")[[1]] - prefixed <- vapply( - tokens, - function(tk) { - any(vapply(vars, function(v) startsWith(tk, v) && substring(tk, nchar(v) + 1) %in% allowed[[v]], logical(1))) - }, - logical(1) - ) - values <- list() - if (all(prefixed)) { - for (tk in tokens) { - v <- vars[vapply(vars, function(v) startsWith(tk, v), logical(1))][1] - values[[v]] <- substring(tk, nchar(v) + 1) - } - } else if (length(tokens) == length(vars)) { - values <- stats::setNames(as.list(tokens), vars) - } else { - values <- stats::setNames(as.list(tokens), vars[lengths(allowed) > 1]) - } - for (v in vars) { - if (is.null(values[[v]])) { - values[[v]] <- allowed[[v]] - } - } - paste(vars, unlist(values[vars]), sep = "=", collapse = ";") -} - -# check that `out` and the emmeans contrasts in `oracle` pair one to one, -# with the same two slopes, the same difference and the same p-value -expect_matches_oracle <- function(out, oracle, allowed) { - oracle <- as.data.frame(oracle) - expect_identical(nrow(out), nrow(oracle)) - oracle_pairs <- strsplit(as.character(oracle$contrast), " - ", fixed = TRUE) - oracle_1 <- vapply(oracle_pairs, function(x) .slope_key(x[1], allowed), character(1)) - oracle_2 <- vapply(oracle_pairs, function(x) .slope_key(x[2], allowed), character(1)) - used <- logical(nrow(oracle)) - for (i in seq_len(nrow(out))) { - key_1 <- .slope_key(as.character(out$Level1[i]), allowed) - key_2 <- .slope_key(as.character(out$Level2[i]), allowed) - same <- which(oracle_1 == key_1 & oracle_2 == key_2) - flipped <- which(oracle_1 == key_2 & oracle_2 == key_1) - expect_identical(length(same) + length(flipped), 1L) - if (length(same) + length(flipped) != 1L) { - next - } - j <- c(same, flipped) - sign <- if (length(same)) 1 else -1 - expect_false(used[j]) - used[j] <- TRUE - expect_equal(out$Difference[i], sign * oracle$estimate[j], tolerance = 1e-5, ignore_attr = TRUE) - expect_equal(out$p[i], oracle$p.value[j], tolerance = 1e-5, ignore_attr = TRUE) - } - expect_true(all(used)) -} - -.oracle <- function(model, filter) { - specs <- if (length(filter$allowed) == 1) ~g else ~ g * a - pairs(emmeans::emtrends(model, specs, "wt", at = filter$at), adjust = "none") +# modelbased orders and orients the pairs differently than emmeans +expect_slope_contrasts <- function(out, at) { + ref <- as.data.frame(pairs( + emmeans::emtrends(m3, ~ g * a, "wt", at = at), + adjust = "none" + )) + expect_equal(sort(abs(out$Difference)), sort(abs(ref$estimate)), tolerance = 1e-5) + expect_equal(sort(out$p), sort(ref$p.value), tolerance = 1e-5) } -probes <- list( - list(model = m3, filter = filters$F1), - list(model = m3, filter = filters$F2), - list(model = m3, filter = filters$F3), - list(model = m6, filter = filters$F1), - list(model = m6, filter = filters$F2), - list(model = m6, filter = filters$F3), - list(model = m2, filter = filters$F4) -) - -test_that("estimate_contrasts, slopes, by-filters, marginaleffects", { - for (probe in probes) { - oracle <- .oracle(probe$model, probe$filter) - for (estimate in c("typical", "average")) { - out <- estimate_contrasts(probe$model, "wt", by = probe$filter$by, estimate = estimate) - expect_matches_oracle(out, oracle, probe$filter$allowed) - } - } -}) +test_that("estimate_contrasts, slopes, by-filters", { + by <- c("g='4'", "a") + expect_slope_contrasts( + estimate_contrasts(m3, "wt", by = by, estimate = "average"), + list(g = "4") + ) + expect_slope_contrasts( + estimate_contrasts(m3, "wt", by = by, backend = "emmeans"), + list(g = "4") + ) + by <- c("g=c('4','6')", "a='1'") + at <- list(g = c("4", "6"), a = "1") + expect_slope_contrasts(estimate_contrasts(m3, "wt", by = by, estimate = "average"), at) + expect_slope_contrasts(estimate_contrasts(m3, "wt", by = by, backend = "emmeans"), at) -test_that("estimate_contrasts, slopes, by-filters, emmeans", { - for (probe in probes) { - oracle <- .oracle(probe$model, probe$filter) - out <- estimate_contrasts(probe$model, "wt", by = probe$filter$by, backend = "emmeans") - expect_matches_oracle(out, oracle, probe$filter$allowed) - } + by <- c("g=c('4','5')", "a") + expect_slope_contrasts( + estimate_contrasts(m3, "wt", by = by, estimate = "average"), + list(g = "4") + ) + expect_slope_contrasts( + estimate_contrasts(m3, "wt", by = by, backend = "emmeans"), + list(g = "4") + ) }) test_that("estimate_contrasts, slopes, by-filters leave fewer than two slopes", { - calls <- list( - list(model = m2, by = "g='4'"), - list(model = m3, by = c("g='4'", "a='1'")) + expect_error( + estimate_contrasts(m2, "wt", by = "g='4'", estimate = "average"), + regexp = "at least two" + ) + expect_error( + estimate_contrasts(m2, "wt", by = "g='4'", backend = "emmeans"), + regexp = "at least two" ) - for (cl in calls) { - expect_error( - estimate_contrasts(cl$model, "wt", by = cl$by, estimate = "typical"), - regexp = "at least two" - ) - expect_error( - estimate_contrasts(cl$model, "wt", by = cl$by, estimate = "average"), - regexp = "at least two" - ) - expect_error( - estimate_contrasts(cl$model, "wt", by = cl$by, backend = "emmeans"), - regexp = "at least two" - ) - } # spaces around "=", and `estimate = "specific"` expect_error( estimate_contrasts(m3, "wt", by = c("g = '4'", "a = '1'"), estimate = "specific"), regexp = "at least two" ) - out <- estimate_contrasts(m3, "wt", by = c("g = '4'", "a"), estimate = "average") - expect_matches_oracle(out, .oracle(m3, filters$F1), filters$F1$allowed) # with `estimate = "average"`, only slopes with data count: for 5 gears, # all cars have am = 1 mt$gear <- factor(mt$gear) @@ -162,39 +68,14 @@ test_that("estimate_contrasts, slopes, by-filters leave fewer than two slopes", estimate_contrasts(m7, "wt", by = c("gear='5'", "a"), estimate = "average"), regexp = "at least two" ) - expect_identical(nrow(estimate_contrasts(m7, "wt", by = c("gear='5'", "a"))), 1L) -}) - - -test_that("estimate_contrasts, slopes, by-filter values that are not levels", { - oracle <- .oracle(m3, filters$F1) - out <- estimate_contrasts(m3, "wt", by = c("g=c('4','5')", "a"), estimate = "average") - expect_matches_oracle(out, oracle, filters$F1$allowed) - out <- estimate_contrasts(m3, "wt", by = c("g=c('4','5')", "a"), backend = "emmeans") - expect_matches_oracle(out, oracle, filters$F1$allowed) }) test_that("estimate_contrasts, slopes, emmeans, numeric by or factors from the formula", { - # emmeans does not filter these, but still returns all contrasts m_factor <- lm(mpg ~ wt * factor(cyl), data = mtcars) out <- estimate_contrasts(m_factor, "wt", by = "cyl=c(4,6)", backend = "emmeans") expect_identical(nrow(out), 3L) - expect_true(all(is.finite(out$Difference))) m_numeric <- lm(mpg ~ wt * hp * g, data = mt) out <- estimate_contrasts(m_numeric, "wt", by = c("g='4'", "hp"), backend = "emmeans") expect_identical(nrow(out), 3L) - expect_true(all(is.finite(out$Difference))) -}) - - -test_that("estimate_contrasts, slopes, no by-filter", { - oracle <- pairs(emmeans::emtrends(m3, ~ g * a, "wt"), adjust = "none") - allowed <- list(g = c("4", "6", "8"), a = c("0", "1")) - for (estimate in c("typical", "average")) { - out <- estimate_contrasts(m3, "wt", by = c("g", "a"), estimate = estimate) - expect_matches_oracle(out, oracle, allowed) - } - out <- estimate_contrasts(m3, "wt", by = c("g", "a"), backend = "emmeans") - expect_matches_oracle(out, oracle, allowed) }) diff --git a/tests/testthat/test-estimate_slopes.R b/tests/testthat/test-estimate_slopes.R index cf3539237..98feacb57 100644 --- a/tests/testthat/test-estimate_slopes.R +++ b/tests/testthat/test-estimate_slopes.R @@ -260,40 +260,18 @@ test_that("estimate_slopes, custom comparison", { test_that("estimate_slopes, p_adjust methods adjust the p-values", { d <- datasets::mtcars - d$safe_var <- factor(d$cyl) - calls <- list( - numeric_by = function(m) { - estimate_slopes( - lm(mpg ~ wt * hp, data = d), - "wt", - by = "hp", - length = 4, - p_adjust = m - ) - }, - factor_by = function(m) { - estimate_slopes(lm(mpg ~ wt * safe_var, data = d), "wt", by = "safe_var", p_adjust = m) - } - ) - methods <- c("holm", "hochberg", "hommel", "bonferroni", "BH", "BY", "fdr") - for (call_name in names(calls)) { - call <- calls[[call_name]] - out0 <- call("none") - p0 <- out0$p - k <- nrow(out0) - expect_gt(k, 1) - for (m in methods) { - out <- call(m) - expect_equal(out$p, stats::p.adjust(p0, m), tolerance = 1e-10, label = paste(call_name, m)) - } - out <- call("sidak") - expect_equal(out$p, 1 - (1 - p0)^k, tolerance = 1e-10, label = paste(call_name, "sidak")) - # slopes are tested against 0, not compared pairwise, so "tukey" is - # replaced by "sidak" - out <- suppressMessages(call("tukey")) - expect_equal(out$p, 1 - (1 - p0)^k, tolerance = 1e-10, label = paste(call_name, "tukey")) - expect_false(isTRUE(all.equal(out$p, p0)), label = paste(call_name, "tukey changed")) - } + d$cyl_f <- factor(d$cyl) + m <- lm(mpg ~ wt * cyl_f, data = d) + p0 <- estimate_slopes(m, "wt", by = "cyl_f", p_adjust = "none")$p + out <- estimate_slopes(m, "wt", by = "cyl_f", p_adjust = "holm") + expect_equal(out$p, stats::p.adjust(p0, "holm"), tolerance = 1e-10) + out <- estimate_slopes(m, "wt", by = "cyl_f", p_adjust = "sidak") + expect_equal(out$p, 1 - (1 - p0)^3, tolerance = 1e-10) + + m <- lm(mpg ~ wt * hp, data = d) + p0 <- estimate_slopes(m, "wt", by = "hp", length = 4, p_adjust = "none")$p + out <- estimate_slopes(m, "wt", by = "hp", length = 4, p_adjust = "holm") + expect_equal(out$p, stats::p.adjust(p0, "holm"), tolerance = 1e-10) }) @@ -303,40 +281,33 @@ test_that("estimate_slopes, p_adjust = 'tukey' uses sidak for slopes", { d$cyl_f <- factor(d$cyl) m <- lm(mpg ~ wt * cyl_f, data = d) - p0 <- estimate_slopes(m, "wt", by = "cyl_f", p_adjust = "none")$p expect_message( out <- estimate_slopes(m, "wt", by = "cyl_f", p_adjust = "tukey"), regexp = "(?s)Tukey.*Sidak", perl = TRUE ) - expect_silent(estimate_slopes(m, "wt", by = "cyl_f", p_adjust = "tukey", verbose = FALSE)) - sidak <- estimate_slopes(m, "wt", by = "cyl_f", p_adjust = "sidak") - - expect_equal(out$p, sidak$p, tolerance = 1e-10) expect_identical(attr(out, "p_adjust"), "sidak") - expect_false(isTRUE(all.equal(out$p, p0))) + # emmeans replaces "tukey" by "sidak" for tests that are not pairwise + em <- suppressMessages(as.data.frame(emmeans::test( + emmeans::emtrends(m, ~cyl_f, var = "wt"), + adjust = "tukey" + ))) + expect_equal(out$p, em$p.value, tolerance = 1e-6) # without p-values, nothing is adjusted, and Sidak is not mentioned expect_no_message( estimate_slopes(m, "wt", by = "cyl_f", p_adjust = "tukey", ci = NA), message = "Sidak" ) - # closed form - expect_equal(out$p, 1 - (1 - p0)^3, tolerance = 1e-6) - # emmeans replaces "tukey" by "sidak" for tests that are not pairwise - em <- suppressMessages(as.data.frame( - emmeans::test(emmeans::emtrends(m, ~cyl_f, var = "wt"), adjust = "tukey") - )) - expect_equal(out$p, em$p.value, tolerance = 1e-6) # the exported getter adjusts the same way - gm_tukey <- suppressMessages( - get_marginaltrends(m, trend = "wt", by = "cyl_f", p_adjust = "tukey") - ) - gm_sidak <- get_marginaltrends(m, trend = "wt", by = "cyl_f", p_adjust = "sidak") - gm_none <- get_marginaltrends(m, trend = "wt", by = "cyl_f", p_adjust = "none") - expect_equal(gm_tukey$p.value, gm_sidak$p.value, tolerance = 1e-10) - expect_false(isTRUE(all.equal(gm_tukey$p.value, gm_none$p.value))) + gm <- suppressMessages(get_marginaltrends( + m, + trend = "wt", + by = "cyl_f", + p_adjust = "tukey" + )) + expect_equal(gm$p.value, em$p.value, tolerance = 1e-6) }) @@ -356,20 +327,15 @@ test_that("estimate_slopes, p_adjust = 'tukey' keeps tukey for pairwise comparis message = "Sidak" ) ref <- estimate_contrasts(m, contrast = "wt", by = "cyl_f", p_adjust = "tukey") - expect_identical(nrow(out), nrow(ref)) - expect_equal(out$Slope, ref$Difference, tolerance = 1e-6) expect_equal(out$p, ref$p, tolerance = 1e-6) # with `post_process`, the rows are also comparisons of slopes - expect_no_message( - out <- estimate_slopes( - m, - "wt", - by = "cyl_f", - post_process = ~pairwise, - p_adjust = "tukey" - ), - message = "Sidak" + out <- estimate_slopes( + m, + "wt", + by = "cyl_f", + post_process = ~pairwise, + p_adjust = "tukey" ) expect_equal(as.numeric(out$p), ref$p, tolerance = 1e-6) @@ -386,14 +352,15 @@ test_that("estimate_slopes, p_adjust with a by-variable named p", { d <- datasets::mtcars d$safe_var <- factor(d$cyl) d$p <- factor(d$cyl) - for (m in c("holm", "sidak")) { - ref <- estimate_slopes(lm(mpg ~ wt * safe_var, data = d), "wt", by = "safe_var", p_adjust = m) - out <- estimate_slopes(lm(mpg ~ wt * p, data = d), "wt", by = "p", p_adjust = m) - expect_true(is.factor(out$p.1), label = m) - expect_identical(levels(out$p.1), levels(d$p), label = m) - expect_identical(as.character(out$p.1), as.character(ref$safe_var), label = m) - expect_equal(out$p, ref$p, tolerance = 1e-10, label = m) - } + ref <- estimate_slopes( + lm(mpg ~ wt * safe_var, data = d), + "wt", + by = "safe_var", + p_adjust = "holm" + ) + out <- estimate_slopes(lm(mpg ~ wt * p, data = d), "wt", by = "p", p_adjust = "holm") + expect_identical(as.character(out$p.1), as.character(ref$safe_var)) + expect_equal(out$p, ref$p, tolerance = 1e-10) }) @@ -404,7 +371,6 @@ test_that("estimate_slopes, p_adjust counts one test per row", { # a single custom hypothesis is one test, so sidak does not change it p0 <- estimate_slopes(m, "wt", by = "cyl_f", comparison = "b2=b1")$p out <- estimate_slopes(m, "wt", by = "cyl_f", comparison = "b2=b1", p_adjust = "sidak") - expect_identical(nrow(out), 1L) expect_equal(out$p, p0, tolerance = 1e-10) # ordinal models have one row per response level and by-level @@ -413,9 +379,7 @@ test_that("estimate_slopes, p_adjust counts one test per row", { m <- MASS::polr(gear_o ~ wt + cyl_f, data = d, Hess = TRUE) p0 <- estimate_slopes(m, "wt", by = "cyl_f")$p out <- estimate_slopes(m, "wt", by = "cyl_f", p_adjust = "sidak") - k <- nrow(out) - expect_identical(k, 9L) - expect_equal(out$p, 1 - (1 - p0)^k, tolerance = 1e-10) + expect_equal(out$p, 1 - (1 - p0)^9, tolerance = 1e-10) }) diff --git a/tests/testthat/test-focal-name-clashes.R b/tests/testthat/test-focal-name-clashes.R index c60addf64..4779d23bc 100644 --- a/tests/testthat/test-focal-name-clashes.R +++ b/tests/testthat/test-focal-name-clashes.R @@ -1,19 +1,7 @@ skip_on_cran() skip_if_not_installed("marginaleffects", minimum_version = "0.29.0") -# Focal variables can have names that insight, parameters or modelbased use -# for their own columns. Each test renames the focal variable and compares the -# output with the output for a neutral name (`safe_var`). - -# names that `insight::standardize_names()` changes, or that modelbased drops -# by name -focal_reported_names <- c( - "type", "p", "CI", "S", "s.value", "predicted", "conf.level", "df_residual", - "df_error", "response", "n", "value", "component", "effect", "weight", - "weights", "sumsq", "meansq", "std.dev", "y.level", "hypothesis", "Median", - "Mean", "MAP", "Dxy", "rho", "r", "tau", "Difference", "Estimate", "F", - "Chi2", "chisq", "Chisq", "W" -) +# focal variables named like columns of insight, marginaleffects or modelbased focal_data <- function(name, numeric = FALSE) { d <- datasets::mtcars @@ -34,579 +22,138 @@ focal_calls <- list( m <- stats::lm(focal_formula(paste0("`", name, "` + hp")), data = d) estimate_means(m, name) }, - means_two = function(name) { - d <- focal_data(name) - m <- stats::lm(focal_formula(paste0("`", name, "` * am_f")), data = d) - estimate_means(m, c("am_f", name)) - }, slopes_by = function(name) { d <- focal_data(name) m <- stats::lm(focal_formula(paste0("`", name, "` * wt")), data = d) estimate_slopes(m, "wt", by = name) }, - contrasts_by = function(name) { + contrasts_by = function(name, ...) { d <- focal_data(name) m <- stats::lm(focal_formula(paste0("`", name, "` * am_f")), data = d) - estimate_contrasts(m, "am_f", by = name) + estimate_contrasts(m, "am_f", by = name, ...) }, - contrasts_by_filter = function(name) { + contrasts = function(name, ...) { d <- focal_data(name) m <- stats::lm(focal_formula(paste0("`", name, "` * am_f")), data = d) - estimate_contrasts(m, "am_f", by = paste0(name, "=c('4', '6')")) - }, - contrasts_by_filter_average = function(name) { - d <- focal_data(name) - m <- stats::lm(focal_formula(paste0("`", name, "` * am_f")), data = d) - estimate_contrasts( - m, - "am_f", - by = paste0(name, "='4'"), - estimate = "average" - ) - }, - contrasts = function(name) { - d <- focal_data(name) - m <- stats::lm(focal_formula(paste0("`", name, "` + hp")), data = d) - estimate_contrasts(m, name) - }, - means_covariate = function(name) { - d <- focal_data(name, numeric = TRUE) - m <- stats::lm(focal_formula(paste0("am_f + `", name, "`")), data = d) - estimate_means(m, "am_f") - }, - contrasts_covariate = function(name) { - d <- focal_data(name, numeric = TRUE) - m <- stats::lm(focal_formula(paste0("am_f + `", name, "`")), data = d) - estimate_contrasts(m, "am_f") - }, - slopes_covariate = function(name) { - d <- focal_data(name, numeric = TRUE) - d$hp_s <- d$hp / 100 - m <- stats::lm(focal_formula(paste0("hp_s * am_f + `", name, "`")), data = d) - estimate_slopes(m, "hp_s", by = "am_f") + estimate_contrasts(m, name, by = "am_f", ...) } ) -# Names that marginaleffects overwrites or uses itself, for each call. These -# calls give an error that asks users to rename the variable. -focal_error_cells <- list( - means_two = c("Predicted", "SE", "CI_low", "CI_high", "df", "rowid", "s.value"), - means_covariate = c("SE", "CI_low", "CI_high", "df"), - contrasts_covariate = c("SE", "CI_low", "CI_high", "df"), - contrasts_by = c("Predicted", "SE", "CI_low", "CI_high", "df", "rowid", "s.value"), - contrasts = c("SE", "CI_low", "CI_high", "df", "rowid"), - slopes_by = c("df", "rowid", "s.value", "predicted", "predicted_hi", "predicted_lo") -) - -# Returns NULL when `out` equals `ref` with "safe_var" replaced by `name`, -# otherwise a short description of the first difference. -focal_difference <- function(ref, out, name) { - if (inherits(out, "error")) { - return(paste("error:", conditionMessage(out))) - } - swap <- function(x) { - if (is.character(x)) gsub("safe_var", name, x, fixed = TRUE) else x - } - ref_names <- swap(colnames(ref)) - if (!identical(ref_names, colnames(out))) { - return(paste0( - "columns [", toString(colnames(out)), "], expected [", - toString(ref_names), "]" - )) - } - for (i in seq_along(ref)) { - r <- ref[[i]] - o <- out[[i]] - if (!identical(class(r), class(o))) { - return(paste0("column ", i, " has class ", class(o)[1])) - } - if (is.factor(r)) { - levels(r) <- swap(levels(r)) - } - r <- swap(r) - if (!isTRUE(all.equal(r, o, check.attributes = FALSE))) { - return(paste0("column ", i, " has other values")) - } - } - for (a in c("focal_terms", "by", "trend", "contrast")) { - if (!identical(swap(attr(ref, a, exact = TRUE)), attr(out, a, exact = TRUE))) { - return(paste0("attribute ", a, " differs")) - } - } - ref_grid <- swap(colnames(attr(ref, "datagrid", exact = TRUE))) - out_grid <- colnames(attr(out, "datagrid", exact = TRUE)) - if (!identical(ref_grid, out_grid)) { - return(paste0("datagrid columns [", toString(out_grid), "]")) - } - NULL -} - -expect_focal_names_match <- function(call_name, names, suffix = "") { - call <- focal_calls[[call_name]] - ref <- call("safe_var") - for (name in names) { - out <- tryCatch(suppressWarnings(call(name)), error = function(e) e) - diff <- focal_difference(ref, out, paste0(name, suffix)) - expect_null(diff, label = paste0(call_name, "(\"", name, "\"): ", diff)) - } +# `out` must equal the output for "safe_var", with that name replaced by +# `new_name` +expect_focal_output <- function(call, name, new_name = name, ...) { + ref <- as.data.frame(call("safe_var", ...)) + out <- as.data.frame(call(name, ...)) + expect_identical(colnames(out), gsub("safe_var", new_name, colnames(ref), fixed = TRUE)) + expect_equal(out, ref, ignore_attr = TRUE) } -test_that("names that marginaleffects overwrites give an informative error", { - n_cells <- 0 - for (call_name in names(focal_error_cells)) { - for (name in focal_error_cells[[call_name]]) { - n_cells <- n_cells + 1 - out <- tryCatch(suppressWarnings(focal_calls[[call_name]](name)), error = function(e) e) - label <- paste0(call_name, "(\"", name, "\")") - expect_true(inherits(out, "error"), label = label) - msg <- if (inherits(out, "error")) conditionMessage(out) else "" - expect_true(grepl(paste0("`", name, "`"), msg, fixed = TRUE), label = label) - expect_true(grepl("rename", msg, fixed = TRUE), label = label) - expect_false(grepl("Sorry, calculating", msg, fixed = TRUE), label = label) - } - } - expect_identical(n_cells, 33) -}) +test_that("focal variables keep names that insight renames", { + expect_focal_output(focal_calls$means, "type") + expect_focal_output(focal_calls$slopes_by, "type") + expect_focal_output(focal_calls$contrasts_by, "type") + expect_focal_output(focal_calls$contrasts_by, "hypothesis") - -test_that("names outside the error cells give the same output as safe names", { - grid_names <- c( - "Predicted", "SE", "CI_low", "CI_high", "df", "rowid", "s.value", - "predicted", "predicted_hi", "predicted_lo", "Parameter", "Level1", - "Level2", "hypothesis" - ) - grid_calls <- c( - "means_two", "means_covariate", "contrasts_by", "contrasts", "slopes_by", - "slopes_covariate" - ) - # a covariate named "rowid" is not part of this check - excluded <- list(means_covariate = "rowid", slopes_covariate = "rowid") - n_combinations <- 0 - for (call_name in grid_calls) { - call <- focal_calls[[call_name]] - ref <- call("safe_var") - skip_names <- c(focal_error_cells[[call_name]], excluded[[call_name]]) - for (name in setdiff(grid_names, skip_names)) { - n_combinations <- n_combinations + 1 - out <- tryCatch(suppressWarnings(call(name)), error = function(e) e) - diff <- focal_difference(ref, out, name) - # the variable may also get the suffix ".1" - if (!is.null(diff)) { - diff <- focal_difference(ref, out, paste0(name, ".1")) - } - expect_null(diff, label = paste0(call_name, "(\"", name, "\"): ", diff)) - } + call <- function(name) { + d <- focal_data(name, numeric = TRUE) + m <- stats::lm(focal_formula(paste0("`", name, "` + hp")), data = d) + newdata <- data.frame(x = c(2.123, 4.567), hp = 150) + colnames(newdata)[1] <- name + estimate_means(m, by = name, newdata = newdata) } - expect_identical(n_combinations, 53) -}) - - -test_that("estimate_means() keeps focal variables named like renamed columns", { - names <- setdiff(focal_reported_names, c("Mean", "s.value")) - expect_focal_names_match("means", names) - expect_focal_names_match("means_two", names) -}) - - -test_that("estimate_slopes() keeps by-variables named like renamed columns", { - names <- setdiff(focal_reported_names, c("p", "s.value", "predicted")) - expect_focal_names_match("slopes_by", names) -}) - - -test_that("estimate_contrasts() keeps by-variables named like renamed columns", { - names <- setdiff( - focal_reported_names, - c("p", "s.value", "hypothesis", "Difference") - ) - expect_focal_names_match("contrasts_by", names) -}) - - -# A by-variable with the name of an output column gets the suffix ".1", and the -# output column keeps its name. The expected output is the output for -# "safe_var", with "safe_var" replaced by ".1". - -# marginaleffects renames columns named "Predicted", "SE", "CI_low" and -# "CI_high" for means and contrasts, so these names are not tested there. - -test_that("estimate_means() suffixes focal variables named like output columns", { - names <- c("Mean", "t") - expect_focal_names_match("means", names, suffix = ".1") - expect_focal_names_match("means_two", names, suffix = ".1") + expect_focal_output(call, "type") }) -test_that("estimate_slopes() suffixes by-variables named like output columns", { - names <- c("Slope", "SE", "CI_low", "CI_high", "t", "p") - expect_focal_names_match("slopes_by", names, suffix = ".1") +test_that("focal variables named like output columns get the suffix '.1'", { + expect_focal_output(focal_calls$means, "t", "t.1") + expect_focal_output(focal_calls$slopes_by, "Slope", "Slope.1") + expect_focal_output(focal_calls$contrasts_by, "p", "p.1") + expect_focal_output(focal_calls$contrasts_by, "Level1", "Level1.1") }) -test_that("estimate_contrasts() suffixes by-variables named like output columns", { - names <- c("Difference", "t", "p") - expect_focal_names_match("contrasts_by", names, suffix = ".1") -}) - - -test_that("estimate_contrasts() suffixes by-variables named Level1 or Level2", { - names <- c("Level1", "Level2") - expect_focal_names_match("contrasts_by", names, suffix = ".1") - for (name in names) { - out <- focal_calls$contrasts_by(name) - expect_false(anyDuplicated(colnames(out)) > 0, label = name) +test_that("estimate_contrasts() filters suffixed by-variables", { + call <- function(name, estimate) { + d <- focal_data(name) + m <- stats::lm(focal_formula(paste0("`", name, "` * am_f")), data = d) + estimate_contrasts(m, "am_f", by = paste0(name, "=c('4', '6')"), estimate = estimate) } -}) - - -test_that("estimate_contrasts() keeps by-variables named hypothesis", { - expect_focal_names_match("contrasts_by", "hypothesis") - out <- focal_calls$contrasts_by("hypothesis") - expect_true(all(c("Level1", "Level2", "hypothesis") %in% colnames(out))) -}) - - -test_that("estimate_contrasts() works for a contrast variable named Parameter", { - expect_focal_names_match("contrasts", "Parameter") + expect_focal_output(call, "p", "p.1", estimate = "typical") + expect_focal_output(call, "p", "p.1", estimate = "average") }) test_that("the suffix skips names that are already taken", { - two_focal <- function(names, call) { - d <- datasets::mtcars - d[[names[1]]] <- factor(d$cyl) - d[[names[2]]] <- factor(d$am) - rhs <- paste0("wt * `", names[1], "` + `", names[2], "`") - m <- stats::lm(focal_formula(rhs), data = d) - call(m, names) - } - slopes <- function(m, names) estimate_slopes(m, "wt", by = names) - means <- function(m, names) estimate_means(m, names) - - ref <- two_focal(c("s1", "s2"), slopes) - out <- two_focal(c("p", "p.1"), slopes) + d <- datasets::mtcars + d$p <- factor(d$cyl) + d$p.1 <- factor(d$am) + out <- estimate_slopes( + stats::lm(mpg ~ wt * p + p.1, data = d), + "wt", + by = c("p", "p.1") + ) expect_identical(colnames(out)[1:2], c("p.2", "p.1")) - expect_identical(as.character(out$p.2), as.character(ref$s1)) - expect_identical(as.character(out$p.1), as.character(ref$s2)) - expect_equal(out$Slope, ref$Slope, tolerance = 1e-10) - expect_equal(out$p, ref$p, tolerance = 1e-10) - - ref <- two_focal(c("s1", "s2"), means) - out <- two_focal(c("t", "t.1"), means) - expect_identical(colnames(out)[1:2], c("t.2", "t.1")) - expect_identical(as.character(out$t.2), as.character(ref$s1)) - expect_identical(as.character(out$t.1), as.character(ref$s2)) - expect_equal(out$Mean, ref$Mean, tolerance = 1e-10) - expect_equal(out$t, ref$t, tolerance = 1e-10) -}) - - -test_that("estimate_contrasts() filters suffixed by-variables", { - names <- c("Difference", "t", "p") - expect_focal_names_match("contrasts_by_filter", names, suffix = ".1") - expect_focal_names_match("contrasts_by_filter_average", names, suffix = ".1") }) test_that("print() and plot() use the suffixed focal column", { skip_if_not_installed("ggplot2") skip_if_not_installed("see") - for (name in c("Mean", "t")) { - out <- focal_calls$means(name) - expect_no_error(utils::capture.output(print(out))) - p <- plot(out) - expect_no_error(ggplot2::ggplot_build(p)) - x_mappings <- unlist(lapply(p$layers, function(l) { - if (is.null(l$mapping$x)) NULL else ggplot2::as_label(l$mapping$x) - })) - expect_true( - length(x_mappings) > 0 && all(x_mappings == paste0(name, ".1")), - label = paste0("x mappings for \"", name, "\": ", toString(x_mappings)) - ) - } + d <- focal_data("t") + attr(d$t, "label") <- "Cylinders" + out <- estimate_means(stats::lm(mpg ~ t + hp, data = d), "t") + expect_no_error(utils::capture.output(print(out))) + p <- plot(out) + expect_identical(ggplot2::as_label(p$layers[[1]]$mapping$x), "t.1") + expect_identical(p$labels$x, "Cylinders") }) -test_that("plot() uses the variable label of a suffixed focal variable", { - skip_if_not_installed("ggplot2") - skip_if_not_installed("see") - for (name in c("safe_var", "Mean", "t")) { - d <- focal_data(name) - attr(d[[name]], "label") <- "Cylinders" - m <- stats::lm(focal_formula(paste0("`", name, "` + hp")), data = d) - p <- plot(estimate_means(m, name)) - expect_identical(p$labels$x, "Cylinders", label = paste0("x label for \"", name, "\"")) - } -}) - - -test_that("focal variable names that worked before still work", { - # names that matched on upstream main b964ea75, computed once from a list - # of 99 candidate names (the reported names and other column, argument and - # output names) - matched_before <- list( - means = c( - "z", "rowidcf", "wts", "predicted_lo", "predicted_hi", "emmean", - "prob", "rate", "lower.CL", "upper.CL", "asymp.LCL", "asymp.UCL", - "Level1", "Level2", "Parameter", "Comparison", "Response", "Threshold", - "Row", "Residuals", "Group", "Level", "Slope", "Coefficient", - "Probability", "Contrast", "Component", "Effects", "Type", "Ratio", - "iter_1", "x", "model", "data", "at", "ci", "trend", "predict", - "comparison", "length", "range", "fun", "i", "out", "means", "grid" - ), - means_two = c( - "z", "rowidcf", "wts", "predicted_lo", "predicted_hi", "emmean", - "prob", "rate", "lower.CL", "upper.CL", "asymp.LCL", "asymp.UCL", - "Level1", "Level2", "Parameter", "Comparison", "Response", "Threshold", - "Row", "Residuals", "Group", "Level", "Slope", "Coefficient", - "Probability", "Contrast", "Component", "Effects", "Type", "Ratio", - "iter_1", "x", "model", "data", "at", "ci", "trend", "predict", - "comparison", "length", "range", "fun", "i", "out", "means", "grid" - ), - slopes_by = c( - "rowidcf", "wts", "emmean", "prob", "rate", "lower.CL", "upper.CL", - "asymp.LCL", "asymp.UCL", "Level1", "Level2", "Response", "Threshold", - "Row", "Residuals", "Group", "Level", "Probability", "Contrast", - "Component", "Effects", "Type", "Ratio", "iter_1", "x", "model", - "data", "at", "ci", "trend", "predict", "comparison", "length", - "range", "fun", "i", "out", "means", "grid" - ), - contrasts_by = c( - "z", "rowidcf", "wts", "predicted_lo", "predicted_hi", "emmean", - "prob", "rate", "lower.CL", "upper.CL", "asymp.LCL", "asymp.UCL", - "Comparison", "Response", "Threshold", "Row", - "Residuals", "Group", "Level", "Probability", "Contrast", "Component", - "Effects", "Type", "Ratio", "iter_1", "x", "model", "data", "at", "ci", - "trend", "predict", "comparison", "length", "range", "fun", "i", "out", - "means", "grid" - ), - contrasts = c( - "z", "F", "Chi2", "chisq", "Chisq", "W", "Median", "Mean", "MAP", - "Dxy", "rho", "r", "tau", "df_residual", "df_error", "std.dev", - "s.value", "component", "effect", "predicted", "response", "y.level", - "conf.level", "n", "type", "weight", "weights", "hypothesis", "value", - "sumsq", "meansq", "Estimate", "rowidcf", "wts", "predicted_lo", - "predicted_hi", "emmean", "prob", "rate", "lower.CL", "upper.CL", - "asymp.LCL", "asymp.UCL", "Predicted", "CI", "Level1", "Level2", - "Comparison", "Response", "Threshold", "Row", "Residuals", "Group", - "Level", "Slope", "Coefficient", "Probability", "Contrast", - "Component", "Effects", "Type", "Ratio", "S", "iter_1", "x", "model", - "data", "at", "ci", "trend", "predict", "comparison", "length", - "range", "fun", "i", "out", "means", "grid" - ) +test_that("names that marginaleffects overwrites give an informative error", { + d <- focal_data("SE") + m <- stats::lm(mpg ~ SE * am_f, data = d) + expect_error(estimate_means(m, c("am_f", "SE")), "Please rename", fixed = TRUE) + d <- focal_data("df", numeric = TRUE) + m <- stats::lm(mpg ~ am_f + df, data = d) + expect_error(estimate_means(m, "am_f"), "Please rename", fixed = TRUE) + d <- focal_data("predicted_hi") + m <- stats::lm(mpg ~ predicted_hi * wt, data = d) + expect_error( + estimate_slopes(m, "wt", by = "predicted_hi"), + "Please rename", + fixed = TRUE ) - for (call_name in names(matched_before)) { - expect_focal_names_match(call_name, matched_before[[call_name]]) - } -}) - - -test_that("estimate_means() keeps numeric focal variables for glm", { - names <- setdiff(focal_reported_names, c("Mean", "s.value")) - call <- function(name) { - d <- focal_data(name, numeric = TRUE) - d$vs <- datasets::mtcars$vs - m <- stats::glm( - stats::as.formula(paste0("vs ~ am_f * `", name, "`")), - family = stats::binomial(), - data = d - ) - estimate_means(m, c("am_f", name)) - } - ref <- suppressWarnings(call("safe_var")) - for (name in names) { - out <- tryCatch(suppressWarnings(call(name)), error = function(e) e) - diff <- focal_difference(ref, out, name) - expect_null(diff, label = paste0("glm(\"", name, "\"): ", diff)) - } -}) - - -test_that("estimate_means() keeps focal variables with custom newdata", { - # the values in `newdata` are not part of the generated data grid - names <- setdiff(focal_reported_names, c("Mean", "s.value")) - call <- function(name) { - d <- focal_data(name, numeric = TRUE) - m <- stats::lm(focal_formula(paste0("`", name, "` + hp")), data = d) - newdata <- data.frame(x = c(2.123, 4.567), hp = 150) - colnames(newdata)[1] <- name - estimate_means(m, by = name, newdata = newdata) - } - ref <- call("safe_var") - for (name in names) { - out <- tryCatch(suppressWarnings(call(name)), error = function(e) e) - diff <- focal_difference(ref, out, name) - expect_null(diff, label = paste0("newdata(\"", name, "\"): ", diff)) - } -}) - - -# Contrasts, ratios, custom hypotheses and slope contrasts, with focal -# variables named like output columns. The output column with that name must -# keep its values and its position. - -clash_names <- c( - "Difference", "Ratio", "Slope", "Mean", "SE", "CI_low", "CI_high", "t", "z", - "df", "p", "Predicted", "Parameter", "Level1", "Statistic", "s.value", - "hypothesis" -) - -clash_data <- function(name) { - d <- datasets::mtcars - d$am_f <- factor(d$am, levels = 0:1, labels = c("auto", "manual")) - d[[name]] <- factor(d$cyl) - d -} - -clash_model <- function(name, rhs) { - d <- clash_data(name) - stats::lm(focal_formula(paste0(rhs, " * `", name, "`")), data = d) -} - -clash_calls <- list( - contrast = function(name) { - estimate_contrasts(clash_model(name, "am_f"), name, by = "am_f") - }, - contrasts_by = function(name) { - estimate_contrasts(clash_model(name, "am_f"), "am_f", by = name) - }, - contrasts_by_ratio = function(name) { - estimate_contrasts( - clash_model(name, "am_f"), - "am_f", - by = name, - comparison = "ratio" - ) - }, - contrasts_by_hypothesis = function(name) { - estimate_contrasts( - clash_model(name, "am_f"), - "am_f", - by = name, - comparison = "b2=b1" - ) - }, - slope_contrasts_by = function(name) { - estimate_contrasts(clash_model(name, "wt"), "wt", by = name) - }, - slopes_by_hypothesis = function(name) { - estimate_slopes(clash_model(name, "wt"), "wt", by = name, comparison = "b2=b1") - } -) - -# names that give the error of `.check_marginaleffects_names()` -clash_error_cells <- list( - contrast = c("SE", "CI_low", "CI_high", "df"), - contrasts_by = c("Predicted", "SE", "CI_low", "CI_high", "df", "s.value"), - contrasts_by_ratio = c("Predicted", "SE", "CI_low", "CI_high", "df", "s.value"), - contrasts_by_hypothesis = c("Predicted", "SE", "CI_low", "CI_high", "df", "s.value"), - slope_contrasts_by = c("df", "s.value"), - slopes_by_hypothesis = c("df", "s.value") -) - -# Returns NULL when `out` equals `ref` with "safe_var" replaced by `name`, -# otherwise a short description of the difference. Label columns (factors and -# character vectors) are compared as character. -clash_difference <- function(ref, out, name) { - if (inherits(out, "error")) { - return(paste("error:", conditionMessage(out))) - } - if (anyDuplicated(colnames(out))) { - return(paste0("duplicated columns [", toString(colnames(out)), "]")) - } - expected <- gsub("safe_var", name, colnames(ref), fixed = TRUE) - if (!identical(expected, colnames(out))) { - return(paste0( - "columns [", toString(colnames(out)), "], expected [", - toString(expected), "]" - )) - } - for (i in seq_along(ref)) { - r <- ref[[i]] - o <- out[[i]] - if (is.factor(r) || is.character(r)) { - r <- gsub("safe_var", name, as.character(r), fixed = TRUE) - o <- if (is.factor(o) || is.character(o)) as.character(o) else o - } - if (!identical(class(r), class(o))) { - return(paste0("column ", colnames(out)[i], " has class ", class(o)[1])) - } - if (!isTRUE(all.equal(r, o, check.attributes = FALSE))) { - return(paste0("column ", colnames(out)[i], " has other values")) - } - } - NULL -} - -test_that("output columns keep their values for focal names like output columns", { - n_errors <- 0 - n_compared <- 0 - for (call_name in names(clash_calls)) { - call <- clash_calls[[call_name]] - ref <- call("safe_var") - for (name in clash_names) { - out <- tryCatch(suppressWarnings(call(name)), error = function(e) e) - label <- paste0(call_name, "(\"", name, "\")") - if (name %in% clash_error_cells[[call_name]]) { - n_errors <- n_errors + 1 - msg <- if (inherits(out, "error")) conditionMessage(out) else "no error" - msg <- gsub("\\s+", " ", msg) - expect_match( - msg, - paste0("Variable `", name, "` has a name that the marginaleffects package uses"), - fixed = TRUE, - label = label - ) - next - } - n_compared <- n_compared + 1 - # the variable has its own name, or a suffix ".1", ".2", ... - for (new_name in c(name, paste0(name, ".", 1:3))) { - diff <- clash_difference(ref, out, new_name) - if (is.null(diff)) { - break - } - } - if (!is.null(diff)) { - diff <- clash_difference(ref, out, name) - } - expect_null(diff, label = paste0(label, ": ", diff)) - } - } - expect_identical(n_errors, 26) - expect_identical(n_compared, 76) }) -test_that("ratio contrasts as formula suffix a by-variable named Ratio", { - call <- function(name) { - estimate_contrasts( - clash_model(name, "am_f"), - "am_f", - by = name, - comparison = ratio ~ pairwise - ) - } - ref <- call("safe_var") - out <- call("Ratio") - expect_null(clash_difference(ref, out, "Ratio.1")) +test_that("output columns of contrasts keep their values for focal names like output columns", { + expect_focal_output(focal_calls$contrasts, "Difference") + expect_focal_output(focal_calls$contrasts_by, "Slope") + expect_focal_output(focal_calls$contrasts_by, "Ratio", "Ratio.1", comparison = "ratio") + expect_focal_output( + focal_calls$contrasts_by, + "Ratio", + "Ratio.1", + comparison = ratio ~ pairwise + ) }) test_that("p_adjust in contrasts of estimate_relation() adjusts the p-values", { - call <- function(name, p_adjust) { - d <- clash_data(name) - m <- stats::lm(focal_formula(paste0("am_f * `", name, "`")), data = d) + d <- focal_data("p.value") + m <- stats::lm(mpg ~ am_f * p.value, data = d) + call <- function(p_adjust) { estimate_contrasts( - estimate_relation(m, by = c("am_f", name)), + estimate_relation(m, by = c("am_f", "p.value")), "am_f", - by = name, + by = "p.value", p_adjust = p_adjust ) } - out0 <- call("p.value", "none") - out <- call("p.value", "holm") + out0 <- call("none") + out <- call("holm") expect_identical(as.character(out$p.value), c("4", "6", "8")) expect_equal(out$p, stats::p.adjust(out0$p, "holm"), tolerance = 1e-10) }) diff --git a/tests/testthat/test-ordinal.R b/tests/testthat/test-ordinal.R index a28649be0..5e3edf021 100644 --- a/tests/testthat/test-ordinal.R +++ b/tests/testthat/test-ordinal.R @@ -233,173 +233,34 @@ 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") - reversed_levels <- c("5", "4", "3", "2", "1") + threshold_names <- paste(text_levels[-5], text_levels[-1], sep = "|") 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 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) - 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, 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", "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, - 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) - } -}) + m <- ordinal::clmm(rating_text ~ temp + contact + (1 | judge), data = wine) - -test_that("estimate_means, backend emmeans, ordinal level names unchanged", { - skip_if_not_installed("emmeans") - skip_if_not_installed("ordinal") - 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 - ), - clmm = list( - ordinal::clmm(rating ~ temp + contact + (1 | judge), data = wine), - c("1", "2", "3", "4", "5") - ) + # emmeans codes the response categories of clmm models as 1, 2, 3, ... + out <- estimate_means(m, "temp", backend = "emmeans") + expected <- as.data.frame(suppressMessages(emmeans::emmeans( + m, + ~ temp | rating_text, + mode = "prob" + ))) + expect_identical( + as.character(out$Response), + text_levels[as.integer(as.character(expected$rating_text))] ) - # 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, - data = housing, - family = glmmTMB::ordinal() - ), - sat_levels - ) - } + expect_identical(levels(out$Response), text_levels) + out <- estimate_means(m, "temp", predict = "cum.prob", backend = "emmeans") + expect_identical(levels(out$Threshold), threshold_names) - 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) - ) - } - } + # `names(m$alpha)` are not the thresholds for `threshold = "symmetric"` + m <- ordinal::clmm( + rating_text ~ temp + (1 | judge), + data = wine, + threshold = "symmetric" + ) + out <- estimate_means(m, "temp", predict = "cum.prob", backend = "emmeans") + expect_identical(levels(out$Threshold), threshold_names) }) @@ -420,162 +281,31 @@ test_that("estimate_means, print bracl", { }) -# contrasts from emmeans, with the contrast split into `Level1` and `Level2`, -# and the values of the pseudo-factor (response categories or thresholds) in -# `group` -.emmeans_ordinal_contrasts <- function(m, specs, mode, group = NULL) { - contrasts <- emmeans::contrast( - suppressMessages(emmeans::emmeans(m, specs, mode = mode)), - "revpairwise", - adjust = "none" - ) - out <- as.data.frame(contrasts) - ci <- as.data.frame(stats::confint(contrasts, adjust = "none")) - contrast_levels <- do.call( - rbind, - strsplit(as.character(out$contrast), " - ", fixed = TRUE) - ) - data.frame( - Level1 = contrast_levels[, 1], - Level2 = contrast_levels[, 2], - group = if (is.null(group)) "" else as.character(out[[group]]), - estimate = out$estimate, - SE = out$SE, - p = out$p.value, - CI_low = ci$asymp.LCL, - CI_high = ci$asymp.UCL, - stringsAsFactors = FALSE - ) -} - - -# rows of `out` in the order of `expected`, matched by levels and the group -# column (`Response` or `Threshold`) -.match_ordinal_contrasts <- function(out, expected, group = NULL) { - out_group <- if (is.null(group)) "" else as.character(out[[group]]) - out_keys <- paste(out$Level1, out$Level2, out_group) - expected_keys <- paste(expected$Level1, expected$Level2, expected$group) - as.data.frame(out)[match(expected_keys, out_keys), ] -} - - -test_that("estimate_contrasts, backend emmeans, ordinal, probability default", { +test_that("estimate_contrasts, backend emmeans, ordinal modes", { skip_if_not_installed("emmeans") - skip_if_not_installed("ordinal") data(housing, package = "MASS") - data(wine, package = "ordinal") - - models <- list( - polr = list( - MASS::polr(Sat ~ Infl + Type + Cont, weights = Freq, data = housing, Hess = TRUE), - "Infl", - "Sat" - ), - clm = list( - ordinal::clm(Sat ~ Infl + Type + Cont, weights = Freq, data = housing), - "Infl", - "Sat" - ), - clmm = list( - ordinal::clmm(rating ~ temp + contact + (1 | judge), data = wine), - "temp", - "rating" - ) - ) - # 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, - data = housing, - family = glmmTMB::ordinal() - ), - "Infl", - "Sat" - ) - } - - for (i in names(models)) { - m <- models[[i]][[1]] - focal <- models[[i]][[2]] - response <- models[[i]][[3]] - expected <- .emmeans_ordinal_contrasts( - m, - stats::as.formula(paste("~", focal, "|", response)), - "prob", - response - ) - - # differences of probabilities within each response category, for the - # default and for predict = "response" - for (predict in list(NULL, "response")) { - out <- estimate_contrasts(m, focal, predict = predict, backend = "emmeans") - expect_identical(nrow(out), nrow(expected), info = i) - expect_true("Response" %in% colnames(out), info = i) - matched <- .match_ordinal_contrasts(out, expected, "Response") - expect_false(anyNA(matched$Difference), info = i) - expect_equal(matched$Difference, expected$estimate, tolerance = 1e-6, info = i) - expect_equal(matched$SE, expected$SE, tolerance = 1e-6, info = i) - expect_equal(matched$p, expected$p, tolerance = 1e-6, info = i) - expect_equal(matched$CI_low, expected$CI_low, tolerance = 1e-6, info = i) - expect_equal(matched$CI_high, expected$CI_high, tolerance = 1e-6, info = i) - } - - # same differences as the marginaleffects backend, for contrasts within one - # response category. marginaleffects does not support clmm, and ordinal - # glmmTMB requires marginaleffects >= 1.0.1 - if (i %in% c("polr", "clm")) { - out <- as.data.frame(estimate_contrasts(m, focal)) - out <- out[as.character(out$Response1) == as.character(out$Response2), ] - out$Response <- out$Response1 - matched <- .match_ordinal_contrasts(out, expected, "Response") - expect_false(anyNA(matched$Difference), info = i) - expect_equal(matched$Difference, expected$estimate, tolerance = 1e-4, info = i) - } + m <- MASS::polr(Sat ~ Infl + Type + Cont, weights = Freq, data = housing, Hess = TRUE) + emmeans_contrasts <- function(specs, mode) { + as.data.frame(emmeans::contrast( + suppressMessages(emmeans::emmeans(m, specs, mode = mode)), + "revpairwise", + adjust = "none" + )) } -}) + out <- estimate_contrasts(m, "Infl", backend = "emmeans") + expected <- emmeans_contrasts(~ Infl | Sat, "prob") + expect_identical(as.character(out$Response), as.character(expected$Sat)) + expect_equal(out$Difference, expected$estimate, tolerance = 1e-6) + expect_equal(out$SE, expected$SE, tolerance = 1e-6) -test_that("estimate_contrasts, backend emmeans, ordinal modes", { - skip_if_not_installed("emmeans") - data(housing, package = "MASS") - m <- MASS::polr(Sat ~ Infl + Type + Cont, weights = Freq, data = housing, Hess = TRUE) + out <- estimate_contrasts(m, "Infl", predict = "cum.prob", backend = "emmeans") + expected <- emmeans_contrasts(~ Infl | cut, "cum.prob") + expect_identical(as.character(out$Threshold), as.character(expected$cut)) + expect_equal(out$Difference, expected$estimate, tolerance = 1e-6) - # the column that holds the pseudo-factor, if any - group_columns <- list( - latent = NULL, - linear.predictor = "Threshold", - cum.prob = "Threshold", - exc.prob = "Threshold", - prob = "Response", - mean.class = NULL - ) - for (mode in names(group_columns)) { - group <- group_columns[[mode]] - pseudo_factor <- switch(mode, prob = "Sat", latent = , mean.class = NULL, "cut") - specs <- if (is.null(pseudo_factor)) { - ~Infl - } else { - stats::as.formula(paste("~ Infl |", pseudo_factor)) - } - expected <- .emmeans_ordinal_contrasts(m, specs, mode, pseudo_factor) - - out <- estimate_contrasts(m, "Infl", predict = mode, backend = "emmeans") - expect_identical( - intersect(c("Response", "Threshold"), colnames(out)), - if (is.null(group)) character(0) else group, - info = mode - ) - expect_identical(nrow(out), nrow(expected), info = mode) - matched <- .match_ordinal_contrasts(out, expected, group) - expect_false(anyNA(matched$Difference), info = mode) - expect_equal(matched$Difference, expected$estimate, tolerance = 1e-6, info = mode) - expect_equal(matched$SE, expected$SE, tolerance = 1e-6, info = mode) - } + out <- estimate_contrasts(m, "Infl", predict = "mean.class", backend = "emmeans") + expect_equal(out$Difference, c(0.5528118, 0.3039427, 0.2488691), tolerance = 1e-6) }) @@ -587,354 +317,90 @@ test_that("estimate_contrasts, backend emmeans, ordinal, numeric by", { m <- MASS::polr(Sat ~ Infl * x, weights = Freq, data = housing, Hess = TRUE) out <- estimate_contrasts(m, "Infl", by = "x", length = 3, backend = "emmeans") - # 3 pairs of levels, 3 values of x, 3 response categories - expect_identical(nrow(out), 27L) - expect_length(unique(out$x), 3) - expect_equal(range(out$x), range(housing$x), tolerance = 1e-3) - expect_identical( - anyDuplicated(as.data.frame(out)[c("Level1", "Level2", "x", "Response")]), - 0L - ) - - expected <- emmeans::contrast( + expected <- as.data.frame(emmeans::contrast( emmeans::emmeans(m, ~ Infl | x + Sat, at = list(x = unique(out$x)), mode = "prob"), "revpairwise", adjust = "none" - ) - ci <- as.data.frame(stats::confint(expected)) - expected <- as.data.frame(expected) - key_out <- paste(out$Level1, out$Level2, format(out$x), out$Response) - key_expected <- paste( - gsub(" - .*", "", expected$contrast), - gsub(".* - ", "", expected$contrast), - format(expected$x), - expected$Sat - ) - matched <- match(key_expected, key_out) - expect_false(anyNA(matched)) - # rows are in the order of emmeans - expect_identical(matched, seq_len(nrow(out))) - expect_equal(out$Difference[matched], expected$estimate, tolerance = 1e-6) - expect_equal(out$SE[matched], expected$SE, tolerance = 1e-6) - expect_equal(out$CI_low[matched], ci$asymp.LCL, tolerance = 1e-6) - expect_equal(out$CI_high[matched], ci$asymp.UCL, tolerance = 1e-6) + )) + # rows in the order of emmeans + expect_identical(nrow(out), 27L) + expect_equal(out$x, expected$x, tolerance = 1e-6) + expect_equal(out$Difference, expected$estimate, tolerance = 1e-6) + expect_equal(out$SE, expected$SE, tolerance = 1e-6) }) -test_that("estimate_contrasts, backend emmeans, clmm shows level names", { +test_that("estimate_contrasts and estimate_slopes, backend emmeans, clmm shows level names", { skip_if_not_installed("emmeans") skip_if_not_installed("ordinal") data(wine, package = "ordinal") - wine$rating_text <- factor(wine$rating, labels = letters[1:5], ordered = TRUE) - all_levels <- list(rating = as.character(1:5), rating_text = letters[1:5]) - - 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) - - out <- estimate_contrasts(m, "temp", backend = "emmeans") - expect_identical(nrow(out), 5L, info = response) - expect_setequal(as.character(out$Response), response_levels) - expect_identical(levels(out$Response), response_levels, info = response) - - for (mode in c("cum.prob", "exc.prob", "linear.predictor")) { - out <- estimate_contrasts(m, "temp", predict = mode, backend = "emmeans") - expect_identical(nrow(out), 4L, info = paste(response, mode)) - expect_setequal(as.character(out$Threshold), threshold_names) - expect_identical( - levels(out$Threshold), - threshold_names, - info = paste(response, mode) - ) - } - } + set.seed(3) + wine$x <- stats::rnorm(nrow(wine)) + as.numeric(wine$temp) + levels(wine$rating) <- letters[1:5] + threshold_names <- c("a|b", "b|c", "c|d", "d|e") + m <- ordinal::clmm(rating ~ temp + x * contact + (1 | judge), data = wine) + + out <- estimate_contrasts(m, "temp", backend = "emmeans") + expect_identical(levels(out$Response), letters[1:5]) + out <- estimate_contrasts(m, "temp", predict = "cum.prob", backend = "emmeans") + expect_identical(levels(out$Threshold), threshold_names) + + out <- estimate_slopes(m, "x", by = "contact", backend = "emmeans") + expect_identical(levels(out$Response), letters[1:5]) + out <- estimate_slopes( + m, + "x", + by = "contact", + predict = "cum.prob", + backend = "emmeans" + ) + expect_identical(levels(out$Threshold), threshold_names) + out <- estimate_contrasts(m, "x", by = "contact", backend = "emmeans") + expect_identical(levels(out$Response), letters[1:5]) }) -# polr fits with a numeric predictor `x` and a factor `g`, for slopes. For the -# second fit, `emmeans::emtrends(mode = "prob")` fails with "Irregular -# reference grid" (emmeans 2.0.4) -.ordinal_slope_models <- function() { - data(housing, package = "MASS") - set.seed(1) - housing$x <- as.numeric(housing$Cont) + stats::rnorm(nrow(housing)) - housing$g <- housing$Infl +# polr fit with a numeric predictor `x` and a factor `g`, for slopes. +# `emmeans::emtrends(mode = "prob")` fails for it with "Irregular reference +# grid" (emmeans 2.0.4) +.ordinal_slope_model <- function() { set.seed(2) n <- 400 - simulated <- data.frame( + d <- data.frame( x = stats::rnorm(n), g = factor(sample(c("a", "b", "c"), n, replace = TRUE)) ) - simulated$y <- cut( - simulated$x + as.numeric(simulated$g) + stats::rlogis(n), + d$y <- cut( + d$x + as.numeric(d$g) + stats::rlogis(n), c(-Inf, 1, 2.5, Inf), labels = c("lo", "mid", "hi"), ordered_result = TRUE ) - list( - housing = MASS::polr(Sat ~ g * x, weights = Freq, data = housing, Hess = TRUE), - simulated = MASS::polr(y ~ x * g, data = simulated, Hess = TRUE) - ) + MASS::polr(y ~ x * g, data = d, Hess = TRUE) } test_that("estimate_slopes, backend emmeans, ordinal", { skip_if_not_installed("emmeans") - models <- .ordinal_slope_models() - - for (i in names(models)) { - m <- models[[i]] - - # slopes of the probability of each response category, for the default - # and for predict = "prob". P(y = k) = P(y <= k) - P(y <= k - 1), so each - # slope is the slope of the cumulative probability at the threshold above - # the category, less the one at the threshold below (0 if there is none) - cumulative <- as.data.frame(suppressMessages(emmeans::emtrends( - m, - ~ g | cut, - var = "x", - mode = "cum.prob" - ))) - response_levels <- m$lev - n_cuts <- length(response_levels) - 1 - cumulative_slope <- function(g, k) { - if (k < 1 || k > n_cuts) { - return(0) - } - cumulative$x.trend[ - as.character(cumulative$g) == g & as.integer(cumulative$cut) == k - ] - } - out_marginaleffects <- as.data.frame(estimate_slopes(m, "x", by = "g")) - for (predict in list(NULL, "prob")) { - out <- as.data.frame(estimate_slopes( - m, - "x", - by = "g", - predict = predict, - backend = "emmeans" - )) - expect_identical(nrow(out), 9L, info = i) - expect_true("Response" %in% colnames(out), info = i) - expected <- mapply( - function(g, response) { - k <- match(response, response_levels) - cumulative_slope(g, k) - cumulative_slope(g, k - 1) - }, - as.character(out$g), - as.character(out$Response) - ) - expect_equal(out$Slope, unname(expected), tolerance = 1e-6, info = i) - # marginaleffects computes the slopes with a smaller step than emtrends(), - # which uses a forward difference with a step of 0.001 times the range - # of `x`. For the simulated fit, single slopes differ by up to 1.2% - # (emmeans 2.0.4, marginaleffects 1.0.0) - matched <- out_marginaleffects[ - match( - paste(out$g, out$Response), - paste(out_marginaleffects$g, out_marginaleffects$Response) - ), - ] - expect_false(anyNA(matched$Slope), info = i) - expect_equal(out$Slope, matched$Slope, tolerance = 1e-2, info = i) - expect_equal(out$SE, matched$SE, tolerance = 1e-2, info = i) - # the probabilities sum to 1, so their slopes sum to 0 - expect_equal( - as.vector(tapply(out$Slope, out$g, sum)), - c(0, 0, 0), - tolerance = 1e-8, - info = i - ) - } - - # modes with slopes from emtrends() - for (mode in c("latent", "linear.predictor", "cum.prob", "exc.prob")) { - specs <- if (mode == "latent") ~g else ~ g | cut - expected <- as.data.frame(suppressMessages(emmeans::emtrends( - m, - specs, - var = "x", - mode = mode - ))) - expected_group <- if (mode == "latent") "" else as.character(expected$cut) - out <- as.data.frame(estimate_slopes( - m, - "x", - by = "g", - predict = mode, - backend = "emmeans" - )) - expect_identical( - intersect(c("Response", "Threshold"), colnames(out)), - if (mode == "latent") character(0) else "Threshold", - info = paste(i, mode) - ) - out_group <- if (mode == "latent") "" else as.character(out$Threshold) - matched <- out[match(paste(expected$g, expected_group), paste(out$g, out_group)), ] - expect_identical(nrow(out), nrow(expected), info = paste(i, mode)) - expect_false(anyNA(matched$Slope), info = paste(i, mode)) - expect_equal( - matched$Slope, - expected$x.trend, - tolerance = 1e-6, - info = paste(i, mode) - ) - expect_equal(matched$SE, expected$SE, tolerance = 1e-6, info = paste(i, mode)) - } - - # the mean class is 1 + the sum of the exceedance probabilities - out <- as.data.frame(estimate_slopes( - m, - "x", - by = "g", - predict = "mean.class", - backend = "emmeans" - )) - out_exceedance <- as.data.frame(estimate_slopes( - m, - "x", - by = "g", - predict = "exc.prob", - backend = "emmeans" - )) - expect_identical(nrow(out), 3L, info = i) - expect_length(intersect(c("Response", "Threshold"), colnames(out)), 0) - sums <- tapply(out_exceedance$Slope, as.character(out_exceedance$g), sum) - expect_equal( - out$Slope, - as.vector(sums[as.character(out$g)]), - tolerance = 1e-6, - info = i - ) - } -}) - - -test_that("estimate_contrasts, backend emmeans, ordinal slopes", { - skip_if_not_installed("emmeans") - models <- .ordinal_slope_models() - - for (i in names(models)) { - m <- models[[i]] - slopes <- as.data.frame(estimate_slopes(m, "x", by = "g", backend = "emmeans")) - out <- as.data.frame(estimate_contrasts(m, "x", by = "g", backend = "emmeans")) - expect_true("Response" %in% colnames(out), info = i) - expect_identical(nrow(out), 9L, info = i) - - # each later level minus each earlier level ("revpairwise") - g_levels <- levels(slopes$g) - pairs <- utils::combn(seq_along(g_levels), 2) - expect_setequal( - paste(out$Level1, out$Level2), - paste(g_levels[pairs[2, ]], g_levels[pairs[1, ]]) - ) - - # differences of the slopes within each response category - slope_of <- function(g, response) { - slopes$Slope[ - as.character(slopes$g) == g & as.character(slopes$Response) == response - ] - } - expected <- mapply( - function(level1, level2, response) { - slope_of(level1, response) - slope_of(level2, response) - }, - as.character(out$Level1), - as.character(out$Level2), - as.character(out$Response) - ) - expect_equal(out$Difference, unname(expected), tolerance = 1e-6, info = i) - } -}) - - -test_that("estimate_slopes, backend emmeans, clmm", { - skip_if_not_installed("emmeans") - skip_if_not_installed("ordinal") - data(wine, package = "ordinal") - set.seed(3) - wine$x <- stats::rnorm(nrow(wine)) + as.numeric(wine$temp) - wine$rating_text <- factor(wine$rating, labels = letters[1:5], ordered = TRUE) - all_levels <- list(rating = as.character(1:5), rating_text = letters[1:5]) - - 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, "~ x * contact + (1 | judge)")) - m <- ordinal::clmm(f, data = wine) - - out <- as.data.frame(estimate_slopes(m, "x", by = "contact", backend = "emmeans")) - expect_identical(nrow(out), 10L, info = response) - expect_identical(levels(out$Response), response_levels, info = response) - - # slopes of the probabilities, compared with a central finite difference - # of the probabilities at the mean of `x` - h <- 1e-4 - probabilities <- function(x) { - as.data.frame(suppressMessages(emmeans::emmeans( - m, - stats::as.formula(paste("~ contact |", response)), - at = list(x = x), - mode = "prob" - ))) - } - upper <- probabilities(mean(wine$x) + h) - lower <- probabilities(mean(wine$x) - h) - finite_difference <- (upper$prob - lower$prob) / (2 * h) - # emmeans codes the response categories of clmm models as 1, 2, 3, ... - upper_levels <- response_levels[as.integer(as.character(upper[[response]]))] - matched <- out[ - match(paste(upper$contact, upper_levels), paste(out$contact, out$Response)), - ] - expect_false(anyNA(matched$Slope), info = response) - expect_equal(matched$Slope, finite_difference, tolerance = 1e-3, info = response) - - for (mode in c("cum.prob", "exc.prob", "linear.predictor")) { - out <- as.data.frame(estimate_slopes( - m, - "x", - by = "contact", - predict = mode, - backend = "emmeans" - )) - expect_identical(nrow(out), 8L, info = paste(response, mode)) - expect_identical( - levels(out$Threshold), - threshold_names, - info = paste(response, mode) - ) - } - - # contrasts of slopes keep the level names - out <- estimate_contrasts(m, "x", by = "contact", backend = "emmeans") - expect_identical(levels(out$Response), response_levels, info = response) - out <- estimate_contrasts( - m, - "x", - by = "contact", - predict = "cum.prob", - backend = "emmeans" - ) - expect_identical(levels(out$Threshold), threshold_names, info = response) - } -}) - - -test_that("estimate_slopes, backend emmeans, ordinal, predict = 'response' and no by", { - skip_if_not_installed("emmeans") - m <- .ordinal_slope_models()$simulated + m <- .ordinal_slope_model() - # predict = "response" returns probabilities, same as the default - out <- estimate_slopes(m, "x", by = "g", predict = "response", backend = "emmeans") - out_default <- estimate_slopes(m, "x", by = "g", backend = "emmeans") - expect_identical(paste(out$g, out$Response), paste(out_default$g, out_default$Response)) - expect_equal(out$Slope, out_default$Slope, tolerance = 1e-10) + # P(y = k) = P(y <= k) - P(y <= k - 1) + cumulative <- as.data.frame(suppressMessages(emmeans::emtrends( + m, + ~ cut | g, + var = "x", + mode = "cum.prob" + ))) + cumulative <- matrix(cumulative$x.trend, nrow = 2) + out <- estimate_slopes(m, "x", by = "g", backend = "emmeans") + expected <- rbind(cumulative[1, ], cumulative[2, ] - cumulative[1, ], -cumulative[2, ]) + expect_equal( + out$Slope, + expected[cbind(as.integer(out$Response), as.integer(out$g))], + tolerance = 1e-6 + ) - # no `by`: slopes of the probabilities, from the cumulative probabilities - out <- as.data.frame(estimate_slopes(m, "x", backend = "emmeans")) - expect_identical(as.character(out$Response), c("lo", "mid", "hi")) + out <- estimate_slopes(m, "x", backend = "emmeans") cumulative <- as.data.frame(suppressMessages(emmeans::emtrends( m, ~cut, @@ -947,50 +413,71 @@ test_that("estimate_slopes, backend emmeans, ordinal, predict = 'response' and n tolerance = 1e-6 ) - # no `by`: the mean class is 1 + the sum of the exceedance probabilities - out <- as.data.frame(estimate_slopes( + for (mode in c("latent", "cum.prob")) { + specs <- if (mode == "latent") ~g else ~ g | cut + expected <- as.data.frame(suppressMessages(emmeans::emtrends( + m, + specs, + var = "x", + mode = mode + ))) + out <- estimate_slopes(m, "x", by = "g", predict = mode, backend = "emmeans") + expect_equal(out$Slope, expected$x.trend, tolerance = 1e-6, info = mode) + } + + # the mean class is 1 + the sum of the exceedance probabilities + out <- estimate_slopes(m, "x", by = "g", predict = "mean.class", backend = "emmeans") + out_exceedance <- estimate_slopes( m, "x", - predict = "mean.class", + by = "g", + predict = "exc.prob", backend = "emmeans" - )) - out_exceedance <- estimate_slopes(m, "x", predict = "exc.prob", backend = "emmeans") - expect_identical(nrow(out), 1L) - expect_equal(out$Slope, sum(out_exceedance$Slope), tolerance = 1e-6) + ) + expect_equal( + out$Slope, + as.vector(tapply(out_exceedance$Slope, out_exceedance$g, sum)), + tolerance = 1e-6 + ) + out <- estimate_slopes(m, "x", predict = "mean.class", backend = "emmeans") + expect_equal(out$Slope, sum(out_exceedance$Slope) / 3, tolerance = 1e-6) }) -test_that("estimate_contrasts and estimate_slopes, backend emmeans, ordinal footer", { +test_that("estimate_contrasts, backend emmeans, ordinal slopes", { skip_if_not_installed("emmeans") - data(housing, package = "MASS") - m <- MASS::polr(Sat ~ Infl + Type + Cont, weights = Freq, data = housing, Hess = TRUE) - expect_output( - print(estimate_contrasts(m, "Infl", backend = "emmeans")), - "Contrasts are on the prob-scale.", - fixed = TRUE + m <- .ordinal_slope_model() + slopes <- estimate_slopes(m, "x", by = "g", backend = "emmeans") + out <- estimate_contrasts(m, "x", by = "g", backend = "emmeans") + slope_of <- function(g, response) { + slopes$Slope[slopes$g == g & slopes$Response == response] + } + expected <- mapply( + function(level1, level2, response) { + slope_of(level1, response) - slope_of(level2, response) + }, + as.character(out$Level1), + as.character(out$Level2), + as.character(out$Response) ) + expect_identical(nrow(out), 9L) + expect_equal(out$Difference, unname(expected), tolerance = 1e-6) +}) + + +test_that("estimate_contrasts and estimate_slopes, backend emmeans, ordinal footer", { + skip_if_not_installed("emmeans") + m <- .ordinal_slope_model() expect_output( - print(estimate_contrasts(m, "Infl", predict = "response", backend = "emmeans")), + print(estimate_contrasts(m, "g", backend = "emmeans")), "Contrasts are on the prob-scale.", fixed = TRUE ) - expect_output( - print(estimate_contrasts(m, "Infl", predict = "cum.prob", backend = "emmeans")), - "Contrasts are on the cum.prob-scale.", - fixed = TRUE - ) - - m <- .ordinal_slope_models()$housing expect_output( print(estimate_slopes(m, "x", by = "g", backend = "emmeans")), "Slopes are on the prob-scale.", fixed = TRUE ) - expect_output( - print(estimate_slopes(m, "x", by = "g", predict = "latent", backend = "emmeans")), - "Slopes are on the latent-scale.", - fixed = TRUE - ) # no scale line for the slopes of other models m <- lm(mpg ~ wt * factor(am), data = mtcars) @@ -1031,227 +518,124 @@ test_that("estimate_contrasts, backend emmeans, ordinal, names of response and p housing$x <- as.numeric(housing$Cont) + stats::rnorm(nrow(housing)) sat_levels <- c("Low", "Medium", "High") - # responses named like emmeans output columns - for (response in c("response", "prob", "rate", "ratio")) { - d <- housing - d[[response]] <- d$Sat - f <- stats::as.formula(paste(response, "~ Infl + Type + Cont")) - m <- MASS::polr(f, weights = Freq, data = d, Hess = TRUE) - out <- estimate_contrasts(m, "Infl", backend = "emmeans") - expect_identical(levels(out$Response), sat_levels, info = response) - expect_type(out$Difference, "double") - - f <- stats::as.formula(paste(response, "~ Infl * x")) - m <- MASS::polr(f, weights = Freq, data = d, Hess = TRUE) - out <- estimate_contrasts(m, "x", by = "Infl", backend = "emmeans") - expect_identical(levels(out$Response), sat_levels, info = response) - expect_type(out$Difference, "double") - } - - # a predictor named `Response` keeps its column, and the response - # categories get a suffix, as for `estimate_means()` - d <- housing - d$Response <- d$Infl - m <- MASS::polr(Sat ~ Response + Type + Cont, weights = Freq, data = d, Hess = TRUE) - out <- estimate_means(m, c("Type", "Response"), backend = "emmeans") - expect_identical(colnames(out)[1:3], c("Type", "Response", "Response.1")) - out <- estimate_contrasts(m, "Type", by = "Response", backend = "emmeans") - expect_identical(nrow(out), 54L) - expect_true(all(c("Response", "Response.1") %in% colnames(out))) - expect_length(unique(paste(out$Response, out$Response.1)), 9) - expect_identical(levels(out$Response.1), sat_levels) -}) - - -test_that("estimate_slopes, backend emmeans, ordinal glmmTMB", { - skip_if_not_installed("emmeans") - # ordinal family in glmmTMB requires glmmTMB >= 1.1.15.2 - skip_if_not_installed("glmmTMB", minimum_version = "1.1.15.2") - data(housing, package = "MASS") - set.seed(1) - housing$x <- as.numeric(housing$Cont) + stats::rnorm(nrow(housing)) - # the same model, fitted with polr and glmmTMB - m_polr <- MASS::polr(Sat ~ Infl * x, weights = Freq, data = housing, Hess = TRUE) - m_glmmTMB <- glmmTMB::glmmTMB( - Sat ~ Infl * x, + # a response named like an emmeans output column + housing$prob <- housing$Sat + m <- MASS::polr(prob ~ Infl + Type + Cont, weights = Freq, data = housing, Hess = TRUE) + out <- estimate_contrasts(m, "Infl", backend = "emmeans") + expect_identical(levels(out$Response), sat_levels) + m <- MASS::polr(prob ~ Infl * x, weights = Freq, data = housing, Hess = TRUE) + out <- estimate_contrasts(m, "x", by = "Infl", backend = "emmeans") + expect_identical(levels(out$Response), sat_levels) + + # a predictor named `Response` keeps its column + housing$Response <- housing$Infl + m <- MASS::polr( + Sat ~ Response + Type + Cont, weights = Freq, data = housing, - family = glmmTMB::ordinal() + Hess = TRUE ) - out_polr <- estimate_slopes(m_polr, "x", by = "Infl", backend = "emmeans") - out <- estimate_slopes(m_glmmTMB, "x", by = "Infl", backend = "emmeans") - expect_identical(nrow(out), 9L) - expect_identical(paste(out$Infl, out$Response), paste(out_polr$Infl, out_polr$Response)) - expect_equal(out$Slope, out_polr$Slope, tolerance = 1e-4) + out <- estimate_contrasts(m, "Type", by = "Response", backend = "emmeans") + expect_identical(nrow(out), 54L) + expect_identical(levels(out$Response.1), sat_levels) }) # `by` filters on the response categories, emmeans backend -------------------- -# the output of a call with a `by` filter on the response categories must be -# the rows of the same call without that filter, for the chosen categories -expect_response_filter_rows <- function(out, out_all, response_levels, info = NULL) { - expected <- out_all[as.character(out_all$Response) %in% response_levels, , drop = FALSE] - expect_identical(nrow(out), nrow(expected), info = info) - expect_gt(nrow(out), 0) - expect_identical(colnames(out), colnames(out_all), info = info) - for (column in colnames(expected)) { - if (is.numeric(expected[[column]])) { - expect_equal( - out[[column]], - expected[[column]], - tolerance = 1e-10, - info = paste(info, column) - ) - } else { - expect_identical( - as.character(out[[column]]), - as.character(expected[[column]]), - info = paste(info, column) - ) - } - } -} - - -ordinal_filter_models <- function() { +ordinal_filter_slope_model <- function() { data(housing, package = "MASS", envir = environment()) - data(wine, package = "ordinal", envir = environment()) - wine_letters <- wine - levels(wine_letters$rating) <- letters[1:5] - list( - polr = list( - model = MASS::polr(Sat ~ Infl + Type + Cont, weights = Freq, data = housing, Hess = TRUE), - factor = "Infl", - response = "Sat", - one = "Low", - two = c("Low", "High") - ), - clm = list( - model = ordinal::clm(Sat ~ Infl + Type + Cont, weights = Freq, data = housing), - factor = "Infl", - response = "Sat", - one = "Medium", - two = c("Low", "High") - ), - clmm = list( - model = ordinal::clmm(rating ~ temp + contact + (1 | judge), data = wine), - factor = "temp", - response = "rating", - one = "2", - two = c("1", "3") - ), - clmm_letters = list( - model = ordinal::clmm(rating ~ temp + contact + (1 | judge), data = wine_letters), - factor = "temp", - response = "rating", - one = "b", - two = c("a", "c") - ) - ) -} - - -response_filter_term <- function(response, response_levels) { - if (length(response_levels) == 1) { - paste0(response, "='", response_levels, "'") - } else { - paste0(response, "=c(", toString(paste0("'", response_levels, "'")), ")") - } + set.seed(1) + housing$x <- stats::rnorm(nrow(housing)) + MASS::polr(Sat ~ x * Infl + Type, data = housing, Hess = TRUE) } -test_that("estimate_means, backend emmeans, ordinal, by filter on the response", { +test_that("estimate_means and estimate_contrasts, backend emmeans, ordinal, by filter on the response", { skip_if_not_installed("emmeans") skip_if_not_installed("ordinal") - models <- ordinal_filter_models() - for (name in names(models)) { - fit <- models[[name]] - out_all <- estimate_means(fit$model, fit$factor, backend = "emmeans") - for (response_levels in list(fit$one, fit$two)) { - term <- response_filter_term(fit$response, response_levels) - out <- estimate_means(fit$model, c(fit$factor, term), backend = "emmeans") - expect_response_filter_rows(out, out_all, response_levels, info = paste(name, term)) - } - } -}) + data(housing, package = "MASS") + m <- MASS::polr(Sat ~ Infl + Type + Cont, weights = Freq, data = housing, Hess = TRUE) + out_all <- as.data.frame(estimate_means(m, "Infl", backend = "emmeans")) + out <- estimate_means(m, c("Infl", "Sat='Low'"), backend = "emmeans") + expect_equal( + as.data.frame(out), + out_all[out_all$Response == "Low", ], + ignore_attr = TRUE + ) + # the response without values keeps all categories, and does not undo a filter + out <- estimate_means(m, c("Infl", "Sat"), backend = "emmeans") + expect_equal(as.data.frame(out), out_all, ignore_attr = TRUE) + out <- estimate_means( + m, + c("Infl", "Sat", "Sat='Low'", "Sat='High'"), + backend = "emmeans" + ) + expect_equal( + as.data.frame(out), + out_all[out_all$Response != "Medium", ], + ignore_attr = TRUE + ) -test_that("estimate_contrasts, backend emmeans, ordinal, by filter on the response", { - skip_if_not_installed("emmeans") - skip_if_not_installed("ordinal") - models <- ordinal_filter_models() - for (name in names(models)) { - fit <- models[[name]] - p_adjusts <- if (name %in% c("polr", "clm")) c("none", "holm") else "none" - for (p_adjust in p_adjusts) { - out_all <- estimate_contrasts( - fit$model, - fit$factor, - backend = "emmeans", - p_adjust = p_adjust - ) - for (response_levels in list(fit$one, fit$two)) { - term <- response_filter_term(fit$response, response_levels) - out <- estimate_contrasts( - fit$model, - fit$factor, - by = term, - backend = "emmeans", - p_adjust = p_adjust - ) - expect_response_filter_rows( - out, - out_all, - response_levels, - info = paste(name, term, p_adjust) - ) - } - } - } + out_all <- as.data.frame(estimate_contrasts( + m, + "Infl", + p_adjust = "holm", + backend = "emmeans" + )) + out <- estimate_contrasts( + m, + "Infl", + by = "Sat=c('Low','High')", + p_adjust = "holm", + backend = "emmeans" + ) + expect_equal( + as.data.frame(out), + out_all[out_all$Response != "Medium", ], + ignore_attr = TRUE + ) + + data(wine, package = "ordinal") + levels(wine$rating) <- letters[1:5] + m <- ordinal::clmm(rating ~ temp + contact + (1 | judge), data = wine) + out_all <- as.data.frame(estimate_means(m, "temp", backend = "emmeans")) + out <- estimate_means(m, c("temp", "rating='b'"), backend = "emmeans") + expect_equal(as.data.frame(out), out_all[out_all$Response == "b", ], ignore_attr = TRUE) }) -test_that("estimate_means, backend emmeans, ordinal, response named in by without values", { +test_that("estimate_slopes, backend emmeans, ordinal, by filter on the response", { skip_if_not_installed("emmeans") - data(housing, package = "MASS") - m <- MASS::polr(Sat ~ Infl + Type + Cont, weights = Freq, data = housing) - out_all <- estimate_means(m, "Infl", backend = "emmeans") - out <- estimate_means(m, c("Infl", "Sat"), backend = "emmeans") - expect_response_filter_rows(out, out_all, levels(housing$Sat)) - expect_identical(nrow(out), 9L) -}) - - -test_that("estimate_slopes, backend emmeans, ordinal, by filter on the response", { - skip_if_not_installed("emmeans") - data(housing, package = "MASS") - set.seed(1) - housing$x <- as.numeric(housing$Cont) + stats::rnorm(nrow(housing)) - m <- MASS::polr(Sat ~ Infl * x, weights = Freq, data = housing, Hess = TRUE) + m <- ordinal_filter_slope_model() - out_all <- estimate_slopes(m, "x", by = "Infl", backend = "emmeans") + out_all <- as.data.frame(estimate_slopes(m, "x", by = "Infl", backend = "emmeans")) out <- estimate_slopes(m, "x", by = c("Infl", "Sat='Low'"), backend = "emmeans") - expect_response_filter_rows(out, out_all, "Low") - expect_identical(nrow(out), 3L) + expect_equal( + as.data.frame(out), + out_all[out_all$Response == "Low", ], + ignore_attr = TRUE + ) - # contrasts of slopes - out_all <- estimate_contrasts(m, "x", by = "Infl", backend = "emmeans") + out_all <- as.data.frame(estimate_contrasts(m, "x", by = "Infl", backend = "emmeans")) out <- estimate_contrasts(m, "x", by = c("Infl", "Sat='Low'"), backend = "emmeans") - expect_response_filter_rows(out, out_all, "Low") - expect_identical(nrow(out), 3L) + expect_equal( + as.data.frame(out), + out_all[out_all$Response == "Low", ], + ignore_attr = TRUE + ) # slopes use `"cum.prob"` internally, but the mode check uses `predict` - err <- expect_error(estimate_slopes( - m, - "x", - by = c("Infl", "Sat='Low'"), - predict = "cum.prob", - backend = "emmeans" - )) - expect_match( - gsub("\\s+", " ", conditionMessage(err)), - "The response `Sat` can only be used in `by` with `predict = \"prob\"`", + expect_error( + estimate_slopes( + m, + "x", + by = c("Infl", "Sat='Low'"), + predict = "cum.prob", + backend = "emmeans" + ), + "The response `Sat` can only be used in `by`", fixed = TRUE ) }) @@ -1259,141 +643,93 @@ test_that("estimate_slopes, backend emmeans, ordinal, by filter on the response" test_that("estimate_means and estimate_contrasts, backend emmeans, ordinal, errors for response filters", { skip_if_not_installed("emmeans") + skip_if_not_installed("ordinal") data(housing, package = "MASS") m <- MASS::polr(Sat ~ Infl + Type + Cont, weights = Freq, data = housing, Hess = TRUE) - # the messages are wrapped at spaces, so we compare them on one line - error_message <- function(expr) { - err <- expect_error(expr) - gsub("\\s+", " ", conditionMessage(err)) - } - - # the response categories can only be filtered for `predict = "prob"` - modes <- c("cum.prob", "exc.prob", "linear.predictor", "latent", "mean.class") - for (mode in modes) { - messages <- c( - error_message(estimate_means( - m, - c("Infl", "Sat='Low'"), - predict = mode, - backend = "emmeans" - )), - error_message(estimate_contrasts( - m, - "Infl", - by = "Sat='Low'", - predict = mode, - backend = "emmeans" - )) - ) - for (msg in messages) { - expect_match(msg, "The response `Sat` can only be used in `by`", fixed = TRUE, info = mode) - expect_match(msg, "`predict = \"prob\"`", fixed = TRUE, info = mode) - } - } - # values that are no response categories - messages <- c( - error_message(estimate_means(m, c("Infl", "Sat='Lowest'"), backend = "emmeans")), - error_message(estimate_contrasts(m, "Infl", by = "Sat='Lowest'", backend = "emmeans")) + expect_error( + estimate_means(m, c("Infl", "Sat='Low'"), predict = "cum.prob", backend = "emmeans"), + "The response `Sat` can only be used in `by`", + fixed = TRUE + ) + expect_error( + estimate_contrasts( + m, + "Infl", + by = "Sat='Low'", + predict = "latent", + backend = "emmeans" + ), + "The response `Sat` can only be used in `by`", + fixed = TRUE + ) + expect_error( + estimate_means(m, c("Infl", "Sat='Lowest'"), backend = "emmeans"), + "The response `Sat` has no category `Lowest`.", + fixed = TRUE + ) + expect_error( + estimate_means(m, c("Infl", "Sat=c('Lowest','Highest','Low')"), backend = "emmeans"), + "has no categories `Lowest` and `Highest`.", + fixed = TRUE + ) + expect_error( + estimate_means(m, c("Infl", "Sat=[fivenum]"), backend = "emmeans"), + "Could not read the response categories in `Sat=[fivenum]`.", + fixed = TRUE ) - for (msg in messages) { - expect_match(msg, "The response `Sat` has no category `Lowest`.", fixed = TRUE) - } - msg <- error_message(estimate_means( - m, - c("Infl", "Sat=c('Lowest','Highest','Low')"), - backend = "emmeans" - )) - expect_match(msg, "has no categories `Lowest` and `Highest`.", fixed = TRUE) - - # values that cannot be read as response categories - msg <- error_message(estimate_means(m, c("Infl", "Sat=[fivenum]"), backend = "emmeans")) - expect_match(msg, "Could not read the response categories in `Sat=[fivenum]`.", fixed = TRUE) # a category of the data that the model does not have - skip_if_not_installed("ordinal") data(wine, package = "ordinal") wine$rating <- factor(wine$rating, levels = 1:6, ordered = TRUE) - m_clm <- ordinal::clm(rating ~ temp, data = wine) - msg <- error_message(estimate_means(m_clm, c("temp", "rating='6'"), backend = "emmeans")) - expect_match(msg, "The response `rating` has no category `6`.", fixed = TRUE) -}) - - -test_that("estimate_means, backend emmeans, ordinal, several by terms on the response", { - skip_if_not_installed("emmeans") - data(housing, package = "MASS") - m <- MASS::polr(Sat ~ Infl + Type + Cont, weights = Freq, data = housing) - out_all <- estimate_means(m, "Infl", backend = "emmeans") - - # the categories of all terms are kept - out <- estimate_means(m, c("Infl", "Sat='Low'", "Sat='High'"), backend = "emmeans") - expect_response_filter_rows(out, out_all, c("Low", "High")) - expect_identical(nrow(out), 6L) - - # the response without values does not undo the filter - out <- estimate_means(m, c("Infl", "Sat", "Sat='Low'"), backend = "emmeans") - expect_response_filter_rows(out, out_all, "Low") - expect_identical(nrow(out), 3L) + m <- ordinal::clm(rating ~ temp, data = wine) + expect_error( + estimate_means(m, c("temp", "rating='6'"), backend = "emmeans"), + "The response `rating` has no category `6`.", + fixed = TRUE + ) }) 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) - ) + m <- MASS::polr(Sat ~ x * Infl + Type, data = d, Hess = TRUE) 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) - } + # one level: the contrasts between the response categories at that level + for (estimate in c("typical", "average")) { + out_all <- as.data.frame(estimate_contrasts(m, "x", by = "Infl", estimate = estimate)) + out <- as.data.frame(estimate_contrasts( + m, + "x", + by = "Infl='Low'", + estimate = estimate + )) + expected <- out_all[out_all$Level1 == "Low" & out_all$Level2 == "Low", ] + expect_identical( + lapply(out[keys], as.character), + lapply(expected[keys], as.character) + ) + expect_equal(out$Difference, expected$Difference, 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)) + m2 <- MASS::polr(Sat ~ x * Infl + Type, data = d2, Hess = TRUE) + out <- estimate_contrasts(m2, "x", by = "Infl='Low'") 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"), + estimate_contrasts(m, "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), + estimate_contrasts(m, "x", by = "Infl='Low'", comparison = ~ pairwise | Infl), "No contrasts to show", fixed = TRUE ) @@ -1408,362 +744,87 @@ test_that("estimate_contrasts, ordinal slopes, by filter that keeps one level", }) -# 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 -# cumulative probabilities, averaged over the levels of `Type` with equal -# weights, as emmeans does. `contrast_fun` maps the 3 x 3 matrix of -# probabilities (rows `Infl`, columns `Sat`) to a vector of contrasts `e`. The -# Jacobian `J` of `e` comes from central finite differences of the parameters, -# so the covariance of `e` is `J V J'`, and the Wald statistic is -# `e' (J V J')^-1 e`. -ordinal_wald_chisq <- function(model, contrast_fun) { - grid <- expand.grid( - Infl = levels(MASS::housing$Infl), - Type = levels(MASS::housing$Type) - ) - design <- stats::model.matrix(~ Infl + Type, grid)[, -1] - n_beta <- length(stats::coef(model)) - probabilities <- function(theta) { - eta <- drop(design %*% theta[seq_len(n_beta)]) - zeta <- theta[-seq_len(n_beta)] - cum_low <- stats::plogis(zeta[1] - eta) - cum_medium <- stats::plogis(zeta[2] - eta) - p <- cbind(Low = cum_low, Medium = cum_medium - cum_low, High = 1 - cum_medium) - out <- t(sapply(levels(MASS::housing$Infl), function(i) { - colMeans(p[grid$Infl == i, , drop = FALSE]) - })) - out - } - theta <- c(stats::coef(model), model$zeta) - e <- contrast_fun(probabilities(theta)) - h <- 1e-6 - J <- sapply(seq_along(theta), function(j) { - up <- down <- theta - up[j] <- up[j] + h - down[j] <- down[j] - h - (contrast_fun(probabilities(up)) - contrast_fun(probabilities(down))) / (2 * h) - }) - V <- J %*% stats::vcov(model) %*% t(J) - drop(t(e) %*% solve(V) %*% e) -} - - test_that("estimate_contrasts, backend emmeans, ordinal, joint tests of the response categories", { 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) - - # reference 1: delta method. Within each `Infl` level, the contrasts are - # Medium - Low and High - Low of the category probabilities - ref_delta <- vapply( - levels(housing$Infl), - function(i) { - ordinal_wald_chisq(m, function(p) { - c(p[i, "Medium"] - p[i, "Low"], p[i, "High"] - p[i, "Low"]) - }) - }, - numeric(1) - ) - - # reference 2: marginaleffects, with the same contrasts of the 36 predicted - # probabilities (4 `Type` x 3 `Infl` x 3 `Sat`). Each column of `H` is a - # contrast at one `Infl` level, with weight 1/4 on the probabilities of - # `Medium` (or `High`) and -1/4 on those of `Low` - nd <- marginaleffects::datagrid( - model = m, - Infl = levels(housing$Infl), - Type = levels(housing$Type) - ) - pr <- marginaleffects::predictions(m, newdata = nd, type = "probs") - ref_me <- vapply( - levels(housing$Infl), - function(i) { - w <- (pr$Infl == i) / 4 - H <- cbind( - w * (pr$group == "Medium") - w * (pr$group == "Low"), - w * (pr$group == "High") - w * (pr$group == "Low") - ) - pc <- marginaleffects::predictions(m, newdata = nd, type = "probs", hypothesis = H) - marginaleffects::hypotheses(pc, joint = TRUE, joint_test = "chisq")$statistic - }, - numeric(1) - ) - expect_equal(ref_delta, ref_me, tolerance = 1e-3, ignore_attr = TRUE) - - for (model in list(m, mc)) { - expect_no_warning({ - out <- estimate_contrasts( - model, - "Sat", - by = "Infl", - comparison = "joint", - backend = "emmeans" - ) - }) - expect_identical(nrow(out), 3L) - expect_identical(as.character(out$Infl), levels(housing$Infl)) - expect_identical(as.character(out$df1), rep("2", 3)) - expect_equal(out$Chisq, ref_delta, tolerance = 1e-3, ignore_attr = TRUE) - expect_equal(out$Chisq, ref_me, tolerance = 1e-3, ignore_attr = TRUE) - } -}) - - -test_that("estimate_contrasts, backend emmeans, ordinal, joint tests of a predictor", { - 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) - - # reference 1: delta method. Within each response category, the contrasts - # are Medium - Low and High - Low of the `Infl` levels - ref_delta <- vapply( - levels(housing$Sat), - function(k) { - ordinal_wald_chisq(m, function(p) { - c(p["Medium", k] - p["Low", k], p["High", k] - p["Low", k]) - }) - }, - numeric(1) - ) - - for (model in list(m, mc)) { - # reference 2: emmeans' joint tests of the pairwise contrasts. For the clm - # fit, emmeans warns with "NaNs produced" when it computes a p-value with - # `pf()`, which does not change these values - ref_emmeans <- suppressWarnings(as.data.frame(emmeans::joint_tests( - emmeans::contrast( - emmeans::emmeans(model, ~ Infl | Sat, mode = "prob"), - "pairwise" - ), - by = "Sat" - ))) - expect_no_warning({ - out <- estimate_contrasts(model, "Infl", comparison = "joint", backend = "emmeans") - }) - expect_identical(nrow(out), 3L) - expect_identical(as.character(out$Response), levels(housing$Sat)) - expect_equal(out$Chisq, ref_delta, tolerance = 1e-3, ignore_attr = TRUE) - expect_equal( - out$Chisq, - ref_emmeans$Chisq[match(levels(housing$Sat), ref_emmeans$Sat)], - tolerance = 1e-3, - ignore_attr = TRUE - ) - } - - # a joint test with another `by` variable gives no warning expect_no_warning({ out <- estimate_contrasts( m, - "Infl", - by = "Type", + "Sat", + by = "Infl", comparison = "joint", backend = "emmeans" ) }) - expect_identical(nrow(out), 12L) + # Wald statistics from `marginaleffects::hypotheses(joint = TRUE)` + expect_equal(out$Chisq, c(45.096, 29.214, 90.344), tolerance = 1e-3) }) -test_that("estimate_contrasts, backend emmeans, ordinal, printed Chisq of joint tests", { +test_that("estimate_contrasts, backend emmeans, ordinal, joint tests of a predictor", { skip_if_not_installed("emmeans") data(housing, package = "MASS") m <- MASS::polr(Sat ~ Infl + Type, weights = Freq, data = housing, Hess = TRUE) - out <- estimate_contrasts(m, "Infl", comparison = "joint", backend = "emmeans") + expected <- as.data.frame(emmeans::joint_tests( + emmeans::contrast(emmeans::emmeans(m, ~ Infl | Sat, mode = "prob"), "pairwise"), + by = "Sat" + )) + expect_no_warning({ + out <- estimate_contrasts(m, "Infl", comparison = "joint", backend = "emmeans") + }) + expect_identical(as.character(out$Response), as.character(expected$Sat)) + expect_equal(out$Chisq, expected$Chisq, tolerance = 1e-3) + + # printed with two decimals, as the F column printed <- capture.output(print(out)) - # the Chisq values are 114.588, 20.014 and 104.514, printed with two - # decimals as the F column - for (chisq in c("114.59", "20.01", "104.51")) { - expect_true(any(grepl(paste0(" ", chisq, " |"), printed, fixed = TRUE)), info = chisq) - } - expect_false(any(grepl("114.588", printed, fixed = TRUE))) + expect_true(any(grepl(" 114.59 |", printed, fixed = TRUE))) + + out <- get_emcontrasts(m, "Infl", comparison = "joint") + expect_identical(class(out), "data.frame") + expect_null(attr(out, "digits")) }) -test_that("estimate_contrasts, backend emmeans, ordinal, joint tests in other modes", { +# `by` as list, data frame or formula, with the response ----------------------- + +test_that("estimate_means, backend emmeans, ordinal, list and data frame filter on the response", { 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) - # emmeans' joint tests of the pairwise `Infl` contrasts, within each - # threshold (`cut`) where the mode has thresholds - specs <- list( - latent = ~Infl, - linear.predictor = ~ Infl | cut, - cum.prob = ~ Infl | cut, - exc.prob = ~ Infl | cut, - mean.class = ~Infl - ) - for (model in list(m, mc)) { - for (mode in names(specs)) { - by <- if (mode %in% c("latent", "mean.class")) NULL else "cut" - ref <- suppressWarnings(as.data.frame(emmeans::joint_tests( - emmeans::contrast( - emmeans::emmeans(model, specs[[mode]], mode = mode), - "pairwise" - ), - by = by - ))) - # for the clm fit, emmeans warns with "NaNs produced" in some of these - # modes, also on the reference above; other warnings are not muffled - out <- withCallingHandlers( - estimate_contrasts( - model, - "Infl", - comparison = "joint", - predict = mode, - backend = "emmeans" - ), - warning = function(w) { - if (identical(conditionMessage(w), "NaNs produced")) { - invokeRestart("muffleWarning") - } - } - ) - expect_identical(nrow(out), nrow(ref), info = mode) - expect_equal(out$Chisq, ref$Chisq, tolerance = 1e-6, ignore_attr = TRUE, info = mode) - } - } - # on the logit scale, the `Infl` effect is the same at each threshold, so - # both thresholds give the latent-scale value - out <- estimate_contrasts( + out_all <- estimate_means(m, "Infl", backend = "emmeans") + + out <- estimate_means( m, - "Infl", - comparison = "joint", - predict = "linear.predictor", + list(Infl = c("Low", "High"), Sat = "Low"), backend = "emmeans" ) - out_latent <- estimate_contrasts( + expect_identical(nrow(out), 2L) + expect_equal(out$Probability, out_all$Probability[c(1, 3)], tolerance = 1e-10) + out <- estimate_means(m, data.frame(Infl = "Low", Sat = "Low"), backend = "emmeans") + expect_equal(out$Probability, out_all$Probability[1], tolerance = 1e-10) + + out <- estimate_contrasts(m, "Infl", by = "Sat", backend = "emmeans") + out_list <- estimate_contrasts( m, "Infl", - comparison = "joint", - predict = "latent", + by = list(Sat = c("Low", "High")), backend = "emmeans" ) - expect_equal(out$Chisq, rep(out_latent$Chisq, 2), tolerance = 1e-6) -}) - - -test_that("get_emcontrasts, joint tests return a data frame", { - skip_if_not_installed("emmeans") - data(housing, package = "MASS") - m <- MASS::polr(Sat ~ Infl + Type, weights = Freq, data = housing, Hess = TRUE) - out <- get_emcontrasts(m, "Infl", comparison = "joint") - expect_identical(class(out), "data.frame") - expect_null(attr(out, "digits")) - expect_no_error(capture.output(print(out))) -}) - - -# `by` as list, data frame or formula, with the response ----------------------- - -ordinal_filter_form_models <- function() { - data(housing, package = "MASS", envir = environment()) - list( - polr = MASS::polr(Sat ~ Infl + Type, weights = Freq, data = housing, Hess = TRUE), - clm = ordinal::clm(Sat ~ Infl + Type, weights = Freq, data = housing) + expect_equal( + out_list$Difference, + out$Difference[out$Response != "Medium"], + tolerance = 1e-10 ) -} - - -ordinal_filter_slope_model <- function() { - data(housing, package = "MASS", envir = environment()) - set.seed(1) - housing$x <- stats::rnorm(nrow(housing)) - MASS::polr(Sat ~ x * Infl + Type, data = housing, Hess = TRUE) -} - - -# the rows of `out` must be the rows of `out_all` with the same values in the -# `keys` columns, with equal numeric columns -expect_matching_rows <- function(out, out_all, keys, columns, info = NULL) { - key <- function(d) do.call(paste, c(lapply(d[keys], as.character), sep = "\r")) - expect_gt(nrow(out), 0) - matched <- match(key(out), key(out_all)) - expect_false(anyNA(matched), info = info) - for (column in columns) { - expect_equal( - out[[column]], - out_all[[column]][matched], - tolerance = 1e-10, - info = paste(info, column) - ) - } -} - - -test_that("estimate_means, backend emmeans, ordinal, list and data frame filter on the response", { - skip_if_not_installed("emmeans") - skip_if_not_installed("ordinal") - models <- ordinal_filter_form_models() - for (name in names(models)) { - m <- models[[name]] - out_all <- estimate_means(m, "Infl", backend = "emmeans") - - out <- estimate_means( - m, - list(Infl = c("Low", "High"), Sat = "Low"), - backend = "emmeans" - ) - expect_identical(nrow(out), 2L, info = name) - expect_identical(as.character(out$Response), c("Low", "Low"), info = name) - expect_identical(as.character(out$Infl), c("Low", "High"), info = name) - expect_matching_rows(out, out_all, c("Infl", "Response"), c("Probability", "SE"), info = name) - - out <- estimate_means(m, data.frame(Infl = "Low", Sat = "Low"), backend = "emmeans") - expect_identical(nrow(out), 1L, info = name) - expect_identical(as.character(out$Response), "Low", info = name) - expect_identical(as.character(out$Infl), "Low", info = name) - expect_matching_rows(out, out_all, c("Infl", "Response"), c("Probability", "SE"), info = name) - } -}) - - -test_that("estimate_contrasts and estimate_slopes, backend emmeans, ordinal, list and data frame filter on the response", { - skip_if_not_installed("emmeans") - skip_if_not_installed("ordinal") - models <- ordinal_filter_form_models() - for (name in names(models)) { - m <- models[[name]] - out_all <- estimate_contrasts(m, "Infl", by = "Sat", backend = "emmeans") - by_forms <- list( - list = list(Sat = c("Low", "High")), - data_frame = data.frame(Sat = c("Low", "High")) - ) - for (form in names(by_forms)) { - out <- estimate_contrasts(m, "Infl", by = by_forms[[form]], backend = "emmeans") - expect_response_filter_rows(out, out_all, c("Low", "High"), info = paste(name, form)) - expect_identical(nrow(out), 6L, info = paste(name, form)) - } - } m2 <- ordinal_filter_slope_model() - out_character <- estimate_slopes(m2, "x", by = "Sat='Low'", backend = "emmeans") - out <- estimate_slopes(m2, "x", by = list(Sat = "Low"), backend = "emmeans") - expect_identical(nrow(out), 1L) - expect_identical(as.character(out$Response), "Low") - expect_equal(out$Slope, out_character$Slope, tolerance = 1e-10) - expect_equal(out$SE, out_character$SE, tolerance = 1e-10) -}) - + out <- estimate_slopes(m2, "x", by = "Sat='Low'", backend = "emmeans") + out_list <- estimate_slopes(m2, "x", by = list(Sat = "Low"), backend = "emmeans") + expect_equal(out_list$Slope, out$Slope, tolerance = 1e-10) -test_that("estimate_means, backend emmeans, ordinal, footer for by forms with the response", { - skip_if_not_installed("emmeans") - data(housing, package = "MASS") - m <- MASS::polr(Sat ~ Infl + Type, weights = Freq, data = housing, Hess = TRUE) footer_lines <- function(out) { grep("^Predictors modulated", insight::trim_ws(utils::capture.output(print(out))), value = TRUE) } - expect_identical( - footer_lines(estimate_means(m, c("Infl", "Sat"), backend = "emmeans")), - "Predictors modulated: Infl" - ) by_forms <- list( formula = ~ Infl + Sat, list = list(Infl = c("Low", "High"), Sat = "Low"), @@ -1773,133 +834,85 @@ test_that("estimate_means, backend emmeans, ordinal, footer for by forms with th out <- estimate_means(m, by_forms[[form]], backend = "emmeans") expect_identical(footer_lines(out), "Predictors modulated: Infl", info = form) } - # the formula has no values, so all response categories are kept - out <- estimate_means(m, ~ Infl + Sat, backend = "emmeans") - expect_identical(nrow(out), 9L) -}) - -test_that("estimate_means and estimate_slopes, backend emmeans, ordinal, errors for list and data frame filters", { - skip_if_not_installed("emmeans") - skip_if_not_installed("ordinal") - error_message <- function(expr) { - err <- expect_error(expr) - gsub("\\s+", " ", conditionMessage(err)) - } - by_forms <- function(value) { - list( - list = list(Infl = "Low", Sat = value), - data_frame = data.frame(Infl = "Low", Sat = value) - ) - } - models <- ordinal_filter_form_models() - m2 <- ordinal_filter_slope_model() - calls <- list( - polr = function(by, ...) estimate_means(models$polr, by, backend = "emmeans", ...), - clm = function(by, ...) estimate_means(models$clm, by, backend = "emmeans", ...), - slopes = function(by, ...) estimate_slopes(m2, "x", by = by, backend = "emmeans", ...) + expect_error( + estimate_means(m, list(Infl = "Low", Sat = "Bogus"), backend = "emmeans"), + "The response `Sat` has no category `Bogus`.", + fixed = TRUE + ) + expect_error( + estimate_slopes( + m2, + "x", + by = data.frame(Sat = "Low"), + predict = "cum.prob", + backend = "emmeans" + ), + "The response `Sat` can only be used in `by`", + fixed = TRUE ) - for (call in names(calls)) { - for (form in c("list", "data_frame")) { - info <- paste(call, form) - msg <- error_message(calls[[call]](by_forms("Bogus")[[form]])) - expect_match(msg, "The response `Sat` has no category `Bogus`.", fixed = TRUE, info = info) - msg <- error_message(calls[[call]](by_forms("Low")[[form]], predict = "cum.prob")) - expect_match(msg, "The response `Sat` can only be used in `by`", fixed = TRUE, info = info) - } - } }) test_that("estimate_contrasts, backend emmeans, ordinal, contrast filter on the response", { skip_if_not_installed("emmeans") - skip_if_not_installed("ordinal") - models <- ordinal_filter_form_models() - for (name in names(models)) { - m <- models[[name]] - out_all <- estimate_contrasts(m, "Sat", by = "Infl", backend = "emmeans") - - # only the contrast of the two named categories, in each level of `Infl` - out <- estimate_contrasts(m, "Sat=c('Low','High')", by = "Infl", backend = "emmeans") - expect_identical(nrow(out), 3L, info = name) - expect_identical(as.character(out$Level1), rep("High", 3), info = name) - expect_identical(as.character(out$Level2), rep("Low", 3), info = name) - expect_identical(as.character(out$Infl), c("Low", "Medium", "High"), info = name) - expect_matching_rows( - out, - out_all, - c("Level1", "Level2", "Infl"), - c("Difference", "SE"), - info = name - ) + data(housing, package = "MASS") + m <- MASS::polr(Sat ~ Infl + Type, weights = Freq, data = housing, Hess = TRUE) + out_all <- estimate_contrasts(m, "Sat", by = "Infl", backend = "emmeans") - # comparisons that are not pairwise use only the named categories - out <- estimate_contrasts( - m, - "Sat=c('Low','Medium')", - by = "Infl", - comparison = "poly", - backend = "emmeans" - ) - expect_identical(nrow(out), 3L, info = name) - expect_identical(as.character(out$Level), rep("linear", 3), info = name) - pairwise <- out_all[out_all$Level1 == "Medium" & out_all$Level2 == "Low", ] - expect_equal( - out$Difference, - pairwise$Difference[match(out$Infl, pairwise$Infl)], - tolerance = 1e-10, - info = name - ) + out <- estimate_contrasts(m, "Sat=c('Low','High')", by = "Infl", backend = "emmeans") + expected <- out_all[out_all$Level1 == "High" & out_all$Level2 == "Low", ] + expect_identical(as.character(out$Infl), as.character(expected$Infl)) + expect_equal(out$Difference, expected$Difference, tolerance = 1e-10) - # one category cannot be contrasted - err <- expect_error(estimate_contrasts(m, "Sat='Low'", by = "Infl", backend = "emmeans")) - expect_match( - gsub("\\s+", " ", conditionMessage(err)), - "Contrasts of the response `Sat` need at least two response categories.", - fixed = TRUE, - info = name - ) + # comparisons that are not pairwise use only the named categories + out <- estimate_contrasts( + m, + "Sat=c('Low','Medium')", + by = "Infl", + comparison = "poly", + backend = "emmeans" + ) + expected <- out_all[out_all$Level1 == "Medium" & out_all$Level2 == "Low", ] + expect_equal(out$Difference, expected$Difference, tolerance = 1e-10) - # values that are no response categories, and modes without categories - err <- expect_error(estimate_contrasts( - m, - "Sat=c('Low','Bogus')", - by = "Infl", - backend = "emmeans" - )) - expect_match( - gsub("\\s+", " ", conditionMessage(err)), - "The response `Sat` has no category `Bogus`.", - fixed = TRUE, - info = name - ) - err <- expect_error(estimate_contrasts( + expect_error( + estimate_contrasts(m, "Sat='Low'", by = "Infl", backend = "emmeans"), + "need at least two response categories", + fixed = TRUE + ) + expect_error( + estimate_contrasts( m, "Sat=c('Low','High')", by = "Infl", predict = "cum.prob", backend = "emmeans" - )) - expect_match( - gsub("\\s+", " ", conditionMessage(err)), - "The response `Sat` can only be used in `contrast` with `predict = \"prob\"`", - fixed = TRUE, - info = name - ) - } + ), + "can only be used in `contrast`", + fixed = TRUE + ) + expect_error( + estimate_contrasts(m, "Sat=c('Low','High')", by = "Sat='Low'", backend = "emmeans"), + "cannot be filtered in `by`", + fixed = TRUE + ) - # the response cannot also be filtered in `by` - err <- expect_error(estimate_contrasts( - models$polr, + out <- estimate_contrasts( + m, "Sat=c('Low','High')", - by = "Sat='Low'", + by = "Infl", + comparison = "joint", backend = "emmeans" - )) - expect_match( - gsub("\\s+", " ", conditionMessage(err)), - "The response `Sat` cannot be filtered in `by` when it is also in `contrast`.", - fixed = TRUE ) + estimated <- emmeans::emmeans(m, specs = c("Infl", "Sat"), mode = "prob") + estimated <- emmeans::regrid(estimated[which(estimated@grid$Sat %in% c("Low", "High"))]) + expected <- as.data.frame(emmeans::test( + emmeans::contrast(estimated, by = "Infl", method = "pairwise"), + joint = TRUE, + by = "Infl" + )) + expect_equal(out$Chisq, expected$Chisq, tolerance = 1e-6) }) @@ -1907,57 +920,32 @@ test_that("estimate_contrasts, backend emmeans, ordinal, contrast filter on the skip_if_not_installed("emmeans") skip_if_not_installed("ordinal") data(wine, package = "ordinal") - # emmeans codes the categories of clmm models as 1, 2, 3, ..., so we use - # names that are not the codes levels(wine$rating) <- letters[1:5] m <- ordinal::clmm(rating ~ temp + contact + (1 | judge), data = wine) out_all <- estimate_contrasts(m, "rating", by = "temp", backend = "emmeans") out <- estimate_contrasts(m, "rating=c('a','c')", by = "temp", backend = "emmeans") - # one contrast of the categories `c` and `a` (codes 3 and 1) per level of `temp` + # emmeans codes the categories of clmm models as 1, 2, 3, ... expected <- out_all[out_all$Level1 == "rating3" & out_all$Level2 == "rating1", ] - expect_identical(nrow(out), 2L) - expect_identical(as.character(out$temp), as.character(expected$temp)) expect_equal(out$Difference, expected$Difference, tolerance = 1e-10) - expect_equal(out$SE, expected$SE, tolerance = 1e-10) }) -test_that("estimate_contrasts, backend emmeans, ordinal, joint test with a contrast filter on the response", { - skip_if_not_installed("emmeans") - data(housing, package = "MASS") - m <- MASS::polr(Sat ~ Infl + Type, weights = Freq, data = housing, Hess = TRUE) - out <- estimate_contrasts( - m, - "Sat=c('Low','High')", - by = "Infl", - comparison = "joint", - backend = "emmeans" - ) - # the joint test of the contrasts of the two named categories in emmeans, - # after `regrid()` as in `get_emcontrasts()` - estimated <- emmeans::emmeans(m, specs = c("Infl", "Sat"), mode = "prob") - estimated <- estimated[which(as.character(estimated@grid$Sat) %in% c("Low", "High"))] - estimated <- emmeans::regrid(estimated) - expected <- as.data.frame(emmeans::test( - emmeans::contrast(estimated, by = "Infl", method = "pairwise"), - joint = TRUE, - by = "Infl" - )) - expect_identical(nrow(out), 3L) - expect_identical(as.character(out$Infl), as.character(expected$Infl)) - expect_equal(out$Chisq, expected$Chisq, tolerance = 1e-6) - expect_equal(out$p, expected$p.value, tolerance = 1e-6) - - # with another term in `contrast`, the joint test fails as without a filter - err <- expect_error(estimate_contrasts( - m, - c("Infl", "Sat=c('Low','High')"), - by = "Type", - comparison = "joint", - backend = "emmeans" - )) - expect_match(conditionMessage(err), "Could not compute joint test", fixed = TRUE) -}) +# the rows of `out` must be the rows of `out_all` with the same values in the +# `keys` columns, with equal numeric columns +expect_matching_rows <- function(out, out_all, keys, columns, info = NULL) { + key <- function(d) do.call(paste, c(lapply(d[keys], as.character), sep = "\r")) + expect_gt(nrow(out), 0) + matched <- match(key(out), key(out_all)) + expect_false(anyNA(matched), info = info) + for (column in columns) { + expect_equal( + out[[column]], + out_all[[column]][matched], + tolerance = 1e-10, + info = paste(info, column) + ) + } +} test_that("estimate_contrasts, backend emmeans, ordinal, slopes with a by filter that keeps one level", { diff --git a/tests/testthat/test-p_adjust.R b/tests/testthat/test-p_adjust.R index 29b2da0c7..bacaf20dd 100644 --- a/tests/testthat/test-p_adjust.R +++ b/tests/testthat/test-p_adjust.R @@ -2,7 +2,6 @@ skip_on_cran() skip_if_not_installed("emmeans") skip_if_not_installed("marginaleffects") -# data and models for the p-value adjustment tests set.seed(12) n <- 200 d_adjust <- data.frame( @@ -17,150 +16,52 @@ d_adjust$y <- 1 + m1 <- lm(y ~ g * h, data = d_adjust) m2 <- lm(y ~ x * g, data = d_adjust) -# p-values of emmeans, in the row order of modelbased. modelbased compares -# "b - a" where emmeans compares "a - b", so rows are matched by the sorted -# pair of levels, and by the levels of `by`, if any -emmeans_p <- function(out, em, by = NULL) { - em <- as.data.frame(summary(em)) - pairs_em <- strsplit(as.character(em$contrast), " - ", fixed = TRUE) - key_em <- vapply(pairs_em, function(i) paste(sort(gsub("[()]", "", i)), collapse = "|"), "") - if (all(c("Level1", "Level2") %in% colnames(out))) { - key_out <- vapply( - seq_len(nrow(out)), - function(i) paste(sort(c(as.character(out$Level1[i]), as.character(out$Level2[i]))), collapse = "|"), - "" - ) - } else { - # compared slopes, labelled like "(b) - (a)" - pairs_out <- strsplit(as.character(out$Parameter), " - ", fixed = TRUE) - key_out <- vapply(pairs_out, function(i) paste(sort(gsub("[()]", "", i)), collapse = "|"), "") - } - if (!is.null(by)) { - key_em <- paste(key_em, em[[by]]) - key_out <- paste(key_out, out[[by]]) - } - em$p.value[match(key_out, key_em)] -} +# p-values of emmeans; modelbased orders and orients the pairs differently +emmeans_p <- function(em) sort(as.data.frame(summary(em))$p.value) # the message can wrap over lines sidak_message <- "(?s)(?=.*Tukey)(?=.*Sidak)" -test_that("p_adjust tukey, pairwise contrasts, marginaleffects", { - out <- estimate_contrasts(m1, "g", p_adjust = "tukey") +test_that("p_adjust tukey and sidak, contrasts, marginaleffects", { em <- suppressMessages(emmeans::emmeans(m1, "g")) - expect_equal(out$p, emmeans_p(out, pairs(em, adjust = "tukey")), tolerance = 1e-5) - expect_identical(attributes(out)$p_adjust, "tukey") + out <- estimate_contrasts(m1, "g", p_adjust = "tukey") + expect_equal(sort(out$p), emmeans_p(pairs(em, adjust = "tukey")), tolerance = 1e-5) + out <- estimate_contrasts(m1, "g", p_adjust = "sidak") + expect_equal(sort(out$p), emmeans_p(pairs(em, adjust = "sidak")), tolerance = 1e-5) - out <- estimate_contrasts(m1, "g", by = "h", p_adjust = "tukey") em <- suppressMessages(emmeans::emmeans(m1, "g", by = "h")) - expect_equal(out$p, emmeans_p(out, pairs(em, adjust = "tukey"), by = "h"), tolerance = 1e-5) -}) - + out <- estimate_contrasts(m1, "g", by = "h", p_adjust = "tukey") + expect_equal(sort(out$p), emmeans_p(pairs(em, adjust = "tukey")), tolerance = 1e-5) -test_that("p_adjust tukey, non-pairwise contrasts use sidak, marginaleffects", { + expect_message( + { + out <- estimate_contrasts(m1, "g", comparison = "reference", p_adjust = "tukey") + }, + regexp = sidak_message, + perl = TRUE + ) em <- suppressMessages(emmeans::emmeans(m1, "g")) - comparisons <- list("reference", ~reference, "sequential", "meandev") - methods <- c("trt.vs.ctrl", "trt.vs.ctrl", "consec", "eff") - for (i in seq_along(comparisons)) { - expect_message( - { - out <- estimate_contrasts( - m1, - "g", - comparison = comparisons[[i]], - p_adjust = "tukey" - ) - }, - regexp = sidak_message, - perl = TRUE - ) - em_p <- as.data.frame(summary(emmeans::contrast(em, methods[i], adjust = "sidak")))$p.value - expect_equal(out$p, em_p, tolerance = 1e-5) - expect_identical(attributes(out)$p_adjust, "sidak") - expect_no_message( - estimate_contrasts( - m1, - "g", - comparison = comparisons[[i]], - p_adjust = "tukey", - verbose = FALSE - ), - message = "Tukey|Sidak" - ) - } - - # the example of the issue - mt <- mtcars - mt$gear <- factor(mt$gear) - m5 <- lm(mpg ~ gear, data = mt) - out <- suppressMessages(estimate_contrasts( - m5, - "gear", - comparison = "reference", - p_adjust = "tukey" - )) em_p <- as.data.frame(summary(emmeans::contrast( - emmeans::emmeans(m5, "gear"), + em, "trt.vs.ctrl", adjust = "sidak" )))$p.value expect_equal(out$p, em_p, tolerance = 1e-5) + expect_identical(attributes(out)$p_adjust, "sidak") }) -test_that("p_adjust sidak, contrasts, family size, marginaleffects", { - out <- estimate_contrasts(m1, "g", p_adjust = "sidak") - em <- suppressMessages(emmeans::emmeans(m1, "g")) - expect_equal(out$p, emmeans_p(out, pairs(em, adjust = "sidak")), tolerance = 1e-5) - - out <- estimate_contrasts(m1, "g", by = "h", p_adjust = "sidak") - em <- suppressMessages(emmeans::emmeans(m1, "g", by = "h")) - expect_equal(out$p, emmeans_p(out, pairs(em, adjust = "sidak"), by = "h"), tolerance = 1e-5) - - em <- suppressMessages(emmeans::emmeans(m1, "g")) - comparisons <- list("reference", ~reference, "sequential", "meandev") - methods <- c("trt.vs.ctrl", "trt.vs.ctrl", "consec", "eff") - for (i in seq_along(comparisons)) { - expect_no_message( - { - out <- estimate_contrasts( - m1, - "g", - comparison = comparisons[[i]], - p_adjust = "sidak" - ) - }, - message = "Sidak" - ) - em_p <- as.data.frame(summary(emmeans::contrast(em, methods[i], adjust = "sidak")))$p.value - expect_equal(out$p, em_p, tolerance = 1e-5) - } -}) - - -test_that("p_adjust tukey and sidak, compared slopes, marginaleffects", { - em <- emmeans::emtrends(m2, "g", var = "x") - for (p_adjust in c("tukey", "sidak")) { - em_pairs <- pairs(em, adjust = p_adjust) - expect_no_message( - { - out <- estimate_slopes(m2, "x", by = "g", comparison = ~pairwise, p_adjust = p_adjust) - }, - message = "Sidak" - ) - expect_equal(out$p, emmeans_p(out, em_pairs), tolerance = 1e-5) - expect_no_message( - { - out <- estimate_contrasts(m2, contrast = "x", by = "g", p_adjust = p_adjust) - }, - message = "Sidak" - ) - expect_equal(out$p, emmeans_p(out, em_pairs), tolerance = 1e-5) - } +test_that("p_adjust tukey, compared slopes, marginaleffects", { + em_p <- emmeans_p(pairs(emmeans::emtrends(m2, "g", var = "x"), adjust = "tukey")) + out <- estimate_slopes(m2, "x", by = "g", comparison = ~pairwise, p_adjust = "tukey") + expect_equal(sort(out$p), em_p, tolerance = 1e-5) + out <- estimate_contrasts(m2, contrast = "x", by = "g", p_adjust = "tukey") + expect_equal(sort(out$p), em_p, tolerance = 1e-5) + out <- estimate_slopes(m2, "x", by = "g", hypo = ~pairwise, p_adjust = "tukey") + expect_equal(sort(out$p), em_p, tolerance = 1e-5) - # post-processing: unadjusted p-values are based on z, so we only check the - # number of slopes in the Tukey adjustment, which uses the model's df + # unadjusted p-values are based on z, so we only check the number of slopes out <- suppressMessages(estimate_slopes( m2, "x", @@ -174,6 +75,25 @@ test_that("p_adjust tukey and sidak, compared slopes, marginaleffects", { stats::ptukey(sqrt(2) * abs(as.vector(out[[z]])), 4, stats::df.residual(m2), lower.tail = FALSE), tolerance = 1e-8 ) + + m3 <- lm(y ~ x * g * h, data = d_adjust) + em_p <- emmeans_p(pairs(emmeans::emtrends(m3, ~ g | h, var = "x"), adjust = "tukey")) + out <- estimate_slopes( + m3, + "x", + by = c("g", "h"), + comparison = ~ pairwise | h, + p_adjust = "tukey" + ) + expect_equal(sort(out$p), em_p, tolerance = 1e-5) + out <- estimate_contrasts( + m3, + contrast = "x", + by = c("g", "h"), + comparison = ~ pairwise | h, + p_adjust = "tukey" + ) + expect_equal(sort(out$p), em_p, tolerance = 1e-5) }) @@ -191,71 +111,10 @@ test_that("p_adjust tukey, slopes tested against a value, marginaleffects", { adjust = "sidak" )))$p.value expect_equal(out$p, em_p, tolerance = 1e-5) - expect_identical(attributes(out)$p_adjust, "sidak") - - # partly matched `hypothesis` argument - em_pairs <- pairs(emmeans::emtrends(m2, "g", var = "x"), adjust = "tukey") - expect_no_message( - { - out <- estimate_slopes(m2, "x", by = "g", hypo = ~pairwise, p_adjust = "tukey") - }, - message = "Sidak" - ) - expect_equal(out$p, emmeans_p(out, em_pairs), tolerance = 1e-5) }) -test_that("p_adjust, slopes, emmeans backend", { - mt <- mtcars - mt$cyl <- factor(mt$cyl) - m4 <- lm(mpg ~ wt * cyl, data = mt) - em <- emmeans::emtrends(m4, "cyl", var = "wt") - for (p_adjust in c("holm", "bonferroni", "sidak", "tukey", "none")) { - out <- estimate_slopes(m4, "wt", by = "cyl", backend = "emmeans", p_adjust = p_adjust) - em_p <- suppressMessages(as.data.frame(summary(emmeans::test(em, adjust = p_adjust)))$p.value) - expect_equal(out$p, em_p, tolerance = 1e-5) - } - p_none <- estimate_slopes(m4, "wt", by = "cyl", backend = "emmeans", p_adjust = "none")$p - p_holm <- estimate_slopes(m4, "wt", by = "cyl", backend = "emmeans", p_adjust = "holm")$p - expect_false(isTRUE(all.equal(p_none, p_holm))) -}) - - -test_that("p_adjust tukey and sidak, slopes compared pairwise within groups", { - m3 <- lm(y ~ x * g * h, data = d_adjust) - em <- emmeans::emtrends(m3, ~ g | h, var = "x") - for (p_adjust in c("tukey", "sidak")) { - expect_no_message( - { - out <- estimate_slopes( - m3, - "x", - by = c("g", "h"), - comparison = ~ pairwise | h, - p_adjust = p_adjust - ) - }, - message = "Sidak" - ) - expect_equal(out$p, emmeans_p(out, pairs(em, adjust = p_adjust), by = "h"), tolerance = 1e-5) - } -}) - - -test_that("p_adjust, families of grouped contrasts and of contrasts of predictions", { - # contrasts of slopes with a group term, Tukey counts the slopes in each group - m3 <- lm(y ~ x * g * h, data = d_adjust) - out <- estimate_contrasts( - m3, - contrast = "x", - by = c("g", "h"), - comparison = ~ pairwise | h, - p_adjust = "tukey" - ) - em <- emmeans::emtrends(m3, "g", by = "h", var = "x") - expect_equal(out$p, emmeans_p(out, pairs(em, adjust = "tukey"), by = "h"), tolerance = 1e-5) - - # contrasts of `estimate_relation()` output, Sidak counts within `by` groups +test_that("p_adjust, contrasts of estimate_relation() output", { out <- estimate_contrasts( estimate_relation(m1, by = c("g", "h")), contrast = "g", @@ -263,7 +122,7 @@ test_that("p_adjust, families of grouped contrasts and of contrasts of predictio p_adjust = "sidak" ) em <- suppressMessages(emmeans::emmeans(m1, "g", by = "h")) - expect_equal(out$p, emmeans_p(out, pairs(em, adjust = "sidak"), by = "h"), tolerance = 1e-5) + expect_equal(sort(out$p), emmeans_p(pairs(em, adjust = "sidak")), tolerance = 1e-5) # missing values of a `by` variable form their own family expect_identical( @@ -273,7 +132,7 @@ test_that("p_adjust, families of grouped contrasts and of contrasts of predictio }) -test_that("p_adjust, no Sidak message when nothing is adjusted, and method names", { +test_that("p_adjust, no Sidak message when nothing is adjusted", { # one contrast or one slope: Sidak does not change the p-value expect_no_message( estimate_contrasts(m1, "g", comparison = "b1 = b2", p_adjust = "tukey"), @@ -294,19 +153,14 @@ test_that("p_adjust, no Sidak message when nothing is adjusted, and method names ), regexp = "not compared pairwise" ) +}) - # emmeans backend: method names in any case, and the `p_adjust` attribute + +test_that("p_adjust, slopes, emmeans backend", { mt <- mtcars mt$cyl <- factor(mt$cyl) m4 <- lm(mpg ~ wt * cyl, data = mt) em <- emmeans::emtrends(m4, "cyl", var = "wt") - out <- estimate_slopes(m4, "wt", by = "cyl", backend = "emmeans", p_adjust = "bh") - expect_equal( - out$p, - as.data.frame(summary(emmeans::test(em, adjust = "BH")))$p.value, - tolerance = 1e-5 - ) - expect_identical(attributes(out)$p_adjust, "BH") out <- suppressMessages(estimate_slopes( m4, "wt", @@ -314,5 +168,14 @@ test_that("p_adjust, no Sidak message when nothing is adjusted, and method names backend = "emmeans", p_adjust = "tukey" )) + em_p <- suppressMessages( + as.data.frame(summary(emmeans::test(em, adjust = "tukey")))$p.value + ) + expect_equal(out$p, em_p, tolerance = 1e-5) expect_identical(attributes(out)$p_adjust, "sidak") + + out <- estimate_slopes(m4, "wt", by = "cyl", backend = "emmeans", p_adjust = "bh") + em_p <- as.data.frame(summary(emmeans::test(em, adjust = "BH")))$p.value + expect_equal(out$p, em_p, tolerance = 1e-5) + expect_identical(attributes(out)$p_adjust, "BH") }) diff --git a/tests/testthat/test-table_footer.R b/tests/testthat/test-table_footer.R index 2c315ff63..0ab4d775e 100644 --- a/tests/testthat/test-table_footer.R +++ b/tests/testthat/test-table_footer.R @@ -22,132 +22,87 @@ test_that("table_footer", { test_that("table_footer, prediction intervals", { pi_line <- "Intervals are prediction intervals." - m_lm <- lm(mpg ~ wt, data = mtcars) - m_glm <- glm(vs ~ wt, data = mtcars, family = binomial()) - - out <- utils::capture.output(estimate_prediction(m_lm, data = "grid", length = 3)) - expect_identical(sum(out == pi_line), 1L) - out <- utils::capture.output(estimate_prediction(m_glm, data = "grid", length = 3)) + m <- lm(mpg ~ wt, data = mtcars) + out <- utils::capture.output(estimate_prediction(m, data = "grid", length = 3)) expect_identical(sum(out == pi_line), 1L) - out <- utils::capture.output(estimate_relation(m_lm, predict = "prediction", length = 3)) + out <- utils::capture.output(estimate_relation(m, predict = "prediction", length = 3)) expect_identical(sum(out == pi_line), 1L) - # no intervals, no line - out <- utils::capture.output(estimate_prediction(m_lm, data = "grid", length = 3, ci = NULL)) + out <- utils::capture.output(estimate_prediction( + m, + data = "grid", + length = 3, + ci = NULL + )) + expect_false(pi_line %in% out) + out <- utils::capture.output(estimate_expectation(m, data = "grid", length = 3)) expect_false(pi_line %in% out) - - # no prediction intervals - for (m in list(m_lm, m_glm)) { - out <- utils::capture.output(estimate_expectation(m, data = "grid", length = 3)) - expect_false(pi_line %in% out) - out <- utils::capture.output(estimate_relation(m, length = 3)) - expect_false(pi_line %in% out) - out <- utils::capture.output(estimate_link(m, length = 3)) - expect_false(pi_line %in% out) - } - - skip_if_not_installed("lme4") - m_lmer <- lme4::lmer(mpg ~ wt + (1 | cyl), data = mtcars) - out <- utils::capture.output(estimate_prediction(m_lmer, data = "grid", length = 3)) - expect_identical(sum(out == pi_line), 1L) }) -test_that("table_footer, link-scale without transform", { +test_that("table_footer, link-scale and transform", { d <- mtcars d$cyl <- factor(d$cyl) g <- glm(am ~ wt + cyl, data = d, family = binomial) - pred_line <- "Predictions are on the link-scale." - contr_line <- "Contrasts are on the link-scale." - - out <- utils::capture.output(estimate_link(g, by = "cyl")) - expect_true(pred_line %in% out) - for (backend in c("marginaleffects", "emmeans")) { - out <- utils::capture.output( - estimate_means(g, "cyl", predict = "link", backend = backend) - ) - expect_true(pred_line %in% out) - out <- utils::capture.output( - estimate_contrasts(g, "cyl", predict = "link", backend = backend) - ) - expect_true(contr_line %in% out) - } -}) - -test_that("table_footer, link-scale with transform", { - d <- mtcars - d$cyl <- factor(d$cyl) - g <- glm(am ~ wt + cyl, data = d, family = binomial) + out <- utils::capture.output(estimate_means( + g, + "cyl", + predict = "link", + backend = "emmeans" + )) + expect_true("Predictions are on the link-scale." %in% out) + out <- utils::capture.output(estimate_contrasts( + g, + "cyl", + predict = "link", + backend = "emmeans" + )) + expect_true("Contrasts are on the link-scale." %in% out) - # a transformation function - out <- utils::capture.output( - estimate_contrasts(g, "cyl", predict = "link", transform = exp) - ) + out <- utils::capture.output(estimate_contrasts( + g, + "cyl", + predict = "link", + transform = exp + )) expect_true("Contrasts are on the transformed link-scale." %in% out) - out <- utils::capture.output(estimate_link(g, by = "cyl", transform = exp)) - expect_true("Predictions are on the transformed link-scale." %in% out) - out <- utils::capture.output( - estimate_means(g, "cyl", predict = "link", transform = exp) - ) - expect_true("Predictions are on the transformed link-scale." %in% out) + out <- utils::capture.output(estimate_contrasts( + g, + "cyl", + predict = "link", + transform = "exp" + )) + expect_true("Contrasts are on the exp-transformed link-scale." %in% out) - # transform = TRUE - out <- utils::capture.output(estimate_link(g, by = "cyl", transform = TRUE)) - expect_false(any(grepl("TRUE-transformed|FALSE-transformed", out))) # the response is not transformed, so `transform = TRUE` changes nothing - out <- utils::capture.output( - estimate_means(g, "cyl", predict = "link", transform = TRUE) - ) + out <- utils::capture.output(estimate_link(g, by = "cyl", transform = TRUE)) expect_true("Predictions are on the link-scale." %in% out) - - # name of the transformation - out <- utils::capture.output( - estimate_contrasts(g, "cyl", predict = "link", transform = "exp") - ) - expect_true("Contrasts are on the exp-transformed link-scale." %in% out) -}) - - -test_that("table_footer, scale of predictions with earlier output as data", { - d <- mtcars - d$cyl <- factor(d$cyl) - g <- glm(am ~ wt + cyl, data = d, family = binomial) - resp_line <- "Predictions are on the response-scale." - - out <- utils::capture.output( - estimate_expectation(g, data = estimate_link(g, by = "cyl")) - ) - expect_true(resp_line %in% out) - out <- utils::capture.output( - estimate_prediction(g, data = estimate_link(g, by = "cyl")) - ) - expect_true(resp_line %in% out) - out <- utils::capture.output( - estimate_link(g, data = estimate_expectation(g, by = "cyl")) - ) + out <- utils::capture.output(estimate_means( + g, + "cyl", + predict = "link", + transform = TRUE + )) expect_true("Predictions are on the link-scale." %in% out) -}) - -test_that("table_footer, link-scale for contrasts of predictions and slopes", { - d <- mtcars - d$cyl <- factor(d$cyl) - g <- glm(am ~ wt + cyl, data = d, family = binomial) + out <- utils::capture.output(estimate_expectation( + g, + data = estimate_link(g, by = "cyl") + )) + expect_true("Predictions are on the response-scale." %in% out) + out <- utils::capture.output(estimate_link( + g, + data = estimate_expectation(g, by = "cyl") + )) + expect_true("Predictions are on the link-scale." %in% out) - # contrasts of an earlier output - out <- utils::capture.output(estimate_contrasts(estimate_link(g, by = "cyl"))) - expect_true("Contrasts are on the link-scale." %in% out) - out <- utils::capture.output( - estimate_contrasts(estimate_link(g, by = "cyl", transform = exp)) - ) + out <- utils::capture.output(estimate_contrasts(estimate_link( + g, + by = "cyl", + transform = exp + ))) expect_true("Contrasts are on the transformed link-scale." %in% out) - - # contrasts of slopes - out <- utils::capture.output( - estimate_contrasts(g, "wt", by = "cyl", predict = "link") - ) - expect_true("Contrasts are on the link-scale." %in% out) out <- utils::capture.output( estimate_contrasts(g, "wt", by = "cyl", predict = "link", transform = exp) ) diff --git a/tests/testthat/test-visualisation_recipe.R b/tests/testthat/test-visualisation_recipe.R index e92ae46de..7871bd981 100644 --- a/tests/testthat/test-visualisation_recipe.R +++ b/tests/testthat/test-visualisation_recipe.R @@ -224,10 +224,6 @@ test_that("visualization_recipe, grid argument modifies facet_grid layer", { r[[which(geoms == "facet_grid")]] } - # default - out <- .grid_layer(visualisation_recipe(x)) - expect_identical(out$scales, "free_x") - # grid argument is used for the facet_grid layer out <- .grid_layer(visualisation_recipe(x, grid = list(scales = "free"))) expect_identical(out$scales, "free") @@ -259,109 +255,65 @@ test_that("visualization_recipe, grid argument modifies facet_grid layer", { ) } -# printed table as a list of character columns -.printed_table <- function(x) { - lapply(as.list(format(x)), as.character) -} - -# the printed table must not depend on the model information -.expect_same_print <- function(x) { - y <- x - attr(y, "model_info") <- NULL - expect_identical(.printed_table(x), .printed_table(y)) -} - test_that("visualization_recipe, emmeans ordinal, color by Response", { skip_if_not_installed("emmeans") skip_if_not_installed("MASS") - skip_if_not_installed("ordinal") - skip_if_not_installed("marginaleffects", minimum_version = "0.29.0") data(housing, package = "MASS") - models <- list( - polr = MASS::polr(Sat ~ Infl + Type, weights = Freq, data = housing, Hess = TRUE), - clm = ordinal::clm(Sat ~ Infl + Type, weights = Freq, data = housing) - ) - for (m in models) { - for (backend in c("emmeans", "marginaleffects")) { - x <- estimate_means(m, by = "Infl", backend = backend) - out <- .recipe_mapping(x) - expect_identical(out$x, "Infl") - expect_identical(out$color, "Response") - expect_null(out$facet) - expect_true("scale_y_continuous" %in% out$geoms) - expect_identical(attributes(x)$focal_terms, "Infl") - if (backend == "emmeans") { - .expect_same_print(x) - } - - x <- estimate_means(m, by = c("Infl", "Type"), backend = backend) - out <- .recipe_mapping(x) - expect_identical(out$x, "Infl") - expect_identical(out$color, "Type") - expect_identical(out$facet, "~Response") - expect_identical(attributes(x)$focal_terms, c("Infl", "Type")) - if (backend == "emmeans") .expect_same_print(x) - } - } - - # printed rows keep their order: the response varies within predictor levels - x <- estimate_means(models$polr, by = c("Infl", "Type"), backend = "emmeans") - printed <- .printed_table(x) - expect_identical( - printed$Infl[1:6], - c("Low", "Low", "Low", "Medium", "Medium", "Medium") - ) - expect_identical(printed$Type[1:6], rep("Tower", 6)) - expect_identical( - printed$Response[1:6], - c("Low", "Medium", "High", "Low", "Medium", "High") - ) - - # mixed models from package ordinal - data(wine, package = "ordinal") - m <- ordinal::clmm(rating ~ temp + contact + (1 | judge), data = wine) - x <- estimate_means(m, by = "temp", backend = "emmeans") + m <- MASS::polr(Sat ~ Infl + Type, weights = Freq, data = housing, Hess = TRUE) + x <- estimate_means(m, by = "Infl", backend = "emmeans") out <- .recipe_mapping(x) - expect_identical(out$x, "temp") + expect_identical(out$x, "Infl") expect_identical(out$color, "Response") - .expect_same_print(x) + expect_true("scale_y_continuous" %in% out$geoms) + expect_identical(attributes(x)$focal_terms, "Infl") + + x <- estimate_means(m, by = c("Infl", "Type"), backend = "emmeans") + out <- .recipe_mapping(x) + expect_identical(out$color, "Type") + expect_identical(out$facet, "~Response") + + # printed rows keep their order: the response varies within predictor levels, + # and the printed table does not depend on the model information + printed <- lapply(as.list(format(x)), as.character) + expect_identical(printed$Response[1:6], rep(c("Low", "Medium", "High"), 2)) + attr(x, "model_info") <- NULL + expect_identical(lapply(as.list(format(x)), as.character), printed) }) test_that("visualization_recipe, emmeans ordinal, color by Threshold", { skip_if_not_installed("emmeans") skip_if_not_installed("MASS") - skip_if_not_installed("ordinal") data(housing, package = "MASS") - models <- list( - polr = MASS::polr(Sat ~ Infl + Type, weights = Freq, data = housing, Hess = TRUE), - clm = ordinal::clm(Sat ~ Infl + Type, weights = Freq, data = housing) - ) - for (m in models) { - for (p in c("cum.prob", "exc.prob", "linear.predictor")) { - x <- estimate_means(m, by = "Infl", backend = "emmeans", predict = p) - out <- .recipe_mapping(x) - expect_identical(out$x, "Infl") - expect_identical(out$color, "Threshold") - expect_null(out$facet) - expect_identical( - "scale_y_continuous" %in% out$geoms, - p %in% c("cum.prob", "exc.prob") - ) - expect_identical(attributes(x)$focal_terms, "Infl") - .expect_same_print(x) - - x <- estimate_means(m, by = c("Infl", "Type"), backend = "emmeans", predict = p) - out <- .recipe_mapping(x) - expect_identical(out$x, "Infl") - expect_identical(out$color, "Type") - expect_identical(out$facet, "~Threshold") - .expect_same_print(x) - } - } + m <- MASS::polr(Sat ~ Infl + Type, weights = Freq, data = housing, Hess = TRUE) + out <- .recipe_mapping(estimate_means( + m, + by = "Infl", + backend = "emmeans", + predict = "cum.prob" + )) + expect_identical(out$color, "Threshold") + expect_true("scale_y_continuous" %in% out$geoms) + out <- .recipe_mapping(estimate_means( + m, + by = c("Infl", "Type"), + backend = "emmeans", + predict = "linear.predictor" + )) + expect_identical(out$color, "Type") + expect_identical(out$facet, "~Threshold") + expect_false("scale_y_continuous" %in% out$geoms) + # no grouping column for the latent scale + out <- .recipe_mapping(estimate_means( + m, + by = "Infl", + backend = "emmeans", + predict = "latent" + )) + expect_null(out$color) # thresholds in model order, not in alphabetical order data(wine, package = "ordinal") @@ -374,102 +326,35 @@ test_that("visualization_recipe, emmeans ordinal, color by Threshold", { m <- MASS::polr(rating ~ temp + contact, data = wine, Hess = TRUE) x <- estimate_means(m, by = "temp", backend = "emmeans", predict = "cum.prob") out <- .recipe_mapping(x) - expect_identical(out$color, "Threshold") expect_identical( levels(out$data$Threshold), c("very low|low", "low|medium", "medium|high", "high|very high") ) - expect_identical( - as.character(unique(out$data$Threshold)), - c("very low|low", "low|medium", "medium|high", "high|very high") - ) - .expect_same_print(x) }) test_that("visualization_recipe, ordinal, focal predictor named Threshold", { skip_if_not_installed("emmeans") skip_if_not_installed("MASS") - skip_if_not_installed("marginaleffects", minimum_version = "0.29.0") data(housing, package = "MASS") d <- housing d$Threshold <- d$Infl m <- MASS::polr(Sat ~ Threshold + Type, weights = Freq, data = d, Hess = TRUE) - for (backend in c("emmeans", "marginaleffects")) { - x <- estimate_means(m, by = "Threshold", backend = backend) - out <- .recipe_mapping(x) - expect_identical(out$x, "Threshold") - expect_identical(out$color, "Response") - expect_null(out$facet) - } + out <- .recipe_mapping(estimate_means(m, by = "Threshold", backend = "emmeans")) + expect_identical(out$x, "Threshold") + expect_identical(out$color, "Response") }) test_that("visualization_recipe, emmeans, percent scale for logistic models", { skip_if_not_installed("emmeans") - skip_if_not_installed("marginaleffects", minimum_version = "0.29.0") - - d <- mtcars - d$cyl <- factor(d$cyl) - m <- glm(am ~ cyl, family = binomial, data = d) - for (backend in c("emmeans", "marginaleffects")) { - x <- estimate_means(m, by = "cyl", backend = backend) - out <- .recipe_mapping(x) - expect_identical(out$geoms, c("pointrange", "labs", "scale_y_continuous")) - } - .expect_same_print(estimate_means(m, by = "cyl", backend = "emmeans")) -}) - - -test_that("visualization_recipe, emmeans, recipes that do not change", { - skip_if_not_installed("emmeans") - skip_if_not_installed("MASS") - - m <- lm(Sepal.Length ~ Species, data = iris) - x <- estimate_means(m, by = "Species", backend = "emmeans") - out <- .recipe_mapping(x, show_data = TRUE) - expect_identical(out$geoms, c("jitter", "pointrange", "labs")) - expect_identical(out$x, "Species") - expect_null(out$color) - .expect_same_print(x) - - m <- lm(log(Sepal.Length) ~ Species, data = iris) - x <- estimate_means(m, by = "Species", backend = "emmeans") - out <- .recipe_mapping(x, show_data = TRUE) - expect_identical(out$geoms, c("jitter", "pointrange", "labs")) - expect_identical(out$x, "Species") - expect_null(out$color) - .expect_same_print(x) d <- mtcars d$cyl <- factor(d$cyl) - m <- glm(carb ~ cyl, family = poisson, data = d) - x <- estimate_means(m, by = "cyl", backend = "emmeans") - out <- .recipe_mapping(x) - expect_identical(out$geoms, c("pointrange", "labs")) - expect_identical(out$x, "cyl") - expect_null(out$color) - .expect_same_print(x) - m <- glm(am ~ cyl, family = binomial, data = d) - x <- estimate_means(m, by = "cyl", backend = "emmeans", predict = "link") - out <- .recipe_mapping(x) - expect_identical(out$geoms, c("pointrange", "labs")) - expect_identical(out$x, "cyl") - expect_null(out$color) - .expect_same_print(x) - - data(housing, package = "MASS") - m <- MASS::polr(Sat ~ Infl + Type, weights = Freq, data = housing, Hess = TRUE) - for (p in c("mean.class", "latent")) { - x <- estimate_means(m, by = "Infl", backend = "emmeans", predict = p) - out <- .recipe_mapping(x) - expect_identical(out$geoms, c("pointrange", "labs")) - expect_identical(out$x, "Infl") - expect_null(out$color) - .expect_same_print(x) - } + out <- .recipe_mapping(estimate_means(m, by = "cyl", backend = "emmeans")) + expect_identical(out$geoms, c("pointrange", "labs", "scale_y_continuous")) }) @@ -519,19 +404,14 @@ test_that("visualization_recipe, informative error when nothing goes on x-axis", test_that("visualization_recipe, emmeans ordinal slopes, color by Response", { skip_if_not_installed("emmeans") skip_if_not_installed("MASS") - skip_if_not_installed("marginaleffects", minimum_version = "0.29.0") data(housing, package = "MASS") set.seed(1) housing$x <- as.numeric(housing$Cont) + stats::rnorm(nrow(housing)) m <- MASS::polr(Sat ~ Infl * x, weights = Freq, data = housing, Hess = TRUE) - for (backend in c("emmeans", "marginaleffects")) { - out <- .recipe_mapping(estimate_slopes(m, "x", by = "Infl", backend = backend)) - expect_identical(out$x, "Infl", info = backend) - expect_identical(out$color, "Response", info = backend) - } - - # thresholds, and no grouping column for the latent scale + out <- .recipe_mapping(estimate_slopes(m, "x", by = "Infl", backend = "emmeans")) + expect_identical(out$x, "Infl") + expect_identical(out$color, "Response") out <- .recipe_mapping(estimate_slopes( m, "x", @@ -540,12 +420,4 @@ test_that("visualization_recipe, emmeans ordinal slopes, color by Response", { backend = "emmeans" )) expect_identical(out$color, "Threshold") - out <- .recipe_mapping(estimate_slopes( - m, - "x", - by = "Infl", - predict = "latent", - backend = "emmeans" - )) - expect_null(out$color) })