Ordination - Part 2

Link to Part 1

##    Callvulg Empenigr Rhodtome Vaccmyrt Vaccviti Pinusylv Descflex Betupube
## 18     0.55    11.13     0.00     0.00    17.80     0.07     0.00        0
## 15     0.67     0.17     0.00     0.35    12.13     0.12     0.00        0
## 24     0.10     1.55     0.00     0.00    13.47     0.25     0.00        0
## 27     0.00    15.13     2.42     5.92    15.97     0.00     3.70        0
## 23     0.00    12.68     0.00     0.00    23.73     0.03     0.00        0
## 19     0.00     8.92     0.00     2.42    10.28     0.12     0.02        0
##    Vacculig Diphcomp Dicrsp Dicrfusc Dicrpoly Hylosple Pleuschr Polypili
## 18     1.60     2.07   0.00     1.62     0.00      0.0     4.67     0.02
## 15     0.00     0.00   0.33    10.92     0.02      0.0    37.75     0.02
## 24     0.00     0.00  23.43     0.00     1.68      0.0    32.92     0.00
## 27     1.12     0.00   0.00     3.63     0.00      6.7    58.07     0.00
## 23     0.00     0.00   0.00     3.42     0.02      0.0    19.42     0.02
## 19     0.00     0.00   0.00     0.32     0.02      0.0    21.03     0.02
##    Polyjuni Polycomm Pohlnuta Ptilcili Barbhatc Cladarbu Cladrang Cladstel
## 18     0.13     0.00     0.13     0.12     0.00    21.73    21.47     3.50
## 15     0.23     0.00     0.03     0.02     0.00    12.05     8.13     0.18
## 24     0.23     0.00     0.32     0.03     0.00     3.58     5.52     0.07
## 27     0.00     0.13     0.02     0.08     0.08     1.42     7.63     2.55
## 23     2.12     0.00     0.17     1.80     0.02     9.08     9.22     0.05
## 19     1.58     0.18     0.07     0.27     0.02     7.23     4.95    22.08
##    Cladunci Cladcocc Cladcorn Cladgrac Cladfimb Cladcris Cladchlo Cladbotr
## 18     0.30     0.18     0.23     0.25     0.25     0.23     0.00     0.00
## 15     2.65     0.13     0.18     0.23     0.25     1.23     0.00     0.00
## 24     8.93     0.00     0.20     0.48     0.00     0.07     0.10     0.02
## 27     0.15     0.00     0.38     0.12     0.10     0.03     0.00     0.02
## 23     0.73     0.08     1.42     0.50     0.17     1.78     0.05     0.05
## 19     0.25     0.10     0.25     0.18     0.10     0.12     0.05     0.02
##    Cladamau Cladsp Cetreric Cetrisla Flavniva Nepharct Stersp Peltapht Icmaeric
## 18     0.08   0.02     0.02     0.00     0.12     0.02   0.62     0.02        0
## 15     0.00   0.00     0.15     0.03     0.00     0.00   0.85     0.00        0
## 24     0.00   0.00     0.78     0.12     0.00     0.00   0.03     0.00        0
## 27     0.00   0.02     0.00     0.00     0.00     0.00   0.00     0.07        0
## 23     0.00   0.00     0.00     0.00     0.02     0.00   1.58     0.33        0
## 19     0.00   0.00     0.00     0.00     0.02     0.00   0.28     0.00        0
##    Cladcerv Claddefo Cladphyl
## 18        0     0.25        0
## 15        0     1.00        0
## 24        0     0.33        0
## 27        0     0.15        0
## 23        0     1.97        0
## 19        0     0.37        0
##       N    P     K    Ca    Mg    S    Al   Fe    Mn   Zn  Mo Baresoil Humdepth
## 18 19.8 42.1 139.9 519.4  90.0 32.3  39.0 40.9  58.1  4.5 0.3     43.9      2.2
## 15 13.4 39.1 167.3 356.7  70.7 35.2  88.1 39.0  52.4  5.4 0.3     23.6      2.2
## 24 20.2 67.7 207.1 973.3 209.1 58.1 138.0 35.4  32.1 16.8 0.8     21.2      2.0
## 27 20.6 60.8 233.7 834.0 127.2 40.7  15.4  4.4 132.0 10.7 0.2     18.7      2.9
## 23 23.8 54.5 180.6 777.0 125.8 39.5  24.2  3.0  50.1  6.6 0.3     46.0      3.0
## 19 22.8 40.9 171.4 691.8 151.4 40.8 104.8 17.6  43.6  9.1 0.4     40.5      3.8
##     pH
## 18 2.7
## 15 2.8
## 24 3.0
## 27 2.8
## 23 2.7
## 19 2.7

Unconstrained Ordination



Visualizing Dissimilarity Matrices

Similarly to a PCA, distance matrices can be expressed in a smaller number of dimensions through the process of ordination. In fact, PCA is a form of ordination that preserves Euclidean distances between points. However, if we want to capture similarity between samples using other distance measures, we need to use an alternative method of ordination, such as PCoA (Principal Coordinates Analysis) or NMDS (Non-Metric Multidimensional Scaling). Here, we’ll discuss use of the latter.

NMDS

NMDS is an increasingly popular unconstrained ordination method to express the non-euclidean relationships between multiple samples. Unlike a PCA, NMDS does not have a closed-form solution, instead relying on computer optimization to express a matrix in a given number of dimensions.

Advantages

  • No eigenvalues! The distance between in reduced space does not depend on the axis. Distance between points in Y = Distance between points in X
  • Different distance measures can be expressed.
  • Often easier to interpret

Disadvantages

  • Results from computer optimization. Can converge on different solutions, so plots look different between runs.
  • Can poorly represent multivariate relationships. Must evaluate model fit using stress and adjust accordingly.


Let’s return to the vegetation patch data, varespec. To generate a matrix of distances between points, we can use the vegdist function from the vegan library. The resulting matrix contains every single pairwise comparison in the data, producing an \(n*n\) square matrix with \(n\choose{2}\) unique values.

# Bray-curtis dissimilarity between samples
distmat <- vegan::vegdist(veg_comm, method = "bray")

# First 5 rows/columns
as.matrix(distmat)[1:5, 1:5]
##           18        15        24        27        23
## 18 0.0000000 0.5310021 0.6680661 0.5621247 0.3747078
## 15 0.5310021 0.0000000 0.3597783 0.4055610 0.3652097
## 24 0.6680661 0.3597783 0.0000000 0.4934947 0.5020306
## 27 0.5621247 0.4055610 0.4934947 0.0000000 0.4286111
## 23 0.3747078 0.3652097 0.5020306 0.4286111 0.0000000

Row names are given here as both column and row names, showing which communities are compared.

Note that this matrix is symmetric – elements in the upper diagonal and lower diagonal are the same, as distance functions are identical no matter the sample order. Elements of the diagonal are all zero, as a community is identical to itself.

To help illustrate the NMDS, we can assign some grouping variable to color the datapopints by. Below, I make a simple group of values above / below the median pH value in the associated environmental dataset.

colors <- c("red", "blue")
ph_group <- as.factor(if_else(veg_envt$pH < median(veg_envt$pH), "pH < 2.9",  "pH > 2.9"))

NMDS is performed in vegan using the metaMDS function. Notice that the orientation of points depends strongly on the distance metric we use:

plot_nmds <- function(my_nmds, title){
  fig <- ordiplot(my_nmds, type = "none", main = title)
  text(fig, "species", col="lightblue", cex=0.9)
  points(fig, "sites", pch=21, bg = colors[as.numeric(ph_group)])
  with(veg_envt, legend("topleft", levels(ph_group), pch = 21,  pt.bg = colors,
  title = "pH"))
}

nmds.eucl <- metaMDS(veg_comm, distance = "euclidean", binary = FALSE, autotransform = FALSE, trace = FALSE)
plot_nmds(nmds.eucl, "Ordination - Euclidean Distance")

nmds.bray <- metaMDS(veg_comm, distance = "bray", binary = FALSE, autotransform = FALSE, trace = FALSE)
plot_nmds(nmds.bray, "Ordination - Bray-Curtis Dissimilarity")

nmds.jacc <- metaMDS(veg_comm, distance = "jaccard", binary = TRUE, autotransform = FALSE, trace = FALSE)
plot_nmds(nmds.jacc, "Ordination - Jaccard Dissimilarity")

  • Are there different patterns, depending on what distance metric is used?

  • What might patterns that depend on metric tell us about the process we observe?



NMDS Stress

How well does this two-dimensional representation of points reflect our observed dissimilarities? From our distance matrix object, we know how similar two points are from one another.

The stress value of an NMDS provides an evaluation of the accuracy of this two-dimensional representation. Stress is a measure based on rank-order discrepancies in distance – whether a point has a neighbor that is closer in low-dimensional distance than another neighbor which is known to be more similar. High stress values indicate that the low-dimensional representation of points has a number of cases where this poor configuration is present.

The stressplot object shows a visualization of this stress, where observed community distances on the X-axis are compared to ordination (xy-coordinate) distances on the Y-axis. Note that the non-linear fit smoothing line has a higher \(R^2\) than the linear one – rank-order differences in points are emphasized in NMDS.

# Stress value
nmds.bray$stress
## [1] 0.1000211
# Visualization of stress
stressplot(nmds.bray)

The NMDS object

If you want to create your plots, the observations and species loadings can be extracted from the NMDS object. Like a PCA, correlations between points and species indicate how well changes in the abundance/occurrence of different species drive differences between samples.

# NMDS points -- X and Y coordinates of the observed community samples
head(nmds.bray$points)
##          MDS1         MDS2
## 18 -0.2198916 -0.204531669
## 15  0.4909467 -0.274533029
## 24  0.8073173  0.001060818
## 27  0.7138549  0.316768770
## 23  0.3001206  0.042184127
## 19  0.1572015  0.203736904
# NMDS species -- the relationships between species loadings and the ordination figure 
head(nmds.bray$species)
##                 MDS1       MDS2
## Callvulg -0.25359926 -0.5770297
## Empenigr  0.07920079  0.2174706
## Rhodtome  0.71930835  0.7284397
## Vaccmyrt  1.01985163  0.5801960
## Vaccviti  0.12354403  0.1944955
## Pinusylv -0.14397878  0.3547154

Adding environmental variables

In some ecological applications, we may only have data on species assemblages and no environmental context. If we’re interested in seeing whether communities are assembling across linear gradients, indirect gradient analysis may be used to capture these patterns.

However, in many cases, samples are accompanied by measurements of 1 or more variables that capture environmental gradients. To visualize how these environmental variables correlate with ordination, envfit will construct permutation-based vectors. Strong correlations have longer arrows, while weaker correlations have smaller ones.

A significance test can be used to determine how ordination structure explains variation in a given environmental factor, though note that this does not provide a strong test of whether these variables drive community dissimilarity observed.

# Environmental Vectors -- projecting environmental variables on an existing NMDS object
env_vec <- envfit(nmds.bray, veg_envt)
env_vec
## 
## ***VECTORS
## 
##             NMDS1    NMDS2     r2 Pr(>r)   
## N         0.01871 -0.99983 0.1885  0.117   
## P         0.66077  0.75059 0.2898  0.031 * 
## K         0.85357  0.52098 0.1830  0.111   
## Ca        0.83355  0.55244 0.3534  0.017 * 
## Mg        0.84813  0.52978 0.2901  0.022 * 
## S         0.41302  0.91072 0.0884  0.380   
## Al       -0.98935  0.14557 0.4889  0.002 **
## Fe       -0.99350  0.11385 0.4081  0.004 **
## Mn        0.99233  0.12359 0.4293  0.004 **
## Zn        0.92740  0.37407 0.1650  0.130   
## Mo       -0.85990 -0.51046 0.0476  0.601   
## Baresoil  0.91937 -0.39340 0.1877  0.120   
## Humdepth  0.98855  0.15088 0.4415  0.002 **
## pH       -0.86064  0.50922 0.2252  0.060 . 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Permutation: free
## Number of permutations: 999
# Adding to previous plot
plot_nmds(nmds.bray, "Ordination - Bray-Curtis Dissimilarity")
plot(env_vec)


Constrained Ordination


The counterpart to linear regression in ordination is known as constrained ordination. Rather than trying to visualize relationships between species and samples, constrained ordination attempts to explain the maximum variation in species composition.

Constrained ordination combines both multiple linear regression and PCA to estimate how environmental variables can be combined to produce environmental gradients that capture changes in species abundance.

From Legendre and Legendre 1998

Constrained ordination produces both a matrix of explained variation that can be visualized through a PCA, as well as a matrix of unexplained variation that may be visualized as a separate PCA.

Types of Constrained Ordination

  • RDA - Redundancy Analysis; Community dissimilarity based on Euclidean distance.
  • CCA - Canonical Correspondence Analysis; Community dissimilarity based on Chi-square distance
  • db-RDA - distance-based Redundancy Analysis: Community dissimilarity based on a chosen distance metric

Fitting the RDA object

vegan will fit all three versions of ordination using rda (RDA), cca (CCA), and capscale or dbrda (db-RDA). Due to similarities in their underlying methods, all objects created by these functions share similar nomenclature, despite terms varying in other sources.

Calling rda

RDA can be called like any linear model. In this case, we’re estimating our community, veg_comm as a function of all our environmental variables ~ ., data = veg_envt

The summary object of an RDA is quite large, but contains a wealth of useful information. Key points to notice are:

  • Variance (Inertia) partitioning: how much of the observed community disimilarity can we explain using our environmental variables?

  • Importance of components: How are axes in constrained variation (RDA) and unconstrained variation (PC) composed? Like a PCA, large eigenvalues describe important axes of variation.

  • Site and species scores: Similar to interpretation of a PCA.

veg_rda <- vegan::rda(veg_comm ~ ., data = veg_envt)
summary(veg_rda)
## 
## Call:
## rda(formula = veg_comm ~ N + P + K + Ca + Mg + S + Al + Fe +      Mn + Zn + Mo + Baresoil + Humdepth + pH, data = veg_envt) 
## 
## Partitioning of variance:
##               Inertia Proportion
## Total          1825.7     1.0000
## Constrained    1459.9     0.7997
## Unconstrained   365.8     0.2003
## 
## Eigenvalues, and their contribution to the variance 
## 
## Importance of components:
##                           RDA1     RDA2      RDA3     RDA4    RDA5     RDA6
## Eigenvalue            820.1042 399.2847 102.56168 47.63169 26.8382 24.04809
## Proportion Explained    0.4492   0.2187   0.05618  0.02609  0.0147  0.01317
## Cumulative Proportion   0.4492   0.6679   0.72409  0.75019  0.7649  0.77806
##                           RDA7      RDA8     RDA9    RDA10    RDA11    RDA12
## Eigenvalue            19.06438 10.166995 4.428786 2.272036 1.535326 0.925528
## Proportion Explained   0.01044  0.005569 0.002426 0.001245 0.000841 0.000507
## Cumulative Proportion  0.78850  0.794069 0.796495 0.797740 0.798581 0.799088
##                           RDA13     RDA14      PC1      PC2      PC3      PC4
## Eigenvalue            0.7155102 0.3118612 186.1917 88.46419 38.18828 18.40207
## Proportion Explained  0.0003919 0.0001708   0.1020  0.04846  0.02092  0.01008
## Cumulative Proportion 0.7994795 0.7996503   0.9016  0.95009  0.97101  0.98109
##                             PC5      PC6      PC7      PC8       PC9
## Eigenvalue            12.839433 10.55197 5.519405 4.521057 1.0922062
## Proportion Explained   0.007033  0.00578 0.003023 0.002476 0.0005983
## Cumulative Proportion  0.988122  0.99390 0.996925 0.999402 1.0000000
## 
## Accumulated constrained eigenvalues
## Importance of components:
##                           RDA1     RDA2      RDA3     RDA4     RDA5     RDA6
## Eigenvalue            820.1042 399.2847 102.56168 47.63169 26.83822 24.04809
## Proportion Explained    0.5618   0.2735   0.07025  0.03263  0.01838  0.01647
## Cumulative Proportion   0.5618   0.8353   0.90551  0.93814  0.95653  0.97300
##                           RDA7      RDA8     RDA9    RDA10    RDA11    RDA12
## Eigenvalue            19.06438 10.166995 4.428786 2.272036 1.535326 0.925528
## Proportion Explained   0.01306  0.006964 0.003034 0.001556 0.001052 0.000634
## Cumulative Proportion  0.98606  0.993021 0.996054 0.997611 0.998662 0.999296
##                           RDA13     RDA14
## Eigenvalue            0.7155102 0.3118612
## Proportion Explained  0.0004901 0.0002136
## Cumulative Proportion 0.9997864 1.0000000
## 
## Scaling 2 for species and site scores
## * Species are scaled proportional to eigenvalues
## * Sites are unscaled: weighted dispersion equal on all dimensions
## * General scaling constant of scores:  14.31485 
## 
## 
## Species scores
## 
##                RDA1       RDA2       RDA3       RDA4       RDA5       RDA6
## Callvulg -0.1516557  0.3244635  8.016e-02  1.0775112 -0.8752695 -0.0895046
## Empenigr -0.1601177 -0.2952352  6.174e-01 -0.2508970  0.6812008 -0.2528771
## Rhodtome -0.1148931 -0.0305967  4.398e-02 -0.0351583  0.0390864  0.0098011
## Vaccmyrt -0.7334271 -0.5778960 -2.229e-01  0.0353005  0.0797716  0.0960784
## Vaccviti -0.1087123 -0.6493169  1.130e+00 -0.6096061  0.6585696 -0.3202957
## Pinusylv  0.0247893 -0.0111517  1.134e-02  0.0123298  0.0077482  0.0094760
## Descflex -0.0989092 -0.0784375 -6.588e-02  0.0055568  0.0349688 -0.0343673
## Betupube -0.0028545  0.0014106  5.602e-03 -0.0046374  0.0010045  0.0016631
## Vacculig -0.0041388  0.2742036 -1.006e-01 -0.0431564  0.1090859 -0.2417334
## Diphcomp -0.0129952  0.0396768  2.102e-02 -0.0203400  0.0366495 -0.0550255
## Dicrsp   -0.3918593 -0.2592629 -1.813e-01 -0.7473832 -0.4388094  0.3455296
## Dicrfusc -0.8586313  0.0483891  6.352e-01  1.6596210  0.7548590  0.4891008
## Dicrpoly -0.0353438 -0.0004213  5.259e-02 -0.1104382 -0.0419125  0.0453031
## Hylosple -0.3502877 -0.4293435 -3.855e-01  0.0477256 -0.0016439 -0.0383548
## Pleuschr -3.6921567 -4.1569004 -1.985e+00  0.2292435  0.1908869 -0.2639675
## Polypili  0.0004240  0.0047989 -3.770e-03 -0.0029230  0.0038193  0.0090228
## Polyjuni -0.0827076 -0.0801868  1.726e-02 -0.0721091  0.1486044 -0.0071481
## Polycomm -0.0069388 -0.0046443  5.123e-03 -0.0079036  0.0027771 -0.0013614
## Pohlnuta  0.0065052 -0.0085815  1.249e-02 -0.0092286 -0.0024448  0.0048080
## Ptilcili -0.1068778  0.0467979  2.492e-01 -0.2085075  0.0406721  0.0464190
## Barbhatc -0.0312925  0.0178933  6.557e-02 -0.0603799  0.0026625  0.0158739
## Cladarbu -0.2793833  2.4862824 -2.205e-01  0.4551938 -0.0672690 -1.2004520
## Cladrang  0.8685685  3.9115659 -2.124e+00 -0.1565390  0.4962173  0.3887768
## Cladstel  8.7125474 -2.1529119 -5.819e-01  0.2637918  0.0983628 -0.1194561
## Cladunci -0.0936083 -0.0939778  4.766e-01  0.0182430 -0.3403172  0.5414676
## Cladcocc  0.0094128  0.0056428  2.086e-03  0.0081529  0.0053307  0.0019166
## Cladcorn -0.0114892 -0.0149133  2.018e-02 -0.0115155  0.0300531 -0.0113060
## Cladgrac -0.0095786  0.0097235  1.134e-02 -0.0096839 -0.0029627 -0.0013762
## Cladfimb  0.0030921  0.0029475  1.458e-02  0.0079253  0.0072732 -0.0086739
## Cladcris -0.0064259 -0.0138225  5.778e-02  0.0211211  0.0079897 -0.0265414
## Cladchlo  0.0125661 -0.0043404  6.770e-03 -0.0094834 -0.0010920 -0.0015495
## Cladbotr -0.0044290  0.0005884  6.002e-03 -0.0038310 -0.0023778  0.0008171
## Cladamau -0.0009165  0.0022457 -1.283e-05 -0.0010579  0.0016960 -0.0008148
## Cladsp    0.0066533 -0.0017477  3.766e-03  0.0034792 -0.0037235 -0.0022040
## Cetreric  0.0049398  0.0023199  5.995e-03 -0.0114634 -0.0358258  0.0174480
## Cetrisla  0.0112714 -0.0044663  1.531e-02 -0.0130466  0.0007059  0.0060393
## Flavniva  0.0508423  0.1169823 -7.238e-02 -0.0293247 -0.3236819 -0.1681646
## Nepharct -0.0414089 -0.0357142 -2.522e-02 -0.0250485  0.0908416 -0.0046090
## Stersp   -0.0741567  0.3160965 -2.008e-01 -0.0944856  0.0999565  0.2977516
## Peltapht -0.0029055 -0.0037065  1.489e-03 -0.0016508  0.0080051 -0.0008611
## Icmaeric -0.0010746  0.0046117 -1.664e-03  0.0003124  0.0012521  0.0037920
## Cladcerv  0.0007786  0.0003074 -7.031e-04 -0.0001470 -0.0011868 -0.0002338
## Claddefo -0.0316284 -0.0149305  1.018e-01  0.0021458 -0.0081164  0.0167849
## Cladphyl  0.0179500 -0.0080406  4.363e-04 -0.0014913  0.0028995  0.0019705
## 
## 
## Site scores (weighted sums of species scores)
## 
##        RDA1    RDA2    RDA3     RDA4    RDA5    RDA6
## 18 -1.08199  2.8291  2.2645 -2.46209  1.9279 -4.8908
## 15 -2.88751 -1.6383 -1.1533  2.41906  0.1465 -1.4756
## 24 -2.70584 -2.0747 -0.1579 -5.39872 -6.5132  4.2994
## 27 -3.48402 -4.5723 -4.1584 -0.89852  4.8101 -2.7348
## 23 -2.10848 -0.2517  3.2354 -3.10308  3.3148 -2.4252
## 19 -0.08590 -1.7570  0.8717 -1.13143 -1.0773 -2.1302
## 22 -2.75642 -1.3568  3.3864 10.30007  5.6539  5.5739
## 16 -2.24200  0.4336  1.4498  6.86770  2.6828  5.9399
## 28 -4.30799 -6.1953 -6.9244  0.34045 -1.1737 -1.9251
## 13 -0.32251  3.0065 -0.1969  4.95192 -9.0286 -3.6287
## 14 -1.79840  0.8953  5.7370  3.17849 -6.4014  2.4628
## 20 -1.82495  0.4712  2.6010 -2.66544 -1.3235  1.5498
## 25 -2.65076 -1.7642  1.5231 -0.66345 -1.2262  2.0652
## 7  -1.01235  6.1917 -2.3113 -1.45940  1.3658 -6.1935
## 5  -0.72750  6.8180 -5.9800 -2.61797  2.0878  6.3410
## 6  -0.03626  5.3562 -1.5201 -0.43468  0.7498 -8.1603
## 3   4.19001  0.9732 -2.5869  0.02795  0.1859  2.2499
## 4   1.43729  1.8085 -0.2818  0.54972 -5.7765 -1.0191
## 2   5.96685 -1.0309 -1.9164 -0.41176  2.1991  3.3238
## 9   6.51120 -3.7132  1.7126 -0.22824  0.6866 -1.1885
## 12  4.52925 -1.7696  0.4311 -1.06298  2.3825 -0.2696
## 10  6.68912 -2.8610  0.9669  0.20437  0.9849 -0.2235
## 11  1.22025  0.4563 -3.7136 -0.79085  0.7115  1.0898
## 21 -0.51110 -0.2546  6.7215 -5.51111  2.6306  1.3695
## 
## 
## Site constraints (linear combinations of constraining variables)
## 
##       RDA1     RDA2     RDA3     RDA4    RDA5     RDA6
## 18 -2.9026  2.89996  2.96998 -1.69026  3.5479 -4.23983
## 15  0.2116 -1.77427  0.73050 -0.09712  0.1127 -3.55952
## 24 -2.1373 -1.16378 -1.29204 -6.06643 -6.0357  2.79874
## 27 -4.1896 -3.90854 -3.86273 -0.62616  1.3639 -2.37415
## 23 -1.9720 -1.25781  3.35488 -0.62268  2.5771 -1.41799
## 19 -0.7834 -2.87861  2.54639 -2.17144 -0.8709 -1.93914
## 22 -2.5321 -0.01799  1.47923  6.59952  2.7141  1.88423
## 16 -1.5311  2.48518  1.16595  3.14926  0.9076  1.99661
## 28 -4.2006 -6.14400 -6.01408  1.72598 -0.5903  0.33411
## 13 -0.8474  2.09226 -0.04308  7.14883 -7.3666 -1.56458
## 14  0.1846 -0.59789  3.61782  1.78651 -0.7940  3.37119
## 20 -1.3270 -0.36980  5.08673 -2.38131 -2.7145  3.49239
## 25 -1.8693 -1.84431 -0.88281 -0.91619  3.7530  0.07439
## 7   1.2482  6.71992 -2.33262 -1.91829  1.2882 -5.29860
## 5  -1.5416  5.81395 -4.25019 -1.34610  1.7026  7.28800
## 6   0.3897  2.15118 -0.91250  1.18848  1.6909 -1.76121
## 3   3.0216  1.79232 -3.75741 -1.37242  0.4210 -1.54716
## 4   0.4539  2.27213 -1.46062 -0.47444 -6.7426 -3.36670
## 2   5.8081 -0.45774 -1.77962  0.39739  0.6294  3.96489
## 9   6.7013 -2.67049  3.02412  0.47160 -0.6095 -0.55045
## 12  1.8277 -0.64637  1.45845  3.52226  1.5646 -0.28305
## 10  5.7568 -3.63420 -1.27682 -1.51524  1.1125 -0.78291
## 11  2.2179 -0.16631 -2.11321 -0.53602  2.0325  2.27426
## 21 -1.9876  1.30519  4.54371 -4.25573  0.3060  1.20646
## 
## 
## Biplot scores for constraining variables
## 
##              RDA1     RDA2     RDA3      RDA4      RDA5      RDA6
## N        -0.37451  0.40005 -0.32028  0.120999 -0.021093  0.248110
## P        -0.07887 -0.67736 -0.18434 -0.212893 -0.346547 -0.006476
## K        -0.20146 -0.54606 -0.07107  0.248513 -0.530396 -0.100114
## Ca       -0.23542 -0.59524  0.04435 -0.126410 -0.106923  0.179456
## Mg       -0.25801 -0.50325  0.05625 -0.250266 -0.347164 -0.037531
## S         0.11042 -0.46983  0.07747  0.053454 -0.725932 -0.094191
## Al        0.62279  0.32064 -0.21908  0.055570 -0.538971 -0.219545
## Fe        0.53855  0.36246 -0.36254 -0.005091 -0.204949 -0.176280
## Mn       -0.53181 -0.55604 -0.10457  0.116995  0.076883 -0.063081
## Zn       -0.18926 -0.46093 -0.08390 -0.214118 -0.557270  0.368514
## Mo        0.08255  0.15128  0.12569 -0.054562 -0.716923 -0.163114
## Baresoil -0.58452  0.02519  0.70509 -0.179725 -0.121501 -0.067668
## Humdepth -0.52825 -0.49892  0.39395  0.105904 -0.005877 -0.228625
## pH        0.48529  0.20916 -0.40764 -0.282358  0.098835  0.168733

Variance Decomposition of RDA

Significance testing can be performed using a permutation-based procedure. Here, random re-samples of the data (without replacement) are used to determine whether the total variation explained differs from random expectation. Like a univariate linear regression, adding additional variables can only increase R^2 – when the number of environmental variables matches or exceeds the number of species, we will achieve a perfect fit.

anova(veg_rda)
## Permutation test for rda under reduced model
## Permutation: free
## Number of permutations: 999
## 
## Model: rda(formula = veg_comm ~ N + P + K + Ca + Mg + S + Al + Fe + Mn + Zn + Mo + Baresoil + Humdepth + pH, data = veg_envt)
##          Df Variance      F Pr(>F)   
## Model    14  1459.89 2.5658  0.006 **
## Residual  9   365.77                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

\(R^2\) = 1459.89 / (1459.89 + 365.77) = 0.80

Constrained RDA Axes

Like a PCA, the eigenvalues of RDA axes represent the fraction of variance that is explained by a given axis. In this case, we see the presence of two large RDA axes, indicating the presence of two strong environmental gradients that were observed to correlate with species abundance.

veg_rda$CCA$eig
##        RDA1        RDA2        RDA3        RDA4        RDA5        RDA6 
## 820.1042107 399.2847431 102.5616781  47.6316940  26.8382218  24.0480875 
##        RDA7        RDA8        RDA9       RDA10       RDA11       RDA12 
##  19.0643756  10.1669954   4.4287860   2.2720357   1.5353257   0.9255277 
##       RDA13       RDA14 
##   0.7155102   0.3118612
screeplot(veg_rda)

Unconstrained PCA Axes (Residuals)

Eigenvalues of the unconstrained variation can also be interesting – large variance across single dimensions reflect strong gradients in the residuals. These can reflect unobserved ecological processes that we may have missed in our sampling.

veg_rda$CA$eig
##        PC1        PC2        PC3        PC4        PC5        PC6        PC7 
## 186.191732  88.464194  38.188277  18.402075  12.839433  10.551972   5.519405 
##        PC8        PC9 
##   4.521057   1.092206

Fitted values

Fitted (estimated) species abundances across sites can be called using predict()

head(predict(veg_rda)[,1:5])
##      Callvulg  Empenigr    Rhodtome   Vaccmyrt Vaccviti
## 18 -0.3011956 12.199953  0.18612486 -2.0077465 19.42163
## 15  2.2815377  9.223692  0.95719513  3.8279863 14.83554
## 24  0.2814485  1.352557 -0.03336222 -0.5692336 13.85235
## 27  0.6930281  8.832771  1.26117287  5.4449661 12.56049
## 23 -1.4321608 11.285113 -0.08854819  3.3501278 18.45637
## 19  0.3111303 11.680112  1.10541064  5.2830738 15.50446

Visualizing the ordination

ordiplot(veg_rda)

vif.cca(veg_rda)
##         N         P         K        Ca        Mg         S        Al        Fe 
##  1.884049  6.184742 11.746873  9.825650  9.595203 18.460717 20.331905  8.849137 
##        Mn        Zn        Mo  Baresoil  Humdepth        pH 
##  5.168278  7.755110  4.575108  2.213443  5.638807  6.910118
step_mod <- ordistep(veg_rda, trace = FALSE)
summary(step_mod)
## 
## Call:
## rda(formula = veg_comm ~ N + K + S + Al + Mo + Baresoil + Humdepth,      data = veg_envt) 
## 
## Partitioning of variance:
##               Inertia Proportion
## Total          1825.7     1.0000
## Constrained    1208.2     0.6618
## Unconstrained   617.4     0.3382
## 
## Eigenvalues, and their contribution to the variance 
## 
## Importance of components:
##                           RDA1     RDA2    RDA3     RDA4      RDA5     RDA6
## Eigenvalue            736.8207 330.2030 90.3746 23.45293 16.177102 9.704509
## Proportion Explained    0.4036   0.1809  0.0495  0.01285  0.008861 0.005316
## Cumulative Proportion   0.4036   0.5845  0.6340  0.64681  0.655669 0.660985
##                            RDA7      PC1       PC2     PC3      PC4     PC5
## Eigenvalue            1.5038193 297.3960 119.00842 95.1253 28.19683 21.3625
## Proportion Explained  0.0008237   0.1629   0.06519  0.0521  0.01544  0.0117
## Cumulative Proportion 0.6618084   0.8247   0.88989  0.9420  0.95744  0.9691
##                            PC6       PC7       PC8     PC9     PC10     PC11
## Eigenvalue            18.93542 10.696216 10.577381 6.66300 5.047257 2.571547
## Proportion Explained   0.01037  0.005859  0.005794 0.00365 0.002765 0.001409
## Cumulative Proportion  0.97952  0.985374  0.991168 0.99482 0.997582 0.998991
##                            PC12      PC13      PC14      PC15      PC16
## Eigenvalue            0.9045629 0.5094595 0.2095761 0.1803777 3.885e-02
## Proportion Explained  0.0004955 0.0002791 0.0001148 0.0000988 2.128e-05
## Cumulative Proportion 0.9994861 0.9997651 0.9998799 0.9999787 1.000e+00
## 
## Accumulated constrained eigenvalues
## Importance of components:
##                           RDA1     RDA2    RDA3     RDA4     RDA5     RDA6
## Eigenvalue            736.8207 330.2030 90.3746 23.45293 16.17710 9.704509
## Proportion Explained    0.6098   0.2733  0.0748  0.01941  0.01339 0.008032
## Cumulative Proportion   0.6098   0.8831  0.9579  0.97733  0.99072 0.998755
##                           RDA7
## Eigenvalue            1.503819
## Proportion Explained  0.001245
## Cumulative Proportion 1.000000
## 
## Scaling 2 for species and site scores
## * Species are scaled proportional to eigenvalues
## * Sites are unscaled: weighted dispersion equal on all dimensions
## * General scaling constant of scores:  14.31485 
## 
## 
## Species scores
## 
##                RDA1       RDA2       RDA3       RDA4       RDA5       RDA6
## Callvulg -0.1862482 -0.0475334  0.3165777  0.7550155 -0.4448616 -5.968e-01
## Empenigr -0.1684226  0.3320329 -0.7136688 -0.6710729 -0.3080404 -1.236e-01
## Rhodtome -0.1094232  0.0063802 -0.0524178 -0.0793187 -0.0724582 -8.933e-04
## Vaccmyrt -0.8763580  0.2377869  0.2393532 -0.2257362 -0.3290690 -5.073e-03
## Vaccviti -0.1079377  0.6673698 -1.2697087 -0.6181152  0.1085505 -6.394e-02
## Pinusylv  0.0220889  0.0167179 -0.0023307 -0.0123700 -0.0049985  1.006e-02
## Descflex -0.0812537  0.0888870  0.0550071 -0.0305238 -0.0324402  1.420e-02
## Betupube -0.0030636 -0.0038765 -0.0067040 -0.0041382 -0.0035766 -6.103e-04
## Vacculig  0.0756432 -0.2085796  0.0414606 -0.0480858 -0.0998604 -3.691e-02
## Diphcomp  0.0015927 -0.0273682 -0.0340156 -0.0142884  0.0056445 -6.349e-03
## Dicrsp   -0.1538544  0.2907850 -0.0933622  0.3260830  0.9480164  2.755e-03
## Dicrfusc -1.1348530  0.1787253  0.2380642 -0.5398900 -0.0434701  2.844e-01
## Dicrpoly -0.0236876 -0.0279748 -0.0872600 -0.0116907  0.0286463 -1.970e-02
## Hylosple -0.4012515  0.2944869  0.3656396  0.0299439 -0.0456992 -1.146e-02
## Pleuschr -3.7896912  3.6026184  1.8609958 -0.1299551 -0.0027511  1.059e-01
## Polypili  0.0013979 -0.0060120  0.0040940 -0.0048090  0.0011612 -4.942e-03
## Polyjuni -0.0795046  0.0712230 -0.0411486 -0.1181828 -0.0172898  8.062e-03
## Polycomm -0.0058036  0.0030963 -0.0085103 -0.0066823 -0.0083713  1.163e-03
## Pohlnuta  0.0056159  0.0083369 -0.0136014 -0.0001664  0.0105031 -4.993e-03
## Ptilcili -0.1138748 -0.1465247 -0.3034056 -0.1599093 -0.1621224 -5.518e-02
## Barbhatc -0.0326116 -0.0475606 -0.0825266 -0.0436306 -0.0451343 -1.090e-02
## Cladarbu  0.1816808 -1.9637553  0.2022528  0.3311706 -0.4169613  6.925e-01
## Cladrang  1.5752765 -3.5618059  1.8933027 -0.5501466  0.2207037 -1.995e-01
## Cladstel  7.9639020  2.5526574  0.5451193 -0.0629022 -0.0698980  9.835e-02
## Cladunci -0.1754203  0.0542227 -0.3437822  0.2493171  0.3788032  2.015e-01
## Cladcocc  0.0079717 -0.0018849  0.0021745 -0.0050123  0.0033119 -1.454e-03
## Cladcorn -0.0118636  0.0168972 -0.0229853 -0.0174099  0.0036693  8.727e-05
## Cladgrac -0.0050918 -0.0078630 -0.0134952  0.0028123  0.0104481  3.178e-03
## Cladfimb -0.0008120 -0.0018688 -0.0109004 -0.0071353 -0.0009692  6.451e-03
## Cladcris -0.0158926  0.0200059 -0.0486115 -0.0063871 -0.0157319 -6.121e-03
## Cladchlo  0.0117633  0.0042669 -0.0103874 -0.0031752  0.0051197 -3.549e-03
## Cladbotr -0.0045652 -0.0027726 -0.0070362 -0.0012758 -0.0042674 -3.475e-03
## Cladamau -0.0003081 -0.0020453 -0.0004938 -0.0012029  0.0010593 -3.696e-04
## Cladsp    0.0051712  0.0033558 -0.0020545  0.0051368  0.0026549  3.774e-03
## Cetreric  0.0066064 -0.0000303 -0.0084706  0.0286618  0.0403751  1.215e-02
## Cetrisla  0.0089261 -0.0001645 -0.0179869 -0.0106784 -0.0033629 -2.336e-03
## Flavniva  0.0556338 -0.1331864  0.0251945  0.4388518 -0.0954164  2.022e-01
## Nepharct -0.0397024  0.0262887  0.0172301 -0.0655570 -0.0078319  1.735e-02
## Stersp   -0.0406558 -0.3204447  0.1787551 -0.1871297  0.2678552 -6.813e-02
## Peltapht -0.0021253  0.0036706 -0.0013506 -0.0042645 -0.0056494 -6.912e-03
## Icmaeric -0.0010427 -0.0045721  0.0020552 -0.0024263  0.0026999 -5.150e-04
## Cladcerv  0.0008093 -0.0003837  0.0004583  0.0019097 -0.0007323  6.752e-04
## Claddefo -0.0525934  0.0030753 -0.0881871  0.0059531  0.0147706 -1.949e-02
## Cladphyl  0.0159442  0.0088529 -0.0001363 -0.0031697  0.0025918 -5.650e-03
## 
## 
## Site scores (weighted sums of species scores)
## 
##       RDA1      RDA2    RDA3    RDA4    RDA5    RDA6
## 18 -0.9194 -3.223737 -2.9503 -1.4559 -2.1692  2.0506
## 15 -3.1401  1.420806  1.8811  1.0533  0.5068  5.6471
## 24 -2.9030  1.821316 -0.6586  5.7286 15.7532 -1.2467
## 27 -3.9000  4.431043  4.1755 -5.4024 -1.4325 -2.7737
## 23 -2.2318 -0.048178 -3.8762 -3.6501 -0.3074 -2.4240
## 19 -0.2216  1.872814 -1.2319  1.1202 -1.9365 -0.1708
## 22 -3.1509  1.270710 -0.6588 -4.2863 -2.6634  5.5540
## 16 -2.4265 -0.710624  0.4608 -2.0013 -1.0429  0.7570
## 28 -4.8853  5.943053  7.3302  0.6026 -2.0322 -0.2788
## 13 -0.1400 -3.160911  1.0299  9.9851 -7.4205 -7.3800
## 14 -1.9449 -1.175762 -4.9440  7.6610  0.5919  9.5219
## 20 -1.9015 -0.872900 -3.0293  0.9084  2.1637 -3.0254
## 25 -2.9061  1.578666 -1.3969  1.5954  6.7358  1.5307
## 7  -0.5429 -6.881358  1.8372 -0.1665 -2.2558  6.9344
## 5  -0.1860 -7.731691  5.7860 -3.5899  5.8774 -5.9809
## 6   0.4105 -5.753799  1.1081  0.9632 -3.2591 12.0074
## 3   4.4830 -0.479605  2.4083 -0.6200  0.3391 -1.7118
## 4   1.6219 -1.787923  0.2096  6.8767 -2.0298  2.5877
## 2   6.1946  1.979733  1.5543 -3.3770  0.2291 -5.2551
## 9   6.5362  5.078624 -2.2586 -0.8314 -0.9887 -0.1602
## 12  4.6256  2.626679 -0.9976 -3.0278 -0.8772 -2.5236
## 10  6.7830  4.162544 -1.4308 -1.3926 -2.3309 -1.4177
## 11  1.3763 -0.350634  3.5617 -1.2491  2.1202 -2.7906
## 21 -0.6314 -0.008867 -7.9099 -5.4440 -3.5710 -9.4515
## 
## 
## Site constraints (linear combinations of constraining variables)
## 
##        RDA1     RDA2    RDA3    RDA4     RDA5    RDA6
## 18 -1.61602 -1.55783 -3.9629 -1.5346  0.99954 -1.0479
## 15  0.19828  1.93338 -1.5782 -1.4833 -1.01753  1.8778
## 24  0.07849  1.82697 -1.0651  4.4723  8.38019 -0.3531
## 27 -2.92772  4.77916  2.6841 -1.1562 -1.33696  0.6669
## 23 -2.38118  1.12766 -3.3177 -0.8641 -0.10998 -1.8853
## 19 -0.07531  3.73237 -3.9184 -0.5513 -4.37213  0.9841
## 22 -3.07228  1.46440  2.0942 -1.2167  0.32560  0.7450
## 16 -2.47906 -2.65543  0.4637 -1.2441 -0.07319  0.1646
## 28 -5.88878  3.16587  6.3220  0.9311 -0.35071 -0.1776
## 13 -0.77497  0.09239  2.3235  5.8269 -3.22383 -5.8622
## 14 -1.10384  0.26872 -2.1652  0.2306  0.30239  1.5966
## 20 -2.85779 -2.42099 -4.5201  3.4472  2.29158 -3.7494
## 25 -1.87885  1.45117  0.6496 -2.7339 -0.19330  0.7504
## 7   3.04490 -5.15424  0.4316 -0.6573 -2.39330 -1.0923
## 5  -1.09410 -6.01759  4.1905 -3.9499  6.24124 -1.1622
## 6   0.87819 -1.90785  1.3027 -0.8287  0.46189  7.1840
## 3   3.29905 -3.79889  3.0045 -0.1746 -1.67110  1.4324
## 4   0.58022 -2.67863  0.5481  9.0965 -2.33801  4.0349
## 2   5.81369  0.87612  1.7836 -0.2944 -0.97659 -2.6130
## 9   5.42077  3.40457 -2.3932  1.2141  3.32990  3.2632
## 12  1.08249  1.68127  0.1228 -1.5954 -0.06757  4.0763
## 10  4.93052  3.90515  0.2444 -2.2494  2.61301 -5.1925
## 11  2.93829 -0.10712  2.4700 -1.6089 -3.87899 -3.0206
## 21 -2.11500 -3.41063 -5.7145 -3.0759 -2.94214 -0.6198
## 
## 
## Biplot scores for constraining variables
## 
##              RDA1    RDA2     RDA3      RDA4       RDA5     RDA6
## N        -0.36178 -0.5000  0.38548 -0.001772 -0.0002297 -0.27184
## K        -0.25110  0.5702  0.13264  0.555365 -0.1749597 -0.46228
## S         0.08164  0.5351 -0.08358  0.756964  0.0924765 -0.34193
## Al        0.67911 -0.2521  0.17347  0.608266 -0.2691127 -0.05236
## Mo        0.09721 -0.1492 -0.15660  0.833802  0.1830340  0.44351
## Baresoil -0.61971 -0.1192 -0.71538  0.099721 -0.1607119 -0.22539
## Humdepth -0.59615  0.4673 -0.34561  0.045876 -0.2712671  0.18862

\(R^2\) = 1050.72 / (1050.72 + 774.94) = 0.58

plot(step_mod)

screeplot(step_mod)

var_comp <- data.frame(inertcomp(veg_rda))
var_comp$frac <- var_comp$CCA / (var_comp$CCA + var_comp$CA)
var_comp$MeanCover <- colMeans(veg_comm)
var_comp <- var_comp[order(var_comp$MeanCover, decreasing = TRUE),]
var_comp
##                   CCA           CA      frac    MeanCover
## Cladstel 7.215844e+02 1.378926e+02 0.8395622 20.279583333
## Cladrang 1.875570e+02 3.476091e+01 0.8436433 16.196250000
## Pleuschr 3.123968e+02 4.455506e+01 0.8751791 15.748750000
## Vaccviti 2.706177e+01 1.175655e+01 0.6971392 11.459583333
## Cladarbu 7.328472e+01 3.913993e+01 0.6518563 10.627083333
## Empenigr 1.217939e+01 1.059943e+01 0.5346805  6.332916667
## Dicrfusc 4.362034e+01 3.927876e+01 0.5261859  4.730000000
## Cladunci 8.083968e+00 1.659875e+01 0.3275153  2.345000000
## Vaccmyrt 1.417196e+01 9.046066e+00 0.6103860  2.112916667
## Callvulg 2.314414e+01 1.733804e+00 0.9303076  1.877916667
## Dicrsp   1.819381e+01 1.133965e+01 0.6160406  1.687500000
## Hylosple 4.615779e+00 1.147427e+00 0.8009048  0.751666667
## Stersp   3.015942e+00 1.308693e+00 0.6973865  0.730000000
## Vacculig 2.005170e+00 8.974036e-01 0.6908248  0.634166667
## Ptilcili 2.469095e+00 1.696077e+00 0.5927955  0.583750000
## Polyjuni 8.213258e-01 1.295774e+00 0.3879486  0.577083333
## Flavniva 3.546769e+00 6.026553e-01 0.8547617  0.493750000
## Claddefo 1.489069e-01 1.106306e-01 0.5737394  0.426250000
## Rhodtome 3.964286e-01 5.363234e-01 0.4250096  0.349583333
## Cladcris 7.568935e-02 1.065916e-01 0.4152345  0.311250000
## Cladcorn 3.135534e-02 4.777437e-02 0.3962525  0.259166667
## Dicrpoly 2.840107e-01 1.760784e-01 0.6172950  0.252500000
## Descflex 3.132998e-01 2.700539e-01 0.5370666  0.233333333
## Nepharct 4.157315e-01 5.680678e-01 0.4225775  0.219166667
## Cladgrac 8.711633e-03 3.826773e-03 0.6947959  0.214166667
## Pinusylv 4.399701e-02 3.110136e-02 0.5858584  0.171250000
## Cladfimb 6.209211e-03 2.690500e-04 0.9584688  0.165000000
## Cetreric 3.004450e-02 1.367724e-02 0.6871753  0.150000000
## Diphcomp 1.044172e-01 8.028714e-02 0.5653208  0.135000000
## Barbhatc 2.229363e-01 1.504592e-01 0.5970514  0.132916667
## Cladcocc 3.778233e-03 4.107093e-03 0.4791473  0.116250000
## Pohlnuta 8.628627e-03 1.683692e-03 0.8367301  0.109166667
## Cetrisla 1.574617e-02 8.297124e-03 0.6549091  0.084583333
## Cladchlo 4.404284e-03 2.392818e-03 0.6479650  0.048333333
## Cladphyl 6.253658e-03 8.477908e-04 0.8806172  0.033333333
## Peltapht 3.917454e-03 2.949212e-03 0.5705031  0.031666667
## Polycomm 4.064233e-03 6.355862e-04 0.8647637  0.029583333
## Polypili 2.552390e-03 1.221342e-03 0.6763569  0.025416667
## Cladsp   1.606891e-03 1.085862e-03 0.5967465  0.021666667
## Cladbotr 1.630399e-03 1.095507e-03 0.5981128  0.019583333
## Betupube 1.484168e-03 1.115651e-03 0.5708737  0.012083333
## Icmaeric 5.381944e-04 7.847230e-05 0.8727476  0.009166667
## Cladamau 1.570736e-04 1.639409e-04 0.4893037  0.005833333
## Cladcerv 1.229109e-04 2.419059e-05 0.8355517  0.004166667
comments powered by Disqus