0
$\begingroup$

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.

enter image description here


Task

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?
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.