This practical will use the PONDS dataset to demonstrate methods of indirect gradient analysis (PCA, CA, and DCA) of species and environmental data. The file pondsenv.csv contains the species data (48 taxa) and pondsenv.csv contains the transformed environmental variables (15 variables) for 30 sites. You will use vegan to analyse these data using a variety of indirect ordination graphical display techniques.
## Loading required package: nlme
## This is mgcv 1.8-40. For overview type 'help("mgcv-package")'.
Principal Components Analysis (PCA) is a common statistical methods and is implemented in R by two functions in the standard installation (prcomp()
and princomp()
), as well as in numerous guises as part of other add-on packages. In this practical you will use the implimentation provided in the vegan package by Jari Oksanen. This is because you will be using vegan for the direct gradient analysis class and the vegan package provides a rich set of analysis and graphical tools that can be applied to a variety of ordination techniques. As such, you only have to learn a single set of commands to run a wide variety of analyses.
Start R and load the vegan package for use. Read in the two data sets:
library("vegan")
## Loading required package: permute
## Loading required package: lattice
## This is vegan 2.6-2
url_env <- "https://bit.ly/pondsenv"
pondsenv <- read.csv(url(url_env))
url_diat <- "https://bit.ly/pondsdiat"
pondsdiat <- read.csv(url(url_diat))
PCA is fitted using the rda()
function, which was designed to implement Redundancy Analysis (RDA). RDA is the constrained form of PCA so running rda()
without any constraints yields PCA. Run a PCA of the Ponds environmental data using rda()
and display the results. The scale
argument to rda()
scales the variables to zero mean and unit standard deviation. This results in a PCA on a correlation rather than covariance matrix. This is appropriate in situations where the variables are measured in different units, as is the case with the hydrochemical data analysed here. It may also be appropriate in situations where species abundance (counts etc.) are being analysed and you wish to focus on explaining variation in all species not just the most abundant ones.
pondspca <- rda(pondsenv, scale = TRUE)
pondspca
## Call: rda(X = pondsenv, scale = TRUE)
##
## Inertia Rank
## Total 15
## Unconstrained 15 15
## Inertia is correlations
##
## Eigenvalues for unconstrained axes:
## PC1 PC2 PC3 PC4 PC5 PC6 PC7 PC8 PC9 PC10 PC11 PC12 PC13
## 4.964 3.094 2.329 1.428 0.772 0.626 0.534 0.398 0.309 0.293 0.128 0.069 0.029
## PC14 PC15
## 0.019 0.010
The ouput shows the results of the PCA, displaying the eigenvalues for each axis. Note that these values differ from the ones report by CANOCO, which scales the total variance (called inertia in vegan) to 1, rda()
does not. To achieve a comparable set of eigenvalues simply divide each eigenvalue by the total inertia. This conveniently, is what the summary()
method for eigenvals()
does:
summary(eigenvals(pondspca))
## Importance of components:
## PC1 PC2 PC3 PC4 PC5 PC6 PC7
## Eigenvalue 4.9639 3.0937 2.3289 1.42778 0.77180 0.62609 0.53372
## Proportion Explained 0.3309 0.2062 0.1553 0.09519 0.05145 0.04174 0.03558
## Cumulative Proportion 0.3309 0.5372 0.6924 0.78763 0.83908 0.88082 0.91640
## PC8 PC9 PC10 PC11 PC12 PC13
## Eigenvalue 0.39793 0.3090 0.29284 0.127767 0.069340 0.028661
## Proportion Explained 0.02653 0.0206 0.01952 0.008518 0.004623 0.001911
## Cumulative Proportion 0.94293 0.9635 0.98305 0.991570 0.996193 0.998103
## PC14 PC15
## Eigenvalue 0.018654 0.009794
## Proportion Explained 0.001244 0.000653
## Cumulative Proportion 0.999347 1.000000
What are the values of \(\lambda_1\) and \(\lambda_1\), the eigenvalues for axes one and two?
How much of the total variance in the environmental data is explained by axes one and two individually and cummulatively?
For a longer print out of the results of the PCA and to display the species and site scores a summary()
method is available. This truncates the output to the first six axes but this can be changed using the axes
argument.
# output not show as it is large!
summary(pondspca, scaling = "species")
The use of the scaling
argument is redundant in this case as the default is scaling = "species"
but it is included here to illustrate a common argument that can be specified by the user.
To assess the likely importance of the axes, it is useful to both plot a scree plot and to compare the sizes of the actual PCA axes with the sizes expected under a random (null) model, such as the broken stick distribution. A scree plot of the results of a PCA can be produced using the screeplot()
method. The broken-stick distribution is simple to calculate and the expected values of the pieces of the broken stick are given by
\[E_j=\frac{1}{n}\sum_{x=j}^n \frac{1}{x}\]
where \(E_j\) is the expected value of the \(j^{th}\) piece, and \(n\) is the number of pieces (axes in this case). We can generate the null hypothesis values using the bstick
argument.
screeplot(pondspca, type = "lines", bstick = TRUE)
How many axes does the scree plot suggest are significant?
Compare the variances of the pieces expected under the broken-stick model with the eigenvalues obtained from the results of the PCA of the Ponds environmental data. The broken-stick distribution represents the null model, so significant axes are those that have an eigenvalue that exceeds the expected value for the \(j^{th}\) piece of the broken-stick distribution.
How many axes does the broken-stick distribution suggest are significant?
Two types of biplot may be used to visualise the results of a PCA; a distance biplot and a correlation biplot. In this analysis of the Ponds environmental data, the focus is on the correlations between the environmental variables (“species” in common parlance, even though environmental data are being analysed!), therefore, a correlation biplot is appropriate. The type of biplot produced is determined by the scaling applied to one or both of the species and site scores. For a correlation biplot, scaling "species"
is used, in which the scores for the \(k^{th}\) species (the \(k^{th}\) eigenvector) are scaled to length \(\sqrt{\lambda_k}\) [^The scores for the \(k^{th}\) are then proportional to the square root of the \(k^{th}\) eigenvalue.], whilst the sites are left unscaled. Biplots are easily obtained using the plot()
method of rda()
.
plot(pondspca, scaling = "species")
Species are traditionally represented by biplot arrows, which we can do for rda()
objects using the biplot()
method
biplot(pondspca, scaling = "species")
What are the main chemical gradients represented by axes one and two?
Are there any outliers sites on axes one or two?
Interpretting ordination diagrams can be improved by enhancing the plot with additional information, commonly in the form of response surfaces. Function ordisurf()
generates response surfaces using a generalized additive model (GAM) to predict a response variable for the ordination configuration. The irregular configuration is then interpolated to a regular grid using linear interpolation routines.
plot(pondspca, scaling = "species", display = "sites")
ordisurf(pondspca ~ TP, data = pondsenv, add = TRUE)
##
## Family: gaussian
## Link function: identity
##
## Formula:
## y ~ s(x1, x2, k = 10, bs = "tp", fx = FALSE)
##
## Estimated degrees of freedom:
## 2.3 total = 3.3
##
## REML score: 10.03019
Firstly, the biplot is redisplayed with the species suppressed (display = "sites"
), then a response surface for the variables TP (Total Phosphorus) produced using ordisurf()
. The argument add = TRUE
is used to add the response surface to the existing plot without clearing the graphics device first.
Make response surface plots for selected variables
Which variables respond non-linearly across the ordination?
Species often show unimodal responses to environmental gradients. PCA assumes a linear response model and as such may not be best suited to the analysis of species data exhibiting unimodal responses as PCA is unlikely to fit the data adequately. Correspondence Analysis (CA) is an indirect ordination technique that assumes an idealised unimodal response in species. The response is idealised in that it assumes that species responses are symmetrical, of equal height and width and are equally spaced.
A CA is performed using function cca()
, also in package vegan
. CA can produced unstable results when there are rare species or odd samples in the data set, therefore, rare species tend to be downweighted. Function downweight()
provides this capability, replicating the behaviour of CANOCO.
pondsca <- cca(downweight(pondsdiat))
pondsca
## Call: cca(X = downweight(pondsdiat))
##
## Inertia Rank
## Total 4.996
## Unconstrained 4.996 29
## Inertia is scaled Chi-square
##
## Eigenvalues for unconstrained axes:
## CA1 CA2 CA3 CA4 CA5 CA6 CA7 CA8
## 0.6728 0.5162 0.4199 0.3619 0.3408 0.3254 0.2913 0.2663
## (Showing 8 of 29 unconstrained eigenvalues)
summary(eigenvals(pondsca))
## Importance of components:
## CA1 CA2 CA3 CA4 CA5 CA6 CA7
## Eigenvalue 0.6728 0.5162 0.41994 0.36191 0.34085 0.32535 0.2913
## Proportion Explained 0.1347 0.1033 0.08406 0.07244 0.06823 0.06513 0.0583
## Cumulative Proportion 0.1347 0.2380 0.32206 0.39450 0.46273 0.52785 0.5862
## CA8 CA9 CA10 CA11 CA12 CA13 CA14
## Eigenvalue 0.26632 0.22195 0.21172 0.20888 0.1659 0.14270 0.13265
## Proportion Explained 0.05331 0.04443 0.04238 0.04181 0.0332 0.02856 0.02655
## Cumulative Proportion 0.63946 0.68389 0.72627 0.76808 0.8013 0.82985 0.85640
## CA15 CA16 CA17 CA18 CA19 CA20 CA21
## Eigenvalue 0.11605 0.10483 0.08433 0.07646 0.06781 0.06392 0.05188
## Proportion Explained 0.02323 0.02098 0.01688 0.01530 0.01357 0.01280 0.01038
## Cumulative Proportion 0.87963 0.90061 0.91749 0.93280 0.94637 0.95917 0.96955
## CA22 CA23 CA24 CA25 CA26 CA27
## Eigenvalue 0.035947 0.034768 0.023473 0.019797 0.014481 0.01244
## Proportion Explained 0.007195 0.006959 0.004699 0.003963 0.002899 0.00249
## Cumulative Proportion 0.976747 0.983707 0.988406 0.992368 0.995267 0.99776
## CA28 CA29
## Eigenvalue 0.009305 0.0019015
## Proportion Explained 0.001863 0.0003806
## Cumulative Proportion 0.999619 1.0000000
What are the values of \(\lambda_{1-4}\), the eigenvalues for axes one to four?
How much of the total variance (inertia) in the species data is explained by axes 1 and 2, individually and combined?
As with PCA, a scree plot is a useful way of visualising the important components for further analysis. Comprison of the eignevalues of the CA with those expected under the broken-stick model is also useful.
screeplot(pondsca, bstick = TRUE, type = "lines")
How many axes are iportant when compared with the null model?
Biplots, are generated by plotting the site and species scores obtained from CA of the species matrix. As with PCA, these scores can be scaled to focus the resulting diagram on relationships among sites (scaling = "sites"
) or species (scaling = "species
) or some compromise of the two (scaling = "symmetric"
). So-called biplot and Hill’s scaling are the two types of scaling that can be applied to three scalings mentioned above. These two types dictate how information on the species data is gleaned from the resulting biplot. With biplot scaling, the biplot rule applies and is most suited for short gradients. Hill’s scaling equalises the average niche breadth for all axes and allows, for long gradients, interpretation via the distance rule. The plot()
method for cca()
produces a biplot of the results of a CA
plot(pondsca)
Are there outliers in the plot (species or samples)?
Is there an arch apparent in the biplot?
It is useful to enhance these plots with further information to aid interpretation. One useful addition is to overlay measures of diversity on to the biplot. Diversity indices can be calculated using the diversity()
function of vegan. One particular measure of diveristy is Hill’s \(\mathrm{N_2}\), which gives the effective number of occurrences of a species across all sites or the effective number of species in an individual plot. Hill’s \(\mathrm{N_2}\) is equivalent to the Inverse Simpson diversity measure we tell diversity()
to use. The third argument is what Jari Oksanen refers to as the MARGIN; rows (1) or columns (2). For the effective number of occurences of a particular species then the calculation should be over the columns, but rows should be used to calculate the effective numbers of species per sample.
spp.n2 <- renyi(t(pondsdiat), scales = 2, hill = TRUE)
site.n2 <- renyi(pondsdiat, scales = 2, hill = TRUE)
Having calculated the relevant diversity information, this can be used to augment the biplot. A simple way to use these data is to scale the plotting symbol according to the Hill’s \(\mathrm{N_2}\) value, using the cex
graphical parameter. The identify
function is used to label the plotting symbols after plotting.
ca.plot <- plot(pondsca, type = "n")
points(pondsca, display = "species", cex = 0.3 * spp.n2)
identify(ca.plot, what = "species", col = "red", ps = 10)
plot(pondsca, type = "n")
points(pondsca, display = "sites", cex = 0.5 * site.n2)
identify(ca.plot, what = "sites", col = "red", ps = 10)
Are outlying species common or rare?
Are outlying samples dominated by a few, rare, species?
oldpar <- par(mfrow = c(1,2), pty = "s")
plot(pondsca, scaling = "species",
main = "Inter-species distance & Biplot scaling")
plot(pondsca, scaling = "sites", hill = TRUE,
main = "Inter-sample distance & Hills scaling")
par(oldpar)
What effect does the choice of scaling have on the ordination plots? Use the code below to display the biplots with two different scalings.
Non-metric multidimensional scaling can be performed using metaMDS()
. Function metaMDS()
in vegan implements helper functions that start the NMDS iterative algorithm at \(k\) randomly selected start points and choose the best model fit (i.e. that reduces the stress the most). metaMDS()
only requires a matrix of data as the function internally calculates the specified dissimilarity matrix for you. We will fit a NMDS model to the Ponds environmental data using Euclidean distances.
library(MASS)
set.seed(123456)
euclid.dis <- vegdist(pondsenv, "euclidean")
nmds.env <- metaMDS(pondsenv, distance = "euclidean", trymax = 50)
## 'comm' has negative data: 'autotransform', 'noshare' and 'wascores' set to FALSE
## Run 0 stress 0.173918
## Run 1 stress 0.2151342
## Run 2 stress 0.1739174
## ... New best solution
## ... Procrustes: rmse 0.0003115365 max resid 0.001443426
## ... Similar to previous best
## Run 3 stress 0.1739187
## ... Procrustes: rmse 0.0006034781 max resid 0.0028008
## ... Similar to previous best
## Run 4 stress 0.1785826
## Run 5 stress 0.2199961
## Run 6 stress 0.188833
## Run 7 stress 0.1888448
## Run 8 stress 0.1888293
## Run 9 stress 0.1785827
## Run 10 stress 0.1739184
## ... Procrustes: rmse 0.001983074 max resid 0.009166062
## ... Similar to previous best
## Run 11 stress 0.1785824
## Run 12 stress 0.1888302
## Run 13 stress 0.2148062
## Run 14 stress 0.1739178
## ... Procrustes: rmse 0.0002252463 max resid 0.001042586
## ... Similar to previous best
## Run 15 stress 0.2400498
## Run 16 stress 0.1739173
## ... New best solution
## ... Procrustes: rmse 8.483775e-05 max resid 0.000388316
## ... Similar to previous best
## Run 17 stress 0.173918
## ... Procrustes: rmse 0.001691557 max resid 0.007818737
## ... Similar to previous best
## Run 18 stress 0.1739181
## ... Procrustes: rmse 0.001738825 max resid 0.008037138
## ... Similar to previous best
## Run 19 stress 0.2147383
## Run 20 stress 0.1739186
## ... Procrustes: rmse 0.001952894 max resid 0.009027558
## ... Similar to previous best
## *** Solution reached
nmds.env
##
## Call:
## metaMDS(comm = pondsenv, distance = "euclidean", trymax = 50)
##
## global Multidimensional Scaling using monoMDS
##
## Data: pondsenv
## Distance: euclidean
##
## Dimensions: 2
## Stress: 0.1739173
## Stress type 1, weak ties
## Two convergent solutions found after 20 tries
## Scaling: centring, PC rotation
## Species: scores missing
NMDS maps the observed distances on to the ordination space in a non-linear fashion. How well this mapping is achieved can be visualised using stressplot()
, which draws a Shepard plot and the fit of the NMDS as a stepped line. stressplot()
also displays two correlation statistics for the goodness of the fit. The correlation based on stress is \(R^2=1-S^2\), and the ``fit-based’’ correlation is the correlation between the fitted values, \(\theta(d)\) and the original distances, \(d\), which is the correlation between the stepped line and the points.
stressplot(nmds.env, euclid.dis)
To draw the ordination diagram for the NMDS model, vegan provides a plot()
method:
plot(nmds.env, type = "text")
## species scores not available
Two ordinations can be very similar but this similarity may be masked as a result of the two ordinations having different scalings, orientations and signs. Procrustes rotation is a good way of of comparing ordination configurations. vegan has function procrustes()
, which performs Procrustes rotation.
pondsenv.pro <- procrustes(nmds.env, pondspca, symmetric = TRUE)
summary(pondsenv.pro)
##
## Call:
## procrustes(X = nmds.env, Y = pondspca, symmetric = TRUE)
##
## Number of objects: 30 Number of dimensions: 2
##
## Procrustes sum of squares:
## 0.245266
## Procrustes root mean squared error:
## 0.09041866
## Quantiles of Procrustes errors:
## Min 1Q Median 3Q Max
## 0.007698021 0.027309880 0.050129855 0.078476015 0.358319320
##
## Rotation matrix:
## [,1] [,2]
## [1,] -0.5064789 0.8622524
## [2,] 0.8622524 0.5064789
##
## Translation of averages:
## [,1] [,2]
## [1,] 1.11458e-18 -6.80119e-19
##
## Scaling of target:
## [1] 0.8687543
The plot()
method for procrustes()
can produce two kinds of plots; 1) the ordination digram showing the comparison between the two configurations, and 2) a residuals plot
par(mfrow = c(1,2))
plot(pondsenv.pro, kind = "1")
plot(pondsenv.pro, kind = "2")
par(mfrow = c(1,1))
The protest()
method allows you to test whether the two configurations are significantly similar to one another by means of a permutation test.
set.seed(123456)
pondsenv.prot <- protest(nmds.env, pondspca)
pondsenv.prot
##
## Call:
## protest(X = nmds.env, Y = pondspca)
##
## Procrustes Sum of Squares (m12 squared): 0.2453
## Correlation in a symmetric Procrustes rotation: 0.8688
## Significance: 0.001
##
## Permutation: free
## Number of permutations: 999
This data set is on bird assemblages in a montane forest in the Velka Fatra Mountains (Slovak Republic).
library("readr")
library("tidyr")
library("dplyr")
library("tibble")
We can load the bird data and process it into two data sets
birds_url <- "https://bit.ly/maed-birds"
birds <- read_csv(birds_url, na = "", skip = 1)
## New names:
## Rows: 43 Columns: 51
## ── Column specification
## ───────────────────────────────────── Delimiter: "," chr
## (4): ...1, ...40, Rocks, Expos dbl (46): AegiCaud, AlauArve,
## AnthTriv, CertFami, ColuPalu, CorvCora, CucuCa... lgl (1):
## ...39
## ℹ Use `spec()` to retrieve the full column specification for
## this data. ℹ Specify the column types or set `show_col_types
## = FALSE` to quiet this message.
## • `` -> `...1`
## • `` -> `...39`
## • `` -> `...40`
birds_spp <- birds %>%
mutate(across(AegiCaud:OriOri, ~ replace_na(.x, 0))) %>%
column_to_rownames("...1")
birds_spp <- birds_spp[, 1:36]
birds_env <- birds %>%
select(c("...1", Altit:Expos)) %>%
mutate(Rocks = factor(Rocks), Expos = factor(Expos)) %>%
column_to_rownames("...1")
The aim of the analysis is to find a good ordination of the species data (choose your ordination method, transformation, dissimilarity, etc according to what features of the species data you want to focus on). Then, having found a good ordination, use the environmental data to try to explain the ordination axes or mapping you have created.
The variables are: