Two-Way PERMANOVA (adonis, vegan-Package) with Customized Contrasts

...say you have a multivariate dataset and a two-way factorial design - you do a PERMANOVA and the aov-table (adonis is using ANOVA or "sum"-contrasts) tells you there is an interaction - how to proceed when you want to go deeper into the analysis?
You could, however somewhat tedious, customize contrasts for the PERMANOVA and check for differences between certain level combinations.







Here's an example with a 2-way factorial design with 4 dependent variables, which represent species abundances. I tested for significant differences between "treatment"-level 1 vs 2 and "treatment"-level 2 vs 3 as well as for simple "impact"-effects at each level of "treatment".
When comparing the default aov-table with the table for the customized contrasts one can see that the residual and the total sum of squares are the same for both models and that the main effects are decomposed in its elements as specified by the contrasts.

library(vegan)
 
# species matrix
spdf <- matrix(NA, 60, 4, dimnames =
               list(1:60, c("sp1", "sp2", "sp3", "sp4")))
spdf <- as.data.frame(spdf)
 
# 1st factor = treatment:
treat <- gl(3, 20, labels = paste("t", 1:3, sep=""))
 
# 2nd factor = impact:
imp <- rep(gl(2, 10, labels = c("yes", "no")), 3)
 
# simulating effect -
# simulation will add similar effects
# and no interactions (see plot):
eff <- sort(rep(1:6, 10))
 
# add noise:
spdf$sp1 = eff + rnorm(60, 0, 0.25)
spdf$sp2 = eff + rnorm(60, 0, 0.25)
spdf$sp3 = eff + rnorm(60, 0, 0.25)
spdf$sp4 = eff + rnorm(60, 0, 0.25)
 
# interaction plot for all species:
par(mfrow=c(2, 2), mar = c(2.5, 2.5, 0.5, 0.5))
for (i in 1:4) {interaction.plot(treat, imp, spdf[, i], lty = c(1, 2), legend = F);
                legend("topleft", bty = "n", cex = 1.2,
                       paste("Species", i, sep = " "));
                legend("bottomright", c("no", "yes"),
                       bty = "n", lty = c(2, 1))}
 
## create a design matrix of the contrasts for "imp"
contrasts(imp) <- c(-1, 1)
Imp <- model.matrix(~ imp)[, -1]
 
## create a design matrix of the contrasts for "treat"
contrasts(treat) <- cbind(c(0,1,0),c(0,0,1))
Treat <- model.matrix(~ treat)[, -1]
 
imp.in.t1 <- Imp * ifelse(treat == "t1", 1, 0)
imp.in.t2 <- Imp * Treat[, 1]
imp.in.t3 <- Imp * Treat[, 2]
 
## specify the orthogonal contrasts for "treat"
contrasts(treat) <- cbind(c(1, -1, 0), c(1, 0, -1))
 
## specify the design matrix of the orthogonal
## contrasts for "treat"
Treat.ortho <- model.matrix(~ treat)[, -1]
 
## create a factor for each of the orthogonal "treat" contrasts
treat1vs2 <- Treat.ortho[, 1]
treat1vs3 <- Treat.ortho[, 2]
 
## do the pm-manova with the full model
fm1 <- adonis(spdf ~ treat * imp, method = "euclidean", perm = 999)
 
## do the pm-manova with the orthogonal contrasts for imp and treat'
## and the interaction contrasts of interest
fm2 <- adonis(spdf ~ treat1vs2 + treat1vs3 +
                     imp.in.t1 + imp.in.t2 + imp.in.t3,
              method = "euclidean", perm = 999)
fm1; fm2
 
## check model with changed effects
eff <- sort(rep(1:6, 10))
eff[treat == "t1" | treat == "t2"] <- 3
 
# add noise:
spdf$sp1 = eff + rnorm(60, 0, 0.25)
spdf$sp2 = eff + rnorm(60, 0, 0.25)
spdf$sp3 = eff + rnorm(60, 0, 0.25)
spdf$sp4 = eff + rnorm(60, 0, 0.25)
 
# interaction plot for all species:
par(mfrow=c(2, 2), mar = c(2.5, 2.5, 0.5, 0.5))
for (i in 1:4) {interaction.plot(treat, imp, spdf[, i], lty = c(1, 2), legend = F);
                legend("topleft", bty = "n", cex = 1.2,
                       paste("Species", i, sep = " "));
                legend("bottomright", c("no", "yes"),
                       bty = "n", lty = c(2, 1))}
 
fm1 <- adonis(spdf ~ treat * imp, method = "euclidean", perm = 999)
fm2 <- adonis(spdf ~ treat1vs2 + treat1vs3 +
                     imp.in.t1 + imp.in.t2 + imp.in.t3,
              method = "euclidean", perm = 999)
fm1; fm2

To cite package ‘vegan’ in publications use:
Jari Oksanen, F. Guillaume Blanchet, Roeland Kindt, Pierre Legendre,
R. B. O'Hara, Gavin L. Simpson, Peter Solymos, M. Henry H. Stevens
and Helene Wagner (2011). vegan: Community Ecology Package. R package
version 1.17-9. http://CRAN.R-project.org/package=vegan