How to perform post-hoc comparison on interaction term with mixed-effects model?

Question

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",            
                        "Site6","Site7","Site8","Site9"),each=12))
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, 
             random=~1|Site/Hydrology/Depth)
anova(mod1)

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):

summary(glht(mod1,linfct=mcp(Sediment="Tukey")))

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, 
            random=~1|Site/Hydrology)
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.

Answer 1

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.

http://cran.r-project.org/web/packages/lsmeans/vignettes/using-lsmeans.pdf

Answer 2

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)
mydata$SHD <-  interaction(mydata$Sediment, mydata$Hydrology, 
                           mydata$Depth)
mod2 <- lme(Variable ~ -1 + SHD, data=mydata, 
                 random=~1|Site/Hydrology/Depth)
summary(glht(mod2, linfct=mcp(SHD="Tukey")))