I have collected data on gas fluxes from plots of soil subjected to 5 different treatments ("D2", "K2", "M", "N", and "O2"), which also possessed variable clay contents. The experiment was laid out in a randomized complete block design, with 4 replications. Within each plot, two separate measurements of flux were performed. 4 blocks x 5 treatments x 2 repeated measurements for 40 total observations

I have built a linear mixed model with nlme, and performed multiple comparisons with Tukey's test, as shown below. The data.frame cum_fluxes can be downloaded from: https://www.dropbox.com/s/e58k4vdnevw2fsm/cum_fluxes

cum_fluxes <- dget("cum_fluxes")
n2o_flux.lme <- lme(cum_flux_n2o_ln ~ treatment * clay + block, random = ~1|subsampling, data = cum_fluxes)
aov.n2o <- anova(n2o_flux.lme)
comp.treat <- glht(n2o_flux.lme, linfct=mcp(treatment="Tukey"))

I have two questions. On the one hand, the resulting ANOVA table looks like this (significant effect of treatment with p = 0.0279):

               numDF denDF   F-value p-value
(Intercept)        1    20 14815.691  <.0001
treatment          4     7     5.281  0.0279
clay               1     7     1.326  0.2874
block              3     7     0.674  0.5951
treatment:clay     4     7     0.810  0.5568

However, there are apparently no significant differences found when considering pair-wise comparisons of the five treatments. Did I make a mistake somewhere, or am I interpreting the results incorrectly?

Second, I get the following warning message when I run the glht() command:

Warning message:
In mcp2matrix(model, linfct = linfct) :
  covariate interactions found -- default contrast might be inappropriate

This makes some sense to me: since I allow for interaction between a covariate (quantitative regressor – i.e., clay content) and a factor (so, non-parallel lines), the difference between two factor levels depends on the value of the covariate. However, I can't figure out how and to what value I should specify the covariate when I compare the factor levels?


1 Answer 1


It's not really possible to assess main effects in a model that includes an interaction term. If you conclude that you don't really need the interaction term, then you need to re-run the model without the interaction term. In order to assess whether you need the interaction term, you could use the AIC. But to compare AICs for mixed models with different fixed effects, you would have to use maximum likelihood (not REML). So you would need to include 'method="ML"' as an argument to lme.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.