# Mixed-model ANCOVA using nlme in R: How do I interpret and report significant interaction effects?

I've got a significant interaction effect and I'm unsure how to interpret and report the results.

In my experiment, I've got plants planted inside fences to protect them from deer (fence plots) and outside fences where they can be eaten (open plots). Within each plot, half the area of gm (garlic mustard) was removed. Before any gm was removed, the initial % cover of gm (ipc.gm) was measured in each subplot.

Variables

height4 = final plant height

height = initial plant height

ipc.gm = initial % cover of gm

fence = factor with 2 levels: Fence, Open

gm = factor with 2 levels: removed, gm

plot/subplot = random effect variable

My Model

model <- lme(height4 ~ height + ipc.gm * fence * gm, random = ~1|plot/subplot)


Analysis

> anova.lme(model, type = "marginal")
numDF denDF   F-value p-value
(Intercept)         1    92  0.821069  0.3672
height              1    92 30.827557  <.0001
ipc.gm              1    92  0.673770  0.4139
fence               1     8 15.483526  0.0043
gm                  1     8  1.066909  0.3319
ipc.gm:fence        1    92  5.021460  0.0274
ipc.gm:gm           1    92  0.023078  0.8796
fence:gm            1     8  1.967297  0.1983
ipc.gm:fence:gm     1    92  2.088226  0.1518


The main effect of fence is significant (P = 0.0043).

> emmeans(model, ~ fence)
NOTE: Results may be misleading due to involvement in interactions
fence   emmean       SE df  lower.CL upper.CL
Fence 24.00245 1.052951  3 20.651490 27.35341
Open  11.67182 1.347438  3  7.383666 15.95996

Results are averaged over the levels of: gm
Degrees-of-freedom method: containment
Confidence level used: 0.95


Plants were taller in Fenced plots (24 cm) than in Open plots (11.7 cm).

The main effect of ipc.gm is not significant (P = 0.4139), but the interaction term ipc.gm:fence is significant (P = 0.0274).

This means that the effect of ipc.gm depended on the level of fence

Plotting height4 ~ ipc.gm, subset by fence level (package: lattice):

xyplot(height4 ~ ipc.gm, group = fence, type = c("p","r"), auto.key = TRUE)


It appears that height4 increases with increasing ipc.gm, but only in fence plots.

My questions are:

1. How do back up my statement "height4 increases with ipc.gm but only in fence plots" with stats? Do I just need to report the slopes of the lines? Or do I need to do a statistical test that show the slope of the line is significantly different from zero?
2. How do I get the slope coefficients of height4 ~ ipc.gm, subset by fence level?
3. What statistical test (if any) do I need to use to show that height4 increases with ipc.gm in Fence plots, and that height4 is unaffected by ipc.gm in Open plots
• I've rolled back your question because you introduced a number of new questions that (I think) muddied the present question after it was answered (and you accepted an answer). Instead, I recommend you post a new question with your additional queries. – mkt - Reinstate Monica Jan 24 at 10:57

You can directly estimate and compare the slopes of those two lines:

emtrends(model, pairwise ~ fence, var = “ipc.gm”)


1) Slope estimates along with confidence intervals will take care of this. It's actually already reasonably clear from the model output (summary & plot) you have shown. You've already done the statistical test above (p = 0.0274 for ipc.gm:fence), so you really don't need any additional tests.

2) This is more an implementation question that is off-topic here. I expect you will get most of the relevant information from coef(lme), though you will need to do some work with the interaction coefficients to extract exactly what you need. The effects package may be helpful for this.

3) This seems to be the same as #1

• Thank you for your answer! I got the slope estimates using emtrends(), but doesn't the p = 0.0274 for ipc.gm:fence just mean that the slopes are not equal? I'd like to be able to say that "Height4 increased with ipc.gm in fenced plots" and of course I will report the positive slope and show the figure, but would I also need to show that the slope is significantly different from zero? – Jay Jan 18 at 0:58
• Doesn’t the second party of that emtrends output give you the difference of the slopes, its SE, and its P value of .0274? – rvl Jan 18 at 2:11
• @Jay The slope estimates along with their CIs will tell you this. If the CIs for open plots include zero, and the CIs for fenced plots do not, this would indicate exactly what you are saying. I am not familiar with emtrends(), so I don't know exactly what information that gives you. – mkt - Reinstate Monica Jan 18 at 7:28
• @mkt Thank you so much, I understand now. I've also updated my post to include the emtrends output and I tried to explain it a bit (for future readers). I am now just wondering if my slope is actually for (height4-height) ~ ipc.gm instead of height4 ~ ipc.gm, any thoughts? – Jay Jan 18 at 17:40
• @Jay It's a bit hard to understand your new question and I would not edit answered questions to add new ones; I suggest instead posting a new question with more detail about height4 and height. – mkt - Reinstate Monica Jan 21 at 9:47

The notion of 'significance' means very little, and the notion of 'insignificance' even less. Pre-specify the contrasts you want to estimate, and estimate them with confidence (compatibility) intervals. Contrasts can include interaction-type contrasts.