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)
As said by @amoeba in the comment, the question is not so much a mixed model question, but more a general question on how to parameterize a regression model with interactions. The full quote from our chapter also provides an answer to your second question (i.e., the why):
A common contrast scheme, which is the default in R, is called
treatment contrasts (i.e., contr.treatment; also called dummy
coding). With treatment contrasts the first factor level serves as the
baseline whereas all other levels are mapped onto exactly one of the
contrast variables with a value of 1. As a consequence, the intercept
corresponds to the mean of the baseline group and not the grand mean.
When fitting models without interactions, this type of contrast has
the advantage that the estimates (i.e., the parameters corresponding
to the contrast variables) indicate whether there is a difference
between the corresponding factor level and the baseline. However, when
including interactions, treatment contrasts lead to results that are
often difficult to interpret. Whereas the highest-order interaction is
unaffected, the lower-order effects (such as main effects) are
estimated at the level of the baseline, ultimately yielding what are
known as simple effects rather than the usually expected lower-order
effects. Importantly, this applies to both the resulting parameter
estimates of the lower order effects as well as their Type III tests.
In other words, a mixed model (or any other regression type model)
that includes interactions with factors using treatment contrasts
produces parameter estimates as well as Type III tests that often do
not correspond to what one wants (e.g., main effects are not what is
commonly understood as a main effect). Therefore we generally
recommend to avoid treatment contrasts for models that include
interactions.
Orthogonal sum-to-zero contrasts are better because they avoid potentially difficult to interpret lower-order effects. That is, for those contrasts all lower order effects are evaluated at the grand mean. For a quick explanation of dummy vs. effect coding difference, please see: http://www.lrdc.pitt.edu/maplelab/slides/Simple_Main_Effects_Fraundorf.pdf
This means for your case, almost all your interpretations are correct with one exception.
- ConditionB - what is the difference in Intercept for Condition B from Condition A, when X is zero.
Hence, if zero is somewhat meaningless for your variable (e.g., it is age and you only observe adult participants) your estimate of Condition (which is now a simple effect of condition at X = 0) becomes meaningless as well.
In general, having interactions with continuous covariates is not trivial and there are at least two books and several papers which extensively discuss this issue. A common solution is centering the covariate on the mean. Whether or not this makes sense depends on your covariate. What I sometimes do when having a variable with a restricted range (e.g., it goes from 0 to 100) is to center on the midpoint of the scale (see e.g., here).
More information on centering can be found in the following references. I recommend you read the first one at least:
- Dalal, D. K., & Zickar, M. J. (2012). Some Common Myths About Centering Predictor Variables in Moderated Multiple Regression and Polynomial Regression. Organizational Research Methods, 15(3), 339-362. doi:10.1177/1094428111430540 [free pdf]
- Cohen, J., Cohen, P., West, S. G., & Aiken, L. S. (2002). Applied Multiple Regression/Correlation Analysis for the Behavioral Sciences. New York: Routledge Academic. [great book]
- Aiken, L. S., & West, S. G. (1991). Multiple regression: testing and interpreting interactions. Newbury Park, Calif.: Sage Publications.
There is also some mixed-model specific discussion on centering, but to me this appears to be mainly relevant for hierarchical structures (i.e., at least two-levels of nesting), e.g.,
- Wang, L., & Maxwell, S. E. (2015). On disaggregating between-person and within-person effects with longitudinal data using multilevel models. Psychological Methods, 20(1), 63–83. https://doi.org/10.1037/met0000030
Potentially also relevant:
- Iacobucci, D., Schneider, M. J., Popovich, D. L., & Bakamitsos, G. A. (2016). Mean centering helps alleviate “micro” but not “macro” multicollinearity. Behavior Research Methods, 48(4), 1308–1317. https://doi.org/10.3758/s13428-015-0624-x
Best Answer
When there is an interaction, there is no unique slope for any of the main effects; those change based on other variables in the interaction.
I have found that the simplest method of interpreting interactions is visually: make graphs of your DV at different levels of your IVs. Since you have two numeric IVs and 3 categorical ones, you could make 6 graphs. For each, plot one numeric IV on the x axis and the predicted value of the DV on the Y axis, with a line for each level of one of the categorical DVs (you will have to make some assumptions about the value of the other IVs - the mean for the continuous and the mode for the categorical may be sensible).
Alternatively, you could make a lattice plot (using the
latticepackage) or use faceting (with theggplot2package).Another way to go is to make a table of the predicted values of the DV for various common combinations of the IVs. (E.g. the quartiles of the two continuous IVs and all the values of the DVs - which would give a table with 3*3*3*2*2 = 108 rows.
As to what the interactions mean - an interaction means that the relationship between one IV and the DV is different at different levels of the other IV.