How to perform a NPMANOVA with R?

NPMANOVA is a non-parametric multivariate test used to compare the multivariate distribution of several groups [1,3]. Its null hypothesis is that the centroids of all groups are equal: thus, a rejection of the null (i.e., a low p-value) means that there is at least one pair of groups with significantly different centroids.

For example, let’s say that you have three groups with \(n_1\), \(n_2\) and \(n_3\) individuals respectively. The individuals are described by continuous variables (e.g., osteometric measurements such as bone lengths or diameters) \(X_1, \cdots, X_p\). The NPMANOVA tests the hypothesis that the three underlying populations have the same “profile” (i.e., the same centroid, or multivariate mean) for the set of variables under consideration.

If this null hypothesis were true, the observed differences among the three centroids should be relatively “small”. “Small” means that they should have the same magnitude that what would be obtained under random re-allocation of individuals to one of the three groups. This is exactly how the test proceeds: it is called permutational MANOVA because it simply checks if the differences actually observed are compatible with the differences we get by permuting the group labels of individuals.

The usual (parametric) MANOVA has relatively restrictive assumptions, that are not always met with the datasets gathered in archaeological sciences. Conversely, the NPMANOVA is a distribution-free method, that makes only one explicit assumption, namely the exchangeability of units under the null hypothesis (i.e., individuals are supposed to be independent, without repeated measures or autocorrelation of any kind).

However, several papers [1,4,3,7] warn that (at least in the case of unbalanced designs, i.e. a different number of individuals in each group), NPMANOVA may lead to spurious results when the groups have a substantially different multivariate dispersion. Indeed, the test statistic may be sensitive to differences in dispersion, so that “significant differences may be caused by different within-group variation (dispersion) instead of different mean values of the groups” (citation taken from the documentation of the R functions adonis()).

Thus, the real (although implicit) assumption is that the groups under comparison have a roughly equal multivariate dispersion, which is similar to the assumption of variance homogeneity in univariate ANOVA. This can be checked either visually with a principal component analysis, or formally with a specific test [2] implemented in the R function vegan::betadisper().

Full mathematical details can be found in the seminal publication by M. J. Anderson [1], and a quick version can be found on Wikipedia. We’ll not get into this in this blog post. But, just as a general idea:

• The reader may be aware of the general result (used in univariate ANOVA) saying that (in a sample divided in groups of individuals) the total variance is equal to the sum of within-group variance and between-groups variance. NPMANOVA also operates a partition of multivariate dispersion: it is thus another technique of partition of sum of squares.
• The test statistic is a “pseudo \(F\) ratio”, which is essentially a ratio of between-groups dispersion over within-group dispersion.
• Let \(F_\text{obs}\) be the actual \(F\) statistic observed in your data. By performing a great number \(B\) of permutations of the data, we can compute \(B\) values of the \(F\) ratio under the null hypothesis. The (empirical) p-value is thus defined by: \(p = \frac{\#\left\{F \geq F_\text{obs}\right\}+1}{B+1}\), i.e. the proportion of cases where we can get, under the null hypothesis, an \(F\) ratio more extreme than \(F_\text{obs}\).

Remark: the adjustment “\(+1\)” is the reason why we often choose a number of permutation \(B = 999\) or \(B = 9999\), so that \(B+1\) can be a “pleasant” number.

This is not always a trivial decision, and this can even be seen as a controversial choice [6].

As a reminder, ANOVA tests for differences between three or more groups for one single continuous response variable. Thus, if you have four variables and want to compare three groups, there are essentially two approaches:

1. You can perform one single MANOVA, comparing the three groups for the considered set of four input variables.
2. You can perform four different ANOVAs or Kruskal-Wallis tests (one per variable). Warning: in this case, using a method of adjustment of p-values might be necessary if you have a great number of response variables.

In fact, both approaches are not mutually exclusive, and you can perfectly do both of them (although this is part of the controversy detailed in [6]). But it is important to note that those two approaches address different research questions: are you interested in differences for each single, isolated outcome variable? or in differences in the whole pattern of multiple dependent variables?

The seminal paper by M. J. Anderson [1] clearly explains the differences between the classical (parametric) MANOVA, and its non-parametric alternative:

• NPMANOVA must be used when the (restrictive) assumptions of parametric MANOVA are not satisfied.
• But NPMANOVA “takes no account of the correlations among variables” (see Figure 4 in [1] for an illustration), so that NPMANOVA is not the appropriate testing framework if you want to know whether the correlation pattern between outcome variables is identical among groups. The classical parametric MANOVA can handle that, but this is not the case of NPMANOVA.

(Reminder: a NPMANOVA only tests whether all centroids are equal. If your p-value is low, you only know that this hypothesis is very unlikely. The negation of “all centroids are equal” is simply “there exists at least one pair of centroids which are not equal”, but you do not know which ones! To investigate this question, one can use post-hoc tests.)

To the best of my knowledge, there are not many R functions to perform post-hoc comparisons between groups after a NPMANOVA. The R package pairwiseAdonis (which can be installed from here) seems to be the only solution available (although creating your own procedure of post-hoc tests is not really complicated if you are comfortable with R programming). Its use is documented on the GitHub repo and is quite straightforward.

Let’s illustrate that on another simple example, making use of the classical iris dataset included in R. We will consider only two response variables: Sepal.Length and Sepal.Width; and we will compare those two outcome variables among three species of iris (50 flowers in each group). First, as usual, one should start with a graphical inspection:

## Load iris dataset:
data(iris)
## Scatterplot for two variables:
library(lattice)
xyplot(Sepal.Length ~ Sepal.Width, groups = Species, data = iris,
       pch = 16, auto.key = TRUE)

A difference in sepal measurements seems to be clear here (and there is no huge difference in within-group dispersions, so that the assumptions of NPMANOVA are met). Let’s verify that with a NPMANOVA:

## Only consider two variables:
Y <- iris[, c("Sepal.Length", "Sepal.Width")]
## Perform the NPMANOVA:
adonis(Y ~ iris$Species, method = "euclidean")

Call:
adonis(formula = Y ~ iris$Species, method = "euclidean") 

Permutation: free
Number of permutations: 999

Terms added sequentially (first to last)

              Df SumsOfSqs MeanSqs F.Model      R2 Pr(>F)
iris$Species   2    74.557  37.279  97.999 0.57143  0.001
Residuals    147    55.918   0.380         0.42857       
Total        149   130.475                 1.00000       
                
iris$Species ***
Residuals       
Total           
---
Signif. codes:  
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Once again, the very low p-value is a strong evidence against the null hypothesis of no difference between groups. Now let’s proceed to the post-hoc analysis:

library(pairwiseAdonis)
pairwise.adonis(Y, iris$Species, sim.method = "euclidean",
                p.adjust.m = "bonferroni")

                    pairs Df SumsOfSqs   F.Model        R2
1    setosa vs versicolor  1   32.4466 102.54269 0.5113260
2     setosa vs virginica  1   67.7210 174.47442 0.6403332
3 versicolor vs virginica  1   11.6680  26.72319 0.2142600
  p.value p.adjusted sig
1   0.001      0.003   *
2   0.001      0.003   *
3   0.001      0.003   *

Even when using the most conservative method ("bonferroni") for the adjustment of p-values, we get for each pair of groups a strong evidence against the hypothesis of no difference.