# interactions involving two continuous predictors with emtrends

Let's say I have a model with two continuous predictors (nitrogen and temperature) and one categorical variable (variety).

dat <- data.frame(blocks=rep(c(1:15),each=2), variety=rep(c("A","B","C"),each=10),soil=runif(30,0,10), nitro = runif(30, 0, 10), temp= rnorm(30, 10, 3));

mod <- lmer(soil ~ variety*nitro*temp + (1|blocks), data=dat)


Let's pretend the variety:nitro:temp interaction is significant. How can I calculate the trend involving these two continuous predictors? For example maybe I believe that the slope between soil and nitrogen will be steeper at increasingly higher temperatures only in variety A. Or, in other words, how can I test if the slope between soil and nitrogen changes as a function of temperature in variety A?

I thought this would give me that information

> emtrends(mod, pairwise~variety, var=c("nitro","temp"))


but it throws an error

Error in .subset2(x, i, exact = exact) : subscript out of bounds


The var argument specifies the variable whose slope you are interested in. Thus it can be only one character string, not a vector; in your case nitro.

The at argument allows you to specify values of other variables for which you wish do to tests and comparisons. at is a named list. In your case, to address the question you ask that would be:

emtrends(mod, pairwise ~ temp, var="nitro", at=list(variety="A", temp=c(20,40))


The above line tests whether the trend of nitro is different between 20° and 40° in variety A. I picked 20° and 40° arbitrarily; a common choice is $$mean-SD$$ and $$mean+SD$$, or the medians of the first and third terciles.

• Thanks @Ous, your answer was helpful. However can't you have a regression slope for each categorical predictor that reflects the relationship between the continuous predictors and the DV? For example, I can imagine a regression line in which the DV is in y axis, nitro in x axis and temp in z axis. Each variety would then get its own slope in this 3D plane. It would then just be a matter of comparing the slope for variety A with the slope for variety B for example Jun 22, 2019 at 23:34
• I think this might be possible to test with the summary() function if the contrasts are done right, but I thought it would be easier with emtrends somehow. Jun 22, 2019 at 23:35
• Actually you would get a 2D plane for the interaction, not a line. The slope of it would be of dimension 2, and I have no idea how you might contrasts such slopes. Imagine you had a 3-way interaction with categorical variables. How would you proceed for post hoc tests using emmeans?
– Ous
Jun 23, 2019 at 0:39
• You are right, it would be a plane not a line. But if you look at the output of summary(mod) you get varietyB:nitro:temp ... 1.75 0.09 I read this like: how different is the nitro:temp interaction for varietyB from the nitro:temp interaction for varietyA (the reference). So here you could potentially calculate whether this interaction is significantly different between the two varieties. Or am I misunderstanding this output? Jun 23, 2019 at 22:59
• Also, the code you provided returns the same t.ratio and p.value regardless of what contrast I specify. Changing the code from at=list(variety="A", temp=c(20,40)) to at=list(variety="A", temp=c(20,30)) results in the same t.ratio and p.value Jun 23, 2019 at 23:01