Unconstrained ordination is a method to arrange sites by their assemblage dissimilarity. Here, we learn how to constract an ordination plot and to visualize an association of environmental variables to the ordination.

Dissimilarity of assemblages and nMDS plot

Non-metric multidimentional scaling (nMDS) constructs an ordination plot showing sites of similar assemblages close together and sites of different assemblages far from each other.

First, load the vegan package require(vegan) and a data file. In the example below, we use a csv file with the first row containing site names. The argument row.names=1 set the site names as row names. Check the data by str(dat) and head(dat) before following calculations.

dat <- read.csv("spcdat.csv", header = TRUE, row.names=1)

Then, get results of nMDS by metaMDS() by setting a dissimilarity by an argument distance=, and show the nMDS plot.
The example below uses the Morisita-Horn dissimilarity.

res1.mds <- metaMDS(dat, distance="horn", autotransform = FALSE, trace = 0)
plot(res1.mds, type="t", display="sites")

Here is a resulting two dimentional nMDS ordination plot.

Original 102 dimentional data, which is the number of species, is now reduced into the two dimentional plot. On this reducing process, we lost some imformation. Degree of this loss is expressed by a stress value.

res1.mds

## 
## Call:
## metaMDS(comm = dat, distance = "horn", autotransform = FALSE,      trace = 0) 
## 
## global Multidimensional Scaling using monoMDS
## 
## Data:     dat 
## Distance: horn 
## 
## Dimensions: 2 
## Stress:     0.1590454 
## Stress type 1, weak ties
## Two convergent solutions found after 20 tries
## Scaling: centring, PC rotation, halfchange scaling 
## Species: expanded scores based on 'dat'

In this example, the stress value is 0.1590454. Generally, stress value < 0.2 is acceptable. Otherwise, increase the dimention of nMDS. To perform nMDS with 3-dimention, set the argument k=3, and run metaMDS(dat, distance="horn", k=3, autotransform = FALSE, trace = 0). The resulting stress value is 0.0948034.

PCoA ordination: a way to dbRDA

In the example above, the function metaMDS proessed a raw assemblage data. This function also deal with a distance matrix between sites.

dis.mh <- vegdist(dat, method="horn")
dis.bc <- vegdist(dat, method="bray")

These lines produce distance matrices based on the Morishita-Horn dissimilarity and the Bray-Curtis dissimilarity. Then function metaMDS can read these distance matrices and produce the same nMDS as before. The followin lines show a nMDS plot by the Bray-Curtis dissimilarity.

res2.mds <- metaMDS(dis.bc, trace=0)
plot(res2.mds, type="t", display="sites")

Using a distance matrix, we perform the Principal Coordinate Analysis (PCoA) by a function cmdscale. Following example shows a plot by the Morisita-Horn dissimilarity by using an argument dis.mh. When you would like to use the Bray-Curtis dissimilarity, use an argument dis.bc. An argument k=2 indicates to calculate up to 2nd axix of the PCoA ordination.

res.pco <- cmdscale(dis.mh, k=2, eig=TRUE)
ordiplot(res.pco, type="t")

PCoA is also called metric multidimentional scaling, which is the base of the distance-based redundancy anaysis (dbRDA) to analyze effects of environmental factors.

We can evaluete differences betwenn the nMDS plot and the PCoA plot by the Procrustes Test protest().

pro <- protest(res1.mds, res.pco)
pro

## 
## Call:
## protest(X = res1.mds, Y = res.pco) 
## 
## Procrustes Sum of Squares (m12 squared):        0.2582 
## Correlation in a symmetric Procrustes rotation: 0.8613 
## Significance:  0.001 
## 
## Permutation: free
## Number of permutations: 999

plot(pro)

Procrustes test shows the probability 0.001, indicating a significant difference between these two plots.

Relationships between the assemblage ordination and environmental factors

We have an expectation that environmental factors may affect any differences in assemblages. So, we would like to know there are any relationships between the assemblage ordination and environmental variables.

First, read an environmental data matrix.

edat <- read.csv("envdat.csv", header=TRUE)

Then, examine linear correlation between the nMDS axses and environmental variables by a function envfit().

ef <- envfit(res2.mds, edat, permu=4999)
ef

## 
## ***VECTORS
## 
##         NMDS1    NMDS2     r2 Pr(>r)  
## PSU   0.29621 -0.95512 0.1732 0.3232  
## PH    0.86206 -0.50680 0.4043 0.0454 *
## ORP3  0.25097  0.96799 0.2325 0.2066  
## ORP5 -0.01330  0.99991 0.3348 0.0926 .
## ORP8 -0.00882  0.99996 0.3121 0.1120  
## IL   -0.22412  0.97456 0.3402 0.0886 .
## Chla -0.67474  0.73806 0.0363 0.7982  
## MdG  -0.72298 -0.69087 0.3608 0.0348 *
## PM    0.16020 -0.98708 0.2993 0.1084  
## TN   -0.04183  0.99912 0.1270 0.4814  
## TOC  -0.91193 -0.41035 0.1318 0.4850  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Permutation: free
## Number of permutations: 4999

ordiplot(res2.mds, type="t", display="sites")
plot(ef, p.max=0.1)

The above result shows that PH, ORP5, IL, and MdG shows a high correlation. On the plot, vectors show the direction and intensity of their increase. The next lines show a surface fitting of these variables.

sgev <- names(ef$vectors$r[ef$vectors$pvals<0.1])
ordiplot(res2.mds, display="sites")
for (i in 1:length(sgev)){
　fmi <- formula(paste("res2.mds~", sgev[i]))
　tmp <- with(edat, ordisurf(fmi, add=TRUE, col=i+1))
}

If you would like to see the contour of each variables, set as with(edat, ordisurf(res2.mds, Ph, add=TRUE, col="red")).

The uneven and curved contous indicate that the variables do not have simple linear relationships.

---

References

• Borcard D, Gille F, Legendre P (2011) Numerical Ecology with R. Springer
• Legendre P, Legendre L (2012) Numerical Ecology. Elsevier
• Magurran AE, McGill BJ (2011) Biological Diversity. Oxford Univ.

Back to "Analyses on Ecological Assemblages by R"