# Emmeans and lmer with quadratic and cubic interaction

I have an experiment where I predict that the two levels of my Cond variable (effect coded as -1 and 1) will have different trends over time on my EDA_cs outcome variable (skin conductance). Specifically, I expect my -1 to be s-shaped (cubic) and my 1 to be u-shaped (quadratic). I am testing this with the following mixed-effect model:

m.EDA = lmer(EDA_cs ~ Cond * Time_sd + Cond*I(Time_sd^2) + Cond *I(Time_sd^3) + (1|Participant) + (1|Stim),data= phys)


All three interactions are significant. When I plot this with sjPlot's package using the following code:

plot_model(m.EDA.H2, type = "pred", terms = c("Time_sd[all]", "Cond"))


I get this beautiful graph confirming the hypothesis (yay):

Now the question is:

How can I further probe that this is difference in shape is significant? Is there some way to use emmeans() at different time_sd points so that I can show where the blue line is higher or lower than the red line? I want something like this code:

emmeans(m.EDA.H2,pairwise~ ~Cond|Time_sd, at = list(Time_sd = c(-1, 0, 1)))


But that, instead of giving the estimated mean difference for the linear Time_sd, gives the estimated mean difference from the cubic interaction.

Does such thing exist?

From the graph etc., I'd suggest using a lot more time points, and separating the means from the comparisons, e.g.,

EMM <- emmeans(m.EDA.H2,~ ~Cond|Time_sd,
at = list(Time_sd = seq(-2, 2, by = 0.1)))
emmip(EMM, Cond ~ Time_sd)   # display the curves
test(contrast(EMM, "consec"), by = NULL, adjust = "none")


The last statement obtains the differences, then summarixes them in a more compact form (rather than a separate table for each time point).

• thank you so much! Commented Oct 9, 2023 at 13:48
• just to understand better: I was wrong, the code I used does not test for the linear but for the cubic right? the only issue is the breaks I chose, because it was only 3 points and all in the same direction it made me think it was testing for a line. Thanks a lot for the solution. Commented Oct 9, 2023 at 13:51
• No. You spec is ~ Cond | time_sd, so you are getting estimated marginal means for Cond and time_sd is your by variable, so you have Cond means at three time_sd values. Commented Oct 16, 2023 at 15:38