Multidimensional Scaling (MDS) in R: cmdscale() & NMDS with vegan

Multidimensional scaling (MDS) takes a table of distances between objects and lays them out in 2D so that close points stay close and far points stay far. In R, base cmdscale() handles exact (metric) distances, and vegan::metaMDS() handles ranks (non-metric) for ecology, surveys, and any non-Euclidean similarity.

What does multidimensional scaling project to 2D?

You have distances, similarities, or a feature table you cannot picture. MDS turns that mess into a 2D scatter plot where geometric distance on the page approximates the dissimilarities you fed in. The classical version, called principal coordinates analysis when applied to a distance matrix, falls out of one call to cmdscale(). The fastest way to see what it does is the built-in eurodist matrix of road distances between 21 European cities.

Read what just happened. We never told cmdscale() about latitude or longitude. We only handed it a table of road distances. The function arranged the cities so geometric distance between any two points on the plot tracks the road distance you would drive between them. Stockholm sits up north on the plot because every distance to Stockholm is large; Madrid sits at one extreme of axis 1 because Madrid is far from everything in central Europe. The signs of the axes are arbitrary, so the result may appear east-west flipped from a real map. That is fine; only relative positions carry meaning.

How do you run classical MDS with cmdscale()?

Now that you have seen what classical MDS produces, look at how it gets there. cmdscale() doubles, centers, and eigendecomposes the distance matrix. The first k eigenvectors give you the coordinates; the matching eigenvalues tell you how much of the dissimilarity each axis recovers. Reading those eigenvalues is how you decide whether 2D was enough.

Use eig = TRUE to get the diagnostics back along with the coordinates. We will scale the four columns of USArrests first so no single variable dominates the Euclidean distance.

The first two eigenvalues swallow most of the variation, and the GOF field reports 0.885, meaning the 2D layout preserves about 88% of the squared input distances. That is solid. If GOF had come back below 0.6, you would either accept that 2D loses real structure or bump k to 3. The eigenvalues drop fast after the second one, so adding a third axis would buy little.

Use the matrix arr_mds$points directly when you want to colour, label, or feed coordinates into another model.

When should you switch from classical to non-metric MDS?

Classical MDS treats your distances as exact numbers. Halve a distance and the layout shifts accordingly. Non-metric MDS, called NMDS, treats your distances as a ranking only. It cares which pair is closest and which is farthest, not by how much. That distinction matters when distances are ordinal (Likert-scale similarities), come from a non-Euclidean metric like Bray-Curtis, or you suspect the absolute numbers are noisy but the ordering is reliable.

Figure 1: Pick metric (cmdscale) when distances are exact and Euclidean; pick non-metric (isoMDS or metaMDS) when only ranks matter or distances come from non-Euclidean metrics like Bray-Curtis.

The base R route to NMDS is MASS::isoMDS(), which implements Kruskal's algorithm. It iteratively rearranges points until the rank order of fitted distances best matches the rank order of input distances. Compare the layout against the classical run from the previous section.

isoMDS() returns coordinates plus a single stress value. Stress is reported as a percentage here (5.6 means 5.6%, well below the "usable" cutoff you will see in the next section). The point arrangement looks similar to the classical version but is not identical: Alaska and Mississippi swap relative positions because NMDS only honoured rank order, not absolute magnitudes.

How do you fit NMDS with vegan::metaMDS()?

MASS::isoMDS() is fine for textbook problems. For real ecology, survey, or community data, the standard tool is vegan::metaMDS(). It runs many random starts, picks the lowest-stress solution, applies a square-root transform when abundances are skewed, and rotates the axes so the first axis spans the longest direction. All of that happens behind one call. The default distance is Bray-Curtis, which is the right answer for community abundance counts.

The dune dataset that ships with vegan is a 20-by-30 matrix of plant cover scores at 20 Dutch dune meadows. Fit NMDS on it directly.

Read the printout top-down. Stress is 0.118, comfortably under the 0.20 "unreliable" line, so the 2D ordination is trustworthy. The line "Best solution was repeated 4 times" means metaMDS found the same near-optimal arrangement four times across its random starts. That is confirmation the result is stable and not a local minimum. metaMDS also automatically applied wisconsin(sqrt(...)) standardisation, which dampens the influence of dominant species.

Now plot the ordination. metaMDS supplies coordinates for both sites (rows) and species (columns), so you can show which species drive which patches of the plot.

Sites near each other share community composition. Species names sit at the centroid of the sites where they are abundant. A species at the far left of axis 1, surrounded by sites also far left, is the species defining that gradient end. That is the entire point of an ordination plot: read it as a community-by-community map.

How do you read stress and Shepard plots to judge MDS fit?

Stress is the single number that tells you whether your MDS is worth interpreting. It measures how badly the rank order of fitted 2D distances disagrees with the rank order of input distances. Lower is better. The standard formula:

$$\text{stress} = \sqrt{\frac{\sum_{i

Where:

• $d_{ij}$ = the input dissimilarity between objects $i$ and $j$
• $\hat d_{ij}$ = the monotone regression fit of input distances on layout distances
• the sum runs over all unordered pairs

Use this rule-of-thumb table from Kruskal:

A Shepard plot pairs every input distance with its fitted layout distance. If MDS fit perfectly, those points would lie on a straight line. Run stressplot() on a metaMDS result to see it.

Read the plot two ways. First, the points should hug the step line; wide vertical scatter at any given x-value means MDS could not honour those rank orders. Second, the printout shows two correlations: a non-metric fit (R² of the monotone fit) and a linear fit. For dune you should see non-metric R² around 0.99 and linear R² around 0.95, both excellent.

How do you choose the right distance metric for MDS?

The distance metric is the most important choice in any MDS workflow. The algorithm only visualises whatever you feed it. Pick the wrong metric, and even a stress of 0.05 means nothing.

Compare three of these on the dune community matrix and see how the choice changes stress.

Bray-Curtis wins on this community data. Jaccard is close because dune's presence patterns track abundance patterns. Euclidean is worst because it treats a column with cover score 9 as nine times more important than a column with cover score 1, which inflates the distances driven by a few common species and hides the rare-species turnover that defines community structure.

Practice Exercises

Exercise 1: MDS of the airquality monthly profiles

Take the built-in airquality data and build an MDS that places each month as a single point. Aggregate to monthly means (drop NAs first), scale the columns, compute Euclidean distance, run cmdscale(), and plot months coloured by season. Save the result to my_mds.

Exercise 2: Bray-Curtis vs Jaccard on lichen communities

The varespec dataset in vegan records cover values for 44 lichen species across 24 sites. Fit two NMDS ordinations: one with Bray-Curtis (counts), one with Jaccard binary (presence/absence). Save them as vare_bray and vare_jac. Print both stress values and decide which metric better fits this data.

Exercise 3: Procrustes correlation between cmdscale and isoMDS

Run both cmdscale() and MASS::isoMDS() on a scaled-Euclidean distance from iris[, 1:4]. Save the layouts as iris_cmd and iris_iso. Use vegan::procrustes() to compute the correlation between the two layouts and print it. Decide whether the choice of metric vs non-metric mattered for this dataset.

Complete Example

Here is a full NMDS pipeline on the dune community data, the same workflow you would run in a real ecology paper. The dune.env companion table records management type for each site, so we can colour the ordination by environmental gradient.

Now you have a fully interpretable plot. Sites cluster by management practice, and the envfit() overlay tells you that soil moisture explains 61% of the layout variance and is highly significant. Soil A1 horizon thickness explains 42% and is also significant. The arrows on the plot point in the direction of increasing values for each environmental variable. That is a publication-grade community ecology figure in fewer than 15 lines of code.

Summary

Figure 2: The five-step MDS pipeline: data, distance, MDS fit, fit diagnostics, plot.

References

1. Cox, T.F. and Cox, M.A.A. Multidimensional Scaling, 2nd Edition. Chapman & Hall/CRC (2001).
2. Borg, I. and Groenen, P.J.F. Modern Multidimensional Scaling: Theory and Applications, 2nd Edition. Springer (2005).
3. Oksanen, J. et al. vegan: Community Ecology Package documentation. Link
4. R Documentation. cmdscale reference. Link
5. Venables, W.N. and Ripley, B.D. Modern Applied Statistics with S, 4th Edition. Springer (2002). Chapter 11 covers MDS via isoMDS().
6. Kruskal, J.B. Nonmetric multidimensional scaling: a numerical method. Psychometrika 29 (1964): 115-129.
7. Legendre, P. and Legendre, L. Numerical Ecology, 3rd English Edition. Elsevier (2012). The reference text for ecological ordination.

Continue Learning

• t-SNE and UMAP in R: sibling dimensionality-reduction methods designed for very-high-dimensional data and tight local neighbourhoods.
• PCA in R: when your distances are exactly Euclidean and linear, PCA is the faster cousin of classical MDS.
• Cluster Analysis in R: natural pairing for NMDS layouts with k-means or hierarchical grouping.