##################################### #Load libraries ##################################### library(BiodiversityR) library(vegan) library(labdsv) ##################################### #Load Data ##################################### community<-read.csv(file="example.csv", row.names='Sites', sep=",", header=TRUE) # Calculate dissimilarity index, you may want to experiment with other indices. distance<-vegdist(community, method="bray") # Do a standard NMDS, k=2 and maximum iterations at 2 nmds<-isoMDS(distance, k=2, tol=0.001, maxit=2) nmds ordiplot(nmds, main=sprintf("k=2 stress=%.3f",nmds$stress)) # Same analysis, greater number of iterations allowed nmds<-isoMDS(distance, k=2, tol=0.001, maxit=50) nmds ordiplot(nmds, main=sprintf("k=2 stress=%.3f",nmds$stress)) # Calculate the weighted average scores for each species # using function wascores. # exapand=TRUE scales scores to match ordination variance, # you may want to try it both ways in your own data. sp_wa_scores<-wascores(nmds$points, community, expand=FALSE) # plot the ordination and then add the species weighted average scores ordiplot(nmds, type="text", main=sprintf("k=2 stress=%.3f",nmds$stress)) text(sp_wa_scores, rownames(sp_wa_scores), col="blue") # Gather data for and then look at the Shepard plot (assess stress) stress_plot_values<-Shepard(distance,nmds$points) plot(stress_plot_values, xlab="Original Bray Curtis Distance", ylab="Distance in NMDS Space", main=sprintf("k=2 stress=%.3f",nmds$stress)) # Sensitivity analysis, assess stress values with variable numbers of axes and iterations # Run through a series of NMDS using from 1 to 50 as the maximum number of iterations # Keep track of stress values, we anticipate stress will drop as we allow # more iterations. # Set k=2 and use the maxit variable to control # of iterations stress_values<-numeric(50) # Calculate dissimilarity index distance<-vegdist(community, method="bray") for (n in 1:50) { stress_values[n] <- isoMDS(distance, k=2, maxit=n) $ stress } plot(stress_values, xlab="Number of iterations", ylab="Stress") # Run through a series of NMDS setting k from 1 to 10 axes stress_values<-numeric(10) # Calculate dissimilarity index distance<-vegdist(community, method="bray") for (n in 1:10) { stress_values[n] <- isoMDS(distance, k=n, maxit=50) $ stress } plot(stress_values, xlab="Number of axes", ylab="Stress") # k=2 nmds and Shepard plot k2_nmds <- isoMDS(distance, k=2, maxit=50) stress_plot_values<-Shepard(distance,k2_nmds$points) plot(stress_plot_values, xlab="Original Bray Curtis Distance", ylab="Distance in NMDS Space", main=sprintf("k=2 stress=%.3f",k2_nmds$stress)) # k=10 nmds and Shepard plot k10_nmds <- isoMDS(distance, k=10, maxit=50) stress_plot_values<-Shepard(distance,k10_nmds$points) plot(stress_plot_values, xlab="Original Bray Curtis Distance", ylab="Distance in NMDS Space", main=sprintf("k=10 stress=%.3f",k10_nmds$stress)) # What about the Shepard plot for just the first two axes of the k=10 configureation? stress_plot_values<-Shepard(distance,k10_nmds$points[,1:2]) plot(stress_plot_values, xlab="Original Bray Curtis Distance", ylab="Distance in First 2 Dimensions of NMDS Space", main=sprintf("k=10 stress=%.3f",k10_nmds$stress)) ##################################### # Function metaMDS ##################################### metanmds<-metaMDS(community,k=2,distance="bray", noshare=0.1, autotransform=FALSE, plot=TRUE) metanmds # This is the closest we have come in any ordination to producing the single linear # gradient in our artificial dataset. ordiplot(metanmds, type="text", main=sprintf("k=2 stress=%.3f",metanmds$stress)) # Goodness of fit values for each site gof<-goodness(metanmds) points(metanmds, display="sites", cex=gof*1000) # Calculate the stress values by hand # Reconstruct community distance as in the nmds, accounting for zeroes comm_dist<-metaMDSredist(metanmds) # Calculate the ordination distance ord_dist<-vegdist(metanmds$points[,1:2],method="euclidean") # Extract the residuals from a monotonic regression (recall nmds works with ranks) dhat<-isoreg(comm_dist, ord_dist) # plot the monotonic regression plot(dhat) # calculate stress with the residuals of this regression and ordination distnaces sqrt(sum(residuals(dhat)^2)/sum(ord_dist^2)) stressplot(metanmds) # Not a true % variance explained, but an estimate from # correlating the original distance matrix with the ordination distance distance<-vegdist(community, method="bray") plot(distance,ord_dist) cor(distance,ord_dist)^2 # To estimate % variance on each axis, do the same but use scores for one # axis at a time ax1_dist<-vegdist(metanmds$points[,1],method="euclidean") plot(distance,ax1_dist) cor(distance,ax1_dist)^2 ax2_dist<-vegdist(metanmds$points[,2],method="euclidean") plot(distance,ax2_dist) cor(distance,ax2_dist)^2