modelbased icon indicating copy to clipboard operation
modelbased copied to clipboard

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

Open profandyfield opened this issue 6 months ago • 19 comments

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.

profandyfield avatar Jun 11 '25 09:06 profandyfield

@mattansb @strengejacke

DominiqueMakowski avatar Jun 11 '25 09:06 DominiqueMakowski

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

strengejacke avatar Jun 11 '25 09:06 strengejacke

@profandyfield how do you interpret this output?

DominiqueMakowski avatar Jun 11 '25 09:06 DominiqueMakowski

@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 avatar Jun 11 '25 10:06 profandyfield

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

profandyfield avatar Jun 11 '25 10:06 profandyfield

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?

mattansb avatar Jun 11 '25 10:06 mattansb

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

DominiqueMakowski avatar Jun 11 '25 10:06 DominiqueMakowski

If you want I can write an example?

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

strengejacke avatar Jun 11 '25 11:06 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

Image

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 avatar Jun 11 '25 11:06 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?

strengejacke avatar Jun 11 '25 12:06 strengejacke

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

mattansb avatar Jun 11 '25 12:06 mattansb

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

strengejacke avatar Jun 11 '25 12:06 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:

> 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 avatar Jun 11 '25 14:06 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() 😆

profandyfield avatar Jun 11 '25 14:06 profandyfield

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

DominiqueMakowski avatar Jun 11 '25 14:06 DominiqueMakowski

From a user perspective, what I would have found intuitive (but maybe not anyone else!) is something like contrast = looks*personality, comparison = list(strategy = "trt.vs.ctrl", looks = "trt.vs.ctrl")), reference = c(2, 2)

So, you set the contrast as being the interaction of interest, you use comparison to set standard contrasts for each variable (or specify one contrast type that is then applied to all variables in the interaction) and you can optionally suppply a vector of reference categories for the variables in the interaction.

profandyfield avatar Jun 11 '25 14:06 profandyfield

Ok, what's the state of this issue? Can it be closed?

strengejacke avatar Jun 30 '25 10:06 strengejacke

From the perspective of my book, I switched the backend. However, I think you should consider making a user-friendly interface to do it with marginaleffects. So, you can probably close this specific issue, but I'd argue you should add it to the 'easystats team feature enhancement discussion`.

profandyfield avatar Jun 30 '25 10:06 profandyfield

Ok, I'll leave it open. To which extent we will be able to support this with our default beackend also depends on how marginaleffects implements this feature.

strengejacke avatar Jun 30 '25 10:06 strengejacke