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

## No comments :

## Post a Comment