Frédéric Santos' notebook
Emacs and R tricks for anthropologists and archaeologists.
07 mai 2020

How to perform a NPMANOVA with R?

Performing a non-parametric multivariate analysis of variance (NPMANOVA), also known as permutational multivariate analysis of variance (PERMANOVA), is relatively easy thanks to the function adonis() implemented in the R package vegan. In this blog post :

• we’ll cover only the one-way NPMANOVA, since more sophisticated extensions are quite rare in archaeological sciences;
• we’ll assume that the reader is reasonably familiar with the classical univariate analysis of variance (ANOVA).

1 General introduction

1.1 Goal

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.

1.2 (Explicit and implicit) assumptions

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

1.3 A (very) quick overview of the maths

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.

2 A simple example in R

Let’s consider a simple example with only two groups and three variables. We’ll use the DSP reference sample [5] to illustrate the method. The DSP reference sample is available in the R package anthrostat. More precisely, we will consider the differences between males and females for three pelvic measurements (DCOX, PUM, SPU) within the population sample of Coimbra.

This example is rather caricatural since the observed difference between the centroids is very clear, and the sample sizes are important (Figure 1), so that one could think that a NPMANOVA is not really mandatory here. But let’s go all the same!

library(anthrostat)
library(scatterplot3d)
data(data_dsp)
## Select a subsample (population sample of Coimbra):
samp <- na.omit(subset(data_dsp, Collection == "Coimbra"))
## Display a 3D scatterplot:
with(samp, scatterplot3d(x = DCOX, y = SPU, z = PUM,
color = as.numeric(Sex), pch = 16,
box = FALSE))

Performing a NPMANOVA with the R function adonis() requires to choose a method, i.e. a dissimilarity that will be used to calculate pairwise distances between individuals. (This is required to compute the $$F$$ ratio.) For most applications in biological anthropology (like here, where we are dealing with continuous variables), the usual euclidean distance should be picked. For more details, one can read the help page of adonis().

Let’s perform the NPMANOVA:

library(vegan)
## Create the Y matrix of variables under comparison:
Y <- samp[, c("DCOX", "SPU", "PUM")]
## Perform a one-way NPMANOVA:
adonis(Y ~ samp$Sex, method = "euclidean", permutations = 999) Call: adonis(formula = Y ~ samp$Sex, permutations = 999, method = "euclidean")

Permutation: free
Number of permutations: 999

Terms added sequentially (first to last)

Df SumsOfSqs MeanSqs F.Model      R2 Pr(>F)
samp$Sex 1 13527 13527.2 137.51 0.41858 0.001 *** Residuals 191 18789 98.4 0.58142 Total 192 32316 1.00000 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Unsurprisingly, the null hypothesis of no differences between males and females is strongly rejected. (Reminder: the p-value is an empirical p-value. With 999 permutations, a p-value of 0.001 is the “most significant” result we can obtain: that means that the observed difference could be replicated by none of the 999 permutation trials.) 3 (Slightly) more advanced topics 3.1 “Should I perform one NPMANOVA or several univariate tests”? 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? 3.2 How to choose between MANOVA and NPMANOVA? 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. Exercise Visit and read this good blog post. In the situation exposed by the author, the differences among the three groups are mostly due to a different correlation pattern between the two response variable. Perform a NPMANOVA on the dataset provided there, and compare the p-value to the p-value obtained with a parametric MANOVA. Do you see how different are those two methods? 3.3 Post-hoc analysis (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",

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