non-Metric Multidimensional Scaling (nMDS)

Author

A. J. Smit

Published

2026/04/06

Material Required for This Chapter

Type	Name	Link
Theory	Numerical Ecology with R	See pages 145-151
Slides	nMDS lecture slides	💾 `BCB743_11_nMDS.pdf`
Data	The Doubs River data	💾 `Doubs.RData`

Tasks to Complete in This Chapter

Task F 1–5

Non-Metric Multidimensional Scaling (nMDS) is a rank-based indirect gradient analysis that uses a distance or dissimilarity matrix as its input (either pre-calculated using vegdist() or constructed internal to the metaMDS() function via the dist argument). Should one supply a ‘raw’ species × site table, the default dissimilarity matrix is Bray-Curtis dissimilarity, but any other dissimilarity index in vegdist() can be specified. Unlike other ordination methods such as Principal Component Analysis (PCA) and Correspondence Analysis (CA), which aim to maximise variance or correspondence between sites, nMDS focuses on representing the pairwise dissimilarities between sites in an ordination space. It does not use the raw distances or dissimilarities directly; instead, these values are replaced with their ranks, which is why the method is termed “non-metric.”

nMDS is the non-metric equivalent of Principal Coordinates Analysis (PCoA); in fact, PCoA is sometimes referred to as metric multidimensional scaling. PCoA and nMDS can both produce ordinations of objects from any distance or dissimilarity matrix. However, nMDS does not preserve the exact distances among objects in an ordination plot. Instead, it tries to represent the ordering (rank) relationships among objects as accurately as possible on a specified number of axes. This nonlinear mapping of dissimilarities onto a low-dimensional ordination space means that the Euclidean distances of points in the ordination space are rank-order similar to the original community dissimilarities.

The ordination space in nMDS is metric, but the regression used to fit the dissimilarities to this space is non-metric. This makes nMDS more robust than the eigen-value methods, especially when the data are not well-represented by a specific distribution as may sometimes be the case for ecological data. As with PCoA, it can handle quantitative, semi-quantitative, qualitative, or mixed variables, so we can flexibly apply it to many ecological problems.

A new concept, not seen in the eigen-approaches, is the idea of ‘stress.’ Stress quantifies the discrepancy between the observed dissimilarities and the distances in the ordination space. Stress is visually presented as the scatter of observed dissimilarities against the expected monotone regression. Lower stress values indicate a better fit of the data to the ordination space. Because rank orders of dissimilarities cannot be exactly preserved by rank-orders of ordination distances in low-dimensional space, some stress cannot be avoided.

nMDS does have some limitations. The rank-based approach means that information about the magnitude of differences between site pairs is lost and this can be a disadvantage when the actual distances are important for interpretation. Also, nMDS can be computationally intensive with large datasets or when trying to minimise stress through numerous iterations.

After performing nMDS, environmental interpretation can be facilitated using vegan’s envfit() and ordisurf() functions, as we have already seen in PCA, CA, and PCoA. As before, they allow for the fitting of environmental variables onto the ordination. This aids in visualising and understanding how the environmental variables influence the species ordination. Again, I must emphasis that this is not the same as doing a formal constrained ordination, which will be discussed in the next section (RDA and CCA). Additionally, if we require a statistical framework to assess the influence of categorical factors on the observed dissimilarities, we can use PERMANOVA (Permutational Multivariate Analysis of Variance) to test for differences that might be attributed to group effects. These ideas will be demonstrated below.

1 Set-Up the Analysis Environment

library(tidyverse)
library(vegan)
library(viridis)

# Files will be referenced using here::here() for absolute paths

2 The Doubs River Data

We continue to use the species data:

load(here::here("data", "BCB743", "NEwR-2ed_code_data", "NEwR2-Data", "Doubs.RData"))
spe <- dplyr::slice(spe, -8)

3 Do the nMDS

1spe_nmds <- metaMDS(spe, distance = "bray", trace = 0)
spe_nmds

1: I use trace = 0 to suppress the output of the iterations.


Call:
metaMDS(comm = spe, distance = "bray", trace = 0) 

global Multidimensional Scaling using monoMDS

Data:     spe 
Distance: bray 

Dimensions: 2 
Stress:     0.07376217 
Stress type 1, weak ties
Best solution was repeated 1 time in 20 tries
The best solution was from try 7 (random start)
Scaling: centring, PC rotation, halfchange scaling 
Species: expanded scores based on 'spe'

As always, reading the help file (accessible as ?vegan::metaMDS) is invaluable, as it is for all other ordination techniques.

There is a summary method available, but it is not particularly useful and I do not display the output here:

summary(spe_nmds)

Although summary(spe_nmds) does not return anything interesting, the species and site scores are nevertheless available directly through the scores() command, and they can be plotted as a layer in ggplot2 if needed:

scores(spe_nmds)

$sites
         NMDS1       NMDS2
1  -1.79011576  0.81435342
2  -1.14368915 -0.15605368
3  -1.00359122 -0.07436490
4  -0.62134248 -0.07823852
5   0.07205006  0.45941345
6  -0.42772071 -0.15265702
7  -0.87067317 -0.22994631
8  -0.01432404 -0.86245214
9  -0.52653413 -0.40369688
10 -1.00859344 -0.37649876
11 -0.97241028 -0.16336226
12 -1.16261896  0.09765316
13 -0.80289405  0.12428915
14 -0.49313700  0.17986162
15 -0.18897791  0.27611467
16  0.08275311  0.12262974
17  0.29913526  0.11272275
18  0.44597391  0.14938363
19  0.76740639  0.28448363
20  0.86245352  0.37752437
21  0.95540773  0.44476962
22  0.75364417 -1.44469525
23  1.13000242 -0.63139570
24  0.85946953 -0.87233669
25  0.93358857  0.12337103
26  0.97535775  0.36339901
27  1.02385650  0.37366465
28  0.79627294  0.55872608
29  1.06925042  0.58333812

$species
          NMDS1       NMDS2
Cogo -0.9122316 -0.10413690
Satr -1.1125981 -0.22548592
Phph -0.7890871 -0.32078845
Babl -0.5341364 -0.28947352
Thth -0.9380679 -0.09671147
Teso -0.5248947  0.16428177
Chna  0.9017546  0.36612148
Pato  0.5172679  0.38527609
Lele  0.3302886  0.28573606
Sqce  0.3586162 -0.15456245
Baba  0.7038968  0.47497292
Albi  0.7282122  0.47029668
Gogo  0.6853089  0.30428569
Eslu  0.6095852  0.42162555
Pefl  0.6172603  0.50364321
Rham  0.9667613  0.58278057
Legi  0.9562887  0.50926512
Scer  0.9659478  0.55191844
Cyca  0.9586891  0.62825289
Titi  0.7242518  0.41801123
Abbr  1.0825655  0.67973327
Icme  1.1282239  0.78234647
Gyce  1.0744386  0.40733462
Ruru  0.7594335  0.15152262
Blbj  1.0952088  0.58040604
Alal  0.9936141  0.02747419
Anan  1.0092010  0.61222867

See Numerical Ecology with R (pp. 145 to 149) for information about the interpretation of a nMDS and the ordination diagrams shown below.

4 Ordination Diagrams

We create the ordination diagrams as before, but new concepts introduced here are stress, Shepard plots, and goodness of fit (Figure 1). The stress indicates the scatter of observed dissimilarities against an expected monotone regression, while a Shepard diagram plots ordination distances against original dissimilarities and adds a monotone or linear fit line to highlight this relationship. The stressplot() function also produces two fit statistics. The goodness-of-fit of the ordination is measured as the $R^{2}$ of either a linear or a non-linear regression of the nMDS distances on the original ones.

opar <- par(no.readonly = TRUE)
par(mfrow = c(2, 2))
stressplot(spe_nmds, main = "Shepard plot")
ordiplot(spe_nmds, type = "t", cex = 1.2,
         main = paste0("nMDS stress = ", round(spe_nmds$stress, 2)))
gof <- goodness(spe_nmds)
plot(spe_nmds, type = "t", cex = 1.2, main = "Goodness of fit")
points(spe_nmds, display = "sites", cex = gof * 200)
# ...bigger bubbles indicate a worse fit
par(opar)

Figure 1: nMDS ordination plots of the Doubs River species data showing a Shepard plot (top, left), the ordination diagram (top, right), and goodness of fit (bottom, right).

A good rule of thumb: stress <0.05 provides an excellent representation in reduced dimensions, <0.1 is great, <0.2 is so-so, and stress <0.3 provides a poor representation.

We can also build ordination plots from scratch to suit specific needs:

pl <- ordiplot(spe_nmds, type = "none", main = "nMDS fish abundances ")
points(pl, "sites", pch = 21, cex = 1.75, col = "grey80", bg = "grey80")
points(pl, "species", pch = 21, col = "turquoise", arrows = TRUE)
text(pl, "species", col = "blue4", cex = 0.9)
text(pl, "sites", col = "red4", cex = 0.9)

nMDS ordination plot of the Doubs River species data assembled from scratch.

Or we can fit response surfaces using ordisurf() and project environmental drivers (Figure 2):

palette(viridis(8))
opar <- par(no.readonly = TRUE)
par(mar = c(4, 4, 0.9, 0.5) + .1, mfrow = c(2, 2))

invisible(ordisurf(spe_nmds ~ Satr, data = spe, bubble = 3,
                   family = quasipoisson, knots = 2, col = 6,
                   display = "sites", main = "Salmo trutta fario"))
abline(h = 0, v = 0, lty = 3)

invisible(ordisurf(spe_nmds ~ Scer, data = spe, bubble = 3,
                   family = quasipoisson, knots = 2, col = 6,
                   display = "sites", main = "Scardinius erythrophthalmus"))
abline(h = 0, v = 0, lty = 3)

invisible(ordisurf(spe_nmds ~ Teso, data = spe, bubble = 3,
                   family = quasipoisson, knots = 2, col = 6,
                   display = "sites", main = "Telestes souffia"))
abline(h = 0, v = 0, lty = 3)

invisible(ordisurf(spe_nmds ~ Cogo, data = spe, bubble = 3,
                   family = quasipoisson, knots = 2, col = 6,
                   display = "sites", main = "Cottus gobio"))
abline(h = 0, v = 0, lty = 3)

env <- env[-8, ] # because we removed the eighth site in the spp data

# A posteriori projection of environmental variables in a CA
# The last plot produced (CA scaling 2) must be active
spe_nmds_env <- envfit(spe_nmds, env, scaling = 2) # Scaling 2 is default
plot(spe_nmds_env)

# Plot significant variables with a different colour
plot(spe_nmds_env, p.max = 0.05, col = "red")
par(opar)

Figure 2: nMDS ordination plots with species response surfaces of the Doubs River species data emphasising four species of fish: A) Satr, B) Scer, C) Teso, and D) Cogo. D) additionally has the environmental vectors projected on the plot, with the significant vectors shown in red.

5 References

Reuse

CC BY-NC-SA 4.0

Citation

BibTeX citation:

@online{smit2026,
  author = {Smit, A. J. and J. Smit, A.},
  title = {Non-Metric {Multidimensional} {Scaling} {(nMDS)}},
  date = {2026-04-06},
  url = {https://tangledbank.netlify.app/BCB743/nMDS.html},
  langid = {en}
}

For attribution, please cite this work as:

Smit AJ, J. Smit A (2026) non-Metric Multidimensional Scaling (nMDS). https://tangledbank.netlify.app/BCB743/nMDS.html.

--- date: last-modified title: "non-Metric Multidimensional Scaling (nMDS)" author: "A. J. Smit" format: html: code-link: false --- ```{r code-brewing-opts, echo=FALSE} knitr::opts_chunk$set( comment = "R>", warning = FALSE, message = FALSE, fig.width = 4.5, fig.height = 2.625, out.width = "75%", fig.asp = NULL, # control via width/height dpi = 300 ) ggplot2::theme_set( ggplot2::theme_minimal(base_size = 8) ) ggplot2::theme_set( ggplot2::theme_bw(base_size = 8) ) ```  ::: callout-tip ## **Material Required for This Chapter** | Type | Name | Link | | :--- | :--- | :--- | | **Theory** | Numerical Ecology with R | See pages 145-151 | | **Slides** | nMDS lecture slides | [💾 `BCB743_11_nMDS.pdf`](../slides/BCB743_11_nMDS.pdf) | | **Data** | The Doubs River data | [💾 `Doubs.RData`](../data/BCB743/NEwR-2ed_code_data/NeWR2-Data/Doubs.RData) | ::: ::: {.callout-important appearance="simple"} ## Tasks to Complete in This Chapter * [Task F 1--5](tasks/Task_F.qmd) ::: Non-Metric Multidimensional Scaling (nMDS) is a rank-based indirect gradient analysis that uses a distance or dissimilarity matrix as its input (either pre-calculated using `vegdist()` or constructed internal to the `metaMDS()` function via the `dist` argument). Should one supply a 'raw' species × site table, the default dissimilarity matrix is Bray-Curtis dissimilarity, but any other dissimilarity index in `vegdist()` can be specified. Unlike other ordination methods such as Principal Component Analysis (PCA) and Correspondence Analysis (CA), which aim to maximise variance or correspondence between sites, nMDS focuses on representing the pairwise dissimilarities between sites in an ordination space. It does not use the raw distances or dissimilarities directly; instead, these values are replaced with their ranks, which is why the method is termed "non-metric." nMDS is the non-metric equivalent of Principal Coordinates Analysis (PCoA); in fact, PCoA is sometimes referred to as metric multidimensional scaling. PCoA and nMDS can both produce ordinations of objects from any distance or dissimilarity matrix. However, nMDS does not preserve the exact distances among objects in an ordination plot. Instead, it tries to represent the ordering (rank) relationships among objects as accurately as possible on a specified number of axes. This nonlinear mapping of dissimilarities onto a low-dimensional ordination space means that the Euclidean distances of points in the ordination space are rank-order similar to the original community dissimilarities. The ordination space in nMDS is metric, but the regression used to fit the dissimilarities to this space is non-metric. This makes nMDS more robust than the eigen-value methods, especially when the data are not well-represented by a specific distribution as may sometimes be the case for ecological data. As with PCoA, it can handle quantitative, semi-quantitative, qualitative, or mixed variables, so we can flexibly apply it to many ecological problems. A new concept, not seen in the eigen-approaches, is the idea of 'stress.' Stress quantifies the discrepancy between the observed dissimilarities and the distances in the ordination space. Stress is visually presented as the scatter of observed dissimilarities against the expected monotone regression. Lower stress values indicate a better fit of the data to the ordination space. Because rank orders of dissimilarities cannot be exactly preserved by rank-orders of ordination distances in low-dimensional space, *some* stress cannot be avoided. nMDS does have some limitations. The rank-based approach means that information about the magnitude of differences between site pairs is lost and this can be a disadvantage when the actual distances are important for interpretation. Also, nMDS can be computationally intensive with large datasets or when trying to minimise stress through numerous iterations. After performing nMDS, environmental interpretation can be facilitated using **vegan**'s `envfit()` and `ordisurf()` functions, as we have already seen in PCA, CA, and PCoA. As before, they allow for the fitting of environmental variables onto the ordination. This aids in visualising and understanding how the environmental variables influence the species ordination. Again, I must emphasis that this is not the same as doing a formal constrained ordination, which will be discussed in the next section (RDA and CCA). Additionally, if we require a statistical framework to assess the influence of categorical factors on the observed dissimilarities, we can use PERMANOVA (Permutational Multivariate Analysis of Variance) to test for differences that might be attributed to group effects. These ideas will be demonstrated below. ## Set-Up the Analysis Environment ```{r code-library-tidyverse, message=FALSE, warning=FALSE} library(tidyverse) library(vegan) library(viridis) # Files will be referenced using here::here() for absolute paths ``` ## The Doubs River Data We continue to use the species data: ```{r code-load-here-here-data} load(here::here("data", "BCB743", "NEwR-2ed_code_data", "NEwR2-Data", "Doubs.RData")) spe <- dplyr::slice(spe, -8) ``` ## Do the nMDS ```{r code-spe-nmds-metamds-spe-distance} spe_nmds <- metaMDS(spe, distance = "bray", trace = 0) # <1> spe_nmds ``` 1. I use `trace = 0` to suppress the output of the iterations. As always, reading the help file (accessible as `?vegan::metaMDS`) is invaluable, as it is for all other ordination techniques. There is a summary method available, but it is not particularly useful and I do not display the output here: ```{r code-summary-spe-nmds} #| eval: false summary(spe_nmds) ``` Although `summary(spe_nmds)` does not return anything interesting, the species and site scores are nevertheless available directly through the `scores()` command, and they can be plotted as a layer in **ggplot2** if needed: ```{r code-scores-spe-nmds} scores(spe_nmds) ``` See *Numerical Ecology with R* (pp. 145 to 149) for information about the interpretation of a nMDS and the ordination diagrams shown below. ## Ordination Diagrams We create the ordination diagrams as before, but new concepts introduced here are **stress**, **Shepard plots**, and **goodness of fit** (@fig-nmds-plots). The stress indicates the scatter of observed dissimilarities against an expected monotone regression, while a Shepard diagram plots ordination distances against original dissimilarities and adds a monotone or linear fit line to highlight this relationship. The `stressplot()` function also produces two fit statistics. The goodness-of-fit of the ordination is measured as the $R^{2}$ of either a linear or a non-linear regression of the nMDS distances on the original ones. ```{r fig-nmds-plots} #| fig-width: 10 #| fig-height: 10.75 #| fig.align: center #| fig.cap: "nMDS ordination plots of the Doubs River species data showing a Shepard plot (top, left), the ordination diagram (top, right), and goodness of fit (bottom, right)." opar <- par(no.readonly = TRUE) par(mfrow = c(2, 2)) stressplot(spe_nmds, main = "Shepard plot") ordiplot(spe_nmds, type = "t", cex = 1.2, main = paste0("nMDS stress = ", round(spe_nmds$stress, 2))) gof <- goodness(spe_nmds) plot(spe_nmds, type = "t", cex = 1.2, main = "Goodness of fit") points(spe_nmds, display = "sites", cex = gof * 200) # ...bigger bubbles indicate a worse fit par(opar) ``` A good rule of thumb: stress <0.05 provides an excellent representation in reduced dimensions, <0.1 is great, <0.2 is so-so, and stress <0.3 provides a poor representation. We can also build ordination plots from scratch to suit specific needs: ```{r code-pca-ordiplot} #| fig-width: 5.625 #| fig-height: 5.625 #| fig.align: center #| fig.cap: "nMDS ordination plot of the Doubs River species data assembled from scratch." pl <- ordiplot(spe_nmds, type = "none", main = "nMDS fish abundances ") points(pl, "sites", pch = 21, cex = 1.75, col = "grey80", bg = "grey80") points(pl, "species", pch = 21, col = "turquoise", arrows = TRUE) text(pl, "species", col = "blue4", cex = 0.9) text(pl, "sites", col = "red4", cex = 0.9) ``` Or we can fit response surfaces using `ordisurf()` and project environmental drivers (@fig-nmds-ordisurf): ```{r fig-nmds-ordisurf} #| fig-width: 10 #| fig-height: 10 #| fig.align: center #| fig.cap: "nMDS ordination plots with species response surfaces of the Doubs River species data emphasising four species of fish: A) Satr, B) Scer, C) Teso, and D) Cogo. D) additionally has the environmental vectors projected on the plot, with the significant vectors shown in red." palette(viridis(8)) opar <- par(no.readonly = TRUE) par(mar = c(4, 4, 0.9, 0.5) + .1, mfrow = c(2, 2)) invisible(ordisurf(spe_nmds ~ Satr, data = spe, bubble = 3, family = quasipoisson, knots = 2, col = 6, display = "sites", main = "Salmo trutta fario")) abline(h = 0, v = 0, lty = 3) invisible(ordisurf(spe_nmds ~ Scer, data = spe, bubble = 3, family = quasipoisson, knots = 2, col = 6, display = "sites", main = "Scardinius erythrophthalmus")) abline(h = 0, v = 0, lty = 3) invisible(ordisurf(spe_nmds ~ Teso, data = spe, bubble = 3, family = quasipoisson, knots = 2, col = 6, display = "sites", main = "Telestes souffia")) abline(h = 0, v = 0, lty = 3) invisible(ordisurf(spe_nmds ~ Cogo, data = spe, bubble = 3, family = quasipoisson, knots = 2, col = 6, display = "sites", main = "Cottus gobio")) abline(h = 0, v = 0, lty = 3) env <- env[-8, ] # because we removed the eighth site in the spp data # A posteriori projection of environmental variables in a CA # The last plot produced (CA scaling 2) must be active spe_nmds_env <- envfit(spe_nmds, env, scaling = 2) # Scaling 2 is default plot(spe_nmds_env) # Plot significant variables with a different colour plot(spe_nmds_env, p.max = 0.05, col = "red") par(opar) ``` ## References ::: {#refs} :::