# Ordination - Part 2

##    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
## 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
## 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
## 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
## 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.

• 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

• 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. 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