Vegan Multivariate Analysis

NMDS

2023-08-13

Ordination with vegan

Data

data(dune)
data(dune.env)

data(varespec)
data(varechem)

Analysis

M1 <- dune
M2 <- dune.env

nmds <- metaMDS(M1 , distance = "bray", k=2)

Run 0 stress 0.1192678 
Run 1 stress 0.1183186 
... New best solution
... Procrustes: rmse 0.02026967  max resid 0.06495895 
Run 2 stress 0.2069997 
Run 3 stress 0.1192679 
Run 4 stress 0.1183186 
... New best solution
... Procrustes: rmse 7.908497e-06  max resid 2.407945e-05 
... Similar to previous best
Run 5 stress 0.1192678 
Run 6 stress 0.1183186 
... Procrustes: rmse 6.91075e-07  max resid 2.078912e-06 
... Similar to previous best
Run 7 stress 0.1183186 
... Procrustes: rmse 1.28847e-05  max resid 4.209671e-05 
... Similar to previous best
Run 8 stress 0.1183186 
... Procrustes: rmse 1.742262e-05  max resid 5.521291e-05 
... Similar to previous best
Run 9 stress 0.1183186 
... New best solution
... Procrustes: rmse 2.518232e-06  max resid 6.11417e-06 
... Similar to previous best
Run 10 stress 0.1939202 
Run 11 stress 0.1192679 
Run 12 stress 0.1183186 
... Procrustes: rmse 1.006626e-05  max resid 3.418378e-05 
... Similar to previous best
Run 13 stress 0.2389616 
Run 14 stress 0.1183186 
... Procrustes: rmse 1.442631e-05  max resid 4.421442e-05 
... Similar to previous best
Run 15 stress 0.1183186 
... Procrustes: rmse 1.08212e-05  max resid 3.343553e-05 
... Similar to previous best
Run 16 stress 0.1183186 
... New best solution
... Procrustes: rmse 2.652172e-06  max resid 8.530837e-06 
... Similar to previous best
Run 17 stress 0.1183186 
... Procrustes: rmse 5.245032e-06  max resid 1.663786e-05 
... Similar to previous best
Run 18 stress 0.1183186 
... Procrustes: rmse 1.597972e-05  max resid 5.013583e-05 
... Similar to previous best
Run 19 stress 0.1183186 
... Procrustes: rmse 6.85407e-06  max resid 2.214136e-05 
... Similar to previous best
Run 20 stress 0.1808911 
*** Best solution repeated 4 times

EnvFit

ord.fit <- envfit(nmds ~ Management, data=dune.env, perm=999)
ord.fit

***FACTORS:

Centroids:
               NMDS1   NMDS2
ManagementBF -0.4534 -0.0102
ManagementHF -0.2636 -0.1282
ManagementNM  0.2958  0.5790
ManagementSF  0.1506 -0.4670

Goodness of fit:
               r2 Pr(>r)  
Management 0.4134  0.011 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Permutation: free
Number of permutations: 999

Plot

co= c("blue", "red", "green", "orange")
shape=c(18,16,19,14)
plot(nmds$points, col=co[M2$Management], pch=shape[M2$Management], cex=1.2, main="community composition", xlab= "axis 1", ylab="axis 2")

ordiellipse(nmds, col= co , groups= M2$Management, label = TRUE)

txt= c("HF","BF","NM","SF")
legend('topleft',legend= c("HF","BF","NM","SF") , pch=c(18,16,19,14), col=co, cex=1, bty= "y")

Testing for differences (Permutations)

#Bootstrapping and testing for differences between the groups
fit <- adonis2(M1 ~ Management, data= M2, permutations=999, method="bray")
fit

Permutation test for adonis under reduced model
Terms added sequentially (first to last)
Permutation: free
Number of permutations: 999

adonis2(formula = M1 ~ Management, data = M2, permutations = 999, method = "bray")
           Df SumOfSqs      R2      F Pr(>F)   
Management  3   1.4686 0.34161 2.7672  0.004 **
Residual   16   2.8304 0.65839                 
Total      19   4.2990 1.00000                 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

#Check assumption of homogeneity of multivariate dispersion
distances_data <- vegdist(M1)
anova(betadisper(distances_data, M2$Management))

Analysis of Variance Table

Response: Distances
          Df  Sum Sq  Mean Sq F value Pr(>F)
Groups     3 0.13831 0.046104  1.9506 0.1622
Residuals 16 0.37816 0.023635               

ggplot2

#Extract the axes scores
datascores <- as.data.frame(scores(nmds$points))  #extract the site scores

#Add/calculate spider diagram
scores <- cbind(as.data.frame(datascores), Management = M2$Management)
centroids <- aggregate(cbind(MDS1, MDS2) ~ Management, data = scores, FUN = mean)
seg <- merge(scores, setNames(centroids, c('Management','oNMDS1','oNMDS2')),
             by = 'Management', sort = FALSE)

#plot
p <- ggplot(data=scores, aes(x = MDS1, y = MDS2, colour = Management)) +
  geom_segment(data = seg,
               mapping = aes(xend = oNMDS1, yend = oNMDS2)) + # add spiders
  geom_point(data = centroids, size = 4) +                    # add centroids
  geom_point() +                                              
  coord_fixed() +                                              
  theme_bw() + 
  theme(legend.position="right",legend.text=element_text(size=10),legend.direction='vertical')

ggplotly(p)