Is it possible to do standard contrasts across an interaction with the default-backend?
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.
@mattansb @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
@profandyfield how do you interpret this output?
@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?'.
@strengejacke thanks - I didn't realise that estimate_contrasts() would accept the interaction = as an argument so that's useful to know.
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?
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
If you want I can write an example?
Yes, if possible, using the function. That should directly work with modelbased as well.
@strengejacke thanks - I didn't realise that
estimate_contrasts()would accept theinteraction =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.
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?
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
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
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.639The 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?'
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() 😆
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
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.
Ok, what's the state of this issue? Can it be closed?
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`.
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.