Online Trick

...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

This is a small example of how custom contrasts can easily be applied with the contrast-package. The package-manual has several useful explanations and the below example was actually grabbed from there.
This example can also be applied to a GLM but I choose to use a LM because the coefficients are more easily interpreted.

set.seed(3)

dat <- data.frame(var1 = as.factor(rep(c("a","b"), each = 8)),
          var2 = as.factor(rep(c("c","d"), 8)),
          response = round(rnorm(16, 10, 4), 0))

mod <- lm(response ~ var1 * var2, data = dat)
summary(mod)
interaction.plot(dat$var2, dat$var1, dat$response)

# the intercept in the first line of summary(mod) tests var1a_var2c == 0
# in the second line of summary(mod) you test var1a_var2c - var1b_var2c == 0,
# which is the effect of var1 within var2c
# in the third line you test var1a_var2c - var1a_var2d == 0,
# which is the effect of var2 within var1a
# in the fourth line you test the interaction, which is:
# (var1a_var2c - var1a_var2d) - (var1b_var2c - var1b_var2d) == 0

# relevel var2, then the intercept is var1a_var2d
# and the second line tests var1a_var2d - var1b_var2d == 0,
# which is the effect of var1 within var2d
dat1 <- dat
dat1$var2 <- relevel(dat$var2, ref = "d")
mod1 <- lm(response ~ var1 * var2, data = dat1)
summary(mod1)

# test this with calling contrast():
require(contrast)
var2_in_var1 <- contrast(mod,
                         list(var2 = levels(dat$var2), var1 = "b"),
                         list(var2 = levels(dat$var2), var1 = "a"))
print(var2_in_var1, X = TRUE)