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

[profandyfield]

Problem

In afex::aov_4() you can use anova_table = list(p_adjust_method = "bonferroni") to adjust p-values in the ANOVA table (obviously you can chose other methods to Bonferroni). However, if you pass the afex model to model_parameters() it reverts to the unadjusted p-values.

goggles_tib <- discovr::goggles
goggles_afx <- afex::aov_4(attractiveness ~ facetype*alcohol + (1|id), data = goggles_tib, anova_table = list(p_adjust_method = "bonferroni"))
#> Contrasts set to contr.sum for the following variables: facetype, alcohol

goggles_afx
#> Anova Table (Type 3 tests, bonferroni-adjusted)
#> 
#> Response: attractiveness
#>             Effect    df  MSE         F  ges p.value
#> 1         facetype 1, 42 1.37 15.58 *** .271   <.001
#> 2          alcohol 2, 42 1.37    6.04 * .223    .015
#> 3 facetype:alcohol 2, 42 1.37   8.51 ** .288    .002
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '+' 0.1 ' ' 1
parameters::model_parameters(goggles_afx, es_type = "omega", ci = 0.95)
#> ANOVA estimation for factorial designs using 'afex'
#> 
#> Parameter        | Sum_Squares | df | Mean_Square |       F |      p
#> --------------------------------------------------------------------
#> (Intercept)      |     1541.33 |  1 |     1541.33 | 1125.84 | < .001
#> facetype         |       21.33 |  1 |       21.33 |   15.58 | < .001
#> alcohol          |       16.54 |  2 |        8.27 |    6.04 | 0.005 
#> facetype:alcohol |       23.29 |  2 |       11.65 |    8.51 | < .001
#> Residuals        |       57.50 | 42 |        1.37 |         |       
#> 
#> Parameter        | Omega2 (partial) | Omega2 95% CI
#> ---------------------------------------------------
#> (Intercept)      |                  |              
#> facetype         |             0.23 |  [0.07, 1.00]
#> alcohol          |             0.17 |  [0.02, 1.00]
#> facetype:alcohol |             0.24 |  [0.06, 1.00]
#> Residuals        |                  |              
#> 
#> Anova Table (Type 3 tests)

^{Created on 2025-05-07 with reprex v2.1.1}

Suggestion

I think ideally, the user should able to use p_adjust within model_parameters() like you can for lm models; that way, the workflow is consistent with model_parameters() being used to extract and augment (for example with effect sizes) information from models. Also the syntax is nicer than anova_table = list(p_adjust_method = "bonferroni"). If that's not possible, then retain p-values from afex to avoid the user thinking they have corrected when in fact they haven't.

[strengejacke]

It's implemented in #1102, but the adjusted p-values slightly differ, because we have one more parameter (Intercept) in our result. @mattansb is that correct, or is the Intercept parameter a bug?

[strengejacke]

goggles_tib <- discovr::goggles
goggles_afx <- afex::aov_4(
  attractiveness ~ facetype * alcohol + (1 | id),
  data = goggles_tib,
  anova_table = list(p_adjust_method = "bonferroni")
)
#> Contrasts set to contr.sum for the following variables: facetype, alcohol
goggles_afx
#> Anova Table (Type 3 tests, bonferroni-adjusted)
#> 
#> Response: attractiveness
#>             Effect    df  MSE         F  ges p.value
#> 1         facetype 1, 42 1.37 15.58 *** .271   <.001
#> 2          alcohol 2, 42 1.37    6.04 * .223    .015
#> 3 facetype:alcohol 2, 42 1.37   8.51 ** .288    .002
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '+' 0.1 ' ' 1


parameters::model_parameters(goggles_afx, ci = 0.95)
#> ANOVA estimation for factorial designs using 'afex'
#> 
#> Parameter        | Sum_Squares | df | Mean_Square |       F |      p
#> --------------------------------------------------------------------
#> (Intercept)      |     1541.33 |  1 |     1541.33 | 1125.84 | < .001
#> facetype         |       21.33 |  1 |       21.33 |   15.58 | < .001
#> alcohol          |       16.54 |  2 |        8.27 |    6.04 | 0.005 
#> facetype:alcohol |       23.29 |  2 |       11.65 |    8.51 | < .001
#> Residuals        |       57.50 | 42 |        1.37 |         |       
#> 
#> Anova Table (Type 3 tests)

parameters::model_parameters(goggles_afx, ci = 0.95, p_adjust = "bonferroni")
#> ANOVA estimation for factorial designs using 'afex'
#> 
#> Parameter        | Sum_Squares | df | Mean_Square |       F |      p
#> --------------------------------------------------------------------
#> (Intercept)      |     1541.33 |  1 |     1541.33 | 1125.84 | < .001
#> facetype         |       21.33 |  1 |       21.33 |   15.58 | 0.001 
#> alcohol          |       16.54 |  2 |        8.27 |    6.04 | 0.020 
#> facetype:alcohol |       23.29 |  2 |       11.65 |    8.51 | 0.003 
#> Residuals        |       57.50 | 42 |        1.37 |         |       
#> 
#> p-value adjustment method: Bonferroni
#> Anova Table (Type 3 tests)

^{Created on 2025-05-07 with reprex v2.1.1}

[mattansb]

I suspect that since the intercept is of little interest to anyone in ANOVA tables, ignoring it makes sense.

But also - why aren't we pulling the p-values directly from the afex table?

[strengejacke]

We do, but the adjustment seems to be done on printing?

goggles_tib <- discovr::goggles
goggles_afx <- afex::aov_4(
  attractiveness ~ facetype * alcohol + (1 | id),
  data = goggles_tib,
  anova_table = list(p_adjust_method = "bonferroni")
)
#> Contrasts set to contr.sum for the following variables: facetype, alcohol
goggles_afx
#> Anova Table (Type 3 tests, bonferroni-adjusted)
#> 
#> Response: attractiveness
#>             Effect    df  MSE         F  ges p.value
#> 1         facetype 1, 42 1.37 15.58 *** .271   <.001
#> 2          alcohol 2, 42 1.37    6.04 * .223    .015
#> 3 facetype:alcohol 2, 42 1.37   8.51 ** .288    .002
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '+' 0.1 ' ' 1
goggles_afx$Anova
#> Anova Table (Type III tests)
#> 
#> Response: dv
#>                   Sum Sq Df   F value    Pr(>F)    
#> (Intercept)      1541.33  1 1125.8435 < 2.2e-16 ***
#> facetype           21.33  1   15.5826 0.0002952 ***
#> alcohol            16.54  2    6.0413 0.0049434 ** 
#> facetype:alcohol   23.29  2    8.5065 0.0007913 ***
#> Residuals          57.50 42                        
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

^{Created on 2025-05-07 with reprex v2.1.1}

And the anova table also contains the intercept, it's just not printed or removed prior to printing. Since I have no experience with anovas, I can't say whether to keep or remove the intercept.

Using our own p_adjust gives us a bit more flexibility anyway, but what to do with the intercept?

[mattansb]

The adjusted p-values are available in $anova_table:

goggles_afx$anova_table
#> Anova Table (Type 3 tests, bonferroni-adjusted)
#> 
#> Response: attractiveness
#>                  num Df den Df    MSE       F     ges    Pr(>F)    
#> facetype              1     42 1.3691 15.5826 0.27061 0.0008857 ***
#> alcohol               2     42 1.3691  6.0413 0.22341 0.0148302 *  
#> facetype:alcohol      2     42 1.3691  8.5065 0.28829 0.0023738 ** 
#> ---
#> Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

($Anova is the car::Anova object, $anova_table is just a data frame of the results.)

[strengejacke]

That table lacks the Sums of Squares. Would it be ok to take that table, and merge matching / missing rows from the other table? But I have no idea how stable this would be. Maybe just take $Anova, remove intercept, and then apply our adjustment?

[mattansb]

I think re-doing the adjustment is fine - it's not computationally expensive, and allows for more flexibility.

I've added automatic detection of the adjutant set (and also made sure it works with mixed designs as well):

goggles_afx <- afex::aov_4(
  attractiveness ~ facetype * alcohol + (1 | id),
  data = discovr::goggles,
  anova_table = list(p_adjust_method = "bonferroni")
)
#> Contrasts set to contr.sum for the following variables: facetype, alcohol

goggles_afx
#> Anova Table (Type 3 tests, bonferroni-adjusted)
#> 
#> Response: attractiveness
#>             Effect    df  MSE         F  ges p.value
#> 1         facetype 1, 42 1.37 15.58 *** .271   <.001
#> 2          alcohol 2, 42 1.37    6.04 * .223    .015
#> 3 facetype:alcohol 2, 42 1.37   8.51 ** .288    .002
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '+' 0.1 ' ' 1

model_parameters(goggles_afx)
#> Error in model_parameters(goggles_afx): could not find function "model_parameters"


a2 <- afex::aov_ez("id", "value", afex::obk.long, between = c("treatment", "gender"), 
                   within = c("phase", "hour"), observed = "gender",
                   anova_table = list(p_adjust_method = "fdr"))
#> Contrasts set to contr.sum for the following variables: treatment, gender

a2
#> Anova Table (Type 3 tests, fdr-adjusted)
#> 
#> Response: value
#>                         Effect          df   MSE        F  ges p.value
#> 1                    treatment       2, 10 22.81     3.94 .198    .205
#> 2                       gender       1, 10 22.81     3.66 .115    .254
#> 3             treatment:gender       2, 10 22.81     2.86 .179    .261
#> 4                        phase 1.60, 15.99  5.02 16.13 ** .151    .002
#> 5              treatment:phase 3.20, 15.99  5.02   4.85 + .097    .063
#> 6                 gender:phase 1.60, 15.99  5.02     0.28 .003    .818
#> 7       treatment:gender:phase 3.20, 15.99  5.02     0.64 .014    .808
#> 8                         hour 1.84, 18.41  3.39 16.69 ** .125    .001
#> 9               treatment:hour 3.68, 18.41  3.39     0.09 .002    .979
#> 10                 gender:hour 1.84, 18.41  3.39     0.45 .004    .808
#> 11       treatment:gender:hour 3.68, 18.41  3.39     0.62 .011    .808
#> 12                  phase:hour 3.60, 35.96  2.67     1.18 .015    .717
#> 13        treatment:phase:hour 7.19, 35.96  2.67     0.35 .009    .979
#> 14           gender:phase:hour 3.60, 35.96  2.67     0.93 .012    .808
#> 15 treatment:gender:phase:hour 7.19, 35.96  2.67     0.74 .019    .808
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '+' 0.1 ' ' 1
#> 
#> Sphericity correction method: GG

parameters::model_parameters(a2)
#> ANOVA estimation for factorial designs using 'afex'
#> 
#> Parameter                   | Sum_Squares | Sum_Squares_Error |   df
#> --------------------------------------------------------------------
#> treatment                   |      179.73 |            228.06 | 2.00
#> gender                      |       83.45 |            228.06 | 1.00
#> treatment:gender            |      130.24 |            228.06 | 2.00
#> phase                       |      129.51 |             80.28 | 1.60
#> treatment:phase             |       77.89 |             80.28 | 3.20
#> gender:phase                |        2.27 |             80.28 | 1.60
#> treatment:gender:phase      |       10.22 |             80.28 | 3.20
#> hour                        |      104.29 |             62.50 | 1.84
#> treatment:hour              |        1.17 |             62.50 | 3.68
#> gender:hour                 |        2.81 |             62.50 | 1.84
#> treatment:gender:hour       |        7.76 |             62.50 | 3.68
#> phase:hour                  |       11.35 |             96.17 | 3.60
#> treatment:phase:hour        |        6.64 |             96.17 | 7.19
#> gender:phase:hour           |        8.96 |             96.17 | 3.60
#> treatment:gender:phase:hour |       14.15 |             96.17 | 7.19
#> 
#> Parameter                   | df (error) | Mean_Square |     F |     p
#> ----------------------------------------------------------------------
#> treatment                   |      10.00 |       22.81 |  3.94 | 0.205
#> gender                      |      10.00 |       22.81 |  3.66 | 0.254
#> treatment:gender            |      10.00 |       22.81 |  2.86 | 0.261
#> phase                       |      15.99 |        5.02 | 16.13 | 0.002
#> treatment:phase             |      15.99 |        5.02 |  4.85 | 0.063
#> gender:phase                |      15.99 |        5.02 |  0.28 | 0.818
#> treatment:gender:phase      |      15.99 |        5.02 |  0.64 | 0.808
#> hour                        |      18.41 |        3.39 | 16.69 | 0.001
#> treatment:hour              |      18.41 |        3.39 |  0.09 | 0.979
#> gender:hour                 |      18.41 |        3.39 |  0.45 | 0.808
#> treatment:gender:hour       |      18.41 |        3.39 |  0.62 | 0.808
#> phase:hour                  |      35.96 |        2.67 |  1.18 | 0.717
#> treatment:phase:hour        |      35.96 |        2.67 |  0.35 | 0.979
#> gender:phase:hour           |      35.96 |        2.67 |  0.93 | 0.808
#> treatment:gender:phase:hour |      35.96 |        2.67 |  0.74 | 0.808
#> 
#> p-value adjustment method: Benjamini & Hochberg (1995)
#> Anova Table (Type 3 tests)

^{Created on 2025-05-08 with reprex v2.1.1}

[profandyfield]

would you like me to test it on a few book examples?

[strengejacke]

Sure, you can install from the PR, if you like.

[profandyfield]

For consistency, would it be a good idea to have p_adjust also work on car::Anova objects?

Example:

goggles_tib <- discovr::goggles
contrasts(goggles_tib$facetype) <- contr.sum(n = 2)
contrasts(goggles_tib$alcohol) <- contr.sum(n = 3)

goggles_lm <- lm(attractiveness ~ facetype*alcohol, data = goggles_tib)

car::Anova(goggles_lm, type = 3) |> 
  parameters::model_parameters(es_type = "omega", ci = 0.95) |> 
  insight::display()
#> Type 3 ANOVAs only give sensible and informative results when covariates
#>   are mean-centered and factors are coded with orthogonal contrasts (such
#>   as those produced by `contr.sum`, `contr.poly`, or `contr.helmert`, but
#>   *not* by the default `contr.treatment`).

Parameter	Sum_Squares	df	Mean_Square	F	p	Omega2 (partial)	Omega2 95% CI
facetype	21.33	1	21.33	15.58	< .001	0.23	(0.07, 1.00)
alcohol	16.54	2	8.27	6.04	0.005	0.17	(0.02, 1.00)
facetype:alcohol	23.29	2	11.65	8.51	< .001	0.24	(0.06, 1.00)
Residuals	57.50	42	1.37

Anova Table (Type 3 tests)

^{Created on 2025-05-08 with reprex v2.1.1}

[mattansb]

Good idea:

goggles_tib <- discovr::goggles
contrasts(goggles_tib$facetype) <- contr.sum(n = 2)
contrasts(goggles_tib$alcohol) <- contr.sum(n = 3)

goggles_lm <- lm(attractiveness ~ facetype*alcohol, data = goggles_tib)

car::Anova(goggles_lm, type = 3) |> 
  parameters::model_parameters(es_type = "omega", ci = 0.95)
#> Type 3 ANOVAs only give sensible and informative results when covariates
#>   are mean-centered and factors are coded with orthogonal contrasts (such
#>   as those produced by `contr.sum`, `contr.poly`, or `contr.helmert`, but
#>   *not* by the default `contr.treatment`).
#> Parameter        | Sum_Squares | df | Mean_Square |     F |      p
#> ------------------------------------------------------------------
#> facetype         |       21.33 |  1 |       21.33 | 15.58 | < .001
#> alcohol          |       16.54 |  2 |        8.27 |  6.04 | 0.005 
#> facetype:alcohol |       23.29 |  2 |       11.65 |  8.51 | < .001
#> Residuals        |       57.50 | 42 |        1.37 |       |       
#> 
#> Parameter        | Omega2 (partial) | Omega2 95% CI
#> ---------------------------------------------------
#> facetype         |             0.23 |  [0.07, 1.00]
#> alcohol          |             0.17 |  [0.02, 1.00]
#> facetype:alcohol |             0.24 |  [0.06, 1.00]
#> Residuals        |                  |              
#> 
#> Anova Table (Type 3 tests)

car::Anova(goggles_lm, type = 3) |> 
  parameters::model_parameters(es_type = "omega", ci = 0.95, p_adjust = "bonferroni")
#> Type 3 ANOVAs only give sensible and informative results when covariates
#>   are mean-centered and factors are coded with orthogonal contrasts (such
#>   as those produced by `contr.sum`, `contr.poly`, or `contr.helmert`, but
#>   *not* by the default `contr.treatment`).
#> Parameter        | Sum_Squares | df | Mean_Square |     F |      p
#> ------------------------------------------------------------------
#> facetype         |       21.33 |  1 |       21.33 | 15.58 | < .001
#> alcohol          |       16.54 |  2 |        8.27 |  6.04 | 0.015 
#> facetype:alcohol |       23.29 |  2 |       11.65 |  8.51 | 0.002 
#> Residuals        |       57.50 | 42 |        1.37 |       |       
#> 
#> Parameter        | Omega2 (partial) | Omega2 95% CI
#> ---------------------------------------------------
#> facetype         |             0.23 |  [0.07, 1.00]
#> alcohol          |             0.17 |  [0.02, 1.00]
#> facetype:alcohol |             0.24 |  [0.06, 1.00]
#> Residuals        |                  |              
#> 
#> Anova Table (Type 3 tests)


Ȁ
#> Error: object 'Ȁ' not found

^{Created on 2025-05-08 with reprex v2.1.1}