Interpretation of interaction effect of mixed linear models

Question

for my master thesis i worked with mixed linear models for the first time. I am trying to get results concerning the effect of tree ring width on various wood anatomical features of two different tree species. I formulated my model in R like this:

> lmer(CD.~scale(MRW)*species+(1|ID)+(1|Year), data= Alles4)

with CD. being an anatomical feature, MRW the tree ring width, ID the individual tree and Year the year of tree ring formation. Species is a categorical variable with only the two tree species.

I got the following results:

Fixed effects:
                     Estimate Std. Error        df t value Pr(>|t|)    
(Intercept)          222.9791     4.9012   37.6322  45.494  < 2e-16 ***
scale(MRW)           -16.2370     0.7431 2373.4519 -21.850  < 2e-16 ***
speciesS             -22.8443     6.7561   34.0283  -3.381  0.00183 ** 
scale(MRW):speciesS    4.1713     0.9364 2352.8566   4.455 8.79e-06 ***

I think I am clear in the understanding of the fixed effects by themselves, however the interpretation of the interaction between these two I am not sure about. I understand that there is a significant interaction, but how do I report it?

At species=1 each step of MRW decreases CD. by 16,2370? So at species=2 each step of MRW decreases CD. by 16,2370-4,1713=12,0657, hence CD. of species 1 is more effected by MRW than CD. of species 2 and this effect is significant?

Answer 1

To sum up my comments, it is always a good idea to (1) write the regression equation with estimated parameters of the regression model in order to interpret the model, and (2) draw an interaction plot.

1. The regression equation tells you that the average response variable, CD., might be predicted as $\hat\beta_0 + \hat\beta_1\text{MRW} + \hat\beta_2\mathbb{I}(\text{Species}=S) + \hat\beta_3\text{MRW}\times \mathbb{I}(\text{Species}=S)$, where $\mathbb{I}(\text{Species}=S)$ is the indicator variable which equals 1 when Species = S, 0 otherwise, and $\hat\beta_j$ are the estimated regression coefficients displayed in your R's output, in the same order. The fact that this is a mixed-effect model doesn't change the interpretation of the "fixed" part of the model, only the predictions. From here, you can put whichever number you want into that equation and you will get the final estimate. Since the interaction is significant, you cannot deduce the difference between the two species from the value of $\hat\beta_2$ only, and details about how to interpret an interaction between two continuous variables using the same notation is available in a related thread. The interpretation remains the same since you are using a binary variable. Look for "interaction plot and ANCOVA" and you will probably find many other related threads. When there is no interaction, the expected value of CD. only depends on the value of Species, which means that if you were to plot the predicted values of CD. as a function of MRW, this would yield two parallel lines. The slightest deviation from this parallelism supposes the idea of an interaction, which might be significant or not. If not, you can summarize your result by saying that there's an overall shift between the two groups, for an amount of $\hat\beta_2$. In your case, the interaction term tells you how MRW and Species relate to the variation of CD., in that for Species S you need to add an extra amount of 4.17 for one SD unit change in MRW. When the interaction is significant, however, you need to tell whether the effects are antagonistic or not (i.e., whether the lines cross or not). In your case again, it seems to be an ordered interaction effect.

2. More than a thousand of words, an interaction plot will help summarizing the results. Below is an example with simulated data, not accounting for random effects. It would be best to "predict" the value of CD. from the mixed-effect model as it will take into account the conditional variance of the random effects. To sum up, the interaction term tells you how much the results depart from the parallel slope assumption, if we were to assume the effect of MRW is constant across species.

---

n <- 100
b0 <- 222.9791
b1 <- -16.2370
b2 <- -22.8443
b3 <- 4.1713

A <- runif(n, 0, 50)
B <- rep(0:1, each = n/2)
e <- rnorm(n, 0, sd=1)

y <- b0 + b1*scale(A) + b2*B + b3*scale(A)*B + e

d <- data.frame(A, B = factor(B, labels = paste0("S", 1:2)), y)

m <- lm(y ~ A * B, data = d)
summary(m)

d$yhat <- predict(m)

library(ggplot2)
library(hrbrthemes)
theme_set(theme_ipsum(base_size = 11))

p <- ggplot(data = d, aes(x = scale(A), y = yhat, color = B))
p <- p + geom_point(aes(x = scale(A), y = y, color = B))
p <- p + geom_line(aes(group = B))
p <- p + guides(col = guide_legend(title = ""))
p <- p + labs(x = "MRW (scaled with unit variance)", y = "CD.", caption = "Simulated data (n=100)")
p <- p + theme(legend.position="top")