adonis in vegan: order of variables or use of strata

Question

I am using the adonis() function in in the vegan package to determine 1) if co-occurring host species vary in their microbial community across multiple sites, and 2) if sites are different. I have examined all posts on CV and SO, and there is no clear answer to how to determine significance of multiple factors using the adonis function.

I first did this, as suggested by https://stackoverflow.com/questions/26768779/vegan-adonis-unbalanced-design-ss-type-ii-or-iii:

where jacc is a dissimilarity matrix using the jaccard metric

adonis <- adonis(jacc ~ Species + Site, data = df_compare)

adonis
Call:
adonis(formula = jacc ~ Species + Site, data = df_compare) 

Permutation: free
Number of permutations: 999

Terms added sequentially (first to last)

          Df SumsOfSqs MeanSqs F.Model      R2 Pr(>F)    
Species    2    0.6055 0.30273  1.7690 0.04981  0.004 ** 
Site       4    2.1378 0.53445  3.1231 0.17587  0.001 ***
Residuals 55    9.4122 0.17113         0.77432           
Total     61   12.1554                 1.00000           
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Then reverse the order:

adonis_2 <- adonis(jacc ~ Site + Species, data = df_compare)

adonis_2

Call:
adonis(formula = jacc ~ Site + Species, data = df_compare) 

Permutation: free
Number of permutations: 999

Terms added sequentially (first to last)

          Df SumsOfSqs MeanSqs F.Model      R2 Pr(>F)    
Site       4    2.4385 0.60962  3.5623 0.20061  0.001 ***
Species    2    0.3048 0.15238  0.8904 0.02507  0.716    
Residuals 55    9.4122 0.17113         0.77432           
Total     61   12.1554                 1.00000           
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

But I don't know how to interpret that, because the order matters, and I'm not really sure if there are differences among species.

After some searching I decided to use strata.

I think this is saying: are co-occurring species different when you only compare species at the same sites.

species_adonis <- adonis(jacc ~ Species, strata = df_compare$Site, data = df_compare)

species_adonis

Call:
adonis(formula = jacc ~ Species, data = df_compare, strata = df_compare$Site) 

Blocks:  strata 
Permutation: free
Number of permutations: 999

Terms added sequentially (first to last)

          Df SumsOfSqs MeanSqs F.Model      R2 Pr(>F)
Species    2    0.6055 0.30273  1.5464 0.04981  0.335
Residuals 59   11.5500 0.19576         0.95019       
Total     61   12.1554                 1.00000  

Then to ask the question about site, I used species in the blocking.

I think this is saying: are sites different when you only comparing the same species

site_adonis <- adonis(jacc ~ Site, strata = df_compare$Species, data = df_compare)

Call:
adonis(formula = jacc ~ Site, data = df_compare, strata = df_compare$Species) 

Blocks:  strata 
Permutation: free
Number of permutations: 999

Terms added sequentially (first to last)

          Df SumsOfSqs MeanSqs F.Model      R2 Pr(>F)    
Site       4    2.4385 0.60962  3.5761 0.20061  0.001 ***
Residuals 57    9.7169 0.17047         0.79939           
Total     61   12.1554                 1.00000           
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

My conclusion is that the microbial community on a given species differs among sites, but that the microbial community does not differ among host species.

Is my approach correct or am I misinterpreting the use of strata (i.e. blocking)?

Or is there a way to somehow average over the tests when I switched the order of the variables?

Answer 1

As you've noted yourself, by running two adonis models with your fixed factors inverted you see that both the variance assigned to each factor, and the P-values differ each time. This occurs in unbalanced designs such as yours, where the degrees of freedom associated with each factor differ.
From the description of your experiment, it looks like a classical case of a nested design where Species is nested in Site. In this case the model you're looking for should look like this:
adonis <- adonis(jacc ~ Site / Species, strata = Site, data = df_compare).
Do note that nestedness should be stated in the model formulation as well as in the strata (see reply by Jari Oksanen).