Is it possible to do standard contrasts across an interaction with the default-backend?

[profandyfield]

I'm trying to reproduce this analysis using modelbased::estimate_contrasts(). It's basically a set of 'out of the box' contrasts fitted across an interaction.

date_tib <- discovr::speed_date

date_afx <- afex::aov_4(date ~ strategy*looks*personality + (looks*personality|id), data = date_tib)
#> Contrasts set to contr.sum for the following variables: strategy

# Desired analysis using emmeans

three_way_emm <- emmeans::emmeans(date_afx, specs = c("strategy", "looks", "personality"), model = "multivariate")

emmeans::contrast(
  three_way_emm,
  interaction = c(strategy = "trt.vs.ctrl",
                  looks = "trt.vs.ctrl",
                  personality = "trt.vs.ctrl"),
  ref = 2,
  adjust = "holm"
  )
#>  strategy_trt.vs.ctrl looks_trt.vs.ctrl personality_trt.vs.ctrl estimate   SE
#>  Normal - Hard to get Low - Average     Low - Average               -4.7 4.07
#>  Normal - Hard to get High - Average    Low - Average              -36.2 4.65
#>  Normal - Hard to get Low - Average     High - Average              18.5 5.41
#>  Normal - Hard to get High - Average    High - Average              -5.8 6.02
#>  df t.ratio p.value
#>  18  -1.155  0.5263
#>  18  -7.789  <.0001
#>  18   3.421  0.0091
#>  18  -0.963  0.5263
#> 
#> P value adjustment: holm method for 4 tests

^{Created on 2025-06-11 with reprex v2.1.1}

I can't see a way to set the contrasts across an interaction in estimate_contrasts() (without manually specifying weights but let's not go there ...). The function documentation doesn't mention an equivalent to interactiuon = that gets passed onto emmeans. But also, highly likley I just don't know how to use your function:) Is there a way? If not, is this something you'd impliment?

Incidentally, I'd prefer to do it without changing the backend to emmeans (because everything else in my book uses the default marginaleffects backend), but if it has to be changed then so be it.

[DominiqueMakowski]

@mattansb @strengejacke

[strengejacke]

The modelbased-way is indeed to switch to emmeans backend, because this is a emmeans-specific syntax. Not sure at the moment how to mimic that with the marginaleffects package?

You could use:

library(easystats)
#> # Attaching packages: easystats 0.7.4.4
#> ✔ bayestestR  0.16.0     ✔ correlation 0.8.7.1 
#> ✔ datawizard  1.1.0.1    ✔ effectsize  1.0.1   
#> ✔ insight     1.3.0.2    ✔ modelbased  0.11.2  
#> ✔ performance 0.14.0.2   ✔ parameters  0.26.0.1
#> ✔ report      0.6.1      ✔ see         0.11.0.3

date_tib <- discovr::speed_date
date_afx <- afex::aov_4(date ~ strategy*looks*personality + (looks*personality|id), data = date_tib)
#> Contrasts set to contr.sum for the following variables: strategy

estimate_contrasts(
  date_afx,
  c("strategy", "looks", "personality"),
  interaction = c(strategy = "trt.vs.ctrl", looks = "trt.vs.ctrl", personality = "trt.vs.ctrl"),
  ref = 2,
  p_adjust = "holm",
  backend = "emmeans"
)
#> Marginal Contrasts Analysis
#> 
#> strategy_trt.vs.ctrl | looks_trt.vs.ctrl | personality_trt.vs.ctrl | Difference
#> -------------------------------------------------------------------------------
#> Normal - Hard to get | High - Average    | High - Average          |      -5.80
#> Normal - Hard to get | High - Average    | Low - Average           |     -36.20
#> Normal - Hard to get | Low - Average     | High - Average          |      18.50
#> Normal - Hard to get | Low - Average     | Low - Average           |      -4.70
#> 
#> strategy_trt.vs.ctrl |           95% CI |   SE | t(18) |      p
#> ---------------------------------------------------------------
#> Normal - Hard to get | [-22.50,  10.90] | 6.02 | -0.96 |  0.526
#> Normal - Hard to get | [-49.09, -23.31] | 4.65 | -7.79 | < .001
#> Normal - Hard to get | [  3.49,  33.51] | 5.41 |  3.42 |  0.009
#> Normal - Hard to get | [-15.99,   6.59] | 4.07 | -1.16 |  0.526
#> 
#> Variable predicted: date
#> Predictors contrasted: strategy, looks, personality
#> p-value adjustment method: Holm (1979)

^{Created on 2025-06-11 with reprex v2.1.1}

[DominiqueMakowski]

@profandyfield how do you interpret this output?

[profandyfield]

@profandyfield how do you interpret this output?

@DominiqueMakowski you swivel your chair round to the shelves in our office, grab a copy of Discovering Statuistics Using SPSS 6th edition and read pages 773-775 😉

Essentially, one interpretation of each contrast is that you're asking 'is the pattern of means across the specified combinations of looks and personality the same in the two strategy groups?'.

[profandyfield]

@strengejacke thanks - I didn't realise that estimate_contrasts() would accept the interaction = as an argument so that's useful to know.

[mattansb]

Interaction contrasts are possible to with marginaleffects, but are still quite easier to do with emmeans.

In marginaleffects you would either need to write a little function to compute those contrasts for you, or you would need to chain hypotheses() |> hypotheses() one after the other. If you want I can write an example?

[DominiqueMakowski]

Essentially, one interpretation of each contrast is that you're asking 'is the pattern of means across the specified combinations of looks and personality the same in the two strategy groups?'.

This sounds similar to what I would say for this though:

> modelbased::estimate_contrasts(date_afx, contrast = "strategy", by = c("looks", "personality"))
Marginal Contrasts Analysis

Level1      | Level2 | looks   | personality | Difference |   SE |           95% CI | t(161) |      p
-----------------------------------------------------------------------------------------------------
Hard to get | Normal | Low     | Low         |      -0.30 | 1.49 | [ -3.25,   2.65] |  -0.20 |  0.841
Hard to get | Normal | Average | Low         |       0.80 | 1.78 | [ -2.71,   4.31] |   0.45 |  0.653
Hard to get | Normal | High    | Low         |      35.50 | 2.04 | [ 31.48,  39.52] |  17.42 | < .001
Hard to get | Normal | Low     | Average     |      -2.90 | 2.42 | [ -7.68,   1.88] |  -1.20 |  0.233
Hard to get | Normal | Average | Average     |       2.90 | 2.34 | [ -1.73,   7.53] |   1.24 |  0.218
Hard to get | Normal | High    | Average     |       1.40 | 2.82 | [ -4.16,   6.96] |   0.50 |  0.620
Hard to get | Normal | Low     | High        |     -29.90 | 2.50 | [-34.83, -24.97] | -11.97 | < .001
Hard to get | Normal | Average | High        |      -5.60 | 3.44 | [-12.40,   1.20] |  -1.63 |  0.106
Hard to get | Normal | High    | High        |      -1.30 | 2.77 | [ -6.76,   4.16] |  -0.47 |  0.639

The difference between strategies for each combination stratum

[strengejacke]

If you want I can write an example?

Yes, if possible, using the function. That should directly work with modelbased as well.

[strengejacke]

@strengejacke thanks - I didn't realise that estimate_contrasts() would accept the interaction = as an argument so that's useful to know.

See docs here (which were updated in the last package release on CRAN):
https://easystats.github.io/modelbased/reference/estimate_contrasts.html

We could add some examples to the docs / elaborate a bit more, so it's more obvious that you can indeed replicate a lot of the original behaviours of the emmeans and marginaleffects package when you use modelbased.

[strengejacke]

Essentially, one interpretation of each contrast is that you're asking 'is the pattern of means across the specified combinations of looks and personality the same in the two strategy groups?'.

This sounds similar to what I would say for this though:

but you don't have the contrasts for looks and personality here, just the levels. But even if you do all pairwise comparisons, it's not the same, not sure exactly what kind of comparison is done here?

[mattansb]

I guess this is how I would do a one-off thing:

marginaleffects::avg_predictions(date_afx, 
  by = c("strategy", "looks", "personality"),

  hypothesis = \(x) {
    x |> 
      dplyr::reframe(
        strategy_contrast = "Normal - Hard to get",
        estimate = estimate[strategy %in% c("Normal")] - estimate[strategy %in% c("Hard to get")],
        
        .by = c(looks, personality)
      ) |>
      dplyr::reframe(
        looks_contrast = c("Low - Average", "High - Average"),
        estimate = estimate[looks %in% c("Low", "High")] - estimate[looks %in% c("Average")],
        
        .by = c(personality, strategy_contrast)
      ) |> 
      dplyr::reframe(
        personality_contrast = c("Low - Average", "High - Average"),
        estimate = estimate[personality %in% c("Low", "High")] - estimate[personality %in% c("Average")],
        
        .by = c(strategy_contrast, looks_contrast)
      ) |> 
      dplyr::mutate(
        term = paste(strategy_contrast, looks_contrast, personality_contrast, sep = " | ")
      )
  })
#> 
#>                                                    Term Estimate Std. Error      z Pr(>|z|)    S 2.5 % 97.5 %
#>  Normal - Hard to get | Low - Average | Low - Average       -4.7       4.07 -1.155    0.248  2.0 -12.7   3.28
#>  Normal - Hard to get | Low - Average | High - Average      18.5       5.41  3.421   <0.001 10.6   7.9  29.10
#>  Normal - Hard to get | High - Average | Low - Average     -36.2       4.65 -7.789   <0.001 47.1 -45.3 -27.09
#>  Normal - Hard to get | High - Average | High - Average     -5.8       6.02 -0.963    0.335  1.6 -17.6   6.00
#> 
#> Type: response

[strengejacke]

Ok, this can also be used in modelbased - though it looks slightly more complicated compared to emmeans:

library(easystats)

date_tib <- discovr::speed_date
date_afx <- afex::aov_4(date ~ strategy*looks*personality + (looks*personality|id), data = date_tib)
#> Contrasts set to contr.sum for the following variables: strategy

estimate_means(date_afx, 
  by = c("strategy", "looks", "personality"),
  hypothesis = \(x) {
    x |> 
      dplyr::reframe(
        strategy_contrast = "Normal - Hard to get",
        estimate = estimate[strategy %in% c("Normal")] - estimate[strategy %in% c("Hard to get")],
        
        .by = c(looks, personality)
      ) |>
      dplyr::reframe(
        looks_contrast = c("Low - Average", "High - Average"),
        estimate = estimate[looks %in% c("Low", "High")] - estimate[looks %in% c("Average")],
        
        .by = c(personality, strategy_contrast)
      ) |> 
      dplyr::reframe(
        personality_contrast = c("Low - Average", "High - Average"),
        estimate = estimate[personality %in% c("Low", "High")] - estimate[personality %in% c("Average")],
        
        .by = c(strategy_contrast, looks_contrast)
      ) |> 
      dplyr::mutate(
        term = paste(strategy_contrast, looks_contrast, personality_contrast, sep = " | ")
      )
  }) |> 
    # remove redundant colums
    data_select(exclude = c("strategy_contrast", "looks_contrast", "personality_contrast")) |> 
    # dont wrap table
    print(table_width = Inf)
#> Warning in summary.Anova.mlm(object$Anova, multivariate = FALSE): HF eps > 1
#> treated as 1
#> Estimated Marginal Means
#> 
#> Parameter                                              | Difference |   SE |           95% CI | t(161) |      p
#> ---------------------------------------------------------------------------------------------------------------
#> Normal - Hard to get | Low - Average | Low - Average   |      -4.70 | 4.07 | [-12.74,   3.34] |  -1.16 |  0.250
#> Normal - Hard to get | Low - Average | High - Average  |      18.50 | 5.41 | [  7.82,  29.18] |   3.42 | < .001
#> Normal - Hard to get | High - Average | Low - Average  |     -36.20 | 4.65 | [-45.38, -27.02] |  -7.79 | < .001
#> Normal - Hard to get | High - Average | High - Average |      -5.80 | 6.02 | [-17.69,   6.09] |  -0.96 |  0.337
#> 
#> Variable predicted: date
#> Predictors modulated: strategy, looks, personality
#> Predictors averaged: id

^{Created on 2025-06-11 with reprex v2.1.1}

[profandyfield]

Essentially, one interpretation of each contrast is that you're asking 'is the pattern of means across the specified combinations of looks and personality the same in the two strategy groups?'.

This sounds similar to what I would say for this though:

> modelbased::estimate_contrasts(date_afx, contrast = "strategy", by = c("looks", "personality"))
Marginal Contrasts Analysis

Level1      | Level2 | looks   | personality | Difference |   SE |           95% CI | t(161) |      p
-----------------------------------------------------------------------------------------------------
Hard to get | Normal | Low     | Low         |      -0.30 | 1.49 | [ -3.25,   2.65] |  -0.20 |  0.841
Hard to get | Normal | Average | Low         |       0.80 | 1.78 | [ -2.71,   4.31] |   0.45 |  0.653
Hard to get | Normal | High    | Low         |      35.50 | 2.04 | [ 31.48,  39.52] |  17.42 | < .001
Hard to get | Normal | Low     | Average     |      -2.90 | 2.42 | [ -7.68,   1.88] |  -1.20 |  0.233
Hard to get | Normal | Average | Average     |       2.90 | 2.34 | [ -1.73,   7.53] |   1.24 |  0.218
Hard to get | Normal | High    | Average     |       1.40 | 2.82 | [ -4.16,   6.96] |   0.50 |  0.620
Hard to get | Normal | Low     | High        |     -29.90 | 2.50 | [-34.83, -24.97] | -11.97 | < .001
Hard to get | Normal | Average | High        |      -5.60 | 3.44 | [-12.40,   1.20] |  -1.63 |  0.106
Hard to get | Normal | High    | High        |      -1.30 | 2.77 | [ -6.76,   4.16] |  -0.47 |  0.639

The difference between strategies for each combination stratum

But this is looking at the effect of strategy in individual combinations of the other variables, so you're never capturing relative effects (and, sure, lots of situations where you wouldn't care about relative effects and the above tells you what you want). The contrasts in my original post capture effects relative to their controls, so it's more akin to doing the above for difference scores. To simplify a bit, let's ignore personality. You'd end up with two contrasts using the method in my original post, that represent:

• hard to get vs normal for (low looks - average looks)
• hard to get vs normal for (high looks - average looks)

If you do your thing above you'd have 3 contrasts

• hard to get vs normal for low looks
• hard to get vs normal for average looks
• hard to get vs normal for high looks

The first method captures effect relative to a control, the second method does not. Neither is better, but they are different. First method answers questions like 'do hard to get and normal differ on the difference score between low and average looks conditions?' the second answers questions like 'do hard to get and normal differ on scores within the average looks condition?'

[profandyfield]

Ok, this can also be used in modelbased - though it looks slightly more complicated compared to emmeans:

library(easystats)

date_tib <- discovr::speed_date
date_afx <- afex::aov_4(date ~ strategylookspersonality + (looks*personality|id), data = date_tib)
#> Contrasts set to contr.sum for the following variables: strategy

estimate_means(date_afx,
by = c("strategy", "looks", "personality"),
hypothesis = (x) {
x |>
dplyr::reframe(
strategy_contrast = "Normal - Hard to get",
estimate = estimate[strategy %in% c("Normal")] - estimate[strategy %in% c("Hard to get")],

    .by = c(looks, personality)
  ) |>
  dplyr::reframe(
    looks_contrast = c("Low - Average", "High - Average"),
    estimate = estimate[looks %in% c("Low", "High")] - estimate[looks %in% c("Average")],
    
    .by = c(personality, strategy_contrast)
  ) |> 
  dplyr::reframe(
    personality_contrast = c("Low - Average", "High - Average"),
    estimate = estimate[personality %in% c("Low", "High")] - estimate[personality %in% c("Average")],
    
    .by = c(strategy_contrast, looks_contrast)
  ) |> 
  dplyr::mutate(
    term = paste(strategy_contrast, looks_contrast, personality_contrast, sep = " | ")
  )

}) |>
# remove redundant colums
data_select(exclude = c("strategy_contrast", "looks_contrast", "personality_contrast")) |>
# dont wrap table
print(table_width = Inf)
#> Warning in summary.Anova.mlm(object$Anova, multivariate = FALSE): HF eps > 1
#> treated as 1
#> Estimated Marginal Means
#>
#> Parameter | Difference | SE | 95% CI | t(161) | p
#> ---------------------------------------------------------------------------------------------------------------
#> Normal - Hard to get | Low - Average | Low - Average | -4.70 | 4.07 | [-12.74, 3.34] | -1.16 | 0.250
#> Normal - Hard to get | Low - Average | High - Average | 18.50 | 5.41 | [ 7.82, 29.18] | 3.42 | < .001
#> Normal - Hard to get | High - Average | Low - Average | -36.20 | 4.65 | [-45.38, -27.02] | -7.79 | < .001
#> Normal - Hard to get | High - Average | High - Average | -5.80 | 6.02 | [-17.69, 6.09] | -0.96 | 0.337
#>
#> Variable predicted: date
#> Predictors modulated: strategy, looks, personality
#> Predictors averaged: id

Created on 2025-06-11 with reprex v2.1.1

for everyone's sanity, I'll just change the backend in the textbook and pass on the interaction argument within estimate_contrast() 😆

[DominiqueMakowski]

I get the idea of differences of differences (interaction of contrasts), but I find the emmeans output quite confusing and hard to correctly interpret to be fair.

Maybe we should think of an API to enable this somewhat natively. Perhaps by allowing multiple arguments in contrast? We could at first throw a warning and automatically switch to emmeans if we don't support that with margianleffects