Interaction term in a linear mixed effect model in R

Question

I am attempting to analyze the effect of two categorical variables (landuse and species) on a continuous variable (carbon) though a linear mixed model analysis. Study sites are included as the random effect in the model (with the random slope and random intercept). Landuse, species (and their interaction) are included as fixed effects.

the model is this -

model1 = lmer(carbon ~ species*landuse + (1+landuse|site), data)

I know that there may be interaction between landuse and species. I know that presence of interaction can change the interpretation of the main effects. I want to know what should I do if there is a significant interaction between landuse and species? In that case, do I study the effect of landuse for each species seprately with the following model -

model.sp1 = lmer(carbon ~ landuse + (1+landuse|site), data.sp1)

and repeat this for all the six species? Do I need any form of corrections (of p-values) due to running multiple tests?

Another question is that if the interaction term is not significant, can I interpret the main effects from model1 or do I run another model (model2) without the interaction effect and interpret the main effects from there?

model2 = lmer(carbon ~ species + landuse + (1+landuse|site), data)

I am fairly new to mixed model and R, so please excuse my naivety!

PS - just to clarify, I have a fairly good idea of what an interaction mean and how to interpret it. I do not know, how to interpret main effects in the presence of an interaction - whether I need to run a seprate analysis to interpret main effects.

Answer 1

Here's what I would do:

First, I would have a look here on how to specify the random term in your model1. I am not quite sure what you are trying to fit. There is also a lot of info on linear mixed effects models here on CV. Click on the lme4-nlme tag, which you also provided. It would also help if you could provide an example dataset, or at least the structure of your data.

Then, you most likely only need one model, which is presumably in the form of:

my_model <- lmer(carbon ~ species + landuse + species : landuse + (1|site), data = mydata)

I specified the random effect to be + (1|site), because you said:

Study sites are included as the random effect in the model.

To get the ANOVA table you can either do:

library(car)
Anova(my_model)

or:

library(afex)
mixed(carbon ~ species + landuse + species : landuse + (1|site), data = mydata)

or instead of running lmer() through the lme4 package, load the lmerTest package and run:

my_model <- lmer(carbon ~ species + landuse + species : landuse + (1|site), data = mydata)
anova(my_model) 

This will give you the ANOVA table you probably need eventually. Make sure to have a look at those functions and their arguments (?Anova, ?mixed, ?lmerTest::anova).

I don't quite understand why would want to exclude species if the interaction is significant and run separate models for all species?!

However, if your main effects are not significant you could consider tossing them out and re-running the model with the interaction only. However, if one or both main effects are significant, I would keep them both in the model and report this together with a potential significant interaction.

In any case, if you have a significant interaction you should focus on interpreting the interaction and not the main effects since their interpretation could now be misleading. The interpretation of the interaction should start by visualizing it. You could do this for example using the emmip() function in the emmeans package:

library(emmeans)
emmip(my_model, landuse ~ species)

Regarding the adjustment of p-values, you only need to do that if you are following up with post-hoc tests.

This could be done with the emmeans() function (also from the emmeans package):

emmeans(my_model, pairwise ~ species : landuse)