Table of Contents
Section: Numerical classification
Evaluation of classification results
The following examples are using the same dataset - Vltava data (please check the dataset description to get familiar with it).
Example 1: Project results of cluster analysis onto the diagram of unconstrained ordination
Combining cluster analysis with unconstrained ordination is a powerful tool, which will help us to understand 1) what is the relationship among individual groups/clusters, 2) whether the clusters are well separated from each other (ie the classification was successful in finding well-defined groups), and 3) possibly what is the relationship of individual clusters with the environment.
First, use Vltava data and calculate cluster analysis. Here we opt for Ward's algorithm combined with the Bray-Curtis distance metric. Note several things: 1) since the raw species composition data are on a percentage scale (cover of species estimated in the field using Braun-Blanquet scale, converted into percentage values), we first log transform them to reduce the importance of species with very high cover (using log (x+1) transformation to ensure that zero after log transformation still stays zero). 2) Since Ward's clustering algorithm requires metric distances, and Bray-Curtis distance is not metric, we need to square-root it to make it metric.
library (vegan) vltava.spe <- read.delim ('https://raw.githubusercontent.com/zdealveindy/anadat-r/master/data/vltava-spe.txt', row.names = 1) spe_log <- log1p (vltava.spe) # log-transform raw species composition data dist_bc <- vegdist (spe_log, method = 'bray') # calculate Bray-Curtis distance matrix dist_bc_sqrt <- sqrt (dist_bc) # square-root BC distance matrix to make distances metric clus_ward <- hclust (dist_bc_sqrt, method = 'ward.D2') # calculate Ward's algorithm (using the correct ''method = 'ward.D2''') clus_ward_cut <- cutree (clus_ward, k = 5) # "cut the tree" - to which groups individual samples belong?
First, draw the cluster dendrogram, and identify individual clusters as rectangles in the dendrogram (function rect.hclust). One important thing to remember: the order in which clusters are numbered in cutree function is not identical with the order of clusters in the dendrogram (i.e. the first rectangular on the left does not need to be cluster 1). See [[en:hier-agglom_examples#example_2ward_cluster_algorithm_applied_on_barro_colorado_island_data|Example 2 in Cluster analysis to see details how to solve this problem (this also explains why we need to create clus_in_dendro object below).
plot (clus_ward, cex = .5) # argument cex reduced the size of the dendrogram leaf labels to make them readable clus_in_dendro <- unique (clus_ward_cut[clus_ward$order]) # make sure to know which box is which cluster! rect.hclust (clus_ward, k = 5, border = clus_in_dendro) legend ('topleft', legend = paste ('Cluster', 1:5), col = 1:5, pch = 22, bty = 'n')
Now, display the same clusters in the unconstrained ordination diagram. We should follow the rule that the distance measure used in ordination should be the same as the distance measure used in cluster analysis, i.e. Bray-Curtis in our example. The only two unconstrained ordination methods which can implement Bray-Curtis distance are non-metric multidimensional scaling (NMDS) and principal coordinate analysis (PCoA) - we choose the first one here. In the case of NMDS, it does not matter whether we use the original Bray-Curtis distances (object dist_bc) or their square rooted version used in Ward's cluster algorithm (dist_bc_sqrt), since NMDS as the first step converts distances into ranks, used for the calculation, and square-rooting is a monotonic transformation, which preserves the rank of the values:
NMDS <- metaMDS (dist_bc) ordiplot (NMDS, type = 'n', main = 'NMDS ordination diagram (Bray-Curtis distances)') points (NMDS, pch = clus_ward_cut+20, bg = clus_ward_cut) legend ('topright', legend = 1:5, pch = (1:5)+20, pt.bg = 1:5)
(adding 20 to the cluster number 1 to 5 allowed using symbols 21 to 25 for displaying the samples in the ordination diagram).
From both the dendrogram and ordination diagram we see that cluster 4 and 5 are rather similar in species composition. In the dendrogram, cluster 4 and 5 merges at a rather low height of 1.5 (note that height here does not represent the Bray-Curtis distance, which is in the range between 0 and 1, but the total error sum of squares, which is the criterion association with Ward's clustering algorithm). In the ordination diagram, their samples are rather close to each other and partially overlapping. In contrast, cluster 2 (the red one) is rather different from other clusters - it branches rather high in the dendrogram, and is quite far from other groups on the ordination diagram.
Also, from the ordination, it seems that the classification results in reasonably distinct clusters, even though there are some samples “ejected” to the other groups, e.g. the lonely sample belonging to cluster 3 (green) in the far right part of the ordination diagram, among samples of cluster 2 (red).
Example 2: Environmental differences between clusters
If we have environmental information about individual samples, it would be useful to quantify and describe environmental differences among individual clusters. In the end, we are ecologists, and this is what we are usually interested in.
Let's import environmental data vltava.env, and use a subset of four variables to draw differences among clusters.
vltava.env <- read.delim ('https://raw.githubusercontent.com/zdealveindy/anadat-r/master/data/vltava-env.txt')
Let's select just two of them (elevation and soil pH) and draw boxplots illustrating the differences:
par (mfrow = c(1,2)) boxplot (vltava.env$ELEVATION ~ clus_ward_cut, notch = TRUE, col = 1:5, xlab = 'Cluster', ylab = 'Elevation [m a.s.l.]', outline = FALSE) boxplot (vltava.env$pH ~ clus_ward_cut, notch = TRUE, col = 1:5, xlab = 'Cluster', ylab = 'Soil pH', outline = FALSE)
Colours of boxplots respect the colours of clusters in Example 1; non-overlapping notches indicate that env. variables among given clusters are likely to be significantly different (check the help file ?boxplot to see details what notch = TRUE means). We also deliberately decided to suppress displaying outliers in the boxplots (outline = FALSE).
Let's take elevation and test the differences between clusters. First, we should do the analysis of variance (ANOVA), which will reveal whether the mean of at least one of the groups is significantly different from the overall mean. If ANOVA is significant, we can go a step further and test pair-wise differences between clusters, using e.g. Tukey HSD comparison.
clus_ward_cut_fact <- as.factor (clus_ward_cut) ANOVA_elev <- aov (vltava.env$ELEVATION ~ clus_ward_cut_fact) summary (ANOVA_elev) # F = 8.832, P < 0.001
Yes - as you see ANOVA is significant, so we can happily resolve pairwise comparison. R contains the function Tukey.HSD, which can do that, but here we use a bit more advanced function HSD.test from the package agricolae (install it if you haven't):
library (agricolae) hsd_elev <- HSD.test(ANOVA_elev, trt = "clus_ward_cut_fact", group = T, console = T) hsd_elev
$statistics MSerror Df Mean CV 548.1152 92 466.1649 5.022227 $parameters test name.t ntr StudentizedRange alpha Tukey clus_ward_cut_fact 5 3.935217 0.05 $means vltava.env$ELEVATION std r Min Max Q25 Q50 Q75 1 464.0625 24.91444 16 420 510 451.25 470.0 476.25 2 434.1429 11.58153 14 405 455 430.00 435.0 438.75 3 474.4545 25.68449 33 430 530 450.00 475.0 495.00 4 464.6429 24.05637 14 430 510 446.25 457.5 482.50 5 477.6500 23.77366 20 430 520 458.75 475.0 495.00 $comparison NULL $groups vltava.env$ELEVATION groups 5 477.6500 a 3 474.4545 a 4 464.6429 a 1 464.0625 a 2 434.1429 b
Important on the output of the function is the last table, which shows differences among clusters, expressed as letters (clusters with the same letter are not significantly different). We can use these letters and display them onto the boxplot (a small added space below the boxes by extending a bit ylim range in boxplot function will help to squeeze them there):
par (mfrow = c(1,2)) boxplot (vltava.env$ELEVATION ~ clus_ward_cut, notch = TRUE, col = 1:5, xlab = 'Cluster', ylab = 'Elevation [m a.s.l.]', outline = FALSE, ylim = c(410, 530)) group_no_elev <- as.numeric (rownames (hsd_elev$groups)) letter_elev <- hsd_elev$groups$groups text (x = group_no_elev, y = 415, labels = letter_elev) ANOVA_pH <- aov (vltava.env$pH ~ clus_ward_cut_fact) summary (ANOVA_pH) # F = 10.92, P < 0.001 hsd_pH <- HSD.test(ANOVA_pH, trt = "clus_ward_cut_fact", group = T, console = T) boxplot (vltava.env$pH ~ clus_ward_cut, notch = TRUE, col = 1:5, xlab = 'Cluster', ylab = 'Soil pH', outline = FALSE, ylim = c(3.1, 6.4)) group_no_pH <- as.numeric (rownames (hsd$groups)) letter_pH <- hsd$groups$groups text (x = group_no_pH, y = 3.2, letter_pH)
We can see that for elevation, plots in cluster 2 are significantly lower than other clusters. When considering soil pH, plots in cluster 3 and 5 are significantly more acid than those in cluster 1 and 2, while cluster 4 stands in between.
Example 3: Fidelity of species to individual clusters
Two common ways how to identify which species are diagnostic to which cluster is the IndVal approach and Φ (phi) coefficient of fidelity. Here we will practice how to calculate Φ, how to test it, how to print results and interpret them. We use the same data as in the exercises above (Vltava data, classified by Ward's clustering algorithm applied on square-rooted Bray-Curtis distances).
We will need library indicspecies written by Miquel de Cáceres to calculate fidelity (install.packages ('indicspecies') if you haven't):
library (indicspecies)
First, we need to transform the species composition data into presences absence (Φ coefficient is calculated on presence-absence data), and then apply the function multipatt (abbreviation for Multi-level pattern analysis) to calculate fidelity. We set the argument fun = 'r.g', where r indicates correlation coefficient (Φ coefficient is closely related to Pearson's r correlation) and g stands for standardizing clusters to the same size (“group size”).
vltava.spe.pa <- decostand (vltava.spe, method = 'pa') phi <- multipatt (vltava.spe.pa, cluster = clus_ward_cut, fun = 'r.g')
The object phi created by function multipatt is a list, and the component $sign contains the table with species indicating to which cluster they are diagnostic. We can filter out the species which are significant (P<0.01 here) and print them:
re <- phi$sign[phi$sign$p.value<=0.01,] re
s.1 s.2 s.3 s.4 s.5 index stat p.value Abiealb23 0 0 0 1 1 15 0.6642737 0.005 Acerpla23 1 1 0 0 0 6 0.4154034 0.005 Alnuglu23 0 1 0 0 0 2 0.8164966 0.005 Fraxexc23 1 1 0 0 0 6 0.3918460 0.005 Pinusyl23 1 0 1 1 0 19 0.4382082 0.005 Prunpad23 0 1 0 0 0 2 0.6740816 0.005 ...
The table contains species abbreviation, and then an indication to which group (s.1 to s.5) is the species diagnostic. The column index indicates which combination of groups the species is diagnostic to (1 means cluster 1, 2 means cluster 2, 6 means cluster 1 & 2, and so on). The column stat is for Φ coefficient (in this case), and p.value for significance.
It may be useful to actually sort the species in the table (otherwise sorted as they appear in the raw data) according to which cluster they are diagnostic to. We want to first list species diagnostic to cluster 1, then cluster 2, and so on, and the cluster 1 and 2, 2 and 3, and so on. Below is using dplyr function arrange, which does exactly this. Note that the sorting is a bit complex: we want to make sure that we first display species diagnostic to only a single cluster (that is why we need to first calculate variable no_diag_clust, indicating what number of clusters the species is diagnostic to; then we sort species in descending order according to their diagnostic values for clusters s.1 to s.5, and finally, also in descending order, according to their Φ value (column stat):
library (tidyverse) no_diag_sp <- rowSums (select (re, s.1, s.2, s.3, s.4, s.5)) re %>% arrange (no_diag_sp, desc (s.1), desc (s.2), desc (s.3), desc (s.4), desc (s.5), desc (stat))
s.1 s.2 s.3 s.4 s.5 index stat p.value Vinchir1 1 0 0 0 0 1 0.6804190 0.005 Euphcyp1 1 0 0 0 0 1 0.6666667 0.005 Hylomax1 1 0 0 0 0 1 0.5693680 0.005 Campper1 1 0 0 0 0 1 0.5479992 0.005 Camprot1 1 0 0 0 0 1 0.5244745 0.005 Genitin1 1 0 0 0 0 1 0.5163978 0.005 Bracpin1 1 0 0 0 0 1 0.5092775 0.005 Polyodo1 1 0 0 0 0 1 0.5092466 0.005 Lychvis1 1 0 0 0 0 1 0.4958468 0.005 Silenut1 1 0 0 0 0 1 0.4958468 0.005 Cytinig1 1 0 0 0 0 1 0.4611528 0.005 Achitan1 1 0 0 0 0 1 0.4588315 0.005 Diancar1 1 0 0 0 0 1 0.4588315 0.005 Pimpsax1 1 0 0 0 0 1 0.4588315 0.005 Trifalp1 1 0 0 0 0 1 0.4588315 0.005 Hypeper1 1 0 0 0 0 1 0.4526211 0.010 Poa.nem1 1 0 0 0 0 1 0.4495051 0.005 Digigra1 1 0 0 0 0 1 0.4246831 0.005 Anthram1 1 0 0 0 0 1 0.4215364 0.005 Geniger1 1 0 0 0 0 1 0.4215364 0.010 Festpal1 1 0 0 0 0 1 0.4010513 0.005 Liguvul23 1 0 0 0 0 1 0.3947710 0.010 Allisen1 1 0 0 0 0 1 0.3947710 0.010 Bromben1 1 0 0 0 0 1 0.3947710 0.005 Festovi1 1 0 0 0 0 1 0.3525932 0.010 Phalaru1 0 1 0 0 0 2 0.9107143 0.005 Lycoeur1 0 1 0 0 0 2 0.9085257 0.005 Stacsyl1 0 1 0 0 0 2 0.8635755 0.005 Festgig1 0 1 0 0 0 2 0.8214286 0.005 Alnuglu23 0 1 0 0 0 2 0.8164966 0.005 Calysep1 0 1 0 0 0 2 0.8014485 0.005 Myospal1 0 1 0 0 0 2 0.7959350 0.005 Stelnem1 0 1 0 0 0 2 0.7777560 0.005 Sympoff1 0 1 0 0 0 2 0.7068371 0.005 Ranurep1 0 1 0 0 0 2 0.6931919 0.005 Aegopod1 0 1 0 0 0 2 0.6770620 0.005 Prunpad23 0 1 0 0 0 2 0.6740816 0.005 Impagla1 0 1 0 0 0 2 0.6530294 0.005 Valeexc1 0 1 0 0 0 2 0.6123724 0.005 Carebri1 0 1 0 0 0 2 0.6089894 0.005 Lysivul1 0 1 0 0 0 2 0.6000616 0.005 Rumeaqu1 0 1 0 0 0 2 0.5566219 0.005 Salifra23 0 1 0 0 0 2 0.5547002 0.005 Bidefro1 0 1 0 0 0 2 0.5547002 0.005 Galiapa1 0 1 0 0 0 2 0.5547002 0.005 Mentx.v1 0 1 0 0 0 2 0.5547002 0.005 Cirsole1 0 1 0 0 0 2 0.5532670 0.005 Angesyl1 0 1 0 0 0 2 0.5476190 0.005 Soladul1 0 1 0 0 0 2 0.5173224 0.005 Artevul1 0 1 0 0 0 2 0.4923660 0.005 Chaehir1 0 1 0 0 0 2 0.4923660 0.005 Myosaqu1 0 1 0 0 0 2 0.4923660 0.005 Humulup1 0 1 0 0 0 2 0.4611148 0.005 Lythsal1 0 1 0 0 0 2 0.4572230 0.010 Urtidio1 0 1 0 0 0 2 0.4340518 0.005 Aconvar1 0 1 0 0 0 2 0.4232074 0.005 Cusceur1 0 1 0 0 0 2 0.4232074 0.010 Poa.pal1 0 1 0 0 0 2 0.4232074 0.005 Lamimac1 0 1 0 0 0 2 0.4220196 0.010 Glycmax1 0 1 0 0 0 2 0.3894008 0.010 Euoneur23 0 1 0 0 0 2 0.3474972 0.010 Vaccmyr1 0 0 1 0 0 3 0.6011300 0.005 Avenfle1 0 0 1 0 0 3 0.4593426 0.005 Betupen23 0 0 1 0 0 3 0.4327338 0.010 Callvul1 0 0 1 0 0 3 0.3791005 0.010 Galespc1 0 0 0 1 0 4 0.3804001 0.005 Lunared1 0 0 0 0 1 5 0.7097582 0.005 Milieff1 0 0 0 0 1 5 0.6437549 0.005 Mercper1 0 0 0 0 1 5 0.5897678 0.005 Prenpur1 0 0 0 0 1 5 0.5897678 0.005 Cardimp1 0 0 0 0 1 5 0.5230789 0.005 Seneova1 0 0 0 0 1 5 0.5158033 0.005 Acerpla23 1 1 0 0 0 6 0.4154034 0.005 Fraxexc23 1 1 0 0 0 6 0.3918460 0.005 Fragmos1 1 1 0 0 0 6 0.3616508 0.010 Hiersab1 1 0 1 0 0 7 0.5596310 0.005 Hierlac1 1 0 1 0 0 7 0.4450991 0.005 Querpet23 1 0 1 0 0 7 0.4330667 0.005 Tilicor23 1 0 0 1 0 8 0.5338429 0.005 Clinvul1 1 0 0 1 0 8 0.4476660 0.005 Pulmobs1 1 0 0 1 0 8 0.3709077 0.005 Scronod1 0 1 0 1 0 11 0.4225365 0.005 Fagusyl23 0 0 1 0 1 14 0.3873290 0.010 Abiealb23 0 0 0 1 1 15 0.6642737 0.005 Dryofil1 0 0 0 1 1 15 0.4911746 0.005 Mycemur1 0 0 0 1 1 15 0.4170078 0.005 Actaspi1 0 0 0 1 1 15 0.4035034 0.005 Dryocar1 0 0 0 1 1 15 0.3994412 0.010 Moehtri1 0 0 0 1 1 15 0.3829024 0.005 Stelhol1 1 1 0 1 0 17 0.4742312 0.005 Pinusyl23 1 0 1 1 0 19 0.4382082 0.005 Luzuluz1 1 0 1 1 0 19 0.3806619 0.005 Impanol1 0 1 0 1 1 24 0.5359159 0.005 Oxalace1 0 1 0 1 1 24 0.5351096 0.010 Asareur1 0 1 0 1 1 24 0.5098884 0.005 Dryodil1 0 1 0 1 1 24 0.3991077 0.010 Galemon1 0 1 0 1 1 24 0.3883148 0.010
This is very similar to the way how vegetation ecologists sort synoptic tables of species when they describe vegetation types distinguished in their study.
Finally, we can also use the function summary to display the output in yet another structured way (using the same significance threshold P < 0.01 as in the table above):
summary (phi, alpha = 0.01)
Multilevel pattern analysis
---------------------------
Association function: r.g
Significance level (alpha): 0.01
Total number of species: 274
Selected number of species: 97
Number of species associated to 1 group: 72
Number of species associated to 2 groups: 17
Number of species associated to 3 groups: 8
Number of species associated to 4 groups: 0
List of species associated to each combination:
Group 1 #sps. 25
stat p.value
Vinchir1 0.680 0.005 **
Euphcyp1 0.667 0.005 **
Hylomax1 0.569 0.005 **
Campper1 0.548 0.005 **
Camprot1 0.524 0.005 **
Genitin1 0.516 0.005 **
Bracpin1 0.509 0.005 **
Polyodo1 0.509 0.005 **
Lychvis1 0.496 0.005 **
Silenut1 0.496 0.005 **
Cytinig1 0.461 0.005 **
Achitan1 0.459 0.005 **
Diancar1 0.459 0.005 **
Pimpsax1 0.459 0.005 **
Trifalp1 0.459 0.005 **
Hypeper1 0.453 0.010 **
Poa.nem1 0.450 0.005 **
Digigra1 0.425 0.005 **
Anthram1 0.422 0.005 **
Geniger1 0.422 0.010 **
Festpal1 0.401 0.005 **
Liguvul23 0.395 0.010 **
Allisen1 0.395 0.010 **
Bromben1 0.395 0.005 **
Festovi1 0.353 0.010 **
Group 2 #sps. 36
stat p.value
Phalaru1 0.911 0.005 **
Lycoeur1 0.909 0.005 **
Stacsyl1 0.864 0.005 **
Festgig1 0.821 0.005 **
Alnuglu23 0.816 0.005 **
Calysep1 0.801 0.005 **
Myospal1 0.796 0.005 **
Stelnem1 0.778 0.005 **
Sympoff1 0.707 0.005 **
Ranurep1 0.693 0.005 **
Aegopod1 0.677 0.005 **
Prunpad23 0.674 0.005 **
Impagla1 0.653 0.005 **
Valeexc1 0.612 0.005 **
Carebri1 0.609 0.005 **
Lysivul1 0.600 0.005 **
Rumeaqu1 0.557 0.005 **
Salifra23 0.555 0.005 **
Bidefro1 0.555 0.005 **
Galiapa1 0.555 0.005 **
Mentx.v1 0.555 0.005 **
Cirsole1 0.553 0.005 **
Angesyl1 0.548 0.005 **
Soladul1 0.517 0.005 **
Artevul1 0.492 0.005 **
Chaehir1 0.492 0.005 **
Myosaqu1 0.492 0.005 **
Humulup1 0.461 0.005 **
Lythsal1 0.457 0.010 **
Urtidio1 0.434 0.005 **
Aconvar1 0.423 0.005 **
Cusceur1 0.423 0.010 **
Poa.pal1 0.423 0.005 **
Lamimac1 0.422 0.010 **
Glycmax1 0.389 0.010 **
Euoneur23 0.347 0.010 **
Group 3 #sps. 4
stat p.value
Vaccmyr1 0.601 0.005 **
Avenfle1 0.459 0.005 **
Betupen23 0.433 0.010 **
Callvul1 0.379 0.010 **
Group 4 #sps. 1
stat p.value
Galespc1 0.38 0.005 **
Group 5 #sps. 6
stat p.value
Lunared1 0.710 0.005 **
Milieff1 0.644 0.005 **
Mercper1 0.590 0.005 **
Prenpur1 0.590 0.005 **
Cardimp1 0.523 0.005 **
Seneova1 0.516 0.005 **
Group 1+2 #sps. 3
stat p.value
Acerpla23 0.415 0.005 **
Fraxexc23 0.392 0.005 **
Fragmos1 0.362 0.010 **
Group 1+3 #sps. 3
stat p.value
Hiersab1 0.560 0.005 **
Hierlac1 0.445 0.005 **
Querpet23 0.433 0.005 **
Group 1+4 #sps. 3
stat p.value
Tilicor23 0.534 0.005 **
Clinvul1 0.448 0.005 **
Pulmobs1 0.371 0.005 **
Group 2+4 #sps. 1
stat p.value
Scronod1 0.423 0.005 **
Group 3+5 #sps. 1
stat p.value
Fagusyl23 0.387 0.01 **
...
...
...
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1