R: Statistically testing NMDS community differences

Situation

I am comparing species communities at 6 different sites. Communities at each site were measured repeatedly over 5 years. I used NMDS plots, based on bray-curtis distance matrices using standardised abundances, to visualise community differences.

The data:

sp.site <- structure(list(Year = structure(c(1L, 1L, 1L, 1L, 1L, 1L, 2L,
2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L, 3L, 4L, 4L, 4L, 4L, 4L,
4L), .Label = c("year1", "year2", "year3", "year4"), class = "factor"),
Site = structure(c(1L, 2L, 3L, 4L, 5L, 6L, 1L, 2L, 3L, 4L,
5L, 6L, 1L, 2L, 3L, 4L, 5L, 6L, 1L, 2L, 3L, 4L, 5L, 6L), .Label = c("Site1",
"Site2", "Site3", "Site4", "Site5", "Site6"), class = "factor"),
Species1 = c(0.95, 1, 0.9, 1.33, 0, 0, 1.13, 0, 1.08, 2.7,
0, 0, 1.44, 0, 0, 2.34, 0, 0, 2.27, 0, 1.99, 2.37, 1.05,
1.5), Species2 = c(0, 0, 0, 0, 0, 0.66, 0, 0, 0, 0, 0, 0,
0, 0, 0, 1.03, 0, 1.18, 0, 0, 0, 0, 0, 2.06), Species3 = c(0,
0, 0, 0, 0, 0, 1.68, 0, 0, 0, 0, 1.28, 2.19, 0, 1.67, 0,
0, 0, 0, 0, 0, 0, 0, 0), Species4 = c(0.72, 0, 0, 0, 0.74,
1.11, 0, 0, 0, 0, 0, 1.06, 0, 1.29, 0, 0, 1.7, 0, 1.5, 0,
0, 0, 1.26, 1.1), Species5 = c(0.72, 1.47, 0.65, 0, 0.53,
0, 0, 1.32, 2.67, 0, 0, 0, 0, 2.45, 1.85, 0, 1.09, 2.22,
0, 1.04, 1.23, 0, 0, 0), Species6 = c(0, 0, 0, 0, 0, 0.62,
2.25, 1.25, 0, 0, 0, 0, 0, 0, 1.07, 0, 0, 1.11, 0, 0, 0,
0, 0, 0), Species7 = c(0, 0, 1.59, 0, 0, 0, 0, 0, 0, 0, 0,
1.06, 2.19, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0), Species8 = c(1.32,
2.65, 3.79, 1.39, 1.92, 1.8, 1.9, 2.39, 3.66, 1.45, 2.56,
1.87, 1.04, 2.45, 2.36, 1.03, 1.7, 0.84, 1.62, 1.37, 3.17,
0, 3.93, 2.95), Species9 = c(1.49, 2.71, 2.78, 0.62, 1.23,
0, 1.68, 2.56, 0, 1.08, 2.21, 1.28, 2.08, 0, 0, 0, 0, 0,
2.19, 1.6, 0, 0, 0, 1.64), Species10 = c(0, 0.56, 0.79, 0.86,
0, 0, 0, 1.28, 1.33, 1.45, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
1.11, 0, 1.04, 1.14), Species11 = c(0, 0.54, 0, 0, 0, 0,
0, 0, 0, 1.32, 0, 0, 1.1, 0, 0, 0, 2.38, 0, 1.1, 1.04, 1.05,
0, 1.65, 0), Species12 = c(1.35, 1.22, 0, 1.14, 1.5, 2.52,
1.64, 1.02, 1.21, 1.39, 2.21, 2.1, 1.04, 0, 1.02, 1, 1.9,
2.54, 0, 1.6, 0.99, 0, 2.03, 2), Species13 = c(1.19, 1.01,
0.72, 1.9, 0.53, 0.62, 1.67, 1.77, 0, 0, 2.22, 0, 2.08, 0,
0, 0, 4.35, 2.22, 0, 1.09, 1.16, 1.78, 1.65, 0), Species14 = c(0,
2.56, 0.78, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 3.69, 0, 0,
0, 0, 0, 0, 0, 0, 0), Species15 = c(0, 1, 0.78, 0, 0, 0.56,
2.25, 1.72, 0, 0, 2.22, 1.16, 2.08, 0, 0, 0, 1.69, 0, 1.02,
2.13, 1.11, 0, 2.09, 0), Species16 = c(0.72, 0.68, 0.56,
1.33, 0, 0, 1.11, 0, 1.08, 2.03, 0, 1.06, 0, 2.21, 1.23,
1.03, 0, 0, 3.85, 1.04, 0, 1.53, 1.02, 1.1), Species17 = c(0,
0, 0, 0, 0.53, 0, 0, 0, 0, 1.32, 0, 0, 0, 0, 1.02, 0, 0,
1.25, 0, 0, 0, 0, 0, 1.64), Species18 = c(1.19, 0.85, 0.6,
0, 1.48, 1.01, 1.67, 0.86, 0, 0, 1.81, 0, 1.04, 1.69, 0,
0, 0, 1.11, 0, 2.01, 1.71, 0, 0, 1.18), Species19 = c(0.78,
1.85, 0.7, 0, 0.62, 0, 1.4, 0, 1.08, 0, 0, 0, 0, 0, 1.28,
0, 1.69, 0, 1.1, 0, 0, 0, 1.03, 1.18), Species20 = c(0.95,
1.85, 0, 0, 0, 0.56, 0, 0, 0, 1.45, 0, 1.23, 1.1, 0, 0, 1,
0, 0, 0, 1.09, 0, 0, 0, 0), Species21 = c(0, 0, 0, 1.15,
0.93, 0, 0, 1.25, 0, 1.39, 0, 1.28, 0, 0, 0, 0, 0, 1.11,
2, 0, 0, 0, 0, 0), Species22 = c(0, 0.85, 0, 0, 1.01, 0.56,
2.67, 1.02, 1.28, 0, 0, 1.28, 0, 0, 0, 0, 1.31, 1.01, 2.19,
0, 0, 2.74, 1.03, 0), Species23 = c(1.23, 0.68, 1.08, 0.86,
0, 1.02, 2.22, 0.86, 0, 2.91, 6.62, 1.28, 1.04, 0, 1.05,
1.69, 0, 3.03, 1.02, 0, 0.99, 0, 0, 1.32), Species24 = c(0,
1.25, 1.39, 0, 0.85, 1.17, 0, 1.32, 2.56, 0, 0, 1.28, 0,
0, 2.04, 0, 2.62, 0, 0, 0, 0.99, 0, 0, 1.64), Species25 = c(0.7,
0, 0, 0.97, 0.74, 0, 0, 0, 1.08, 1.32, 0, 0, 0, 2.21, 0,
0, 0, 0, 0, 0, 0, 1.05, 0, 2.06), Species26 = c(0.95, 2.38,
1.38, 0.97, 0.63, 0.67, 1.49, 1.95, 1.08, 1.22, 0, 1.23,
1.28, 2.21, 2.22, 0, 3.33, 1.46, 0, 0, 1.11, 0, 4.96, 1.5
), Species27 = c(1.9, 0, 0, 0, 0.85, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 1.28, 0, 0, 0, 0, 0, 0), Species28 = c(1.19,
0, 0.87, 0, 0, 0, 0, 1.32, 1.33, 0, 0, 0, 2.08, 1.29, 1.82,
1, 0, 0, 1.22, 1.09, 0, 0, 0, 1.14), Species29 = c(0.95,
1.04, 0, 0, 0, 0, 0, 0, 0, 0, 2.21, 0, 1.04, 0, 0, 0, 0,
1.52, 0, 6.03, 0, 0, 0, 0), Species30 = c(0.85, 0, 0.67,
0, 0, 0, 1.9, 1.9, 0, 0, 0, 1.23, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0), Species31 = c(0.72, 0.56, 1.19, 0, 0, 1.85,
0, 0.95, 0, 0, 0, 1.23, 0, 0, 0, 0, 0, 0, 1.5, 0, 0, 0, 0,
4.93), Species32 = c(0.71, 0.93, 0, 0, 0, 0.6, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 3.01, 1.09, 0, 0, 0, 0), Species33 = c(1.06,
1.33, 0.84, 0, 1.11, 0.6, 1.12, 1.56, 0, 1.45, 2.21, 0, 0,
2.21, 2.2, 1.69, 0, 0, 1.5, 0, 1.63, 2.11, 0, 1.18), Species34 = c(0.72,
0.89, 0.62, 0.63, 0.74, 1.21, 2.09, 1.32, 0, 1.38, 2.64,
1.06, 2.08, 0, 0, 0, 1.33, 0, 0, 1.37, 0, 2.37, 2.03, 2.02
), Species35 = c(0, 0.74, 0, 0, 0.74, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1.03, 0), Species36 = c(1,
1.63, 1, 0.86, 0.53, 0, 1.12, 1.32, 0, 0, 4.44, 0, 2.08,
0, 0, 0, 0, 0, 1.1, 2.03, 0, 2.72, 0, 0), Species37 = c(0,
1.25, 1.15, 0, 2.02, 0.92, 1.12, 0, 3.7, 1.08, 2.08, 0, 0,
0, 0, 1.52, 1.31, 0.96, 1.1, 1.06, 2.62, 2.37, 1.05, 1.44
), Species38 = c(0, 0.55, 0, 0, 0, 0, 1.9, 0, 1.08, 0, 0,
0, 1.04, 0, 0, 0, 1.09, 0, 0, 2.03, 0, 1.75, 0, 0), Species39 = c(0.7,
0.54, 0, 0.86, 0, 0, 0, 0, 0, 0, 0, 1.23, 0, 0, 1.23, 0,
0, 0, 0, 0, 0, 2.72, 0, 0), Species40 = c(0.95, 1.08, 1.39,
0, 0, 0, 1.11, 1.32, 1.33, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1.02,
0, 1.19, 0, 1.02, 0)), class = "data.frame", row.names = c(NA,
-24L))


The plots. Every point represents communities sampled in one year, and I used two different methods to add an ellipse to the data.

There seem to be differences in species communities, and I tried to test this statistically.

What I've done so far

1. General comparison

I used envfit in the vegan package, using the strata argument to restrict permutations by the repeated measures (Year).

envfit(data_dist ~ Site, data = sp.site,
strata = sp.site$Year, permutations = 9999)  If the test results are correct, the sites are significantly different.  Goodness of fit: r2 Pr(>r) Site 0.3261 0.0041 **  1. Pairwise comparison I wanted to know which sites are different, and which aren't, so I ran the test for all site comparisons. library(vegan) library(dplyr) data_dist = as.matrix((vegan::vegdist(sp.site[, -c(1:2)], "bray"))) site_combs <- combn(unique(sp.site$Site), 2)

df <- data.frame(Site1 = site_combs[1,], Site2 = site_combs[2,], envfit.p_value = NA)

for(i in 1:length(rownames(df))){
temp <- sp.site %>%
dplyr::filter(Site == df$$Site1[i] | Site == df$$Site2[i])
df$$envfit.p_value[i] <- as.numeric(vegan::envfit(data_dist ~ sp.site$$Site, strata = sp.site$$Year, perm = 999)$$factors$pvals) }  According to this, every single site combination is significantly different!  > df Site1 Site2 envfit.p_value 1 Site1 Site2 0.0032 2 Site1 Site3 0.0047 3 Site1 Site4 0.0039 4 Site1 Site5 0.0045 5 Site1 Site6 0.0031 6 Site2 Site3 0.0034 7 Site2 Site4 0.0041 8 Site2 Site5 0.0048 9 Site2 Site6 0.0032 10 Site3 Site4 0.0050 11 Site3 Site5 0.0027 12 Site3 Site6 0.0042 13 Site4 Site5 0.0043 14 Site4 Site6 0.0051 15 Site5 Site6 0.0038  Questions 1. Are those results reliable? I find it hard to believe that e.g. Site5 and Site6 are significantly different. Looking at the plot with the 95% confidence intervals, I'm even sceptical that there are statistically clear differences at all. 2. Are there better ways to compare the communities, accounting for the repeated measures design? 1 Answer The first issue is that there is an error in your pairwise comparison function - you are just repeating the "general comparison" test 15 times. The p-values are slightly different each time due to the random nature of the permutations. Your function should look like this: for(i in 1:length(rownames(df))){ temp <- sp.site %>% dplyr::filter(Site == df$$Site1[i] | Site == df$$Site2[i]) temp_data_dist <- as.matrix((vegan::vegdist(temp[, -c(1:2)], "bray"))) df$$envfit.p_value[i] <- as.numeric(vegan::envfit(temp_data_dist ~ temp$$Site, strata = temp$$Year, perm = 999)$$factors$pvals)
}


This will show that none of your pairwise comparisons have significant p-values - likely due to small sample size.

I'm not an expert in envfit, so may be incorrect, but I believe it is testing for how well the centroids of your factor (Site) fit the ordination you made - is that what you want? You are constraining it to only permute sites within the Year groupings - I'm not sure what this leaves it to test...? There will always be differences between two sampled communities - what do you mean when you say you want to test this statistically? Is there a threshold above which these differences are "significant" or "non-significant"?

I would recommend having a look at the vegan function "adonis" rather than "envfit", and also thinking about what it is you want to test with your permutations. For example, it would be possible to test whether your sites show a greater similarity to themselves when they are re-sampled in subsequent years than they do to the other sites. Adonis would do this by calculating distances between the same Site across multiple Years, and then seeing if this is significantly smaller than if you calculated distances to the other Sites instead. I can't really think of any other test you could do with this dataset.