Pairwise PERMANOVA models in R using vegan::adonis

pairwise_permanova.R

#' Run pairwise comparisons by a grouping variable using `vegan::adonis`
#'
#' @param sp_matrix Community data matrix with columns for each species and rows for each site / observation
#' @param group_var Vector of categorical variable of length equal to the number of rows in sp_matrix.
#' @param dist dissimilarity index. See `vegan::vegdist` for options. Defaults to "bray"
#' @param adj P-value adjustment. Defaults to "fdr". See `stats::p.adjust` for options.
#' @param perm Number of permutations per model. Defaults to 10,000.
#'
#' @return A list with data frame containing contrasts and other elements for p-value adjustment and permutations.
#' @export
#'
#' @examples
#' 
#' library("vegan")
#' 
#' # Load species matrix and covariates
#' data("dune", package = "vegan")
#' data("dune.env", package = "vegan")
#' 
#' # pairwise comparisons among land use categories
#' pairwise_permanova(dune, dune.env$Use)
#' 
pairwise_permanova <- function(sp_matrix, group_var, dist = "bray", adj = "fdr", perm = 10000) {
  
  require(vegan)
  
  ## list contrasts
  group_var <- as.character(group_var)
  groups <- as.data.frame(t(combn(unique(group_var), m = 2)))
  
  contrasts <- data.frame(
    group1 = groups$V1, group2 = groups$V2,
    R2 = NA, F_value = NA, df1 = NA, df2 = NA, p_value = NA
  )
  
  for (i in seq(nrow(contrasts))) {
    sp_subset <- group_var == contrasts$group1[i] | group_var == contrasts$group2[i] 
    contrast_matrix <- sp_matrix[sp_subset,]
    
    ## fit contrast using adonis
    fit <- vegan::adonis2(
      contrast_matrix ~ group_var[sp_subset],
      method = dist, 
      perm = perm
    )
    
    contrasts$R2[i] <- round(fit$R2[1], digits = 3)
    contrasts$F_value[i] <- round(fit[["F"]][1], digits = 3)
    contrasts$df1[i] <- fit$Df[1]
    contrasts$df2[i] <- fit$Df[2]
    contrasts$p_value[i] <- fit$`Pr(>F)`[1]
  }
  
  ## adjust p-values for multiple comparisons
  contrasts$p_value <- round(p.adjust(contrasts$p_value, method = adj), digits = 3)
  
  return(list(
    contrasts = contrasts, 
    "p-value adjustment" = adj, 
    permutations = perm
  ))
}

Hi, I am new here and perhaps I might be talking non-sense, but I am in search of a pairwise permanova function to do a post hoc on the output of Adonis use with the adonis funtion from the vegan package. I've tried to use pairwise_permanova, and pairwise.adonis but they are not in the newest versions of R. Do you know any other function that could do that? Thank you

Hi! Sorry that this isn't very well documented, but I wrote this function to do exactly that - instead of using the output from adonis, the function takes as an input a matrix with species as columns, a row for each site, and abundances in each cell - the same input that adonis takes. The other main argument is a grouping variable, which would be a vector of the same length as the number of rows in the species matrix - if there are multiple grouping variables, you'd just paste them together. It's a pretty simple function, all it does is loop through the possible combinations of the grouping variable, subset the data, and run adonis on the subset, then at the end it adjusts p-values for multiple comparisons.

I use this function in this script from a repository for a recent paper. A minimal example:

## Load vegan library and pairwise_permanova function from R script
library("vegan")
source("pairwise_permanova.R") 

## Load species matrix and covariates
data("dune", package = "vegan")
data("dune.env", package = "vegan")

## Run pairwise permanova
groups <- as.character(dune.env$Use) # Apparently, my function doesn't like factors
pairwise_permanova(dune, groups)

I hadn't heard of the pairwise.adonis function you mentioned but to be honest it looks like nicer code than mine, although they both basically do the same thing.

Hope this helps!

Thank you very much for your reply.

Unfortunately, R is not finding your function even after installing vegan and RVAideMemoire

See the error I get below:

source("pairwise_permanova.R") 

Error in file(filename, "r", encoding = encoding) :
cannot open the connection
In addition: Warning message:
In file(filename, "r", encoding = encoding) :
cannot open file 'pairwise_permanova.R': No such file or directory

Because the function is one I wrote, it isn't included in base R or the vegan package - to read it in, press "download zip" above the gist and copy the "pairwise_permanova.R" file into your working directory, or copy the code above into a file and save it as "pairwise_permanova.R". The source function then simply runs the code in the script, which defines the function. vegan has to be installed for the function to work, but RVAideMemoire isn't needed. Cheers!

Thank you for not giving up on me :)
I did exactly that but keep receiving errors.

'adonis' will be deprecated: use 'adonis2' instead
Error in vegdist(lhs, method = method, ...) : input data must be numeric

My data have three columns: first is a character with 5 groups with 6 levels, then I have two columns with d13C and d15N values which are numeric.

I suspect that some functions are deprecated as for the Adonis, I did substitute for adonis2, but the method error message still appears... is "adj" not valid?

No worries, happy to help! Typically the data you feed into these models would be a numeric matrix as a response, with variables (e.g. species) as columns, observations (e.g. different sites) as rows, and some measure of abundance or prevalence in the cells. Covariates do not go in the response matrix, and the method isn't equipped to deal with character data, so I imagine the issue is that you've got a character covariate / predictor in the community data matrix. So, you might do something like this:

sp_matrix <- as.matrix(<data>[,c("d13C", "d15N")])
pairwise_permanova(sp_matrix, <data>$<group_name>)

Where <data> is whatever your data frame is named, and <group_name> is for the categorical predictor that you want to do pairwise comparisons between the levels of.

Of course, I don't actually know what your data look like so I can't say if this is actually what you're looking for. Also, if you only have two variables (d13C and d15N) you might just want to run two regressions, one with d13C as a response and the other with d15N as a response - typically this sort of analysis is used for high-dimensional data (e.g. counts of many different species).

Thank you very much, indeed my data is organized differently. I have many other analysis to do right now, but once I finish I will do like you said and get back to you, have a great week!

Hi Maurice, just to let you know it worked :) the problem was really that, that I was using also non-numeric data, thank you!

Glad to hear it, and you're welcome!

Hi again, I am now trying Permanova pairwise with the data I wanted to use since the beginning. When I told you it worked I was using other data and it did run smoothly, but now I have some error that I am trying to figure out, see below

Error in round(fit$R2[1], digits = 3) :
non-numeric argument to mathematical function
In addition: Warning message:
In vegdist(as.matrix(lhs), method = method, ...) :
results may be meaningless because data have negative entries
in method “bray”

I do have negative number (d13C) but I did use Permanova in the past with Primer with the same negative numbers and it worked so I am not sure why I can't do it now. Any clue? Thank you in advance!

That is the head of my data

I'm guessing the issue is because the R^2 for one or more of the pairwise models is NA or NULL, probably because the dissimilarity index used is inappropriate for the data. Typically the Bray-Curtis dissimilarity is used for species abundance values, which are strictly non-negative - hence the warning message:

results may be meaningless because data have negative entries in method “bray”

you might try dist = "euclidean" in your call to pairwise_permanova (and method = "euclidean" for the full adonis / adonis2 model) - but I'm not sure what type of distance / dissimilarity index is typical for your kind of data so there might be a better option.

Hi Maurice, thank you once more, the Euclidean worked, indeed is the Euclidean that is used, but I am still receiving the R² error
Error in round(fit$R2[1], digits = 3) :
non-numeric argument to mathematical function

So just to make sure it works on my end, I tested this with dist = "euclidean" on some simulated data, and it worked just fine - maybe the problem is something in your data? Do you have any NA, inf, or character values in either of those columns? Are you passing a separate vector for the group_var argument and are those all non-NA?

Ok, so I was able to get the same error as you by replacing adonis with adonis2 inside the pairwise_permanova function - is that what you did? If so the error is happening just because the output of adonis and adonis2 are formatted differently. I've made a few changes to the gist above so that the function uses adonis2 - give that a try :)

my data does not have missing values or characters in the columns. Yes that is correct. I had replaced adonis by adonis2 because it wasn't allowing me to run with adonis. Now it is working. smoothly, many thanks Maurice!!

Great to hear, and you're welcome!