mismatch between {emmeans}\s and model_parameters()'s adjusted p-values

Fixes #1103

Author: strengejacke

DESCRIPTION

@@ -1,7 +1,7 @@
 Type: Package
 Package: parameters
 Title: Processing of Model Parameters
-Version: 0.28.3.3
+Version: 0.28.3.4
 Authors@R:
     c(person(given = "Daniel",
              family = "Lüdecke",

NEWS.md

@@ -17,6 +17,8 @@
 * Fixed issue where wrong (non-robust) standard errors were calculated for
   `coxph` and `svycoxph` objects.
 
+* Fixes issues with Tukey-p-value adjustment for *emmeans* objects.
+
 # parameters 0.28.3
 
 * fixed bug in `standardize_info(<fixest>)` that was preventing

R/format_p_adjust.R

@@ -14,7 +14,8 @@
 format_p_adjust <- function(method) {
   method <- tolower(method)
 
-  switch(method,
+  switch(
+    method,
     holm = "Holm (1979)",
     hochberg = "Hochberg (1988)",
     hommel = "Hommel (1988)",
@@ -68,13 +69,19 @@ format_p_adjust <- function(method) {
       } else if (tolower(p_adjust) == "sidak") {
         # sidak adjustment
         params$p <- 1 - (1 - params$p)^(nrow(params) / rank_adjust)
-      }  else if (tolower(p_adjust) == "sup-t") {
+      } else if (tolower(p_adjust) == "sup-t") {
         # sup-t adjustment
         params <- .p_adjust_supt(model, params)
       }
 
-      if (isTRUE(all(old_p_vals == params$p)) && !identical(p_adjust, "none") && verbose) {
-        insight::format_warning(paste0("Could not apply ", p_adjust, "-adjustment to p-values. Either something went wrong, or the non-adjusted p-values were already very large.")) # nolint
+      if (
+        isTRUE(all(old_p_vals == params$p)) && !identical(p_adjust, "none") && verbose
+      ) {
+        insight::format_warning(paste0(
+          "Could not apply ",
+          p_adjust,
+          "-adjustment to p-values. Either something went wrong, or the non-adjusted p-values were already very large."
+        )) # nolint
       }
     } else if (verbose) {
       insight::format_alert(paste0("`p_adjust` must be one of ", toString(all_methods)))
@@ -93,7 +100,10 @@ format_p_adjust <- function(method) {
       by_vars <- model@misc$by.vars
       if (!is.null(by_vars) && by_vars %in% colnames(params)) {
         correction <- insight::n_unique(params[[by_vars]])
+      } else {
+        correction <- prod(vapply(model@model.info$xlev, length, numeric(1)))
       }
+
       correction
     },
     error = function(e) {
@@ -106,22 +116,26 @@ format_p_adjust <- function(method) {
 # tukey adjustment -----
 
 .p_adjust_tukey <- function(params, stat_column, rank_adjust = 1, verbose = TRUE) {
-  df_column <- colnames(params)[stats::na.omit(match(c("df", "df_error"), colnames(params)))][1]
+  df_column <- colnames(params)[stats::na.omit(match(
+    c("df", "df_error"),
+    colnames(params)
+  ))][1]
   if (!is.na(df_column) && length(stat_column)) {
     params$p <- suppressWarnings(stats::ptukey(
       sqrt(2) * abs(params[[stat_column]]),
-      nmeans = nrow(params) / rank_adjust,
+      nmeans = rank_adjust,
       df = params[[df_column]],
       lower.tail = FALSE
     ))
     # for specific contrasts, ptukey might fail, and the tukey-adjustement
     # could just be simple p-value calculation
     if (all(is.na(params$p))) {
-      params$p <- 2 * stats::pt(
-        abs(params[[stat_column]]),
-        df = params[[df_column]],
-        lower.tail = FALSE
-      )
+      params$p <- 2 *
+        stats::pt(
+          abs(params[[stat_column]]),
+          df = params[[df_column]],
+          lower.tail = FALSE
+        )
       verbose <- FALSE
     }
   }
@@ -132,7 +146,10 @@ format_p_adjust <- function(method) {
 # scheffe adjustment -----
 
 .p_adjust_scheffe <- function(model, params, stat_column, rank_adjust = 1) {
-  df_column <- colnames(params)[stats::na.omit(match(c("df", "df_error"), colnames(params)))][1]
+  df_column <- colnames(params)[stats::na.omit(match(
+    c("df", "df_error"),
+    colnames(params)
+  ))][1]
   if (!is.na(df_column) && length(stat_column)) {
     # 1st try
     scheffe_ranks <- try(qr(model@linfct)$rank, silent = TRUE)
@@ -146,7 +163,8 @@ format_p_adjust <- function(method) {
       scheffe_ranks <- nrow(params)
     }
     scheffe_ranks <- scheffe_ranks / rank_adjust
-    params$p <- stats::pf(params[[stat_column]]^2 / scheffe_ranks,
+    params$p <- stats::pf(
+      params[[stat_column]]^2 / scheffe_ranks,
       df1 = scheffe_ranks,
       df2 = params[[df_column]],
       lower.tail = FALSE
@@ -182,7 +200,9 @@ format_p_adjust <- function(method) {
   # get correlation matrix, based on the covariance matrix
   vc <- .safe(stats::cov2cor(insight::get_varcov(model, component = component)))
   if (is.null(vc)) {
-    insight::format_warning("Could not calculate covariance matrix for `sup-t` adjustment.")
+    insight::format_warning(
+      "Could not calculate covariance matrix for `sup-t` adjustment."
+    )
     return(params)
   }
   # get confidence interval level, or set default
@@ -197,18 +217,30 @@ format_p_adjust <- function(method) {
   }
   # calculate updated confidence interval level, based on simultaenous
   # confidence intervals (https://onlinelibrary.wiley.com/doi/10.1002/jae.2656)
-  crit <- mvtnorm::qmvt(ci_level, df = params[[df_column]][1], tail = "both.tails", corr = vc)$quantile
+  crit <- mvtnorm::qmvt(
+    ci_level,
+    df = params[[df_column]][1],
+    tail = "both.tails",
+    corr = vc
+  )$quantile
   # update confidence intervals
   params$CI_low <- params$Coefficient - crit * params$SE
   params$CI_high <- params$Coefficient + crit * params$SE
   # udpate p-values
   for (i in 1:nrow(params)) {
-    params$p[i] <- 1 - mvtnorm::pmvt(
-      lower = rep(-abs(stats::qt(params$p[i] / 2, df = params[[df_column]][i])), nrow(vc)),
-      upper = rep(abs(stats::qt(params$p[i] / 2, df = params[[df_column]][i])), nrow(vc)),
-      corr = vc,
-      df = params[[df_column]][i]
-    )
+    params$p[i] <- 1 -
+      mvtnorm::pmvt(
+        lower = rep(
+          -abs(stats::qt(params$p[i] / 2, df = params[[df_column]][i])),
+          nrow(vc)
+        ),
+        upper = rep(
+          abs(stats::qt(params$p[i] / 2, df = params[[df_column]][i])),
+          nrow(vc)
+        ),
+        corr = vc,
+        df = params[[df_column]][i]
+      )
   }
   params
 }

tests/testthat/test-p_adjust.R

@@ -29,7 +29,12 @@ test_that("model_parameters, p-adjust after keep/drop", {
   )
 
   expect_message(
-    mp <- model_parameters(model, include_info = TRUE, keep = c("wt", "hp"), p_adjust = "bonferroni"),
+    mp <- model_parameters(
+      model,
+      include_info = TRUE,
+      keep = c("wt", "hp"),
+      p_adjust = "bonferroni"
+    ),
     "more than 1 element"
   )
   expect_equal(
@@ -40,7 +45,12 @@ test_that("model_parameters, p-adjust after keep/drop", {
   )
 
   expect_message(
-    mp <- model_parameters(model, include_info = TRUE, keep = c("cyl", "gear"), p_adjust = "bonferroni"),
+    mp <- model_parameters(
+      model,
+      include_info = TRUE,
+      keep = c("cyl", "gear"),
+      p_adjust = "bonferroni"
+    ),
     "more than 1 element"
   )
   expect_equal(
@@ -58,9 +68,26 @@ test_that("model_parameters, emmeans, p-adjust", {
   mp <- model_parameters(m)
   expect_equal(mp$p, as.data.frame(m)$p.value, tolerance = 1e-4)
 
-  m <- pairs(emmeans::emmeans(aov(Sepal.Width ~ Species, data = iris), ~Species), adjust = "scheffe")
+  m <- pairs(
+    emmeans::emmeans(aov(Sepal.Width ~ Species, data = iris), ~Species),
+    adjust = "scheffe"
+  )
   mp <- model_parameters(m, p_adjust = "scheffe")
   expect_equal(mp$p, as.data.frame(m)$p.value, tolerance = 1e-4)
+
+  m <- pairs(
+    emmeans::emmeans(aov(Sepal.Width ~ Species, data = iris), ~Species),
+    adjust = "tukey"
+  )
+  mp <- model_parameters(m, p_adjust = "tukey")
+  expect_equal(mp$p, as.data.frame(m)$p.value, tolerance = 1e-4)
+
+  data(warpbreaks)
+  warp.lm <- lm(breaks ~ wool * tension, data = warpbreaks)
+
+  myc <- pairs(emmeans::emmeans(warp.lm, ~ wool + tension))
+  mp <- model_parameters(myc, p_adjust = "tukey")
+  expect_equal(mp$p, as.data.frame(myc)$p.value, tolerance = 1e-4)
 })
 
 
@@ -84,7 +111,11 @@ test_that("model_parameters, simultaenous confidence intervals", {
   set.seed(123)
   mp <- model_parameters(model, p_adjust = "sup-t")
   expect_equal(mp$p, c(0, 0.27904, 0, 0, NA, NA), tolerance = 1e-3)
-  expect_equal(mp$CI_low, c(67.70195, -0.48377, -4.66802, -7.51974, 8.42651, 11.50991), tolerance = 1e-3)
+  expect_equal(
+    mp$CI_low,
+    c(67.70195, -0.48377, -4.66802, -7.51974, 8.42651, 11.50991),
+    tolerance = 1e-3
+  )
 
   skip_if_not_installed("glmmTMB")
   data("Salamanders", package = "glmmTMB")
@@ -99,18 +130,49 @@ test_that("model_parameters, simultaenous confidence intervals", {
   expect_equal(
     mp$p,
     c(
-      0.56769, 0.57466, 0.98029, 0.83123, 0.22681, 0.06271, 0.99876,
-      0.00068, 0.61786, 0.95269, 0.81296, 0.60973, 0.97504, 0.80566,
-      0.81871, 0.00024, NA, NA
+      0.56769,
+      0.57466,
+      0.98029,
+      0.83123,
+      0.22681,
+      0.06271,
+      0.99876,
+      0.00068,
+      0.61786,
+      0.95269,
+      0.81296,
+      0.60973,
+      0.97504,
+      0.80566,
+      0.81871,
+      0.00024,
+      NA,
+      NA
     ),
     tolerance = 1e-3
   )
   expect_equal(
     mp$CI_low,
     c(
-      -1.6935, -2.6841, -0.45839, -1.30244, -0.14911, -0.01933, -0.76516,
-      0.4494, -0.77256, -2.41501, -3.08449, -0.87083, -2.50859, -2.91223,
-      -8.38616, -4.18299, 0.92944, 0.16671
+      -1.6935,
+      -2.6841,
+      -0.45839,
+      -1.30244,
+      -0.14911,
+      -0.01933,
+      -0.76516,
+      0.4494,
+      -0.77256,
+      -2.41501,
+      -3.08449,
+      -0.87083,
+      -2.50859,
+      -2.91223,
+      -8.38616,
+      -4.18299,
+      0.92944,
+      0.16671
     ),
-    tolerance = 1e-3)
+    tolerance = 1e-3
+  )
 })