Adjusted p values from afex::aov_4() not retained by model_parameters() (#1102)

* Adjusted p values from `afex::aov_4()` not retained by `model_parameters()`
Fixes #1101

* add test

* fix test

* automatically detect p-adjustment in afex

* fix

* fix., add test

* update and add test, news

* p_adjust for all anova stuff

#1101

* docs

* styler

---------

Co-authored-by: Mattan S. Ben-Shachar <[email protected]>

Author: strengejacke

DESCRIPTION

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

NEWS.md

@@ -6,6 +6,11 @@
   `brmsfit` and `stanreg` gets an additional `"grouplevel"` option, to return
   the group-level estimates for random effects.
 
+* `model_parameters()` for Anova-objects gains a `p_adjust` argument, to apply
+  p-adjustment where possible. Furthermore, for models from package *afex*, where
+  p-adjustment was applied during model-fitting, the correct p-values are now
+  returned (before, unadjusted p-values were returned in some cases).
+
 # parameters 0.25.0
 
 ## Changes

R/extract_parameters_anova.R

@@ -1,5 +1,8 @@
 #' @keywords internal
-.extract_parameters_anova <- function(model, test = "multivariate") {
+.extract_parameters_anova <- function(model,
+                                      test = "multivariate",
+                                      p_adjust = NULL,
+                                      verbose = TRUE) {
   # Processing
   if (inherits(model, "manova")) {
     parameters <- .extract_anova_manova(model)
@@ -19,6 +22,12 @@
     parameters <- .extract_anova_aov_svyglm(model)
   }
 
+  # remove intercept
+  intercepts <- parameters$Parameter %in% c("Intercept", "(Intercept)")
+  if (any(intercepts)) {
+    parameters <- parameters[!intercepts, ]
+  }
+
   # Rename
 
   # p-values
@@ -79,6 +88,11 @@
   )
   parameters <- parameters[col_order[col_order %in% names(parameters)]]
 
+  # ==== adjust p-values?
+  if (!is.null(p_adjust)) {
+    parameters <- .p_adjust(parameters, p_adjust, model, verbose)
+  }
+
   insight::text_remove_backticks(parameters, verbose = FALSE)
 }
 

R/methods_aov.R

@@ -102,6 +102,7 @@ model_parameters.aov <- function(model,
                                  df_error = NULL,
                                  ci = NULL,
                                  alternative = NULL,
+                                 p_adjust = NULL,
                                  test = NULL,
                                  power = FALSE,
                                  es_type = NULL,
@@ -138,7 +139,7 @@ model_parameters.aov <- function(model,
   }
 
   # extract standard parameters
-  params <- .extract_parameters_anova(model, test)
+  params <- .extract_parameters_anova(model, test, p_adjust = p_adjust, verbose = verbose)
 
   # add effect sizes, if available
   params <- .effectsizes_for_aov(
@@ -251,16 +252,19 @@ model_parameters.afex_aov <- function(model,
                                       type = NULL,
                                       keep = NULL,
                                       drop = NULL,
+                                      p_adjust = NULL,
                                       verbose = TRUE,
                                       ...) {
   if (inherits(model$Anova, "Anova.mlm")) {
     params <- model$anova_table
     with_df_and_p <- summary(model$Anova)$univariate.tests
     params$`Sum Sq` <- with_df_and_p[-1, 1]
     params$`Error SS` <- with_df_and_p[-1, 3]
-    out <- .extract_parameters_anova(params, test = NULL)
+    out <- .extract_parameters_anova(params, test = NULL, p_adjust = NULL, verbose)
+    p_adjust <- .extract_p_adjust_afex(model, p_adjust)
   } else {
-    out <- .extract_parameters_anova(model$Anova, test = NULL)
+    p_adjust <- .extract_p_adjust_afex(model, p_adjust)
+    out <- .extract_parameters_anova(model$Anova, test = NULL, p_adjust, verbose)
   }
 
   out <- .effectsizes_for_aov(
@@ -273,7 +277,15 @@ model_parameters.afex_aov <- function(model,
   )
 
   # add attributes
-  out <- .add_anova_attributes(out, model, ci, test = NULL, alternative = NULL, ...)
+  out <- .add_anova_attributes(
+    out,
+    model,
+    ci,
+    test = NULL,
+    alternative = NULL,
+    p_adjust = p_adjust,
+    ...
+  )
 
   # filter parameters
   if (!is.null(keep) || !is.null(drop)) {
@@ -561,3 +573,17 @@ model_parameters.seqanova.svyglm <- model_parameters.aov
 
   data[, col_order]
 }
+
+
+#' @keywords internal
+.extract_p_adjust_afex <- function(model, p_adjust) {
+  if (is.null(p_adjust) && inherits(model, "afex_aov")) {
+    p_adjust <- attr(model$anova_table, "p_adjust_method")
+
+    if (p_adjust == "none") {
+      p_adjust <- NULL
+    }
+  }
+
+  p_adjust
+}

R/utils_model_parameters.R

@@ -388,13 +388,14 @@
 
 
 #' @keywords internal
-.add_anova_attributes <- function(params, model, ci, test = NULL, alternative = NULL, ...) {
+.add_anova_attributes <- function(params, model, ci, test = NULL, alternative = NULL, p_adjust = NULL, ...) {
   dot.arguments <- lapply(match.call(expand.dots = FALSE)$`...`, function(x) x) # nolint
 
   attr(params, "ci") <- ci
   attr(params, "model_class") <- class(model)
   attr(params, "anova_type") <- .anova_type(model)
   attr(params, "text_alternative") <- .anova_alternative(params, alternative)
+  attr(params, "p_adjust") <- p_adjust
 
   if (inherits(model, "Anova.mlm") && !identical(test, "univariate")) {
     attr(params, "anova_test") <- model$test

man/model_parameters.aov.Rd

@@ -12,16 +12,76 @@ test_that("afex_aov", {
   mp1 <- model_parameters(m_between, verbose = FALSE)
   mp2 <- model_parameters(m_within, verbose = FALSE)
 
-  expect_equal(c(nrow(mp1), ncol(mp1)), c(5, 7))
-  expect_equal(mp1$Sum_Squares, c(450.62069, 11.98202, 5.56322, 8.68275, 15.2037), tolerance = 1e-3)
+  expect_equal(c(nrow(mp1), ncol(mp1)), c(4, 7))
+  expect_equal(mp1$Sum_Squares, c(11.98202, 5.56322, 8.68275, 15.2037), tolerance = 1e-3)
   expect_equal(c(nrow(mp2), ncol(mp2)), c(3, 9))
   expect_equal(mp2$Sum_Squares, c(167.5, 106.29167, 11.08333), tolerance = 1e-3)
-  expect_equal(
-    colnames(mp1),
+  expect_named(
+    mp1,
     c("Parameter", "Sum_Squares", "df", "Mean_Square", "F", "p", "Method")
   )
-  expect_equal(
-    colnames(mp2),
-    c("Parameter", "Sum_Squares", "Sum_Squares_Error", "df", "df_error", "Mean_Square", "F", "p", "Method")
+  expect_named(
+    mp2,
+    c(
+      "Parameter", "Sum_Squares", "Sum_Squares_Error", "df", "df_error",
+      "Mean_Square", "F", "p", "Method"
+    )
+  )
+})
+
+
+test_that("afex_aov, p-adjustement", {
+  skip_if_not_installed("afex")
+  data(laptop_urry, package = "afex")
+  afx <- afex::aov_4(
+    overall ~ condition * talk + (1 | pid),
+    data = laptop_urry,
+    anova_table = list(p_adjust_method = "bonferroni")
+  )
+  out1 <- model_parameters(afx, ci = 0.95)
+  out2 <- model_parameters(afx, ci = 0.95, p_adjust = "bonferroni")
+
+  expect_identical(dim(out1), c(4L, 7L))
+  expect_equal(out1$Sum_Squares, c(115.01087, 6703.72241, 1944.0391, 29101.23396), tolerance = 1e-3)
+  expect_named(
+    out1,
+    c("Parameter", "Sum_Squares", "df", "Mean_Square", "F", "p", "Method")
+  )
+  expect_equal(out1$p, c(1, 0, 0.2157, NA), tolerance = 1e-3)
+  expect_equal(out2$p, c(1, 0, 0.2157, NA), tolerance = 1e-3)
+
+  afx <- afex::aov_4(
+    overall ~ condition * talk + (1 | pid),
+    data = laptop_urry
+  )
+  out3 <- model_parameters(afx, ci = 0.95)
+  out4 <- model_parameters(afx, ci = 0.95, p_adjust = "bonferroni")
+  expect_equal(out3$p, c(0.4714, 0, 0.0719, NA), tolerance = 1e-3)
+  expect_equal(out4$p, c(1, 0, 0.2157, NA), tolerance = 1e-3)
+})
+
+
+test_that("afex_aov_ez, p-adjustement", {
+  skip_if_not_installed("afex")
+  data(obk.long, package = "afex")
+  a2 <- afex::aov_ez(
+    "id",
+    "value",
+    data = obk.long,
+    between = c("treatment", "gender"),
+    within = c("phase", "hour"),
+    observed = "gender",
+    anova_table = list(p_adjust_method = "fdr")
+  )
+
+  out <- model_parameters(a2)
+  expect_equal(a2$anova_table$`Pr(>F)`, out$p, tolerance = 1e-4)
+  expect_identical(dim(out), c(15L, 9L))
+  expect_named(
+    out,
+    c(
+      "Parameter", "Sum_Squares", "Sum_Squares_Error", "df", "df_error",
+      "Mean_Square", "F", "p", "Method"
+    )
   )
 })

tests/testthat/test-model_parameters.afex_aov.R

@@ -68,20 +68,20 @@ test_that("print-model_parameters", {
 })
 
 
-test_that("model_parameters_Anova.mlm", {
+test_that("model_parameters_Anova.mlm-1", {
   skip_if_not_installed("car")
 
   m <- lm(cbind(hp, mpg) ~ factor(cyl) * am, data = mtcars)
   a <- car::Anova(m, type = 3, test.statistic = "Pillai")
   mp <- model_parameters(a, verbose = FALSE)
 
   expect_named(mp, c("Parameter", "df", "Statistic", "df_num", "df_error", "F", "p"))
-  expect_equal(mp[["F"]], c(158.2578, 6.60593, 3.71327, 3.28975), tolerance = 1e-3)
-  expect_equal(mp$Statistic, c(0.9268, 0.67387, 0.22903, 0.4039), tolerance = 1e-3)
+  expect_equal(mp[["F"]], c(6.60593, 3.71327, 3.28975), tolerance = 1e-3)
+  expect_equal(mp$Statistic, c(0.67387, 0.22903, 0.4039), tolerance = 1e-3)
 })
 
 
-test_that("model_parameters_Anova.mlm", {
+test_that("model_parameters_Anova.mlm-2", {
   skip_if_not_installed("MASS")
   skip_if_not_installed("car")
   data(housing, package = "MASS")