An ordination is a data reduction technique that uses significantly more complicated math than the above examples. The basic premise is that an ordination fits a linearly uncorrelated components through mulitvariate data to maximize the most variation explained. Similar to fitting a line-of-best-fit through a cloud of data points. There are different approaches both among ordinations and within ordinations, but conceptually this is what is occurring. Ordinations can generally be separated into three categories: unconstrained (indirect gradient analysis), constrained (direct gradient analysis), and NMDS. This website does an excellent job breaking down all the different examples of each with a short description of what is occuring.
Technically, NMDS is a type of unconstrained ordination. However, I have intentionally left it in a third category because it does something very different. NMDS calculates rank order and is a non-parametric version of other ordinations. Similar to how Kruskal-Wallis is the rank order equivalent of ANOVA. NMDS is hugely flexible and versatile, but is often criticised for not matching a distribution (again, like how Kruskal-Wallis does not follow a dsitribution) and the output may be suggested as correct, when it is not the optimal representation of data. The ecological community continues to debate this and for this exercise I will focus on principal component analysis (PCA) and redundancy analysis (RDA). However, for those interested in a great tutorial on NMDS check here.
PCA is a indirect gradient analysis. Indirect, in that we are trying to determine a principal component (line of best fit) through the data that explains the most variation but not attribute that component to anything. We use the dune
dataset within vegan as an example.
library(vegan)
library(LearnCommAnalysis)
## Load data
data(dune)
## conduct pca
pca <- rda(dune)
## see output variables of PCA
summary(pca)
##
## Call:
## rda(X = dune)
##
## Partitioning of variance:
## Inertia Proportion
## Total 84.12 1
## Unconstrained 84.12 1
##
## Eigenvalues, and their contribution to the variance
##
## Importance of components:
## PC1 PC2 PC3 PC4 PC5 PC6
## Eigenvalue 24.7953 18.1466 7.62913 7.15277 5.6950 4.33331
## Proportion Explained 0.2947 0.2157 0.09069 0.08503 0.0677 0.05151
## Cumulative Proportion 0.2947 0.5105 0.60115 0.68618 0.7539 0.80539
## PC7 PC8 PC9 PC10 PC11 PC12
## Eigenvalue 3.19936 2.78186 2.4820 1.85377 1.74712 1.31358
## Proportion Explained 0.03803 0.03307 0.0295 0.02204 0.02077 0.01561
## Cumulative Proportion 0.84342 0.87649 0.9060 0.92803 0.94880 0.96441
## PC13 PC14 PC15 PC16 PC17 PC18
## Eigenvalue 0.99051 0.637794 0.550827 0.350584 0.199556 0.148798
## Proportion Explained 0.01177 0.007582 0.006548 0.004167 0.002372 0.001769
## Cumulative Proportion 0.97619 0.983768 0.990316 0.994483 0.996855 0.998624
## PC19
## Eigenvalue 0.115753
## Proportion Explained 0.001376
## 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: 6.322924
##
##
## Species scores
##
## PC1 PC2 PC3 PC4 PC5 PC6
## Achimill -0.603786 0.12392 0.008464 0.159574 0.40871 0.127857
## Agrostol 1.373953 -0.96401 0.166905 0.266466 -0.08765 0.047368
## Airaprae 0.023415 0.25078 -0.194768 -0.326043 0.05574 -0.079619
## Alopgeni 0.531234 -1.42784 -0.505241 -0.042885 -0.44293 0.278566
## Anthodor -0.559138 0.56761 -0.476205 0.015781 0.34408 -0.135783
## Bellpere -0.333560 -0.18881 0.140638 -0.084177 0.12541 0.134771
## Bromhord -0.523468 -0.19656 0.164222 0.005671 0.38612 0.257634
## Chenalbu 0.017494 -0.05462 -0.055349 -0.010582 0.02664 0.016405
## Cirsarve 0.002398 -0.10237 0.063716 -0.048735 -0.03212 -0.036055
## Comapalu 0.168933 0.10522 0.063625 0.052352 0.13056 0.108129
## Eleopalu 1.278257 0.21782 0.469213 0.667986 0.20877 0.189927
## Elymrepe -0.450692 -0.80310 0.340783 -0.243514 0.25145 -0.691715
## Empenigr 0.014054 0.10956 -0.099378 -0.161788 -0.02289 -0.001195
## Hyporadi -0.014612 0.42079 -0.223096 -0.535685 -0.10309 -0.025185
## Juncarti 0.679423 -0.07604 0.243642 0.310903 -0.08877 -0.248736
## Juncbufo 0.065583 -0.45959 -0.548944 -0.018900 -0.10571 -0.087168
## Lolipere -1.455985 -0.39306 1.013109 0.170493 -0.52430 0.114397
## Planlanc -0.913938 0.55455 -0.244341 0.617631 -0.12494 -0.098047
## Poaprat -0.899147 -0.55712 0.542805 -0.042467 -0.27815 -0.026353
## Poatriv -0.756003 -1.56056 -0.480385 0.351099 0.36641 0.044066
## Ranuflam 0.625121 0.06099 0.124760 0.233953 0.13645 0.087328
## Rumeacet -0.582581 0.06663 -0.574256 0.775879 -0.08772 -0.361433
## Sagiproc 0.156823 -0.42388 -0.331722 -0.454322 -0.43262 0.037181
## Salirepe 0.293607 0.45555 -0.023780 -0.196209 -0.20176 -0.097569
## Scorautu -0.453771 0.39268 -0.212281 -0.382424 -0.27635 0.395164
## Trifprat -0.417853 0.16572 -0.234524 0.570030 -0.09646 -0.128045
## Trifrepe -0.581801 -0.02115 -0.167299 0.196535 0.18714 0.928758
## Vicilath -0.106710 0.11571 0.092827 -0.055592 -0.15433 0.129733
## Bracruta 0.148626 0.47690 -0.168758 0.509177 -0.96307 0.029481
## Callcusp 0.538513 0.17963 0.175086 0.238876 0.25531 0.169209
##
##
## Site scores (weighted sums of species scores)
##
## PC1 PC2 PC3 PC4 PC5 PC6
## 1 -0.85678 -0.1724 2.6079 -1.1296 0.45074 -2.49113
## 2 -1.64477 -1.2299 0.8867 -0.9859 2.03463 1.81057
## 3 -0.44010 -2.3827 0.9297 -0.4601 -1.02783 -0.05183
## 4 0.04795 -2.0463 1.2737 -0.9742 -0.64210 -0.72074
## 5 -1.62445 0.2900 -1.5927 1.5398 1.86008 -2.21191
## 6 -1.97427 1.0802 -1.1501 3.3534 -1.52026 0.03127
## 7 -1.79263 0.3220 -0.2200 1.4714 0.01245 -0.42583
## 8 0.88980 -1.0905 0.9250 0.5165 -1.08897 0.94777
## 9 0.00904 -1.6570 -0.4661 -0.2826 -0.10821 -2.16570
## 10 -1.91463 0.4940 0.7058 0.2676 1.36985 2.62386
## 11 -1.04110 1.2081 1.4203 -0.9566 -2.71745 1.11200
## 12 1.01822 -1.4598 -3.2509 -0.3247 -1.75331 1.01550
## 13 0.69939 -2.1837 -2.2128 -0.4231 1.06502 0.65585
## 14 1.49047 0.9772 0.5447 0.2733 2.38875 2.47896
## 15 1.88644 1.1261 0.7271 0.7732 0.22113 -0.31750
## 16 2.84848 -0.2081 0.7041 2.1012 0.29311 -0.08124
## 17 0.04666 1.7279 -0.9135 -1.6663 1.80070 -1.55572
## 18 -0.26936 1.7157 0.1648 -0.5770 -2.10498 0.33880
## 19 0.28094 2.1901 -1.9865 -3.2341 -0.45760 -0.02389
## 20 2.34069 1.2991 0.9029 0.7178 -0.07573 -0.96909
## plot PCA
biplot(pca)
We used the function rda
from vegan because it calculates a PCA when no constraining factor is specified. We use this instead of the base R version of PCA princomp
because it is less versatile and because the output from the rda
function has more options that can be used in other analyses.
Here PCA1 and PCA2 represents our indirect gradients that determine these species compositions. We are unsure what the mechanism behind these components are and that can be better understood in the constrained ordinations later. First, we check our ordination to make sure PCA1 and PCA2 are the best representations of our data. From the summary we can see that 29.5% and 21.6% of the variation is explained in those two axes, whereas the third only explains 9%. This suggests the first two axes are the best representation of our data. Finally, we look at the plot and use the clustering of species together to determine trends. Species that are parallel are correlated, species that are perpendicular are uncorrelated, and species that are opposite are negatively correlated.
If we are working with species abundance data, it is commonly recommended to use Hellinger transformation. Hellinger has a few advantages including not giving rare species high weight and low sensitivity to double zeros. The double zero problem is when two sites mutually are without a certain species and thus are sorted as more similar because of this share dissimilarity. The work of Legendre does a great job of explaining when transformations of ordinations are required.
Let’s transform our data and present it again.
## Load data
dune.hell <- decostand(dune, method="hellinger")
## conduct pca
pca <- rda(dune.hell)
## see output variables of PCA
summary(pca)
##
## Call:
## rda(X = dune.hell)
##
## Partitioning of variance:
## Inertia Proportion
## Total 0.5559 1
## Unconstrained 0.5559 1
##
## Eigenvalues, and their contribution to the variance
##
## Importance of components:
## PC1 PC2 PC3 PC4 PC5 PC6
## Eigenvalue 0.1699 0.1111 0.05634 0.04397 0.04097 0.03019
## Proportion Explained 0.3057 0.1998 0.10135 0.07910 0.07371 0.05431
## Cumulative Proportion 0.3057 0.5054 0.60679 0.68590 0.75960 0.81391
## PC7 PC8 PC9 PC10 PC11 PC12
## Eigenvalue 0.02104 0.01794 0.01379 0.01090 0.009991 0.00912
## Proportion Explained 0.03785 0.03226 0.02481 0.01961 0.017972 0.01641
## Cumulative Proportion 0.85177 0.88403 0.90885 0.92846 0.946428 0.96283
## PC13 PC14 PC15 PC16 PC17
## Eigenvalue 0.006376 0.005166 0.002942 0.002745 0.001651
## Proportion Explained 0.011469 0.009293 0.005292 0.004939 0.002969
## Cumulative Proportion 0.974303 0.983597 0.988888 0.993827 0.996796
## PC18 PC19
## Eigenvalue 0.001138 0.0006432
## Proportion Explained 0.002047 0.0011570
## Cumulative Proportion 0.998843 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: 1.80276
##
##
## Species scores
##
## PC1 PC2 PC3 PC4 PC5 PC6
## Achimill -0.2241355 0.057048 -0.068674 -0.173281 -0.0498964 -0.0006022
## Agrostol 0.4142446 -0.180068 0.020078 -0.020756 0.0029054 0.0193453
## Airaprae -0.0350694 0.154644 0.116317 -0.039999 -0.1384374 -0.0376930
## Alopgeni 0.1386550 -0.318444 0.186419 0.001329 0.0030171 -0.0215368
## Anthodor -0.1962910 0.238579 0.095768 -0.179434 -0.0380470 -0.0471885
## Bellpere -0.1148332 -0.054380 -0.057640 0.064957 0.0010453 0.0896219
## Bromhord -0.1410948 -0.053134 -0.049586 -0.057365 0.0210053 0.0938974
## Chenalbu 0.0121539 -0.025919 0.036102 -0.015353 -0.0038179 0.0173353
## Cirsarve -0.0026599 -0.031975 0.004872 0.021813 -0.0182940 0.0064877
## Comapalu 0.1100812 0.048635 -0.061981 -0.028257 0.0102828 0.1064664
## Eleopalu 0.3796078 0.069099 -0.185862 -0.064110 0.0144632 0.0035959
## Elymrepe -0.1456860 -0.248314 -0.094249 0.010195 -0.1482516 -0.0563573
## Empenigr 0.0004895 0.057213 0.063829 0.026314 -0.0393467 -0.0066958
## Hyporadi -0.0559119 0.197829 0.139068 0.036814 -0.1427093 -0.0326104
## Juncarti 0.2478729 -0.012497 -0.092252 0.005525 0.0265396 -0.1569375
## Juncbufo 0.0190348 -0.121457 0.179346 -0.051491 0.0931930 -0.0211675
## Lolipere -0.3366259 -0.173707 -0.206664 0.154192 0.0092155 -0.0325721
## Planlanc -0.2547353 0.197554 -0.034482 -0.039980 0.1469090 -0.0256698
## Poaprat -0.2808972 -0.142814 -0.082463 0.084918 -0.0502270 -0.0278725
## Poatriv -0.1317430 -0.363370 0.059442 -0.137436 0.0581366 -0.0136469
## Ranuflam 0.2783344 0.024399 -0.081013 -0.052662 0.0009058 0.0274683
## Rumeacet -0.1079283 -0.024094 0.054027 -0.122507 0.2108960 -0.0801798
## Sagiproc 0.0409392 -0.082141 0.246528 0.127747 -0.0148906 -0.0086761
## Salirepe 0.0651846 0.156446 0.007210 0.136628 -0.0166452 -0.0581183
## Scorautu -0.0562239 0.158686 0.094988 0.095244 0.0265794 0.1029567
## Trifprat -0.1001443 0.028383 -0.024352 -0.095019 0.1389927 -0.0417313
## Trifrepe -0.0471194 -0.006999 0.052812 0.045812 0.1395303 0.2537096
## Vicilath -0.0541878 0.052774 -0.023104 0.114929 0.0381791 0.0343033
## Bracruta 0.0852752 0.107955 -0.003783 0.181534 0.2198174 -0.1398609
## Callcusp 0.2045115 0.057120 -0.104072 -0.065216 -0.0146101 0.0579472
##
##
## Site scores (weighted sums of species scores)
##
## PC1 PC2 PC3 PC4 PC5 PC6
## 1 -0.423423 -0.42754 -6.533e-01 0.01583 -0.761705 -0.463469
## 2 -0.356333 -0.36816 -1.991e-01 -0.08370 -0.355682 0.630146
## 3 -0.080900 -0.55215 -5.753e-02 0.27275 -0.217514 0.016724
## 4 -0.041005 -0.49292 7.510e-02 0.33626 -0.282018 0.100013
## 5 -0.458003 0.05220 -1.432e-01 -0.54110 0.351891 -0.117334
## 6 -0.389342 0.20337 -1.022e-01 -0.30549 0.758789 -0.245696
## 7 -0.451813 0.06864 -6.835e-02 -0.41821 0.585536 -0.138732
## 8 0.295113 -0.29375 -6.362e-02 0.19828 0.010903 -0.185034
## 9 0.035973 -0.49264 2.714e-01 0.07945 0.086166 -0.419257
## 10 -0.453094 0.13639 -2.391e-01 -0.12201 0.144355 0.462946
## 11 -0.273214 0.29633 3.931e-05 0.87663 0.126762 0.097096
## 12 0.246622 -0.33173 9.205e-01 -0.03527 0.493526 -0.017653
## 13 0.226908 -0.48390 6.740e-01 -0.28663 -0.071279 0.323642
## 14 0.570818 0.25288 -3.165e-01 -0.32869 -0.085439 1.196310
## 15 0.654415 0.28845 -3.733e-01 0.01035 0.196968 0.002255
## 16 0.720030 -0.10689 -2.727e-01 -0.35992 0.007984 -0.504518
## 17 -0.317465 0.75272 3.395e-01 -0.64284 -0.803236 -0.262492
## 18 -0.201131 0.39818 -2.007e-01 0.89875 0.365682 0.086314
## 19 0.006263 0.73205 8.167e-01 0.33669 -0.503443 -0.085674
## 20 0.689579 0.36848 -4.077e-01 0.09888 -0.048247 -0.475587
## plot PCA
biplot(pca)
We have improved the amount of variation explained and centralized the biplot. Overall, this is an improvement and we know the transformation has address underlying issues with the analysis.
The next thing we want to do is determine whether we have captured the shape of the gradient. A PCA assumes a linear gradient, whereas a correspondance analysis (CA) assumes a unimodal relationship. Over larger gradients it is assumed that the species data will resemble a unimodal relationship. We can check this using a detrended correspondance analysis. A detrended correspondance analysis (DCA) can be used as its own form of analysis but it is another topic of debate like NMDS. However, DCA tends to be argued against more often than in favour. For the purposes of determining gradient shape, it is perfect.
## conduct DCA
dca <- decorana(dune.hell)
dca
##
## Call:
## decorana(veg = dune.hell)
##
## Detrended correspondence analysis with 26 segments.
## Rescaling of axes with 4 iterations.
##
## DCA1 DCA2 DCA3 DCA4
## Eigenvalues 0.5249 0.3042 0.14047 0.08709
## Decorana values 0.5444 0.2830 0.08498 0.04289
## Axis lengths 3.5605 2.9874 1.50175 1.13223
The output of our dca suggests that our axis lengths for the first axes is 3.56. This correspond to the standard deviation of that components and suggest there could be unimodal distribution in the data. The approximate evaluation is as follows: PCA & RDA <3 and 4< CA & CCA
Our number is in between suggesting that both linear and unimodal approximations of the data would be acceptable.
Illustration of the rule how to select whether to use linear ordination methods (like PCA or RDA) or unimodal (CA, DCA or CCA) on the data. Upper diagram shows simulated community structured by a single environmental gradient, with number of species response curves. The diagram below shows the relationship between the length of the gradient sampled in the simulated community on the //x/-axis (in arbitrary units), and length of the first DCA ordination axies (in units of S.D.). The dataset which according to DCA is rather homogeneous (< 3 S.D.) has environmental gradient up to 2000 units long; the longer gradient results into heterogeneous dataset for which linear methods are not suitable. OP of figure & caption
There are two common uses of an ordination. The first is to observe trends in the similarities among sites and species. The second is to utilize the data reduction technique to correlate the gradient among species and sites to other gradients.
We are going to load the varechem
and varespec
datasets from vegan
. We then extract the scores from the ordination and compare them to environmental variables.
data(varechem)
data(varespec)
## transform data
species <- decostand(varespec, method="hellinger")
## create ordination
pca <- rda(species)
## extract the scores from the ordination
score.vals <- scores(pca)$sites
## compare
cor(score.vals, varechem)
## N P K Ca Mg S
## PC1 0.1880740 0.2204900 0.2784250 0.3521916 0.3175540 0.03184488
## PC2 0.3959471 -0.4827799 -0.3678565 -0.4714881 -0.3989006 -0.34226350
## Al Fe Mn Zn Mo Baresoil
## PC1 -0.6791178 -0.6232759 0.5944998 0.2836699 -0.1433314 0.50544429
## PC2 0.1877119 0.2188943 -0.3322152 -0.3267015 0.1709478 0.05213326
## Humdepth pH
## PC1 0.570643 -0.50455541
## PC2 -0.376158 0.04775287
Finally let’s do a contrained ordination, RDA. A constrained ordination is a direct graident analysis because you suggesting what the composition of the components are that describe the data. Instead of the PC1 representing the x-axis, we will have RDA1 that is a combination of the predictor variables that we include. We can take a look at what this means.
## load soil chem variables
data(varechem)
## conduct RDA
rda <- rda(species, varechem)
## output
rda
## Call: rda(X = species, Y = varechem)
##
## Inertia Proportion Rank
## Total 0.3647 1.0000
## Constrained 0.2551 0.6995 14
## Unconstrained 0.1096 0.3005 9
## Inertia is variance
##
## Eigenvalues for constrained axes:
## RDA1 RDA2 RDA3 RDA4 RDA5 RDA6 RDA7 RDA8 RDA9
## 0.11145 0.06171 0.02137 0.01843 0.01370 0.00787 0.00677 0.00460 0.00246
## RDA10 RDA11 RDA12 RDA13 RDA14
## 0.00206 0.00154 0.00128 0.00124 0.00064
##
## Eigenvalues for unconstrained axes:
## PC1 PC2 PC3 PC4 PC5 PC6 PC7 PC8 PC9
## 0.04830 0.02129 0.01089 0.01025 0.00738 0.00559 0.00304 0.00194 0.00090
anova(rda)
## Permutation test for rda under reduced model
## Permutation: free
## Number of permutations: 999
##
## Model: rda(X = species, Y = varechem)
## Df Variance F Pr(>F)
## Model 14 0.25511 1.4968 0.051 .
## Residual 9 0.10957
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## plot output
plot(rda)
There are two quick ways to determine if the RDA accurately describes the data. The key word in RDA is redundancy, implying that your gradient is essentially described twice. Once by the species data and once by the environmental variables. The more similar those describe the same trend, the greater redundancy and the better the ordination. Here, the eigenvalues for the constrained axes are about 40% of the unconstrained values that suggests there is some degree of redundancy, which is good. When we conduct a permutation test on the RDA we find no significance.
The p value being close to 0.05 and the relatively high level of redundancy suggest that these nutrient variables are not accurately explaining all the variation in our dataset. When we look at the plot, we can see that RDA1 is comprised mostly of the metals. RDA2 instead appears to be described by other nutrients such as nitrogen and phosphorus. As it stands now though, we cannot present this dataset because of the poor explanation of variation and high p.value.
A cluster analysis groups together different sites or plots, typically by dissimiliarity. One of the more widely accepted methods is using the squared version of Ward’s method. Ward’s method uses a minimum variance criterion that minimizes the total within-cluster variance (see for more Murtagh and Legendre 2014).
We again calculate dissimiliarity indices using the Bray-Curtis method. On our initial plotting, it can be seen that are two general catgorizations of the sites. We can cut this dendrogram into separate groups either arbitrary based on the number of groups or by the height of the dendrogram.
## calculate distances
dis <- vegdist(species, method="bray")
## cluster analysis
clus <- hclust(dis, "ward.D2")
## plot cluster
plot(clus)
## cut the tree based on a specific height
grp <- cutree(clus, h=0.7)
## cut the tree into four clusters
grp <- cutree(clus, 4)
We can plot the dendrogram to be more aesthetically pleasing using the cluster
and factoextra
packages. Although this is not necessary. Here we can clearly highlight the clustered groups that we generated.
library(cluster)
library(factoextra)
fviz_dend(clus, k = 4, # Cut in four groups
cex = 0.7, cex.axis=2, # label size
color_labels_by_k = TRUE, # color labels by groups
rect = TRUE # Add rectangle around groups
)
The next question is do these cluster map on to a treatment or measured variate. We can compare the clustered groups against different environmental variables. Let’s compare the measured nitrogen against the groups that were generated. We can also common parametric analyses, such as ANOVAs, to determine if there are statistical differences among the groups.
## Compare environmental variables to cluster groups
boxplot(varechem$N ~ grp)
## Statistically compare groups
m1 <- aov(varechem$N ~ grp)
summary(m1)
## Df Sum Sq Mean Sq F value Pr(>F)
## grp 1 118.2 118.24 4.449 0.0465 *
## Residuals 22 584.6 26.57
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Finally, if these clusters make sense ecologically, are statistically different, and answering a research question, it makes sense to plot these clusters on the ordination. These ordinations can be added using the ordiellipse
function and can be used to determine the species that belong in the clusters.
## conduct ordination
ord <- rda(species)
## calculate priority
spp.priority <- colSums(species)
## plot ordination with clusters
plot(ord, type="n")
ordiellipse(ord, grp, lty = 2, col = "grey80", draw="polygon", alpha=150)
orditorp(ord, display = "species", cex = 0.7, col = "darkorange3", priority=spp.priority, air=0.8)
orditorp(ord, display = "sites", cex = 0.7, col = "darkslateblue", air=0.1)
Indicator species analysis explores statistical significance between particular groups. It is a particularly relevant follow up for ordinations to determine the key species that are driving differences in composition between groups. It can also identify unique species assembalges that are present in a group.
library(indicspecies)
data(dune.env)
## Compare species assemblages among manure levels
indval <- multipatt(dune, dune.env$Manure, control=how(nperm=999))
summary(indval)
##
## Multilevel pattern analysis
## ---------------------------
##
## Association function: IndVal.g
## Significance level (alpha): 0.05
##
## Total number of species: 30
## Selected number of species: 3
## Number of species associated to 1 group: 1
## Number of species associated to 2 groups: 0
## Number of species associated to 3 groups: 0
## Number of species associated to 4 groups: 2
##
## List of species associated to each combination:
##
## Group 1 #sps. 1
## stat p.value
## Vicilath 0.756 0.033 *
##
## Group 1+2+3+4 #sps. 2
## stat p.value
## Poatriv 0.964 0.001 ***
## Lolipere 0.878 0.013 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Compare species assemblages among management strategies
indval2 <- multipatt(dune, dune.env$Management, control=how(nperm=999))
summary(indval2)
##
## Multilevel pattern analysis
## ---------------------------
##
## Association function: IndVal.g
## Significance level (alpha): 0.05
##
## Total number of species: 30
## Selected number of species: 6
## Number of species associated to 1 group: 3
## Number of species associated to 2 groups: 0
## Number of species associated to 3 groups: 3
##
## List of species associated to each combination:
##
## Group BF #sps. 1
## stat p.value
## Vicilath 0.756 0.037 *
##
## Group HF #sps. 2
## stat p.value
## Rumeacet 0.851 0.004 **
## Trifprat 0.775 0.005 **
##
## Group BF+HF+SF #sps. 3
## stat p.value
## Poatriv 0.964 0.002 **
## Poaprat 0.896 0.027 *
## Lolipere 0.875 0.014 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Here we can see that the addition of manure results in the presence of two additional species that when manure is missing. Similarly, the different management strategies each form unique assemblages.
Try an example on your own from one of the datasets we looked at earlier. The multivar
dataset has measured the community composition of different shrub and open microsites along an elevational gradient over multiple years. The research question you are to explore is how does the shrub effect changes in community composition along this elevation gradient.
head(multivar)
## year Region Site Elevation Rep Microsite S.barbatus B.rubens
## 1 2016 Mojave 2.5 1100 1 Shrub 1 12
## 2 2016 Mojave 2.5 1100 1 Open 3 9
## 3 2016 Mojave 2.5 1100 2 Shrub 0 15
## 4 2016 Mojave 2.5 1100 2 Open 0 3
## 5 2016 Mojave 2.5 1100 3 Shrub 2 3
## 6 2016 Mojave 2.5 1100 3 Open 0 2
## E.cicutarium B.hordeaceus P.secunda D.kelloggii P.ovata C.barbigera
## 1 0 0 0 0 0 0
## 2 0 0 0 0 0 0
## 3 0 0 0 0 0 0
## 4 0 0 0 0 0 0
## 5 0 0 0 0 0 1
## 6 0 0 0 0 0 0
## L.gracilis C.exserta A.wrangelianus D.capitatum Pectocarya.spp.
## 1 0 0 0 0 3
## 2 0 0 0 0 2
## 3 0 0 0 0 1
## 4 0 0 0 0 0
## 5 0 0 0 0 5
## 6 0 0 0 0 0
## E.penduliflora A.tessellata A.adscensionis C.fremontii E.minutiflora
## 1 0 0 0 0 0
## 2 0 0 0 0 0
## 3 0 0 0 0 0
## 4 0 0 0 0 0
## 5 0 0 0 0 0
## 6 0 0 0 0 0
## P.membranaceum L.lasiocarpa M.glabrata R.neomexicana E.eremicum C.rigida
## 1 0 0 0 0 0 0
## 2 0 0 0 0 0 0
## 3 0 0 0 0 0 0
## 4 0 1 0 0 0 0
## 5 0 0 0 0 0 0
## 6 0 0 0 0 0 0
## P.crenulata E.wallacei C.brevipes C.claviformis Erinoginum.spp.
## 1 0 0 0 0 0
## 2 0 0 0 0 0
## 3 0 0 0 0 0
## 4 0 0 0 0 2
## 5 0 0 0 0 0
## 6 0 0 0 0 0
## M.bellioides E.exilis C.lasiophyllus L.setosissima A.tuberosa B.nigria
## 1 0 0 0 0 0 0
## 2 1 0 0 0 0 0
## 3 2 0 0 0 0 0
## 4 0 0 0 0 0 0
## 5 0 0 0 0 0 0
## 6 0 0 0 0 0 0
## N.demissum A.lentiginosus A.grandiflora H.vulgare C.intermedia
## 1 0 0 0 0 0
## 2 0 0 0 0 0
## 3 0 0 0 0 0
## 4 0 0 0 0 0
## 5 0 0 0 0 0
## 6 0 0 0 0 0