From 146d5ae31a0e24e8bd14fddb7dd6474fdbf1b717 Mon Sep 17 00:00:00 2001 From: Jeffrey Girard Date: Mon, 19 Oct 2026 13:42:32 -0500 Subject: [PATCH 1/8] Test joint tests of ordinal models against a delta-method Wald test --- tests/testthat/test-ordinal.R | 160 ++++++++++++++++++++++++++++++++++ 1 file changed, 160 insertions(+) diff --git a/tests/testthat/test-ordinal.R b/tests/testthat/test-ordinal.R index 71b5d20af..c41896162 100644 --- a/tests/testthat/test-ordinal.R +++ b/tests/testthat/test-ordinal.R @@ -1336,3 +1336,163 @@ test_that("estimate_means, backend emmeans, ordinal, several by terms on the res expect_response_filter_rows(out, out_all, "Low") expect_identical(nrow(out), 3L) }) + + +# 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` gives + # weight 1/4 to the probabilities of one category at one `Infl` level + 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) + } +}) + + +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" in a step that 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", + comparison = "joint", + backend = "emmeans" + ) + }) + expect_identical(nrow(out), 12L) +}) From cac451823d1bbc32c25164036270d1b4654c9b67 Mon Sep 17 00:00:00 2001 From: Jeffrey Girard Date: Mon, 19 Oct 2026 13:43:20 -0500 Subject: [PATCH 2/8] Regrid ordinal emmeans grids before joint tests --- R/get_emcontrasts.R | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/R/get_emcontrasts.R b/R/get_emcontrasts.R index 4dc665a17..3941ed7da 100644 --- a/R/get_emcontrasts.R +++ b/R/get_emcontrasts.R @@ -137,6 +137,13 @@ get_emcontrasts <- function( specs <- .add_emmeans_pseudo_factor(my_args$emmeans_specs, ordinal$pseudo_factor) fun_args <- c(fun_args, list(specs = specs, at = my_args$emmeans_at)) estimated <- suppressMessages(do.call(emmeans::emmeans, fun_args)) + # Joint tests use the linear functions of the emmeans grid. For contrasts + # between the response categories of ordinal models, these gave values + # about 50 times the Wald statistic of the contrasts (emmeans 2.0.4). + # After `regrid()`, the joint tests equal the Wald statistic. + if (!is.null(ordinal$mode) && identical(comparison, "joint")) { + estimated <- emmeans::regrid(estimated) + } # Find by variables emm_by <- specs[!specs %in% my_args$contrast] if (length(emm_by) == 0) { From c7ba7ffc5154db2ae49288bc8537d8d1c03a7315 Mon Sep 17 00:00:00 2001 From: Jeffrey Girard Date: Mon, 19 Oct 2026 13:46:51 -0500 Subject: [PATCH 3/8] Print the Chisq column of emmeans joint tests with two decimals --- R/get_jointtest.R | 6 ++++++ tests/testthat/test-ordinal.R | 15 +++++++++++++++ 2 files changed, 21 insertions(+) diff --git a/R/get_jointtest.R b/R/get_jointtest.R index 3908feb3d..8e0c85819 100644 --- a/R/get_jointtest.R +++ b/R/get_jointtest.R @@ -153,6 +153,12 @@ result$df1 <- insight::format_value(result$df1, protect_integers = TRUE) result$df2 <- insight::format_value(result$df2, protect_integers = TRUE) result$`F.ratio` <- insight::format_value(result$`F.ratio`, protect_integers = TRUE) + # "insight::format_table()" uses the attribute "digits" for the "Chisq" + # column. For objects of class "summary_emm", `as.data.frame()` sets it to + # `getOption("digits")` when it is missing, so we return a data frame + # without it + class(result) <- "data.frame" + attr(result, "digits") <- NULL # rename statistic column colnames(result)[1] <- "Contrast" diff --git a/tests/testthat/test-ordinal.R b/tests/testthat/test-ordinal.R index c41896162..2479886a0 100644 --- a/tests/testthat/test-ordinal.R +++ b/tests/testthat/test-ordinal.R @@ -1496,3 +1496,18 @@ test_that("estimate_contrasts, backend emmeans, ordinal, joint tests of a predic }) expect_identical(nrow(out), 12L) }) + + +test_that("estimate_contrasts, backend emmeans, ordinal, printed Chisq of joint tests", { + 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") + 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))) +}) From 09252ddd591db52c124b3a40ca4d6566fb632058 Mon Sep 17 00:00:00 2001 From: Jeffrey Girard Date: Mon, 19 Oct 2026 13:47:48 -0500 Subject: [PATCH 4/8] NEWS entry for ordinal joint tests; bump version --- DESCRIPTION | 2 +- NEWS.md | 9 +++++++++ 2 files changed, 10 insertions(+), 1 deletion(-) diff --git a/DESCRIPTION b/DESCRIPTION index b1c5da55f..72c0fb47f 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,7 +1,7 @@ Type: Package Package: modelbased Title: Estimation of Model-Based Predictions, Contrasts and Means -Version: 0.18.0.2 +Version: 0.18.0.3 Authors@R: c(person(given = "Dominique", family = "Makowski", diff --git a/NEWS.md b/NEWS.md index cdd4a67c5..7e5cd8675 100644 --- a/NEWS.md +++ b/NEWS.md @@ -43,6 +43,15 @@ replaced by `"sidak"` with a message, because the slopes are tested against 0. +* Joint tests with `comparison = "joint"` and `backend = "emmeans"` now give + the correct statistics for contrasts between the response categories of + ordinal models, such as `MASS::polr()` and `ordinal::clm()`. An example is + `estimate_contrasts(m, "Sat", by = "Infl", comparison = "joint")`. Before, + the F and Chi-squared statistics were about 50 times too large. Joint tests + of ordinal models also no longer warn with "NaNs produced". The `Chisq` + column of joint tests with `backend = "emmeans"` now prints with two + decimals. + # modelbased 0.18.0 ## Breaking Changes From 983f202f72871b1459e96b0727fccc4444c71ae8 Mon Sep 17 00:00:00 2001 From: Jeffrey Girard Date: Mon, 19 Oct 2026 13:50:30 -0500 Subject: [PATCH 5/8] Correct the factor in the NEWS entry and comments; add backend to the NEWS example --- NEWS.md | 8 +++++--- R/get_emcontrasts.R | 2 +- R/get_jointtest.R | 9 +++++---- tests/testthat/test-ordinal.R | 9 +++++---- 4 files changed, 16 insertions(+), 12 deletions(-) diff --git a/NEWS.md b/NEWS.md index 7e5cd8675..b37607cc3 100644 --- a/NEWS.md +++ b/NEWS.md @@ -46,9 +46,11 @@ * Joint tests with `comparison = "joint"` and `backend = "emmeans"` now give the correct statistics for contrasts between the response categories of ordinal models, such as `MASS::polr()` and `ordinal::clm()`. An example is - `estimate_contrasts(m, "Sat", by = "Infl", comparison = "joint")`. Before, - the F and Chi-squared statistics were about 50 times too large. Joint tests - of ordinal models also no longer warn with "NaNs produced". The `Chisq` + `estimate_contrasts(m, "Sat", by = "Infl", comparison = "joint", + backend = "emmeans")`. Before, the F and Chi-squared statistics were too + large, for example by a factor of 4 to 52 for + `polr(Sat ~ Infl + Type, data = MASS::housing)`. Joint tests of + ordinal models also no longer warn with "NaNs produced". The `Chisq` column of joint tests with `backend = "emmeans"` now prints with two decimals. diff --git a/R/get_emcontrasts.R b/R/get_emcontrasts.R index 3941ed7da..5b14d9cdb 100644 --- a/R/get_emcontrasts.R +++ b/R/get_emcontrasts.R @@ -139,7 +139,7 @@ get_emcontrasts <- function( estimated <- suppressMessages(do.call(emmeans::emmeans, fun_args)) # Joint tests use the linear functions of the emmeans grid. For contrasts # between the response categories of ordinal models, these gave values - # about 50 times the Wald statistic of the contrasts (emmeans 2.0.4). + # 4 to 52 times the Wald statistic of the contrasts (emmeans 2.0.4). # After `regrid()`, the joint tests equal the Wald statistic. if (!is.null(ordinal$mode) && identical(comparison, "joint")) { estimated <- emmeans::regrid(estimated) diff --git a/R/get_jointtest.R b/R/get_jointtest.R index 8e0c85819..8bba93ff1 100644 --- a/R/get_jointtest.R +++ b/R/get_jointtest.R @@ -153,10 +153,11 @@ result$df1 <- insight::format_value(result$df1, protect_integers = TRUE) result$df2 <- insight::format_value(result$df2, protect_integers = TRUE) result$`F.ratio` <- insight::format_value(result$`F.ratio`, protect_integers = TRUE) - # "insight::format_table()" uses the attribute "digits" for the "Chisq" - # column. For objects of class "summary_emm", `as.data.frame()` sets it to - # `getOption("digits")` when it is missing, so we return a data frame - # without it + # "insight::format_table()" uses the attribute "digits" as the default + # number of digits for numeric columns, and "Chisq" is the only one not + # formatted above. For objects of class "summary_emm", `as.data.frame()` + # sets it to `getOption("digits")` when it is missing, so we return a data + # frame without it class(result) <- "data.frame" attr(result, "digits") <- NULL diff --git a/tests/testthat/test-ordinal.R b/tests/testthat/test-ordinal.R index 2479886a0..7eabe267a 100644 --- a/tests/testthat/test-ordinal.R +++ b/tests/testthat/test-ordinal.R @@ -1399,8 +1399,9 @@ test_that("estimate_contrasts, backend emmeans, ordinal, joint tests of the resp ) # reference 2: marginaleffects, with the same contrasts of the 36 predicted - # probabilities (4 `Type` x 3 `Infl` x 3 `Sat`). Each column of `H` gives - # weight 1/4 to the probabilities of one category at one `Infl` level + # 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), @@ -1461,8 +1462,8 @@ test_that("estimate_contrasts, backend emmeans, ordinal, joint tests of a predic for (model in list(m, mc)) { # reference 2: emmeans' joint tests of the pairwise contrasts. For the clm - # fit, emmeans warns with "NaNs produced" in a step that does not change - # these values + # 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"), From c264d8dd57ec89506ebe112c7b20022c9e7455df Mon Sep 17 00:00:00 2001 From: Jeffrey Girard Date: Mon, 19 Oct 2026 14:01:47 -0500 Subject: [PATCH 6/8] Regrid only probability grids before ordinal joint tests; test the other modes --- NEWS.md | 9 ++--- R/get_emcontrasts.R | 6 ++-- R/get_jointtest.R | 4 +-- tests/testthat/test-ordinal.R | 67 +++++++++++++++++++++++++++++++++++ 4 files changed, 78 insertions(+), 8 deletions(-) diff --git a/NEWS.md b/NEWS.md index b37607cc3..cbc71e93f 100644 --- a/NEWS.md +++ b/NEWS.md @@ -49,10 +49,11 @@ `estimate_contrasts(m, "Sat", by = "Infl", comparison = "joint", backend = "emmeans")`. Before, the F and Chi-squared statistics were too large, for example by a factor of 4 to 52 for - `polr(Sat ~ Infl + Type, data = MASS::housing)`. Joint tests of - ordinal models also no longer warn with "NaNs produced". The `Chisq` - column of joint tests with `backend = "emmeans"` now prints with two - decimals. + `polr(Sat ~ Infl + Type, weights = Freq, data = MASS::housing)`. Joint + tests of ordinal models also no longer warn with "NaNs produced". The + `Chisq` column of joint tests with `backend = "emmeans"` now prints with + two decimals, and `get_emcontrasts()` with `comparison = "joint"` now + returns a data frame that prints. # modelbased 0.18.0 diff --git a/R/get_emcontrasts.R b/R/get_emcontrasts.R index 5b14d9cdb..64d4f311e 100644 --- a/R/get_emcontrasts.R +++ b/R/get_emcontrasts.R @@ -140,8 +140,10 @@ get_emcontrasts <- function( # Joint tests use the linear functions of the emmeans grid. For contrasts # between the response categories of ordinal models, these gave values # 4 to 52 times the Wald statistic of the contrasts (emmeans 2.0.4). - # After `regrid()`, the joint tests equal the Wald statistic. - if (!is.null(ordinal$mode) && identical(comparison, "joint")) { + # After `regrid()`, the joint tests equal the Wald statistic. Only the + # "prob" mode is regridded: for "linear.predictor", `regrid()` would + # back-transform the grid to cumulative probabilities. + if (identical(ordinal$mode, "prob") && identical(comparison, "joint")) { estimated <- emmeans::regrid(estimated) } # Find by variables diff --git a/R/get_jointtest.R b/R/get_jointtest.R index 8bba93ff1..b7c64a8d2 100644 --- a/R/get_jointtest.R +++ b/R/get_jointtest.R @@ -154,8 +154,8 @@ result$df2 <- insight::format_value(result$df2, protect_integers = TRUE) result$`F.ratio` <- insight::format_value(result$`F.ratio`, protect_integers = TRUE) # "insight::format_table()" uses the attribute "digits" as the default - # number of digits for numeric columns, and "Chisq" is the only one not - # formatted above. For objects of class "summary_emm", `as.data.frame()` + # number of digits for numeric columns that are not formatted above, such + # as "Chisq". For objects of class "summary_emm", `as.data.frame()` # sets it to `getOption("digits")` when it is missing, so we return a data # frame without it class(result) <- "data.frame" diff --git a/tests/testthat/test-ordinal.R b/tests/testthat/test-ordinal.R index 7eabe267a..98cf64915 100644 --- a/tests/testthat/test-ordinal.R +++ b/tests/testthat/test-ordinal.R @@ -1437,6 +1437,7 @@ test_that("estimate_contrasts, backend emmeans, ordinal, joint tests of the resp 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) } }) @@ -1512,3 +1513,69 @@ test_that("estimate_contrasts, backend emmeans, ordinal, printed Chisq of joint } expect_false(any(grepl("114.588", printed, fixed = TRUE))) }) + + +test_that("estimate_contrasts, backend emmeans, ordinal, joint tests in other modes", { + 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 + ))) + out <- suppressWarnings(estimate_contrasts( + model, + "Infl", + comparison = "joint", + predict = mode, + backend = "emmeans" + )) + 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( + m, + "Infl", + comparison = "joint", + predict = "linear.predictor", + backend = "emmeans" + ) + out_latent <- estimate_contrasts( + m, + "Infl", + comparison = "joint", + predict = "latent", + 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))) +}) From 8f355a192f6a5caab85560c3c9fabe8568acbf04 Mon Sep 17 00:00:00 2001 From: Jeffrey Girard Date: Mon, 19 Oct 2026 14:12:02 -0500 Subject: [PATCH 7/8] Narrow the NEWS entry for joint tests; check row counts and warnings in the mode test --- NEWS.md | 12 ++++++------ R/get_emcontrasts.R | 6 ++++-- tests/testthat/test-ordinal.R | 24 +++++++++++++++++------- 3 files changed, 27 insertions(+), 15 deletions(-) diff --git a/NEWS.md b/NEWS.md index cbc71e93f..319f63f05 100644 --- a/NEWS.md +++ b/NEWS.md @@ -47,13 +47,13 @@ the correct statistics for contrasts between the response categories of ordinal models, such as `MASS::polr()` and `ordinal::clm()`. An example is `estimate_contrasts(m, "Sat", by = "Infl", comparison = "joint", - backend = "emmeans")`. Before, the F and Chi-squared statistics were too - large, for example by a factor of 4 to 52 for + backend = "emmeans")`. Before, the F and Chi-squared statistics were wrong. + For example, they were 4 to 52 times too large for `polr(Sat ~ Infl + Type, weights = Freq, data = MASS::housing)`. Joint - tests of ordinal models also no longer warn with "NaNs produced". The - `Chisq` column of joint tests with `backend = "emmeans"` now prints with - two decimals, and `get_emcontrasts()` with `comparison = "joint"` now - returns a data frame that prints. + tests of probabilities (the default `predict = "prob"`) also no longer warn + with "NaNs produced". The `Chisq` column of joint tests with + `backend = "emmeans"` now prints with two decimals, and `get_emcontrasts()` + with `comparison = "joint"` now returns a data frame that prints. # modelbased 0.18.0 diff --git a/R/get_emcontrasts.R b/R/get_emcontrasts.R index 64d4f311e..ff2b77b1d 100644 --- a/R/get_emcontrasts.R +++ b/R/get_emcontrasts.R @@ -138,8 +138,10 @@ get_emcontrasts <- function( fun_args <- c(fun_args, list(specs = specs, at = my_args$emmeans_at)) estimated <- suppressMessages(do.call(emmeans::emmeans, fun_args)) # Joint tests use the linear functions of the emmeans grid. For contrasts - # between the response categories of ordinal models, these gave values - # 4 to 52 times the Wald statistic of the contrasts (emmeans 2.0.4). + # between the response categories of ordinal models, these gave wrong + # values, for example 4 to 52 times the Wald statistic of the contrasts + # for `polr(Sat ~ Infl + Type, weights = Freq, data = MASS::housing)` + # (emmeans 2.0.4). # After `regrid()`, the joint tests equal the Wald statistic. Only the # "prob" mode is regridded: for "linear.predictor", `regrid()` would # back-transform the grid to cumulative probabilities. diff --git a/tests/testthat/test-ordinal.R b/tests/testthat/test-ordinal.R index 98cf64915..72598ef53 100644 --- a/tests/testthat/test-ordinal.R +++ b/tests/testthat/test-ordinal.R @@ -1540,13 +1540,23 @@ test_that("estimate_contrasts, backend emmeans, ordinal, joint tests in other mo ), by = by ))) - out <- suppressWarnings(estimate_contrasts( - model, - "Infl", - comparison = "joint", - predict = mode, - backend = "emmeans" - )) + # 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) } } From 73493a2f9a23a104c6ab807c06160e05f4f724c5 Mon Sep 17 00:00:00 2001 From: Jeffrey Girard Date: Tue, 20 Oct 2026 15:17:22 -0500 Subject: [PATCH 8/8] Restore snapshot files that a local test run removed --- tests/testthat/_snaps/estimate_contrasts.md | 918 ++++++++++++++++++ .../estimate_contrasts_counterfactual.md | 34 + .../_snaps/estimate_contrasts_methods.md | 149 +++ tests/testthat/_snaps/joint_test.md | 14 + 4 files changed, 1115 insertions(+) create mode 100644 tests/testthat/_snaps/estimate_contrasts.md create mode 100644 tests/testthat/_snaps/estimate_contrasts_counterfactual.md create mode 100644 tests/testthat/_snaps/estimate_contrasts_methods.md create mode 100644 tests/testthat/_snaps/joint_test.md diff --git a/tests/testthat/_snaps/estimate_contrasts.md b/tests/testthat/_snaps/estimate_contrasts.md new file mode 100644 index 000000000..736383b01 --- /dev/null +++ b/tests/testthat/_snaps/estimate_contrasts.md @@ -0,0 +1,918 @@ +# estimate_contrasts - Frequentist, Three factors 1 + + Code + print(estimate_contrasts(model, contrast = c("vs", "am"), by = "gear", backend = "marginaleffects"), + zap_small = TRUE, table_width = Inf) + Output + Marginal Contrasts Analysis + + Level1 | Level2 | gear | Difference | SE | 95% CI | t(25) | p + ------------------------------------------------------------------------------------ + vs 0, am 1 | vs 0, am 0 | 3 | 6.98 | 2.33 | [ 2.17, 11.79] | 2.99 | 0.006 + vs 1, am 0 | vs 0, am 0 | 3 | 5.28 | 2.33 | [ 0.47, 10.09] | 2.26 | 0.033 + vs 1, am 1 | vs 0, am 0 | 3 | 12.27 | 3.30 | [ 5.47, 19.07] | 3.71 | 0.001 + vs 1, am 0 | vs 0, am 1 | 3 | -1.70 | 3.30 | [-8.50, 5.10] | -0.51 | 0.611 + vs 1, am 1 | vs 0, am 1 | 3 | 5.28 | 2.33 | [ 0.47, 10.09] | 2.26 | 0.033 + vs 1, am 1 | vs 1, am 0 | 3 | 6.98 | 2.33 | [ 2.17, 11.79] | 2.99 | 0.006 + vs 0, am 1 | vs 0, am 0 | 4 | 6.98 | 2.33 | [ 2.17, 11.79] | 2.99 | 0.006 + vs 1, am 0 | vs 0, am 0 | 4 | 7.03 | 2.95 | [ 0.95, 13.12] | 2.38 | 0.025 + vs 1, am 1 | vs 0, am 0 | 4 | 14.02 | 4.31 | [ 5.15, 22.88] | 3.26 | 0.003 + vs 1, am 0 | vs 0, am 1 | 4 | 0.05 | 3.13 | [-6.40, 6.50] | 0.02 | 0.987 + vs 1, am 1 | vs 0, am 1 | 4 | 7.03 | 2.95 | [ 0.95, 13.12] | 2.38 | 0.025 + vs 1, am 1 | vs 1, am 0 | 4 | 6.98 | 2.33 | [ 2.17, 11.79] | 2.99 | 0.006 + vs 0, am 1 | vs 0, am 0 | 5 | 6.98 | 2.33 | [ 2.17, 11.79] | 2.99 | 0.006 + vs 1, am 0 | vs 0, am 0 | 5 | 11.27 | 4.04 | [ 2.95, 19.60] | 2.79 | 0.010 + vs 1, am 1 | vs 0, am 0 | 5 | 18.26 | 4.67 | [ 8.64, 27.88] | 3.91 | < .001 + vs 1, am 0 | vs 0, am 1 | 5 | 4.29 | 4.67 | [-5.33, 13.91] | 0.92 | 0.367 + vs 1, am 1 | vs 0, am 1 | 5 | 11.27 | 4.04 | [ 2.95, 19.60] | 2.79 | 0.010 + vs 1, am 1 | vs 1, am 0 | 5 | 6.98 | 2.33 | [ 2.17, 11.79] | 2.99 | 0.006 + + Variable predicted: mpg + Predictors contrasted: vs, am + p-values are uncorrected. + +# estimate_contrasts - marginaleffects, comparisons, validate against predict + + Code + print(estimate_contrasts(m, c("time", "coffee"), backend = "marginaleffects", + p_adjust = "none"), zap_small = TRUE, table_width = Inf) + Output + Marginal Contrasts Analysis + + Level1 | Level2 | Difference | SE | 95% CI | t(113) | p + --------------------------------------------------------------------------------------------- + morning, control | morning, coffee | -5.78 | 1.99 | [-9.73, -1.83] | -2.90 | 0.004 + noon, coffee | morning, coffee | -1.93 | 1.99 | [-5.88, 2.02] | -0.97 | 0.336 + noon, control | morning, coffee | 0.00 | 1.99 | [-3.95, 3.95] | 0.00 | > .999 + afternoon, coffee | morning, coffee | 1.93 | 1.99 | [-2.02, 5.88] | 0.97 | 0.336 + afternoon, control | morning, coffee | 0.00 | 1.99 | [-3.95, 3.95] | 0.00 | > .999 + noon, coffee | morning, control | 3.86 | 1.99 | [-0.09, 7.81] | 1.93 | 0.056 + noon, control | morning, control | 5.78 | 1.99 | [ 1.83, 9.73] | 2.90 | 0.004 + afternoon, coffee | morning, control | 7.71 | 1.99 | [ 3.76, 11.66] | 3.87 | < .001 + afternoon, control | morning, control | 5.78 | 1.99 | [ 1.83, 9.73] | 2.90 | 0.004 + noon, control | noon, coffee | 1.93 | 1.99 | [-2.02, 5.88] | 0.97 | 0.336 + afternoon, coffee | noon, coffee | 3.86 | 1.99 | [-0.09, 7.81] | 1.93 | 0.056 + afternoon, control | noon, coffee | 1.93 | 1.99 | [-2.02, 5.88] | 0.97 | 0.336 + afternoon, coffee | noon, control | 1.93 | 1.99 | [-2.02, 5.88] | 0.97 | 0.336 + afternoon, control | noon, control | 0.00 | 1.99 | [-3.95, 3.95] | 0.00 | > .999 + afternoon, control | afternoon, coffee | -1.93 | 1.99 | [-5.88, 2.02] | -0.97 | 0.336 + + Variable predicted: alertness + Predictors contrasted: time, coffee + Predictors averaged: sex + p-values are uncorrected. + +--- + + Code + print(estimate_contrasts(m, c("time", "coffee"), backend = "marginaleffects", + p_adjust = "none", comparison = ratio ~ reference | coffee), zap_small = TRUE, + table_width = Inf) + Output + Marginal Contrasts Analysis + + Level1 | Level2 | coffee | Ratio | SE | 95% CI | t(113) | p + ---------------------------------------------------------------------------- + noon | morning | coffee | 0.89 | 0.11 | [0.67, 1.11] | -1.02 | 0.309 + afternoon | morning | coffee | 1.11 | 0.12 | [0.87, 1.36] | 0.91 | 0.363 + noon | morning | control | 1.51 | 0.22 | [1.06, 1.95] | 2.27 | 0.025 + afternoon | morning | control | 1.51 | 0.22 | [1.06, 1.95] | 2.27 | 0.025 + + Variable predicted: alertness + Predictors contrasted: time + Predictors averaged: sex + p-values are uncorrected. + +--- + + Code + print(estimate_contrasts(m, c("time", "coffee"), backend = "marginaleffects", + p_adjust = "none", comparison = "(b2-b1)=(b4-b3)"), zap_small = TRUE, + table_width = Inf) + Output + Marginal Contrasts Analysis + + Parameter | Difference | SE | 95% CI | t(113) | p + ---------------------------------------------------------------- + b2-b1=b4-b3 | 5.78 | 2.82 | [0.20, 11.37] | 2.05 | 0.043 + + Variable predicted: alertness + Predictors contrasted: time, coffee + Predictors averaged: sex + p-values are uncorrected. + Parameters: + b2 = time [noon], coffee [coffee] + b1 = time [morning], coffee [coffee] + b4 = time [morning], coffee [control] + b3 = time [afternoon], coffee [coffee] + +--- + + Code + print(estimate_contrasts(m, c("time", "coffee"), backend = "marginaleffects", + p_adjust = "none", comparison = "b5=b3"), zap_small = TRUE, table_width = Inf) + Output + Marginal Contrasts Analysis + + Parameter | Difference | SE | 95% CI | t(113) | p + -------------------------------------------------------------- + b5=b3 | -1.93 | 1.99 | [-5.88, 2.02] | -0.97 | 0.336 + + Variable predicted: alertness + Predictors contrasted: time, coffee + Predictors averaged: sex + p-values are uncorrected. + Parameters: + b5 = time [noon], coffee [control] + b3 = time [afternoon], coffee [coffee] + +--- + + Code + estimate_contrasts(fit, c("e42dep", "c172code"), comparison = "b6-b3=0", + backend = "marginaleffects") + Output + Marginal Contrasts Analysis + + Parameter | Difference | SE | 95% CI | t(800) | p + -------------------------------------------------------------- + b6-b3=0 | -0.34 | 0.56 | [-1.44, 0.77] | -0.60 | 0.552 + + Variable predicted: neg_c_7 + Predictors contrasted: e42dep, c172code + Predictors averaged: c12hour (42), barthtot (65), c161sex + p-values are uncorrected. + Parameters: + b6 = e42dep [slightly dependent], c172code [intermediate level of education] + b3 = e42dep [moderately dependent], c172code [low level of education] + +# estimate_contrasts - marginaleffects vs emmeans + + Code + estimate_contrasts(model, backend = "emmeans") + Message + No variable was specified for contrast estimation. Selecting `contrast = + "Species"`. + Output + Marginal Contrasts Analysis + + Level1 | Level2 | Difference | 95% CI | SE | z | p + ----------------------------------------------------------------------------- + versicolor | setosa | 0.68 | [ 0.54, 0.82] | 0.07 | 9.27 | < .001 + virginica | setosa | 0.50 | [ 0.33, 0.67] | 0.08 | 5.90 | < .001 + virginica | versicolor | -0.18 | [-0.35, -0.01] | 0.08 | -2.12 | 0.034 + + Variable predicted: y + Predictors contrasted: Species + p-values are uncorrected. + Contrasts are on the response-scale. + +--- + + Code + estimate_contrasts(model, backend = "marginaleffects") + Message + We selected `contrast=c("Species")`. + Output + Marginal Contrasts Analysis + + Level1 | Level2 | Difference | SE | 95% CI | z | p + ----------------------------------------------------------------------------- + versicolor | setosa | 0.68 | 0.07 | [ 0.54, 0.82] | 9.27 | < .001 + virginica | setosa | 0.50 | 0.08 | [ 0.33, 0.67] | 5.90 | < .001 + virginica | versicolor | -0.18 | 0.08 | [-0.35, -0.01] | -2.12 | 0.034 + + Variable predicted: y + Predictors contrasted: Species + p-values are uncorrected. + Contrasts are on the response-scale. + +--- + + Code + estimate_contrasts(model, backend = "emmeans", p_adjust = "holm") + Message + No variable was specified for contrast estimation. Selecting `contrast = + "Species"`. + Output + Marginal Contrasts Analysis + + Level1 | Level2 | Difference | 95% CI | SE | z | p + ---------------------------------------------------------------------------- + versicolor | setosa | 0.68 | [ 0.50, 0.86] | 0.07 | 9.27 | < .001 + virginica | setosa | 0.50 | [ 0.30, 0.70] | 0.08 | 5.90 | < .001 + virginica | versicolor | -0.18 | [-0.38, 0.02] | 0.08 | -2.12 | 0.034 + + Variable predicted: y + Predictors contrasted: Species + p-value adjustment method: Holm (1979) + Contrasts are on the response-scale. + +--- + + Code + estimate_contrasts(model, backend = "marginaleffects", p_adjust = "holm") + Message + We selected `contrast=c("Species")`. + Output + Marginal Contrasts Analysis + + Level1 | Level2 | Difference | SE | 95% CI | z | p + ----------------------------------------------------------------------------- + versicolor | setosa | 0.68 | 0.07 | [ 0.54, 0.82] | 9.27 | < .001 + virginica | setosa | 0.50 | 0.08 | [ 0.33, 0.67] | 5.90 | < .001 + virginica | versicolor | -0.18 | 0.08 | [-0.35, -0.01] | -2.12 | 0.034 + + Variable predicted: y + Predictors contrasted: Species + p-value adjustment method: Holm (1979) + Contrasts are on the response-scale. + +# estimate_contrasts - filtering works + + Code + print(estimate_contrasts(fit, "c172code", backend = "marginaleffects"), + table_width = Inf, zap_small = TRUE) + Output + Marginal Contrasts Analysis + + Level1 | Level2 | Difference | SE | 95% CI | t(827) | p + -------------------------------------------------------------------- + mid | low | -0.09 | 0.34 | [-0.76, 0.58] | -0.25 | 0.802 + high | low | 0.61 | 0.43 | [-0.24, 1.45] | 1.40 | 0.162 + high | mid | 0.69 | 0.36 | [-0.02, 1.40] | 1.92 | 0.055 + + Variable predicted: neg_c_7 + Predictors contrasted: c172code + Predictors averaged: e16sex, c161sex + p-values are uncorrected. + +--- + + Code + print(estimate_contrasts(fit, c("c161sex", "c172code"), backend = "marginaleffects"), + table_width = Inf, zap_small = TRUE) + Output + Marginal Contrasts Analysis + + Level1 | Level2 | Difference | SE | 95% CI | t(825) | p + ------------------------------------------------------------------------------- + Male, mid | Male, low | 0.29 | 0.71 | [-1.10, 1.68] | 0.41 | 0.684 + Male, high | Male, low | 0.69 | 0.83 | [-0.95, 2.32] | 0.82 | 0.410 + Female, low | Male, low | 1.05 | 0.69 | [-0.31, 2.40] | 1.51 | 0.131 + Female, mid | Male, low | 0.85 | 0.64 | [-0.40, 2.10] | 1.33 | 0.183 + Female, high | Male, low | 1.65 | 0.71 | [ 0.24, 3.05] | 2.30 | 0.022 + Male, high | Male, mid | 0.40 | 0.68 | [-0.94, 1.73] | 0.59 | 0.558 + Female, low | Male, mid | 0.76 | 0.50 | [-0.22, 1.73] | 1.52 | 0.129 + Female, mid | Male, mid | 0.56 | 0.42 | [-0.26, 1.38] | 1.35 | 0.178 + Female, high | Male, mid | 1.36 | 0.53 | [ 0.32, 2.39] | 2.58 | 0.010 + Female, low | Male, high | 0.36 | 0.66 | [-0.95, 1.66] | 0.54 | 0.589 + Female, mid | Male, high | 0.16 | 0.61 | [-1.03, 1.35] | 0.27 | 0.789 + Female, high | Male, high | 0.96 | 0.69 | [-0.39, 2.30] | 1.40 | 0.163 + Female, mid | Female, low | -0.20 | 0.39 | [-0.96, 0.57] | -0.51 | 0.613 + Female, high | Female, low | 0.60 | 0.51 | [-0.39, 1.59] | 1.18 | 0.236 + Female, high | Female, mid | 0.80 | 0.43 | [-0.04, 1.63] | 1.87 | 0.061 + + Variable predicted: neg_c_7 + Predictors contrasted: c161sex, c172code + Predictors averaged: e16sex + p-values are uncorrected. + +--- + + Code + print(estimate_contrasts(fit, "c161sex", "c172code", backend = "marginaleffects"), + table_width = Inf, zap_small = TRUE) + Output + Marginal Contrasts Analysis + + Level1 | Level2 | c172code | Difference | SE | 95% CI | t(825) | p + ------------------------------------------------------------------------------- + Female | Male | low | 1.05 | 0.69 | [-0.31, 2.40] | 1.51 | 0.131 + Female | Male | mid | 0.56 | 0.42 | [-0.26, 1.38] | 1.35 | 0.178 + Female | Male | high | 0.96 | 0.69 | [-0.39, 2.30] | 1.40 | 0.163 + + Variable predicted: neg_c_7 + Predictors contrasted: c161sex + Predictors averaged: e16sex + p-values are uncorrected. + +--- + + Code + print(estimate_slopes(fit, "barthtot", backend = "marginaleffects"), + table_width = Inf, zap_small = TRUE) + Output + Estimated Marginal Effects + + Slope | SE | 95% CI | t(810) | p + ----------------------------------------------- + -0.05 | 0.00 | [-0.06, -0.05] | -12.77 | < .001 + + Marginal effects estimated for barthtot + Type of slope was dY/dX + +--- + + Code + print(estimate_slopes(fit, "barthtot", by = "c172code", backend = "marginaleffects"), + table_width = Inf, zap_small = TRUE) + Output + Estimated Marginal Effects + + c172code | Slope | SE | 95% CI | t(808) | p + ---------------------------------------------------------- + low | -0.06 | 0.01 | [-0.08, -0.05] | -7.08 | < .001 + mid | -0.05 | 0.01 | [-0.06, -0.04] | -9.82 | < .001 + high | -0.05 | 0.01 | [-0.07, -0.03] | -4.53 | < .001 + + Marginal effects estimated for barthtot + Type of slope was dY/dX + +--- + + Code + print(estimate_contrasts(fit, "barthtot", "c172code", backend = "marginaleffects"), + table_width = Inf, zap_small = TRUE) + Output + Marginal Contrasts Analysis + + Level1 | Level2 | Difference | SE | 95% CI | t(808) | p + -------------------------------------------------------------------- + mid | low | 0.01 | 0.01 | [-0.01, 0.03] | 1.17 | 0.243 + high | low | 0.02 | 0.01 | [-0.01, 0.04] | 1.10 | 0.270 + high | mid | 0.00 | 0.01 | [-0.02, 0.03] | 0.27 | 0.786 + + Variable predicted: neg_c_7 + Predictors contrasted: barthtot + Predictors averaged: e16sex, barthtot (65) + p-values are uncorrected. + +--- + + Code + print(estimate_contrasts(fit, "barthtot", c("c172code", "e16sex"), backend = "marginaleffects"), + table_width = Inf, zap_small = TRUE) + Output + Marginal Contrasts Analysis + + Level1 | Level2 | Difference | SE | 95% CI | t(803) | p + ------------------------------------------------------------------------------- + low, female | low, male | -0.01 | 0.02 | [-0.04, 0.03] | -0.35 | 0.729 + mid, male | low, male | -0.01 | 0.02 | [-0.04, 0.02] | -0.46 | 0.648 + mid, female | low, male | 0.02 | 0.01 | [-0.01, 0.05] | 1.24 | 0.217 + high, male | low, male | 0.00 | 0.02 | [-0.05, 0.04] | -0.16 | 0.876 + high, female | low, male | 0.02 | 0.02 | [-0.01, 0.06] | 1.19 | 0.236 + mid, male | low, female | 0.00 | 0.02 | [-0.03, 0.03] | -0.08 | 0.940 + mid, female | low, female | 0.02 | 0.01 | [ 0.00, 0.05] | 1.76 | 0.079 + high, male | low, female | 0.00 | 0.02 | [-0.04, 0.05] | 0.12 | 0.908 + high, female | low, female | 0.03 | 0.02 | [-0.01, 0.06] | 1.58 | 0.115 + mid, female | mid, male | 0.03 | 0.01 | [ 0.00, 0.05] | 2.23 | 0.026 + high, male | mid, male | 0.00 | 0.02 | [-0.04, 0.05] | 0.18 | 0.859 + high, female | mid, male | 0.03 | 0.02 | [ 0.00, 0.06] | 1.83 | 0.068 + high, male | mid, female | -0.02 | 0.02 | [-0.06, 0.02] | -1.08 | 0.280 + high, female | mid, female | 0.00 | 0.01 | [-0.02, 0.03] | 0.26 | 0.792 + high, female | high, male | 0.03 | 0.02 | [-0.02, 0.07] | 1.11 | 0.269 + + Variable predicted: neg_c_7 + Predictors contrasted: barthtot + Predictors averaged: barthtot (65) + p-values are uncorrected. + +# estimate_contrasts - simple contrasts and with - in levels works + + Code + print(estimate_contrasts(model, "Species", backend = "marginaleffects"), + table_width = Inf) + Output + Marginal Contrasts Analysis + + Level1 | Level2 | Difference | SE | 95% CI | t(146) | p + ---------------------------------------------------------------------------- + versicolor | setosa | 1.46 | 0.11 | [1.24, 1.68] | 13.01 | < .001 + virginica | setosa | 1.95 | 0.10 | [1.75, 2.14] | 19.47 | < .001 + virginica | versicolor | 0.49 | 0.09 | [0.31, 0.67] | 5.41 | < .001 + + Variable predicted: Sepal.Length + Predictors contrasted: Species + Predictors averaged: Sepal.Width (3.1) + p-values are uncorrected. + +--- + + Code + print(estimate_contrasts(m, c("time", "coffee"), backend = "marginaleffects"), + zap_small = TRUE, table_width = Inf) + Output + Marginal Contrasts Analysis + + Level1 | Level2 | Difference | SE | 95% CI | t(113) | p + --------------------------------------------------------------------------------------------- + morning, control | morning, coffee | -5.78 | 1.99 | [-9.73, -1.83] | -2.90 | 0.004 + noon, coffee | morning, coffee | -1.93 | 1.99 | [-5.88, 2.02] | -0.97 | 0.336 + noon, control | morning, coffee | 0.00 | 1.99 | [-3.95, 3.95] | 0.00 | > .999 + afternoon, coffee | morning, coffee | 1.93 | 1.99 | [-2.02, 5.88] | 0.97 | 0.336 + afternoon, control | morning, coffee | 0.00 | 1.99 | [-3.95, 3.95] | 0.00 | > .999 + noon, coffee | morning, control | 3.86 | 1.99 | [-0.09, 7.81] | 1.93 | 0.056 + noon, control | morning, control | 5.78 | 1.99 | [ 1.83, 9.73] | 2.90 | 0.004 + afternoon, coffee | morning, control | 7.71 | 1.99 | [ 3.76, 11.66] | 3.87 | < .001 + afternoon, control | morning, control | 5.78 | 1.99 | [ 1.83, 9.73] | 2.90 | 0.004 + noon, control | noon, coffee | 1.93 | 1.99 | [-2.02, 5.88] | 0.97 | 0.336 + afternoon, coffee | noon, coffee | 3.86 | 1.99 | [-0.09, 7.81] | 1.93 | 0.056 + afternoon, control | noon, coffee | 1.93 | 1.99 | [-2.02, 5.88] | 0.97 | 0.336 + afternoon, coffee | noon, control | 1.93 | 1.99 | [-2.02, 5.88] | 0.97 | 0.336 + afternoon, control | noon, control | 0.00 | 1.99 | [-3.95, 3.95] | 0.00 | > .999 + afternoon, control | afternoon, coffee | -1.93 | 1.99 | [-5.88, 2.02] | -0.97 | 0.336 + + Variable predicted: alertness + Predictors contrasted: time, coffee + Predictors averaged: sex + p-values are uncorrected. + +--- + + Code + print(estimate_contrasts(m, contrast = "time", by = "coffee", backend = "marginaleffects"), + zap_small = TRUE, table_width = Inf) + Output + Marginal Contrasts Analysis + + Level1 | Level2 | coffee | Difference | SE | 95% CI | t(113) | p + ----------------------------------------------------------------------------------- + noon | morning | coffee | -1.93 | 1.99 | [-5.88, 2.02] | -0.97 | 0.336 + afternoon | morning | coffee | 1.93 | 1.99 | [-2.02, 5.88] | 0.97 | 0.336 + afternoon | noon | coffee | 3.86 | 1.99 | [-0.09, 7.81] | 1.93 | 0.056 + noon | morning | control | 5.78 | 1.99 | [ 1.83, 9.73] | 2.90 | 0.004 + afternoon | morning | control | 5.78 | 1.99 | [ 1.83, 9.73] | 2.90 | 0.004 + afternoon | noon | control | 0.00 | 1.99 | [-3.95, 3.95] | 0.00 | > .999 + + Variable predicted: alertness + Predictors contrasted: time + Predictors averaged: sex + p-values are uncorrected. + +--- + + Code + print(estimate_contrasts(model, contrast = c("mined", "spp"), backend = "marginaleffects"), + zap_small = TRUE, table_width = Inf) + Output + Marginal Contrasts Analysis + + Level1 | Level2 | Difference | SE | 95% CI | z | p + ----------------------------------------------------------------------------- + yes, PR | yes, GP | 0.06 | 0.05 | [-0.04, 0.16] | 1.10 | 0.271 + yes, DM | yes, GP | 0.32 | 0.11 | [ 0.10, 0.54] | 2.89 | 0.004 + yes, EC-A | yes, GP | 0.04 | 0.05 | [-0.05, 0.13] | 0.80 | 0.421 + yes, EC-L | yes, GP | 0.17 | 0.08 | [ 0.02, 0.32] | 2.19 | 0.028 + yes, DES-L | yes, GP | 0.43 | 0.14 | [ 0.17, 0.70] | 3.19 | 0.001 + yes, DF | yes, GP | 0.43 | 0.14 | [ 0.17, 0.70] | 3.19 | 0.001 + no, GP | yes, GP | 2.59 | 0.55 | [ 1.50, 3.67] | 4.67 | < .001 + no, PR | yes, GP | 0.51 | 0.16 | [ 0.20, 0.82] | 3.21 | 0.001 + no, DM | yes, GP | 2.86 | 0.61 | [ 1.68, 4.05] | 4.73 | < .001 + no, EC-A | yes, GP | 1.10 | 0.27 | [ 0.57, 1.64] | 4.03 | < .001 + no, EC-L | yes, GP | 4.67 | 0.94 | [ 2.82, 6.52] | 4.95 | < .001 + no, DES-L | yes, GP | 4.62 | 0.93 | [ 2.79, 6.45] | 4.95 | < .001 + no, DF | yes, GP | 2.24 | 0.49 | [ 1.28, 3.20] | 4.58 | < .001 + yes, DM | yes, PR | 0.26 | 0.11 | [ 0.05, 0.48] | 2.43 | 0.015 + yes, EC-A | yes, PR | -0.02 | 0.06 | [-0.13, 0.09] | -0.33 | 0.740 + yes, EC-L | yes, PR | 0.11 | 0.08 | [-0.04, 0.27] | 1.43 | 0.154 + yes, DES-L | yes, PR | 0.38 | 0.13 | [ 0.12, 0.64] | 2.86 | 0.004 + yes, DF | yes, PR | 0.38 | 0.13 | [ 0.12, 0.64] | 2.86 | 0.004 + no, GP | yes, PR | 2.53 | 0.56 | [ 1.44, 3.62] | 4.54 | < .001 + no, PR | yes, PR | 0.45 | 0.16 | [ 0.13, 0.77] | 2.75 | 0.006 + no, DM | yes, PR | 2.80 | 0.61 | [ 1.61, 4.00] | 4.61 | < .001 + no, EC-A | yes, PR | 1.05 | 0.28 | [ 0.50, 1.59] | 3.76 | < .001 + no, EC-L | yes, PR | 4.61 | 0.95 | [ 2.76, 6.47] | 4.87 | < .001 + no, DES-L | yes, PR | 4.56 | 0.94 | [ 2.73, 6.40] | 4.87 | < .001 + no, DF | yes, PR | 2.18 | 0.49 | [ 1.22, 3.15] | 4.44 | < .001 + yes, EC-A | yes, DM | -0.28 | 0.11 | [-0.50, -0.07] | -2.59 | 0.010 + yes, EC-L | yes, DM | -0.15 | 0.11 | [-0.36, 0.06] | -1.39 | 0.164 + yes, DES-L | yes, DM | 0.11 | 0.13 | [-0.14, 0.36] | 0.89 | 0.375 + yes, DF | yes, DM | 0.11 | 0.13 | [-0.14, 0.36] | 0.89 | 0.375 + no, GP | yes, DM | 2.27 | 0.58 | [ 1.14, 3.40] | 3.93 | < .001 + no, PR | yes, DM | 0.19 | 0.20 | [-0.21, 0.58] | 0.93 | 0.352 + no, DM | yes, DM | 2.54 | 0.63 | [ 1.31, 3.77] | 4.04 | < .001 + no, EC-A | yes, DM | 0.78 | 0.31 | [ 0.18, 1.38] | 2.56 | 0.011 + no, EC-L | yes, DM | 4.35 | 0.96 | [ 2.46, 6.24] | 4.51 | < .001 + no, DES-L | yes, DM | 4.30 | 0.95 | [ 2.43, 6.17] | 4.51 | < .001 + no, DF | yes, DM | 1.92 | 0.51 | [ 0.91, 2.93] | 3.74 | < .001 + yes, EC-L | yes, EC-A | 0.13 | 0.08 | [-0.02, 0.29] | 1.68 | 0.093 + yes, DES-L | yes, EC-A | 0.40 | 0.13 | [ 0.14, 0.66] | 2.97 | 0.003 + yes, DF | yes, EC-A | 0.40 | 0.13 | [ 0.14, 0.66] | 2.97 | 0.003 + no, GP | yes, EC-A | 2.55 | 0.56 | [ 1.46, 3.64] | 4.58 | < .001 + no, PR | yes, EC-A | 0.47 | 0.16 | [ 0.15, 0.79] | 2.90 | 0.004 + no, DM | yes, EC-A | 2.82 | 0.61 | [ 1.63, 4.01] | 4.65 | < .001 + no, EC-A | yes, EC-A | 1.06 | 0.28 | [ 0.52, 1.61] | 3.85 | < .001 + no, EC-L | yes, EC-A | 4.63 | 0.95 | [ 2.78, 6.48] | 4.90 | < .001 + no, DES-L | yes, EC-A | 4.58 | 0.94 | [ 2.75, 6.42] | 4.90 | < .001 + no, DF | yes, EC-A | 2.20 | 0.49 | [ 1.24, 3.17] | 4.49 | < .001 + yes, DES-L | yes, EC-L | 0.26 | 0.13 | [ 0.02, 0.51] | 2.08 | 0.037 + yes, DF | yes, EC-L | 0.26 | 0.13 | [ 0.02, 0.51] | 2.08 | 0.037 + no, GP | yes, EC-L | 2.42 | 0.57 | [ 1.31, 3.53] | 4.28 | < .001 + no, PR | yes, EC-L | 0.34 | 0.18 | [-0.01, 0.69] | 1.89 | 0.058 + no, DM | yes, EC-L | 2.69 | 0.62 | [ 1.48, 3.90] | 4.37 | < .001 + no, EC-A | yes, EC-L | 0.93 | 0.29 | [ 0.37, 1.50] | 3.23 | 0.001 + no, EC-L | yes, EC-L | 4.50 | 0.95 | [ 2.63, 6.37] | 4.72 | < .001 + no, DES-L | yes, EC-L | 4.45 | 0.94 | [ 2.60, 6.30] | 4.71 | < .001 + no, DF | yes, EC-L | 2.07 | 0.50 | [ 1.09, 3.05] | 4.14 | < .001 + yes, DF | yes, DES-L | 0.00 | 0.13 | [-0.26, 0.26] | 0.00 | > .999 + no, GP | yes, DES-L | 2.15 | 0.59 | [ 1.00, 3.31] | 3.67 | < .001 + no, PR | yes, DES-L | 0.07 | 0.22 | [-0.36, 0.50] | 0.33 | 0.738 + no, DM | yes, DES-L | 2.43 | 0.64 | [ 1.18, 3.68] | 3.81 | < .001 + no, EC-A | yes, DES-L | 0.67 | 0.32 | [ 0.04, 1.29] | 2.09 | 0.037 + no, EC-L | yes, DES-L | 4.24 | 0.97 | [ 2.33, 6.14] | 4.36 | < .001 + no, DES-L | yes, DES-L | 4.19 | 0.96 | [ 2.30, 6.07] | 4.35 | < .001 + no, DF | yes, DES-L | 1.81 | 0.52 | [ 0.78, 2.83] | 3.45 | < .001 + no, GP | yes, DF | 2.15 | 0.59 | [ 1.00, 3.31] | 3.67 | < .001 + no, PR | yes, DF | 0.07 | 0.22 | [-0.36, 0.50] | 0.33 | 0.738 + no, DM | yes, DF | 2.43 | 0.64 | [ 1.18, 3.68] | 3.81 | < .001 + no, EC-A | yes, DF | 0.67 | 0.32 | [ 0.04, 1.29] | 2.09 | 0.037 + no, EC-L | yes, DF | 4.24 | 0.97 | [ 2.33, 6.14] | 4.36 | < .001 + no, DES-L | yes, DF | 4.19 | 0.96 | [ 2.30, 6.07] | 4.35 | < .001 + no, DF | yes, DF | 1.81 | 0.52 | [ 0.78, 2.83] | 3.45 | < .001 + no, PR | no, GP | -2.08 | 0.48 | [-3.02, -1.14] | -4.35 | < .001 + no, DM | no, GP | 0.27 | 0.37 | [-0.46, 1.00] | 0.73 | 0.466 + no, EC-A | no, GP | -1.49 | 0.41 | [-2.29, -0.68] | -3.61 | < .001 + no, EC-L | no, GP | 2.08 | 0.58 | [ 0.95, 3.21] | 3.61 | < .001 + no, DES-L | no, GP | 2.03 | 0.57 | [ 0.92, 3.15] | 3.57 | < .001 + no, DF | no, GP | -0.35 | 0.35 | [-1.04, 0.35] | -0.98 | 0.328 + no, DM | no, PR | 2.35 | 0.53 | [ 1.32, 3.39] | 4.47 | < .001 + no, EC-A | no, PR | 0.59 | 0.23 | [ 0.14, 1.05] | 2.56 | 0.011 + no, EC-L | no, PR | 4.16 | 0.86 | [ 2.49, 5.84] | 4.87 | < .001 + no, DES-L | no, PR | 4.11 | 0.85 | [ 2.45, 5.77] | 4.86 | < .001 + no, DF | no, PR | 1.73 | 0.42 | [ 0.92, 2.55] | 4.15 | < .001 + no, EC-A | no, DM | -1.76 | 0.46 | [-2.65, -0.87] | -3.86 | < .001 + no, EC-L | no, DM | 1.81 | 0.55 | [ 0.73, 2.89] | 3.29 | < .001 + no, DES-L | no, DM | 1.76 | 0.54 | [ 0.70, 2.82] | 3.24 | 0.001 + no, DF | no, DM | -0.62 | 0.38 | [-1.36, 0.12] | -1.65 | 0.100 + no, EC-L | no, EC-A | 3.57 | 0.77 | [ 2.07, 5.07] | 4.66 | < .001 + no, DES-L | no, EC-A | 3.52 | 0.76 | [ 2.03, 5.00] | 4.65 | < .001 + no, DF | no, EC-A | 1.14 | 0.36 | [ 0.43, 1.85] | 3.16 | 0.002 + no, DES-L | no, EC-L | -0.05 | 0.48 | [-0.99, 0.89] | -0.10 | 0.918 + no, DF | no, EC-L | -2.43 | 0.61 | [-3.63, -1.22] | -3.95 | < .001 + no, DF | no, DES-L | -2.38 | 0.61 | [-3.57, -1.19] | -3.92 | < .001 + + Variable predicted: count + Predictors contrasted: mined, spp + Predictors averaged: cover (8.7e-11), site + p-values are uncorrected. + Contrasts are on the response-scale. + +--- + + Code + print(estimate_contrasts(model, contrast = "mined", by = "spp", backend = "marginaleffects"), + zap_small = TRUE, table_width = Inf) + Output + Marginal Contrasts Analysis + + Level1 | Level2 | spp | Difference | SE | 95% CI | z | p + -------------------------------------------------------------------------- + no | yes | GP | 2.59 | 0.55 | [1.50, 3.67] | 4.67 | < .001 + no | yes | PR | 0.45 | 0.16 | [0.13, 0.77] | 2.75 | 0.006 + no | yes | DM | 2.54 | 0.63 | [1.31, 3.77] | 4.04 | < .001 + no | yes | EC-A | 1.06 | 0.28 | [0.52, 1.61] | 3.85 | < .001 + no | yes | EC-L | 4.50 | 0.95 | [2.63, 6.37] | 4.72 | < .001 + no | yes | DES-L | 4.19 | 0.96 | [2.30, 6.07] | 4.35 | < .001 + no | yes | DF | 1.81 | 0.52 | [0.78, 2.83] | 3.45 | < .001 + + Variable predicted: count + Predictors contrasted: mined + Predictors averaged: cover (8.7e-11), site + p-values are uncorrected. + Contrasts are on the response-scale. + +# estimate_contrasts - examples from docs work as intendec + + Code + print(estimate_contrasts(model, contrast = "Petal.Width", by = "Species"), + zap_small = TRUE, table_width = Inf) + Output + Marginal Contrasts Analysis + + Level1 | Level2 | Difference | SE | 95% CI | t(144) | p + ---------------------------------------------------------------------------- + versicolor | setosa | 0.22 | 0.46 | [-0.70, 1.13] | 0.47 | 0.640 + virginica | setosa | -0.21 | 0.44 | [-1.07, 0.66] | -0.47 | 0.638 + virginica | versicolor | -0.42 | 0.27 | [-0.95, 0.11] | -1.58 | 0.116 + + Variable predicted: Sepal.Width + Predictors contrasted: Petal.Width + Predictors averaged: Petal.Width (1.2) + p-values are uncorrected. + +--- + + Code + print(estimate_contrasts(model, contrast = c("Species", "Petal.Width"), length = 2), + zap_small = TRUE, table_width = Inf) + Output + Marginal Contrasts Analysis + + Level1 | Level2 | Difference | SE | 95% CI | t(144) | p + ---------------------------------------------------------------------------------------- + setosa, 2.5 | setosa, 0.1 | 2.01 | 0.98 | [ 0.08, 3.94] | 2.06 | 0.041 + versicolor, 0.1 | setosa, 0.1 | -1.83 | 0.28 | [-2.38, -1.28] | -6.55 | < .001 + versicolor, 2.5 | setosa, 0.1 | 0.70 | 0.27 | [ 0.17, 1.23] | 2.61 | 0.010 + virginica, 0.1 | setosa, 0.1 | -1.55 | 0.31 | [-2.17, -0.93] | -4.95 | < .001 + virginica, 2.5 | setosa, 0.1 | -0.03 | 0.11 | [-0.25, 0.19] | -0.29 | 0.773 + versicolor, 0.1 | setosa, 2.5 | -3.84 | 0.96 | [-5.73, -1.95] | -4.01 | < .001 + versicolor, 2.5 | setosa, 2.5 | -1.31 | 0.95 | [-3.19, 0.58] | -1.37 | 0.172 + virginica, 0.1 | setosa, 2.5 | -3.56 | 0.97 | [-5.47, -1.65] | -3.68 | < .001 + virginica, 2.5 | setosa, 2.5 | -2.04 | 0.92 | [-3.86, -0.22] | -2.21 | 0.028 + versicolor, 2.5 | versicolor, 0.1 | 2.53 | 0.52 | [ 1.50, 3.56] | 4.86 | < .001 + virginica, 0.1 | versicolor, 0.1 | 0.28 | 0.41 | [-0.52, 1.08] | 0.69 | 0.492 + virginica, 2.5 | versicolor, 0.1 | 1.80 | 0.28 | [ 1.24, 2.35] | 6.35 | < .001 + virginica, 0.1 | versicolor, 2.5 | -2.25 | 0.40 | [-3.04, -1.46] | -5.64 | < .001 + virginica, 2.5 | versicolor, 2.5 | -0.73 | 0.27 | [-1.27, -0.20] | -2.70 | 0.008 + virginica, 2.5 | virginica, 0.1 | 1.52 | 0.37 | [ 0.77, 2.26] | 4.04 | < .001 + + Variable predicted: Sepal.Width + Predictors contrasted: Species, Petal.Width + p-values are uncorrected. + +--- + + Code + print(estimate_contrasts(model, contrast = c("Species", "Petal.Width=c(1, 2)")), + zap_small = TRUE, table_width = Inf) + Output + Marginal Contrasts Analysis + + Level1 | Level2 | Difference | SE | 95% CI | t(144) | p + ------------------------------------------------------------------------------------ + setosa, 2 | setosa, 1 | 0.84 | 0.41 | [ 0.03, 1.64] | 2.06 | 0.041 + versicolor, 1 | setosa, 1 | -1.63 | 0.32 | [-2.27, -1.00] | -5.09 | < .001 + versicolor, 2 | setosa, 1 | -0.58 | 0.35 | [-1.26, 0.10] | -1.68 | 0.096 + virginica, 1 | setosa, 1 | -1.73 | 0.35 | [-2.43, -1.04] | -4.93 | < .001 + virginica, 2 | setosa, 1 | -1.10 | 0.31 | [-1.72, -0.48] | -3.52 | < .001 + versicolor, 1 | setosa, 2 | -2.47 | 0.72 | [-3.89, -1.05] | -3.43 | < .001 + versicolor, 2 | setosa, 2 | -1.42 | 0.73 | [-2.86, 0.03] | -1.94 | 0.055 + virginica, 1 | setosa, 2 | -2.57 | 0.73 | [-4.02, -1.12] | -3.50 | < .001 + virginica, 2 | setosa, 2 | -1.94 | 0.72 | [-3.35, -0.52] | -2.71 | 0.008 + versicolor, 2 | versicolor, 1 | 1.05 | 0.22 | [ 0.62, 1.48] | 4.86 | < .001 + virginica, 1 | versicolor, 1 | -0.10 | 0.19 | [-0.47, 0.27] | -0.54 | 0.589 + virginica, 2 | versicolor, 1 | 0.53 | 0.09 | [ 0.35, 0.71] | 5.72 | < .001 + virginica, 1 | versicolor, 2 | -1.15 | 0.23 | [-1.60, -0.71] | -5.13 | < .001 + virginica, 2 | versicolor, 2 | -0.52 | 0.16 | [-0.84, -0.21] | -3.31 | 0.001 + virginica, 2 | virginica, 1 | 0.63 | 0.16 | [ 0.32, 0.94] | 4.04 | < .001 + + Variable predicted: Sepal.Width + Predictors contrasted: Species, Petal.Width=c(1, 2) + p-values are uncorrected. + +--- + + Code + print(estimate_contrasts(model, by = "Petal.Width", length = 4), zap_small = TRUE, + table_width = Inf) + Message + We selected `contrast=c("Species")`. + Output + Marginal Contrasts Analysis + + Level1 | Level2 | Petal.Width | Difference | SE | 95% CI | t(144) | p + -------------------------------------------------------------------------------------------- + versicolor | setosa | 0.10 | -1.83 | 0.28 | [-2.38, -1.28] | -6.55 | < .001 + virginica | setosa | 0.10 | -1.55 | 0.31 | [-2.17, -0.93] | -4.95 | < .001 + virginica | versicolor | 0.10 | 0.28 | 0.41 | [-0.52, 1.08] | 0.69 | 0.492 + versicolor | setosa | 0.90 | -1.65 | 0.29 | [-2.22, -1.08] | -5.74 | < .001 + virginica | setosa | 0.90 | -1.71 | 0.32 | [-2.35, -1.07] | -5.28 | < .001 + virginica | versicolor | 0.90 | -0.06 | 0.21 | [-0.47, 0.35] | -0.28 | 0.780 + versicolor | setosa | 1.70 | -1.48 | 0.60 | [-2.67, -0.29] | -2.47 | 0.015 + virginica | setosa | 1.70 | -1.88 | 0.60 | [-3.06, -0.70] | -3.14 | 0.002 + virginica | versicolor | 1.70 | -0.40 | 0.11 | [-0.62, -0.17] | -3.50 | < .001 + versicolor | setosa | 2.50 | -1.31 | 0.95 | [-3.19, 0.58] | -1.37 | 0.172 + virginica | setosa | 2.50 | -2.04 | 0.92 | [-3.86, -0.22] | -2.21 | 0.028 + virginica | versicolor | 2.50 | -0.73 | 0.27 | [-1.27, -0.20] | -2.70 | 0.008 + + Variable predicted: Sepal.Width + Predictors contrasted: Species + p-values are uncorrected. + +# estimate_contrasts - test all combinations of contrast and by, with filtering + + Code + print(estimate_contrasts(model2, c("grp", "time", "x")), zap_small = TRUE, + table_width = Inf) + Output + Marginal Contrasts Analysis + + Level1 | Level2 | Difference | SE | 95% CI | t(992) | p + -------------------------------------------------------------------------------------- + control, 1, b | control, 1, a | 0.15 | 0.13 | [-0.10, 0.40] | 1.20 | 0.231 + control, 2, a | control, 1, a | 0.03 | 0.13 | [-0.23, 0.28] | 0.19 | 0.847 + control, 2, b | control, 1, a | 0.09 | 0.13 | [-0.16, 0.34] | 0.73 | 0.465 + treatment, 1, a | control, 1, a | 0.19 | 0.13 | [-0.06, 0.45] | 1.49 | 0.137 + treatment, 1, b | control, 1, a | -0.03 | 0.13 | [-0.28, 0.22] | -0.22 | 0.824 + treatment, 2, a | control, 1, a | 0.02 | 0.13 | [-0.23, 0.27] | 0.19 | 0.850 + treatment, 2, b | control, 1, a | 0.05 | 0.12 | [-0.19, 0.28] | 0.37 | 0.712 + control, 2, a | control, 1, b | -0.12 | 0.13 | [-0.38, 0.13] | -0.97 | 0.332 + control, 2, b | control, 1, b | -0.06 | 0.12 | [-0.30, 0.19] | -0.46 | 0.648 + treatment, 1, a | control, 1, b | 0.04 | 0.13 | [-0.20, 0.29] | 0.34 | 0.730 + treatment, 1, b | control, 1, b | -0.18 | 0.12 | [-0.42, 0.06] | -1.45 | 0.147 + treatment, 2, a | control, 1, b | -0.13 | 0.12 | [-0.37, 0.12] | -1.02 | 0.307 + treatment, 2, b | control, 1, b | -0.10 | 0.12 | [-0.34, 0.13] | -0.89 | 0.376 + control, 2, b | control, 2, a | 0.07 | 0.13 | [-0.19, 0.32] | 0.52 | 0.604 + treatment, 1, a | control, 2, a | 0.17 | 0.13 | [-0.09, 0.43] | 1.26 | 0.207 + treatment, 1, b | control, 2, a | -0.05 | 0.13 | [-0.31, 0.20] | -0.41 | 0.680 + treatment, 2, a | control, 2, a | 0.00 | 0.13 | [-0.26, 0.25] | -0.01 | 0.991 + treatment, 2, b | control, 2, a | 0.02 | 0.13 | [-0.23, 0.27] | 0.16 | 0.876 + treatment, 1, a | control, 2, b | 0.10 | 0.13 | [-0.15, 0.35] | 0.78 | 0.437 + treatment, 1, b | control, 2, b | -0.12 | 0.13 | [-0.37, 0.12] | -0.97 | 0.333 + treatment, 2, a | control, 2, b | -0.07 | 0.13 | [-0.32, 0.18] | -0.55 | 0.582 + treatment, 2, b | control, 2, b | -0.05 | 0.12 | [-0.29, 0.19] | -0.40 | 0.691 + treatment, 1, b | treatment, 1, a | -0.22 | 0.13 | [-0.47, 0.03] | -1.73 | 0.083 + treatment, 2, a | treatment, 1, a | -0.17 | 0.13 | [-0.42, 0.08] | -1.32 | 0.187 + treatment, 2, b | treatment, 1, a | -0.15 | 0.12 | [-0.39, 0.09] | -1.20 | 0.230 + treatment, 2, a | treatment, 1, b | 0.05 | 0.12 | [-0.19, 0.30] | 0.42 | 0.676 + treatment, 2, b | treatment, 1, b | 0.07 | 0.12 | [-0.16, 0.31] | 0.61 | 0.541 + treatment, 2, b | treatment, 2, a | 0.02 | 0.12 | [-0.21, 0.26] | 0.18 | 0.861 + + Variable predicted: score + Predictors contrasted: grp, time, x + p-values are uncorrected. + +--- + + Code + print(estimate_contrasts(model2, c("grp", "time"), by = "x"), zap_small = TRUE, + table_width = Inf) + Output + Marginal Contrasts Analysis + + Level1 | Level2 | x | Difference | SE | 95% CI | t(992) | p + ------------------------------------------------------------------------------------ + control, 2 | control, 1 | a | 0.03 | 0.13 | [-0.23, 0.28] | 0.19 | 0.847 + treatment, 1 | control, 1 | a | 0.19 | 0.13 | [-0.06, 0.45] | 1.49 | 0.137 + treatment, 2 | control, 1 | a | 0.02 | 0.13 | [-0.23, 0.27] | 0.19 | 0.850 + treatment, 1 | control, 2 | a | 0.17 | 0.13 | [-0.09, 0.43] | 1.26 | 0.207 + treatment, 2 | control, 2 | a | 0.00 | 0.13 | [-0.26, 0.25] | -0.01 | 0.991 + treatment, 2 | treatment, 1 | a | -0.17 | 0.13 | [-0.42, 0.08] | -1.32 | 0.187 + control, 2 | control, 1 | b | -0.06 | 0.12 | [-0.30, 0.19] | -0.46 | 0.648 + treatment, 1 | control, 1 | b | -0.18 | 0.12 | [-0.42, 0.06] | -1.45 | 0.147 + treatment, 2 | control, 1 | b | -0.10 | 0.12 | [-0.34, 0.13] | -0.89 | 0.376 + treatment, 1 | control, 2 | b | -0.12 | 0.13 | [-0.37, 0.12] | -0.97 | 0.333 + treatment, 2 | control, 2 | b | -0.05 | 0.12 | [-0.29, 0.19] | -0.40 | 0.691 + treatment, 2 | treatment, 1 | b | 0.07 | 0.12 | [-0.16, 0.31] | 0.61 | 0.541 + + Variable predicted: score + Predictors contrasted: grp, time + p-values are uncorrected. + +--- + + Code + print(estimate_contrasts(model2, "grp", by = c("time", "x")), zap_small = TRUE, + table_width = Inf) + Output + Marginal Contrasts Analysis + + Level1 | Level2 | time | x | Difference | SE | 95% CI | t(992) | p + ----------------------------------------------------------------------------------- + treatment | control | 1 | a | 0.19 | 0.13 | [-0.06, 0.45] | 1.49 | 0.137 + treatment | control | 2 | a | 0.00 | 0.13 | [-0.26, 0.25] | -0.01 | 0.991 + treatment | control | 1 | b | -0.18 | 0.12 | [-0.42, 0.06] | -1.45 | 0.147 + treatment | control | 2 | b | -0.05 | 0.12 | [-0.29, 0.19] | -0.40 | 0.691 + + Variable predicted: score + Predictors contrasted: grp + p-values are uncorrected. + +--- + + Code + print(estimate_contrasts(model2, "grp", by = "time"), zap_small = TRUE, + table_width = Inf) + Output + Marginal Contrasts Analysis + + Level1 | Level2 | time | Difference | SE | 95% CI | t(992) | p + ------------------------------------------------------------------------------- + treatment | control | 1 | 0.01 | 0.09 | [-0.17, 0.18] | 0.09 | 0.931 + treatment | control | 2 | -0.02 | 0.09 | [-0.20, 0.15] | -0.28 | 0.780 + + Variable predicted: score + Predictors contrasted: grp + Predictors averaged: x + p-values are uncorrected. + +--- + + Code + print(estimate_contrasts(model2, c("grp", "time", "x='a'")), zap_small = TRUE, + table_width = Inf) + Output + Marginal Contrasts Analysis + + Level1 | Level2 | Difference | SE | 95% CI | t(992) | p + -------------------------------------------------------------------------------- + control, 2 | control, 1 | 0.03 | 0.13 | [-0.23, 0.28] | 0.19 | 0.847 + treatment, 1 | control, 1 | 0.19 | 0.13 | [-0.06, 0.45] | 1.49 | 0.137 + treatment, 2 | control, 1 | 0.02 | 0.13 | [-0.23, 0.27] | 0.19 | 0.850 + treatment, 1 | control, 2 | 0.17 | 0.13 | [-0.09, 0.43] | 1.26 | 0.207 + treatment, 2 | control, 2 | 0.00 | 0.13 | [-0.26, 0.25] | -0.01 | 0.991 + treatment, 2 | treatment, 1 | -0.17 | 0.13 | [-0.42, 0.08] | -1.32 | 0.187 + + Variable predicted: score + Predictors contrasted: grp, time, x='a' + p-values are uncorrected. + +--- + + Code + print(estimate_contrasts(model2, c("grp", "time=1"), by = "x"), zap_small = TRUE, + table_width = Inf) + Output + Marginal Contrasts Analysis + + Level1 | Level2 | x | Difference | SE | 95% CI | t(992) | p + ---------------------------------------------------------------------------- + treatment | control | a | 0.19 | 0.13 | [-0.06, 0.45] | 1.49 | 0.137 + treatment | control | b | -0.18 | 0.12 | [-0.42, 0.06] | -1.45 | 0.147 + + Variable predicted: score + Predictors contrasted: grp, time=1 + p-values are uncorrected. + +--- + + Code + print(estimate_contrasts(model2, "grp", by = c("time", "x='a'")), zap_small = TRUE, + table_width = Inf) + Output + Marginal Contrasts Analysis + + Level1 | Level2 | time | x | Difference | SE | 95% CI | t(992) | p + ----------------------------------------------------------------------------------- + treatment | control | 1 | a | 0.19 | 0.13 | [-0.06, 0.45] | 1.49 | 0.137 + treatment | control | 2 | a | 0.00 | 0.13 | [-0.26, 0.25] | -0.01 | 0.991 + + Variable predicted: score + Predictors contrasted: grp + p-values are uncorrected. + +--- + + Code + print(estimate_contrasts(model2, "time=c(1,2)", by = "grp"), zap_small = TRUE, + table_width = Inf) + Output + Marginal Contrasts Analysis + + Level1 | Level2 | grp | Difference | SE | 95% CI | t(994) | p + -------------------------------------------------------------------------------- + 2 | 1 | control | 0.14 | 0.11 | [-0.08, 0.36] | 1.24 | 0.216 + 2 | 1 | treatment | -0.07 | 0.11 | [-0.28, 0.14] | -0.63 | 0.529 + + Variable predicted: score + Predictors contrasted: time=c(1,2) + p-values are uncorrected. + +--- + + Code + print(estimate_contrasts(model2, c("grp", "time=2")), zap_small = TRUE, + table_width = Inf) + Output + Marginal Contrasts Analysis + + Level1 | Level2 | Difference | SE | 95% CI | t(994) | p + ------------------------------------------------------------------------ + treatment | control | -0.15 | 0.11 | [-0.36, 0.05] | -1.47 | 0.143 + + Variable predicted: score + Predictors contrasted: grp, time=2 + p-values are uncorrected. + +# estimate_contrast, slopes with emmeans-2 + + Code + print(out, table_width = Inf) + Output + Marginal Contrasts Analysis + + Level1 | Level2 | Difference | 95% CI | p + ---------------------------------------------------------- + 0, 0.725 | 0, -2.012 | 1.12 | [0.37, 3.43] | 0.839 + 0, 3.463 | 0, -2.012 | 1.26 | [0.14, 11.76] | 0.839 + 1, -2.012 | 0, -2.012 | 1.00 | [0.35, 2.82] | 0.999 + 1, 0.725 | 0, -2.012 | 1.12 | [0.23, 5.43] | 0.887 + 1, 3.463 | 0, -2.012 | 1.26 | [0.10, 15.76] | 0.858 + 0, 3.463 | 0, 0.725 | 1.12 | [0.37, 3.43] | 0.839 + 1, -2.012 | 0, 0.725 | 0.89 | [0.20, 3.87] | 0.877 + 1, 0.725 | 0, 0.725 | 1.00 | [0.35, 2.82] | 0.999 + 1, 3.463 | 0, 0.725 | 1.12 | [0.23, 5.43] | 0.887 + 1, -2.012 | 0, 3.463 | 0.79 | [0.07, 8.71] | 0.849 + 1, 0.725 | 0, 3.463 | 0.89 | [0.20, 3.87] | 0.877 + 1, 3.463 | 0, 3.463 | 1.00 | [0.35, 2.82] | 0.999 + 1, 0.725 | 1, -2.012 | 1.12 | [0.37, 3.43] | 0.839 + 1, 3.463 | 1, -2.012 | 1.26 | [0.14, 11.76] | 0.839 + 1, 3.463 | 1, 0.725 | 1.12 | [0.37, 3.43] | 0.839 + + Variable predicted: outcome + Predictors contrasted: var_binom, var_cont + p-values are uncorrected. + Contrasts are on the exp-transformed link-scale. + diff --git a/tests/testthat/_snaps/estimate_contrasts_counterfactual.md b/tests/testthat/_snaps/estimate_contrasts_counterfactual.md new file mode 100644 index 000000000..3bee94140 --- /dev/null +++ b/tests/testthat/_snaps/estimate_contrasts_counterfactual.md @@ -0,0 +1,34 @@ +# estimate_contrast, counterfactual, snapshots, Level-columns + + Code + print(out, table_width = Inf) + Output + Counterfactual Contrasts Analysis (G-computation) + + Level1 | Level2 | Difference | SE | 95% CI | z | p + ------------------------------------------------------------------ + tx | control | 5.65 | 1.95 | [1.84, 9.47] | 2.91 | 0.004 + + Variable predicted: QoL + Predictors contrasted: treatment + Predictors averaged: time, education, hospital (0.95), age (0.22), phq4 (-0.076), ID + p-values are uncorrected. + +--- + + Code + print(out, table_width = Inf) + Output + Counterfactual Contrasts Analysis (G-computation) + + Level1 | Level2 | time | Difference | SE | 95% CI | z | p + ----------------------------------------------------------------------------- + tx | control | pre | -1.40 | 2.26 | [-5.82, 3.03] | -0.62 | 0.535 + tx | control | post | 3.32 | 2.26 | [-1.10, 7.75] | 1.47 | 0.141 + tx | control | end | 15.00 | 2.26 | [10.57, 19.42] | 6.64 | < .001 + + Variable predicted: QoL + Predictors contrasted: treatment + Predictors averaged: education, hospital (0.95), age (0.22), phq4 (-0.076), ID + p-values are uncorrected. + diff --git a/tests/testthat/_snaps/estimate_contrasts_methods.md b/tests/testthat/_snaps/estimate_contrasts_methods.md new file mode 100644 index 000000000..75ab54781 --- /dev/null +++ b/tests/testthat/_snaps/estimate_contrasts_methods.md @@ -0,0 +1,149 @@ +# estimate_contrasts - Random Effects Levels, pairwise + + Code + print(estimate_contrasts(estim, contrast = c("gender", "employed", "age")), + zap_small = TRUE, table_width = Inf) + Output + Model-based Contrasts Analysis + + Level1 | Level2 | Difference | SE | 95% CI | Statistic | p + ------------------------------------------------------------------------------------------------- + Male, no, -40 | Female, no, -40 | 0.07 | 1.11 | [-2.11, 2.25] | 0.06 | 0.951 + Male, no, -40 | Male, yes, -40 | -0.88 | 1.18 | [-3.18, 1.43] | -0.75 | 0.455 + Female, no, -40 | Male, yes, -40 | -0.95 | 1.04 | [-2.99, 1.10] | -0.91 | 0.364 + Male, no, -40 | Female, yes, -40 | -0.79 | 1.07 | [-2.89, 1.30] | -0.74 | 0.457 + Female, no, -40 | Female, yes, -40 | -0.86 | 0.92 | [-2.67, 0.94] | -0.94 | 0.348 + Male, yes, -40 | Female, yes, -40 | 0.08 | 0.99 | [-1.86, 2.03] | 0.08 | 0.932 + Male, no, -40 | Male, no, 41-64 | 0.30 | 1.10 | [-1.84, 2.45] | 0.28 | 0.781 + Female, no, -40 | Male, no, 41-64 | 0.24 | 0.95 | [-1.63, 2.10] | 0.25 | 0.805 + Male, yes, -40 | Male, no, 41-64 | 1.18 | 1.02 | [-0.83, 3.19] | 1.15 | 0.248 + Female, yes, -40 | Male, no, 41-64 | 1.10 | 0.90 | [-0.66, 2.86] | 1.22 | 0.221 + Male, no, -40 | Female, no, 41-64 | 1.49 | 0.94 | [-0.35, 3.32] | 1.59 | 0.112 + Female, no, -40 | Female, no, 41-64 | 1.42 | 0.76 | [-0.08, 2.91] | 1.86 | 0.063 + Male, yes, -40 | Female, no, 41-64 | 2.37 | 0.85 | [ 0.70, 4.04] | 2.78 | 0.005 + Female, yes, -40 | Female, no, 41-64 | 2.28 | 0.69 | [ 0.92, 3.64] | 3.29 | 0.001 + Male, no, 41-64 | Female, no, 41-64 | 1.18 | 0.74 | [-0.26, 2.63] | 1.60 | 0.109 + Male, no, -40 | Male, yes, 41-64 | -0.23 | 1.04 | [-2.26, 1.81] | -0.22 | 0.827 + Female, no, -40 | Male, yes, 41-64 | -0.30 | 0.88 | [-2.03, 1.44] | -0.33 | 0.738 + Male, yes, -40 | Male, yes, 41-64 | 0.65 | 0.96 | [-1.23, 2.54] | 0.68 | 0.498 + Female, yes, -40 | Male, yes, 41-64 | 0.57 | 0.83 | [-1.05, 2.19] | 0.69 | 0.492 + Male, no, 41-64 | Male, yes, 41-64 | -0.53 | 0.86 | [-2.22, 1.16] | -0.62 | 0.538 + Female, no, 41-64 | Male, yes, 41-64 | -1.71 | 0.65 | [-2.98, -0.44] | -2.65 | 0.008 + Male, no, -40 | Female, yes, 41-64 | 1.15 | 0.94 | [-0.70, 2.99] | 1.22 | 0.224 + Female, no, -40 | Female, yes, 41-64 | 1.08 | 0.77 | [-0.43, 2.59] | 1.40 | 0.162 + Male, yes, -40 | Female, yes, 41-64 | 2.02 | 0.86 | [ 0.34, 3.71] | 2.36 | 0.018 + Female, yes, -40 | Female, yes, 41-64 | 1.94 | 0.70 | [ 0.56, 3.32] | 2.76 | 0.006 + Male, no, 41-64 | Female, yes, 41-64 | 0.84 | 0.74 | [-0.62, 2.30] | 1.13 | 0.258 + Female, no, 41-64 | Female, yes, 41-64 | -0.34 | 0.48 | [-1.28, 0.60] | -0.71 | 0.476 + Male, yes, 41-64 | Female, yes, 41-64 | 1.37 | 0.66 | [ 0.09, 2.66] | 2.09 | 0.036 + Male, no, -40 | Male, no, 65++ | 0.83 | 1.08 | [-1.29, 2.94] | 0.77 | 0.443 + Female, no, -40 | Male, no, 65++ | 0.76 | 0.93 | [-1.07, 2.58] | 0.82 | 0.415 + Male, yes, -40 | Male, no, 65++ | 1.71 | 1.01 | [-0.26, 3.68] | 1.70 | 0.090 + Female, yes, -40 | Male, no, 65++ | 1.62 | 0.88 | [-0.10, 3.34] | 1.85 | 0.064 + Male, no, 41-64 | Male, no, 65++ | 0.52 | 0.91 | [-1.26, 2.31] | 0.57 | 0.565 + Female, no, 41-64 | Male, no, 65++ | -0.66 | 0.71 | [-2.05, 0.73] | -0.93 | 0.354 + Male, yes, 41-64 | Male, no, 65++ | 1.06 | 0.84 | [-0.59, 2.70] | 1.26 | 0.209 + Female, yes, 41-64 | Male, no, 65++ | -0.32 | 0.72 | [-1.72, 1.09] | -0.44 | 0.658 + Male, no, -40 | Female, no, 65++ | 1.71 | 0.98 | [-0.21, 3.63] | 1.75 | 0.081 + Female, no, -40 | Female, no, 65++ | 1.64 | 0.81 | [ 0.04, 3.23] | 2.01 | 0.044 + Male, yes, -40 | Female, no, 65++ | 2.59 | 0.90 | [ 0.83, 4.35] | 2.88 | 0.004 + Female, yes, -40 | Female, no, 65++ | 2.50 | 0.75 | [ 1.03, 3.97] | 3.34 | < .001 + Male, no, 41-64 | Female, no, 65++ | 1.40 | 0.79 | [-0.14, 2.95] | 1.78 | 0.076 + Female, no, 41-64 | Female, no, 65++ | 0.22 | 0.55 | [-0.85, 1.29] | 0.40 | 0.686 + Male, yes, 41-64 | Female, no, 65++ | 1.94 | 0.71 | [ 0.55, 3.32] | 2.74 | 0.006 + Female, yes, 41-64 | Female, no, 65++ | 0.56 | 0.56 | [-0.53, 1.65] | 1.01 | 0.312 + Male, no, 65++ | Female, no, 65++ | 0.88 | 0.77 | [-0.62, 2.38] | 1.15 | 0.250 + Male, no, -40 | Male, yes, 65++ | -0.10 | 1.39 | [-2.83, 2.63] | -0.07 | 0.942 + Female, no, -40 | Male, yes, 65++ | -0.17 | 1.28 | [-2.69, 2.35] | -0.13 | 0.895 + Male, yes, -40 | Male, yes, 65++ | 0.78 | 1.34 | [-1.85, 3.40] | 0.58 | 0.561 + Female, yes, -40 | Male, yes, 65++ | 0.69 | 1.25 | [-1.75, 3.13] | 0.56 | 0.577 + Male, no, 41-64 | Male, yes, 65++ | -0.41 | 1.27 | [-2.89, 2.08] | -0.32 | 0.750 + Female, no, 41-64 | Male, yes, 65++ | -1.59 | 1.13 | [-3.81, 0.64] | -1.40 | 0.162 + Male, yes, 41-64 | Male, yes, 65++ | 0.13 | 1.22 | [-2.26, 2.52] | 0.10 | 0.917 + Female, yes, 41-64 | Male, yes, 65++ | -1.25 | 1.14 | [-3.48, 0.99] | -1.09 | 0.274 + Male, no, 65++ | Male, yes, 65++ | -0.93 | 1.25 | [-3.39, 1.53] | -0.74 | 0.459 + Female, no, 65++ | Male, yes, 65++ | -1.81 | 1.17 | [-4.10, 0.48] | -1.55 | 0.122 + Male, no, -40 | Female, yes, 65++ | 0.39 | 1.36 | [-2.27, 3.05] | 0.29 | 0.775 + Female, no, -40 | Female, yes, 65++ | 0.32 | 1.25 | [-2.12, 2.76] | 0.26 | 0.798 + Male, yes, -40 | Female, yes, 65++ | 1.27 | 1.30 | [-1.28, 3.82] | 0.97 | 0.330 + Female, yes, -40 | Female, yes, 65++ | 1.18 | 1.20 | [-1.18, 3.54] | 0.98 | 0.326 + Male, no, 41-64 | Female, yes, 65++ | 0.08 | 1.23 | [-2.33, 2.49] | 0.07 | 0.946 + Female, no, 41-64 | Female, yes, 65++ | -1.10 | 1.09 | [-3.23, 1.04] | -1.01 | 0.313 + Male, yes, 41-64 | Female, yes, 65++ | 0.62 | 1.18 | [-1.69, 2.92] | 0.52 | 0.601 + Female, yes, 41-64 | Female, yes, 65++ | -0.76 | 1.09 | [-2.90, 1.39] | -0.69 | 0.489 + Male, no, 65++ | Female, yes, 65++ | -0.44 | 1.21 | [-2.82, 1.94] | -0.36 | 0.717 + Female, no, 65++ | Female, yes, 65++ | -1.32 | 1.13 | [-3.53, 0.89] | -1.17 | 0.241 + Male, yes, 65++ | Female, yes, 65++ | 0.49 | 1.50 | [-2.45, 3.43] | 0.33 | 0.745 + + Variable predicted: qol + Predictors contrasted: gender, employed, age + +--- + + Code + print(estimate_contrasts(estim, contrast = c("gender", "employed"), by = "age"), + zap_small = TRUE, table_width = Inf) + Output + Model-based Contrasts Analysis + + Level1 | Level2 | age | Difference | SE | 95% CI | Statistic | p + ----------------------------------------------------------------------------------------- + Male, no | Female, no | -40 | 0.07 | 1.11 | [-2.11, 2.25] | 0.06 | 0.951 + Male, no | Male, yes | -40 | -0.88 | 1.18 | [-3.18, 1.43] | -0.75 | 0.455 + Female, no | Male, yes | -40 | -0.95 | 1.04 | [-2.99, 1.10] | -0.91 | 0.364 + Male, no | Female, yes | -40 | -0.79 | 1.07 | [-2.89, 1.30] | -0.74 | 0.457 + Female, no | Female, yes | -40 | -0.86 | 0.92 | [-2.67, 0.94] | -0.94 | 0.348 + Male, yes | Female, yes | -40 | 0.08 | 0.99 | [-1.86, 2.03] | 0.08 | 0.932 + Male, no | Female, no | 41-64 | 1.18 | 0.74 | [-0.26, 2.63] | 1.60 | 0.109 + Male, no | Male, yes | 41-64 | -0.53 | 0.86 | [-2.22, 1.16] | -0.62 | 0.538 + Female, no | Male, yes | 41-64 | -1.71 | 0.65 | [-2.98, -0.44] | -2.65 | 0.008 + Male, no | Female, yes | 41-64 | 0.84 | 0.74 | [-0.62, 2.30] | 1.13 | 0.258 + Female, no | Female, yes | 41-64 | -0.34 | 0.48 | [-1.28, 0.60] | -0.71 | 0.476 + Male, yes | Female, yes | 41-64 | 1.37 | 0.66 | [ 0.09, 2.66] | 2.09 | 0.036 + Male, no | Female, no | 65+ | 0.88 | 0.77 | [-0.62, 2.38] | 1.15 | 0.250 + Male, no | Male, yes | 65+ | -0.93 | 1.25 | [-3.39, 1.53] | -0.74 | 0.459 + Female, no | Male, yes | 65+ | -1.81 | 1.17 | [-4.10, 0.48] | -1.55 | 0.122 + Male, no | Female, yes | 65+ | -0.44 | 1.21 | [-2.82, 1.94] | -0.36 | 0.717 + Female, no | Female, yes | 65+ | -1.32 | 1.13 | [-3.53, 0.89] | -1.17 | 0.241 + Male, yes | Female, yes | 65+ | 0.49 | 1.50 | [-2.45, 3.43] | 0.33 | 0.745 + + Variable predicted: qol + Predictors contrasted: gender, employed + +--- + + Code + print(estimate_contrasts(estim, contrast = "employed", by = c("age", "gender")), + zap_small = TRUE, table_width = Inf) + Output + Model-based Contrasts Analysis + + Level1 | Level2 | age | gender | Difference | SE | 95% CI | Statistic | p + ---------------------------------------------------------------------------------------- + no | yes | -40 | Female | -0.86 | 0.92 | [-2.67, 0.94] | -0.94 | 0.348 + no | yes | -40 | Male | -0.88 | 1.18 | [-3.18, 1.43] | -0.75 | 0.455 + no | yes | 41-64 | Female | -0.34 | 0.48 | [-1.28, 0.60] | -0.71 | 0.476 + no | yes | 41-64 | Male | -0.53 | 0.86 | [-2.22, 1.16] | -0.62 | 0.538 + no | yes | 65+ | Female | -1.32 | 1.13 | [-3.53, 0.89] | -1.17 | 0.241 + no | yes | 65+ | Male | -0.93 | 1.25 | [-3.39, 1.53] | -0.74 | 0.459 + + Variable predicted: qol + Predictors contrasted: employed + +# estimate_contrasts - Random Effects Levels, interaction + + Code + print(estimate_contrasts(estim, contrast = c("age", "employed"), comparison = "interaction"), + zap_small = TRUE, table_width = Inf) + Output + Model-based Contrasts Analysis + + age | employed | Difference | SE | 95% CI | Statistic | p + ------------------------------------------------------------------------------- + -40-41-64 | no and yes | 0 | 2.21 | [-4.34, 4.34] | 0 | > .999 + -40-65+ | no and yes | 0 | 2.21 | [-4.34, 4.34] | 0 | > .999 + 41-64-65+ | no and yes | 0 | 2.21 | [-4.34, 4.34] | 0 | > .999 + + Variable predicted: qol + Predictors contrasted: age, employed + diff --git a/tests/testthat/_snaps/joint_test.md b/tests/testthat/_snaps/joint_test.md new file mode 100644 index 000000000..2567a8d60 --- /dev/null +++ b/tests/testthat/_snaps/joint_test.md @@ -0,0 +1,14 @@ +# estimate_contrasts - joint test, works with select printing + + Code + estimate_contrasts(m, contrast = "time", by = "coffee", comparison = "joint") + Output + Marginal Joint Test + + Contrast | coffee | df1 | df2 | F | p + --------------------------------------------- + time | coffee | 2 | 114 | 1.76 | 0.176 + time | control | 2 | 114 | 5.29 | 0.006 + + p-values are uncorrected. +