I’m working on a data set in order to evaluate the impact of drying on sediment microbial activities. The objective is to determine if the impact of drying varies among sediment types and/or depth within the sediment.

The experimental design is as follows:

  • The first factor Sediment corresponds to three types of sediment (coded Sed1, Sed2, Sed3). For each type of Sediment, sampling was carried out on three sites (3 sites for Sed1, 3 sites for Sed2, 3 sites for Sed3).
  • Site is coded : Site1, Site2, ..., Site9.
  • The next factor is Hydrology: within each site, sampling is carried out in a dry plot and in a wet plot (coded Dry/Wet).

Within each of the previous plot, sampling is carried out at two Depths (D1, D2) in triplicate.

There are a total of n = 108 samples = 3 Sediment * 3 Sites * 2 Hydrology * 2 Depths * 3 replicates.

I use the lme() function in R (nlme package) as follows :

Sediment <- as.factor(rep(c("Sed1","Sed2","Sed3"),each=36))
Site <- as.factor(rep(c("Site1","Site2","Site3","Site4","Site5",            
Hydrology <- as.factor(rep(rep(c("Dry","Wet"),each=6),9))
Depth <- as.factor(rep(rep(c("D1","D2"),each=3),18))
Variable <- rnorm(108)

mydata <- data.frame(Sediment,Site,Hydrology,Depth,Variable)

mod1 <- lme(Variable ~ Sediment*Hydrology*Depth, data=mydata, 

I would like to run a post-hoc comparison to test whether a term is significant or not.

I'm able to do it for a simple main effect (e.g., Sediment):


But the glht() function doesn't work for interaction terms.

I found that the following could work for a 2-way anova :

mod1 <- lme(Variable~Sediment*Hydrology, data=mydata, 
mydata$SH <- interaction(mydata$Sediment, mydata$Hydrology)
mod2 <- lme(Variable ~ -1 + SH, data=mydata, random=~1|Site/Hydrology)
summary(glht(mod2, linfct=mcp(SH="Tukey")))

Is it possible to use the same approach in the case of a 3-way anova? Any help on the way to do post-hoc comparison on interaction terms in this case would be much appreciated.


I found the package "lsmeans" quite useful especially when there is a x*z*v interaction. However, the package is available only for newer versions of R.


| cite | improve this answer | |

Do you mean you want to do all pairwise comparisons for the three factors?

mod1<-lme(Variable~Sediment*Hydrology*Depth, data=mydata, random=~1|Site/Hydrology/Depth)
mod2<-lme(Variable~-1+SHD, data=mydata, random=~1|Site/Hydrology/Depth)
| cite | improve this answer | |
  • $\begingroup$ I want to do pairwise comparisons for the two-way interaction terms, for example "SedimentHydrology". If the anova table indicates that the SedimentHydrology interaction is significant for the full model mod1<-lme(Variable~Sediment*Hydrology*Depth, data=mydata, random=~1|Site/Hydrology/Depth) I don't think that it is correct to do post-hoc comparisons (on the Sediment*Hydrology interaction) the following way : mydata$SH<-interaction(mydata$Sediment,mydata$Hydrology) mod2<-lme(Variable~-1+SH, data=mydata, random=~1|Site/Hydrology/Depth) summary(glht(mod2,linfct=mcp(SH="Tukey"))) $\endgroup$ – John Smith Feb 22 '12 at 8:30

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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