Script

IntroOrdinations.R

Video lecture

Video on YouTube

Packages and data set-up

pacman::p_load(tidyverse, vegan) 
data("varechem") 
chem_d <- select(varechem, N:Mo)
as_tibble(chem_d)

## # A tibble: 24 x 11
##        N     P     K    Ca    Mg     S    Al    Fe    Mn    Zn    Mo
##    <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
##  1  19.8  42.1  140.  519.  90    32.3  39    40.9  58.1   4.5   0.3
##  2  13.4  39.1  167.  357.  70.7  35.2  88.1  39    52.4   5.4   0.3
##  3  20.2  67.7  207.  973. 209.   58.1 138    35.4  32.1  16.8   0.8
##  4  20.6  60.8  234.  834  127.   40.7  15.4   4.4 132    10.7   0.2
##  5  23.8  54.5  181.  777  126.   39.5  24.2   3    50.1   6.6   0.3
##  6  22.8  40.9  171.  692. 151.   40.8 105.   17.6  43.6   9.1   0.4
##  7  26.6  36.7  171.  739.  94.9  33.8  20.7   2.5  77.6   7.4   0.3
##  8  24.2  31    138.  395.  45.3  27.1  74.2   9.8  24.4   5.2   0.3
##  9  29.8  73.5  260   749. 105.   42.5  17.9   2.4 107.    9.3   0.3
## 10  28.1  40.5  314.  541. 119.   60.2 330.  110.   61.7   9.1   0.5
## # ... with 14 more rows

Fit a cluster diagram that we’ll end up using later:

chem_m <- vegdist(chem_d, method="euclidean")
chem_clust <- hclust(chem_m, method="average")

Fitting an ordination

A common ordination method is Principal Components Analysis, or PCA. All ordinations are based on an underlying distance matrix, and PCA is no different–it is based on the Euclidean distance measure. So while many ordination methods allow the user to specify the distance measure, when Euclidean is specified, the result is a PCA no matter what function produced it.

We use here vegan::capscale to perform our PCA because capscale offers the most flexibility for conducting a variety of ordinations in \({\bf\textsf{R}}\), even ones beyond what we’ll cover in this course. The base stats package in \({\bf\textsf{R}}\) provides a few functions to conduct a PCA², but they are limited to PCA specifically, and it is an objective of this course to introduce you to the vegan package.

capscale returns easy-to-interpret ordination results, including PCA if one sets the arguments correctly. These settings dial back a lot of the full functionality of capscale, which is actually a flexible and powerful function.

Setting 1: Perform unconstrained ordination

capscale is designed to perform constrained ordinations aka direct gradient analysis. I am not a fan of constrained ordinations and do not cover them here, but it is relevant because capscale uses formula notation Y ~ x in which community data \(\mathbb{Y}\) are fit against one or more constraining variables \(x\).

Because PCA is an unconstrained ordination we use Y ~ 1 and capscale forgets about constraints.

Setting 2: Perform Euclidean-based ordination

capscale can perform an ordination based on any of the \(15+\) distance measures listed in ?vegdist, and any you create yourself with designdist. Since PCA is an ordination based on the Euclidean distance, we simply set method = "euc":

chem_pca <- capscale(chem_d ~ 1, method = "euc")

Biplot

Dissimilarities are calculated in the distance matrix–if one wants to know exactly how similiar or dissimilar two sites are, one can look it up in the distance matrix. But matrices are unwieldy and a series of numbers does little to inform general questions about trends in the data across all the sites.

At its core, an ordination is a tool to identify where variation occurs in a multivariate dataset and reduce that variation into just a few composite variables that let one see the variation in fewer dimensions, ideally few enough that they fit on a screen or page. To see this, we need a way to graph the ordination.

The basic graph of an ordination result is called a biplot, but the bi part doesn’t refer to the two axes plotted, \(\textsf{x}\) and \(\textsf{y}\). Rather, biplot refers to the two components of information in the plot: site scores, and species scores.

• Site scores are the main result of the ordination.
  • Numerically, site scores represent the location of each sampled point, or site, in ordination space.
  • Conceptually, they represent the dissimilarity among the sites. Samples that are more similar are closer in the ordination, and samples that are more dissimilar are farther apart in the ordination.
  • Because each row in a data frame of multivariate data represents one sample point, the length of the site scores equals that of the wide-format data frame.
• Species scores, when returned by an ordination,³, help one interpret which columns of the multivariate data correspond to variation in the ordination.
  • They are called species scores because in conventional community ecology, the columns in multivariate datasets are the presence/absence or abundance in plant or animal species. But the columns can just as easily be trait data, etc.
  • Scores further from the origin denote species with relatively greater contributions to overall variation.
  • Especially at the periphery of an ordination–where most of the variation is apparent–sites closer to a species score likely tended to have a higher abundance of that species than, say, a site on the opposite side of the ordination.

By default, species scores are plotted as the column names, which is fairly easy to interpret, but site scores are plotted as row numbers. While base stats does have a biplot function…

biplot(chem_pca, col=c("blue", "red")) 

… graphics::plot can also handle vegan objects:

plot(chem_pca)

plot also offers more control, such as over whether one or both type of scores are shown…

plot(chem_pca, display = "sites")

plot(chem_pca, display = "species")

… and can even allow a blank plot (type='n') that can be filled with sequential calls to text:

plot(chem_pca, type="n") # Empty plot
  text(chem_pca, display = "species", 
       col="red", cex=1.5)
  text(chem_pca, display = "sites", 
       col="blue", cex=1.5)

Important note on axis labels

All axes labels are printed with \(\textsf{MDS}\), which stands for Multi-dimensional Scaling. MDS refers to the general approach to unconstrained ordination employed by capscale, and \(\textsf{MDS}\) would be returned for any distance measure. Many times, PCAs are plotted as \(\textsf{PCA2 x PCA1}\), and one might consider changing the axis labels to ensure readers remember you are reporting a PCA, even though as savvy vegan users we know an MDS fit with the Euclidean distance measure is a PCA.

Interpreting scores

Ordinations are about visualizing patterns in multivariate data, but two things are important for ordinations to be useful:

• Insight into what accounts for the patterns
• Confidence that the picture is meaningful

Both of these are addressed by analyzing the site and species scores, and both begin by retrieving them with vegan::scores and the display= argument:

Site scores (display = "sites"); note rownames as row numbers from original dataset:

scores(chem_pca, display = "sites", choices = c(1:10)) %>%
  head %>% round(1) %>% pander::pander()

Species scores (display = "species"); note rownames as column names from original dataset:

scores(chem_pca, display = "species", choices = c(1:10)) %>%
  round(2) %>% pander::pander() 

In metric ordinations such as those returned by capscale, there are \(n-1\) dimensions, with \(n=\) the number of species (columns). Only the first few are of importance, so we narrow it here to 10 via choices=.

summary(chem_pca)$cont %>% 
  as.data.frame %>%       
    round(2) %>%
    pander::pander() 

screeplot(chem_pca, type="lines")

Ordination vs. clustering

Ordination projections and cluster diagrams are two ways to view the relationships in the distance matrix. While the cluster diagram might have the best fidelity in terms of representing the actual distances between pairs of sites, it is difficult to compare relative differences between more than two sites or see how groups relate to each other, because the information is conveyed by how high up the branch occurs. As two-dimensional projections of points arranged by an algorithm in multi-dimensional space, ordinations are often an abstraction of the distance matrix that better shows relative relationships. If we make sure the ordination is a good representation of the variation in the distance matrix, we can live with the abstraction of the projection.

While you will rarely do this, we can see how the relationships in the cluster diagram are arranged in the ordination. This exercise also introduces the important group of ordi* functions vegan provides for showing groups and clusters in the ordination, which visualize the groups and gradients we test post-hoc for statistical correlations with the latent variable (and thus interpret as possible sources of variation in the multivariate data).

First, we use an old-school, base graphics way to put plots side-by-side.⁴ Then we use vegan::ordicluster to connect the site scores in the ordination chem_pca by the dendrogram stored in the cluster diagram object, chem_clust.

par(mgp=c(4, 1, 0), 
    mar=c(6, 6, 1, 1), 
     mfrow=c(1,2), 
    las=1, cex.lab=1, cex.axis=1) 
plot(chem_clust, las=1, 
      labels=chem_d$name, 
      xlab="Sample", 
     ylab="Euclidean distance")
plot(chem_pca, display = "sites", las=1) 
ordicluster(chem_pca, chem_clust)

In the cluster diagram, we see that sites \(24\) and \(11\) are paired together but off on their own; this is apparent in the ordination projection, as well. Likewise sites \(2\), \(13\), and \(4\) share a branch in the dendrogram, and group closely together in the ordination.

The utility of the ordination comes in the visualization of very dissimilar sites, which are less obvious in the cluster diagram. Take, for example, how sites \(3\) and \(7\) are on the left of the primary axis of variation (\(\textsf{MDS1}\)), while \(24\) and \(11\) are on the far right. This is to be determined as these site groups to be maximally dissimilar, but the clusters appear side-by-side in the dendrogram! One must look up, and see thatalthough the clusters are printed side-by-side, they actually only connect at a very high level (large Euclidean distance, high dissimilarity).

Scaled data

We have applied scaling functions to our data a couple of times before, most recently prior to calculating the distance matrix for cluster analysis. Scaling data can have major implications on results, especially for multivariate analyses.

There are arguments for and against scaling, but for our purposes here, I generally advise scaling data prior to multivariate analysis. It has its greatest effect–and I would often argue, is most necessary–when the values across columns in the data occur on different scales.

Recall our initial ordination on unscaled data:

plot(chem_pca)

summary(chem_pca)$species %>%  
  as.data.frame %>%   
  select(MDS1) %>%             
  round(2) %>%
  t() %>% pander::pander() 

The plot and primary axis scores indicate that the Calcium \(\text{Ca}\) gradient has the strongest influence on variation in the dataset (by an order of magnitude; \(31.1~vs.~5.59\). One might reasonably conclude, then, that calcium is the most important of the nutrients measured.

But let’s look back at the original data:

head(chem_d) %>% pander::pander() 

The \(\text{Ca}\) column has much greater values than most other columns: min \(\text{Ca} = 188.5\); max \(\text{N} = 33.1\), max \(\text{Mo} = 1.1\). Curiously, these lower-range elements cluster in the middle of the ordination. What if the apparent importance of calcium is simply an artefact of having the largest values, not the widest gradient of relative values, which is what the ordination is looking for?

chem_s <- scale(chem_d)
chem_pca2 <- capscale(chem_s ~ 1, method = "euc")
plot(chem_pca2)

summary(chem_pca2)$species %>%  
  as.data.frame %>%   
  select(MDS1) %>%             
  round(2) %>%
  t() %>% pander::pander()

summary(chem_pca2)$cont %>% 
  as.data.frame %>%  
    .[1:6] %>%  round(2) %>% pander::pander() 

screeplot(chem_pca2, type="lines") 
  abline(a=1, b=0, lty=2)