## PCA# we use the same dataset of flower characteristics of three species of irispairs(iris, col = iris$Species)
head(iris)## Sepal.Length Sepal.Width Petal.Length Petal.Width Species## 1 5.1 3.5 1.4 0.2 setosa## 2 4.9 3.0 1.4 0.2 setosa## 3 4.7 3.2 1.3 0.2 setosa## 4 4.6 3.1 1.5 0.2 setosa## 5 5.0 3.6 1.4 0.2 setosa## 6 5.4 3.9 1.7 0.4 setosapca =prcomp(iris[, 1:4], scale = T) # always set scale = T# when data is very skewed --> better transform e.g. logsummary(pca)## Importance of components:## PC1 PC2 PC3 PC4## Standard deviation 1.7084 0.9560 0.38309 0.14393## Proportion of Variance 0.7296 0.2285 0.03669 0.00518## Cumulative Proportion 0.7296 0.9581 0.99482 1.00000# standard deviation^2 is variance!!!# cum prop of PC2 is the variance that is visualized in a biplot# plot the result# absolute variance of each componentplot(pca) # see row1 of the summary(pca): (sd)^2 = variance
# rel variance of each componentbarplot(summary(pca)$importance[2, ], ylab="proportion of variance explained") # displays % of variance explained by PCs
# Biplotbiplot(pca) # displays PC1 and PC2 AND rotation (vectors) of the different variables AND observations## distance-based approach: NMDSlibrary(vegan)## Loading required package: permute## Loading required package: lattice## This is vegan 2.6-4
?vegan# community dataset for plants in dunes (included in vegan package):data("dune")str(dune) # display structure of the dataset## 'data.frame': 20 obs. of 30 variables:## $ Achimill: num 1 3 0 0 2 2 2 0 0 4 ...## $ Agrostol: num 0 0 4 8 0 0 0 4 3 0 ...## $ Airaprae: num 0 0 0 0 0 0 0 0 0 0 ...## $ Alopgeni: num 0 2 7 2 0 0 0 5 3 0 ...## $ Anthodor: num 0 0 0 0 4 3 2 0 0 4 ...## $ Bellpere: num 0 3 2 2 2 0 0 0 0 2 ...## $ Bromhord: num 0 4 0 3 2 0 2 0 0 4 ...## $ Chenalbu: num 0 0 0 0 0 0 0 0 0 0 ...## $ Cirsarve: num 0 0 0 2 0 0 0 0 0 0 ...## $ Comapalu: num 0 0 0 0 0 0 0 0 0 0 ...## $ Eleopalu: num 0 0 0 0 0 0 0 4 0 0 ...## $ Elymrepe: num 4 4 4 4 4 0 0 0 6 0 ...## $ Empenigr: num 0 0 0 0 0 0 0 0 0 0 ...## $ Hyporadi: num 0 0 0 0 0 0 0 0 0 0 ...## $ Juncarti: num 0 0 0 0 0 0 0 4 4 0 ...## $ Juncbufo: num 0 0 0 0 0 0 2 0 4 0 ...## $ Lolipere: num 7 5 6 5 2 6 6 4 2 6 ...## $ Planlanc: num 0 0 0 0 5 5 5 0 0 3 ...## $ Poaprat : num 4 4 5 4 2 3 4 4 4 4 ...## $ Poatriv : num 2 7 6 5 6 4 5 4 5 4 ...## $ Ranuflam: num 0 0 0 0 0 0 0 2 0 0 ...## $ Rumeacet: num 0 0 0 0 5 6 3 0 2 0 ...## $ Sagiproc: num 0 0 0 5 0 0 0 2 2 0 ...## $ Salirepe: num 0 0 0 0 0 0 0 0 0 0 ...## $ Scorautu: num 0 5 2 2 3 3 3 3 2 3 ...## $ Trifprat: num 0 0 0 0 2 5 2 0 0 0 ...## $ Trifrepe: num 0 5 2 1 2 5 2 2 3 6 ...## $ Vicilath: num 0 0 0 0 0 0 0 0 0 1 ...## $ Bracruta: num 0 0 2 2 2 6 2 2 2 2 ...## $ Callcusp: num 0 0 0 0 0 0 0 0 0 0 ...?dunesummary(dune) # display summary of the dataset (summary statistics for each variable)## Achimill Agrostol Airaprae Alopgeni Anthodor ## Min. :0.0 Min. :0.0 Min. :0.00 Min. :0.00 Min. :0.00 ## 1st Qu.:0.0 1st Qu.:0.0 1st Qu.:0.00 1st Qu.:0.00 1st Qu.:0.00 ## Median :0.0 Median :1.5 Median :0.00 Median :0.00 Median :0.00 ## Mean :0.8 Mean :2.4 Mean :0.25 Mean :1.80 Mean :1.05 ## 3rd Qu.:2.0 3rd Qu.:4.0 3rd Qu.:0.00 3rd Qu.:3.25 3rd Qu.:2.25 ## Max. :4.0 Max. :8.0 Max. :3.00 Max. :8.00 Max. :4.00 ## Bellpere Bromhord Chenalbu Cirsarve Comapalu ## Min. :0.00 Min. :0.00 Min. :0.00 Min. :0.0 Min. :0.0 ## 1st Qu.:0.00 1st Qu.:0.00 1st Qu.:0.00 1st Qu.:0.0 1st Qu.:0.0 ## Median :0.00 Median :0.00 Median :0.00 Median :0.0 Median :0.0 ## Mean :0.65 Mean :0.75 Mean :0.05 Mean :0.1 Mean :0.2 ## 3rd Qu.:2.00 3rd Qu.:0.50 3rd Qu.:0.00 3rd Qu.:0.0 3rd Qu.:0.0 ## Max. :3.00 Max. :4.00 Max. :1.00 Max. :2.0 Max. :2.0 ## Eleopalu Elymrepe Empenigr Hyporadi Juncarti ## Min. :0.00 Min. :0.0 Min. :0.0 Min. :0.00 Min. :0.00 ## 1st Qu.:0.00 1st Qu.:0.0 1st Qu.:0.0 1st Qu.:0.00 1st Qu.:0.00 ## Median :0.00 Median :0.0 Median :0.0 Median :0.00 Median :0.00 ## Mean :1.25 Mean :1.3 Mean :0.1 Mean :0.45 Mean :0.90 ## 3rd Qu.:1.00 3rd Qu.:4.0 3rd Qu.:0.0 3rd Qu.:0.00 3rd Qu.:0.75 ## Max. :8.00 Max. :6.0 Max. :2.0 Max. :5.00 Max. :4.00 ## Juncbufo Lolipere Planlanc Poaprat Poatriv ## Min. :0.00 Min. :0.0 Min. :0.0 Min. :0.0 Min. :0.00 ## 1st Qu.:0.00 1st Qu.:0.0 1st Qu.:0.0 1st Qu.:0.0 1st Qu.:0.00 ## Median :0.00 Median :2.0 Median :0.0 Median :3.0 Median :4.00 ## Mean :0.65 Mean :2.9 Mean :1.3 Mean :2.4 Mean :3.15 ## 3rd Qu.:0.00 3rd Qu.:6.0 3rd Qu.:3.0 3rd Qu.:4.0 3rd Qu.:5.00 ## Max. :4.00 Max. :7.0 Max. :5.0 Max. :5.0 Max. :9.00 ## Ranuflam Rumeacet Sagiproc Salirepe Scorautu ## Min. :0.0 Min. :0.0 Min. :0 Min. :0.00 Min. :0.0 ## 1st Qu.:0.0 1st Qu.:0.0 1st Qu.:0 1st Qu.:0.00 1st Qu.:2.0 ## Median :0.0 Median :0.0 Median :0 Median :0.00 Median :2.0 ## Mean :0.7 Mean :0.9 Mean :1 Mean :0.55 Mean :2.7 ## 3rd Qu.:2.0 3rd Qu.:0.5 3rd Qu.:2 3rd Qu.:0.00 3rd Qu.:3.0 ## Max. :4.0 Max. :6.0 Max. :5 Max. :5.00 Max. :6.0 ## Trifprat Trifrepe Vicilath Bracruta Callcusp ## Min. :0.00 Min. :0.00 Min. :0.0 Min. :0.00 Min. :0.0 ## 1st Qu.:0.00 1st Qu.:1.00 1st Qu.:0.0 1st Qu.:1.50 1st Qu.:0.0 ## Median :0.00 Median :2.00 Median :0.0 Median :2.00 Median :0.0 ## Mean :0.45 Mean :2.35 Mean :0.2 Mean :2.45 Mean :0.5 ## 3rd Qu.:0.00 3rd Qu.:3.00 3rd Qu.:0.0 3rd Qu.:4.00 3rd Qu.:0.0 ## Max. :5.00 Max. :6.00 Max. :2.0 Max. :6.00 Max. :4.0NMDS =metaMDS(dune)## Run 0 stress 0.1192678 ## Run 1 stress 0.1886532 ## Run 2 stress 0.1192678 ## ... Procrustes: rmse 1.258026e-05 max resid 3.164268e-05 ## ... Similar to previous best## Run 3 stress 0.1183186 ## ... New best solution## ... Procrustes: rmse 0.02027014 max resid 0.06496123 ## Run 4 stress 0.1183186 ## ... Procrustes: rmse 4.230384e-06 max resid 1.475206e-05 ## ... Similar to previous best## Run 5 stress 0.1192678 ## Run 6 stress 0.1809578 ## Run 7 stress 0.1889638 ## Run 8 stress 0.1808911 ## Run 9 stress 0.3680059 ## Run 10 stress 0.1183186 ## ... Procrustes: rmse 1.159455e-05 max resid 3.442588e-05 ## ... Similar to previous best## Run 11 stress 0.1183186 ## ... Procrustes: rmse 4.643337e-06 max resid 1.388314e-05 ## ... Similar to previous best## Run 12 stress 0.2361935 ## Run 13 stress 0.1192678 ## Run 14 stress 0.1192679 ## Run 15 stress 0.1183186 ## ... Procrustes: rmse 1.447333e-05 max resid 4.510486e-05 ## ... Similar to previous best## Run 16 stress 0.1192678 ## Run 17 stress 0.1192678 ## Run 18 stress 0.1192679 ## Run 19 stress 0.1808911 ## Run 20 stress 0.1808911 ## *** Best solution repeated 4 times# algorithm is iterativeNMDS # gives information on NMDS: distance measure, stress (should be low)## ## Call:## metaMDS(comm = dune) ## ## global Multidimensional Scaling using monoMDS## ## Data: dune ## Distance: bray ## ## Dimensions: 2 ## Stress: 0.1183186 ## Stress type 1, weak ties## Best solution was repeated 4 times in 20 tries## The best solution was from try 3 (random start)## Scaling: centring, PC rotation, halfchange scaling ## Species: expanded scores based on 'dune'# stress of >= 0.2 = be suspicious, stress >=0.3 indicates that ordination is arbitrary# increase Dimensions if k is too highordiplot(NMDS, type ="t") #"t" = text
# if we have time:# distance measure can be changed (default is Bray-Curtis): see ?vegdist # some recommendations thereNMDS2 =metaMDS(dune, distance="euclidean")## Run 0 stress 0.1174523 ## Run 1 stress 0.1174523 ## ... Procrustes: rmse 3.580249e-06 max resid 1.150625e-05 ## ... Similar to previous best## Run 2 stress 0.1174523 ## ... Procrustes: rmse 2.054645e-06 max resid 3.818306e-06 ## ... Similar to previous best## Run 3 stress 0.1174523 ## ... Procrustes: rmse 2.128171e-06 max resid 5.656827e-06 ## ... Similar to previous best## Run 4 stress 0.1174523 ## ... Procrustes: rmse 2.022272e-06 max resid 6.26572e-06 ## ... Similar to previous best## Run 5 stress 0.1174523 ## ... Procrustes: rmse 2.488733e-06 max resid 8.218612e-06 ## ... Similar to previous best## Run 6 stress 0.1177339 ## ... Procrustes: rmse 0.01706434 max resid 0.05525565 ## Run 7 stress 0.1174523 ## ... New best solution## ... Procrustes: rmse 1.353662e-06 max resid 2.894393e-06 ## ... Similar to previous best## Run 8 stress 0.1174523 ## ... Procrustes: rmse 1.40202e-06 max resid 4.031046e-06 ## ... Similar to previous best## Run 9 stress 0.1174523 ## ... Procrustes: rmse 2.805574e-06 max resid 7.787809e-06 ## ... Similar to previous best## Run 10 stress 0.1174523 ## ... Procrustes: rmse 1.140966e-06 max resid 2.993735e-06 ## ... Similar to previous best## Run 11 stress 0.1177339 ## ... Procrustes: rmse 0.0170643 max resid 0.05525551 ## Run 12 stress 0.1174523 ## ... Procrustes: rmse 1.964847e-06 max resid 5.765074e-06 ## ... Similar to previous best## Run 13 stress 0.1174523 ## ... Procrustes: rmse 1.231259e-06 max resid 2.345061e-06 ## ... Similar to previous best## Run 14 stress 0.1174523 ## ... Procrustes: rmse 1.055035e-06 max resid 2.452641e-06 ## ... Similar to previous best## Run 15 stress 0.1177339 ## ... Procrustes: rmse 0.01706422 max resid 0.05525614 ## Run 16 stress 0.1177339 ## ... Procrustes: rmse 0.01706442 max resid 0.05525521 ## Run 17 stress 0.1177339 ## ... Procrustes: rmse 0.01706456 max resid 0.05525899 ## Run 18 stress 0.1177339 ## ... Procrustes: rmse 0.01706413 max resid 0.05525604 ## Run 19 stress 0.1177339 ## ... Procrustes: rmse 0.01706403 max resid 0.05525579 ## Run 20 stress 0.1174523 ## ... Procrustes: rmse 2.007087e-06 max resid 5.652376e-06 ## ... Similar to previous best## *** Best solution repeated 8 timesNMDS2 ## ## Call:## metaMDS(comm = dune, distance = "euclidean") ## ## global Multidimensional Scaling using monoMDS## ## Data: dune ## Distance: euclidean ## ## Dimensions: 2 ## Stress: 0.1174523 ## Stress type 1, weak ties## Best solution was repeated 8 times in 20 tries## The best solution was from try 7 (random start)## Scaling: centring, PC rotation ## Species: expanded scores based on 'dune'ordiplot(NMDS2, type ="t")
# hierarchical clustering: library(cluster) # clusteringlibrary(ape) # phylogenetic analyses, here to get pretty dendrogram# example for distance matrixdist(iris[1:3, 1:4]) # creates a distance matrix (comparison of all possible sample pairs), ## 1 2## 2 0.5385165 ## 3 0.5099020 0.3000000# default method ="euclidean", but can be changed# get Hierarchical Clusteringhc =hclust(dist(iris[, 1:4]))plot(hc)
# for colors use package apeplot(as.phylo(hc), tip.color =as.numeric(iris$Species)) # as.phylo converts object type "hcclust" into object type "phylo"
# change plotting type:plot(as.phylo(hc), tip.color =as.numeric(iris$Species), type ="fan")
In this exercise you will practice analyzing multivariate datasets using:
ordination methods which explore higher-order correlations among variables and similarities between observations
clustering methods which aggregate individual observations into clusters.
10.4.0.1 Useful functions for multivariate data analyses
pairs() - create pairplot (plot all variables against each other) prcomp() - calculate PCA plot(pca_object) - plot variance explained by PC axes biplot() - plot PCA metaMDS() - calculate NMDS ordiplot(nmds_object) - plot NMDS rda() - calculate RDA plot(rda_object) - plot RDA dist() - create distance matrix hclust() - perform hierarchical clustering; apply to distance matrix plot(as.phylo(hc_object)) - plot dendrogram from hclust output agnes() - another clustering algorithm kmeans() - perform kmeans non-hierarchical clustering; access cluster assignments using kmeans_object$cluster
10.4.0.2 Background on the dataset
Marine coastal ecosystems harbor the most productive and diverse communities on Earth. However, they are extremely vulnerable to climate change and human activities (such as landclaim, pollution, recreation purposes, …) As a consequence, the performance of this ecosystem has decreased considerably. To better understand these ecosystems, the dutch governmental institute (RIKZ: Rijksinstituut voor Kust en Zee) started a research project on sandy beaches to investigate the relationship between environmental factors and how they might affect benthic fauna. (from the LMU and Zuur, Ieno, Smith (2007), Chapter 12.8-12.10 and 27)
Question
We now want to use ordination methods to explore
the variability in environmental conditions
the variability in species richness
We also want to use clustering methods to define 3 environmental types and hierarchically cluster the samples with respect to their species richness.
Dataset
Read in the dataset as follows:
dat =read.table(file ="http://evol.bio.lmu.de/_statgen/Multivariate/11SS/RIKZGroups.txt", header = T)head(dat)
You already know the functions str() and summary() to get an overview of a dataset and to see which variables the dataset contains.
columns 2:5 is species richness within four larger faunal groups
columns 7:16 are different environmental variables
Not important for us are:
column 6 (week), a time stamp
column 17 (sorting1), variable from the observational design
Let’s get into the analysis!
10.4.1 Unconstrained ordination (PCA)
Carry out the following analyses:
Make two PCAs, one for the environmental and one for the species richness data (see columns above).
Name the results pca_environment and pca_species.
Create a biplot for each PCA.
Create a barplot for the proportion of variance explained by each component.
For example, the result for species richness should look like this:
You need the following functions and the package vegan:
prcomp()
biplot()
barplot()
summary()
library(vegan)
Hints
Don’t forget to scale the variables in the PCA.
In order to get the proportion of the explained variance, have a look at the summary of your analysis. str() shows you what the summary contains. You want to specifically look at the importance, here the second row contains the proportion of variance of all components. So what we want to plot is: summary(pca_species)$importance[2,].
In which multivariate dataset do the first and second components explain more variance?
What are the two samples that have the highest score on the respective third PCA axes?
From looking at the biplot: Which environmental variables would you choose, if you have only the resource to measure three variables?
NAP, angle 1 and chalk
salinity, penetrability and grain size
grain size, chalk and exposure
Describe the following correlations:
salinity and humus
angle1 and exposure
Polychaeta and Mollusca
Solution
## 1.) Conduct the principle component analyses (PCAs):# a) PCA of environmental datapca_environment =prcomp(dat[,7:16], scale = T)summary(pca_environment)## Importance of components:## PC1 PC2 PC3 PC4 PC5 PC6 PC7## Standard deviation 1.8476 1.3300 1.0919 1.0297 0.99521 0.7746 0.72784## Proportion of Variance 0.3414 0.1769 0.1192 0.1060 0.09904 0.0600 0.05298## Cumulative Proportion 0.3414 0.5182 0.6375 0.7435 0.84254 0.9025 0.95551## PC8 PC9 PC10## Standard deviation 0.52299 0.34727 0.22541## Proportion of Variance 0.02735 0.01206 0.00508## Cumulative Proportion 0.98286 0.99492 1.00000op <-par(mfrow =c(1,2))biplot(pca_environment, main ="PCA") # plot the results of the PCA as a rotation matrixbarplot(summary(pca_environment)$importance[2,], # get the importance measuremain ="Components of environment",ylab ="Proportion of variance explained")
# b) PCA of species richness datapca_species =prcomp(dat[,2:5], scale = T)summary(pca_species)## Importance of components:## PC1 PC2 PC3 PC4## Standard deviation 1.1177 1.0251 1.0095 0.8251## Proportion of Variance 0.3123 0.2627 0.2548 0.1702## Cumulative Proportion 0.3123 0.5750 0.8298 1.0000biplot(pca_species, main ="PCA") # plot the results of the PCA as a rotation matrixbarplot(summary(pca_species)$importance[2,], # get the importance measuremain ="Components of species",ylab ="Proportion of variance explained")
par(op)# From the *summary()* output we can see that the first and second components explain more variance in the species PCA.## 2.) What are the two samples that have the highest score on the third PCA axis?# -> order the samples by their PC3 coordinate:order(pca_environment$x[,'PC3'], decreasing = T) # -> 24 is highest## [1] 24 21 23 28 25 16 22 40 36 19 8 12 6 30 9 7 39 32 44 1 15 4 31 13 14## [26] 5 10 35 37 45 27 20 18 34 11 17 26 29 38 3 41 33 2 42 43order(pca_species$x[,'PC3'], decreasing = T) # -> 7 is highest## [1] 7 8 6 1 37 5 13 3 35 4 14 38 11 42 29 39 28 43 2 45 26 34 17 27 36## [26] 18 20 15 23 30 41 21 19 31 12 16 33 25 32 24 40 44 22 10 9## 3.) Which environmental variables would you choose, if you have only the resource to measure three variables? # From looking at the biplot, we choose 3 variables that describe a lot of variation (i.e. have a large length in the biplot) and have little collinearity. # -> For example, an appropriate choice would be salinity, penetrability and grain size.## 4.) We can get information on the correlations of variables by looking at their representation in the biplot: # a) salinity and humus: same direction -> positively correlated# b) angle1 and exposure: opposite directions -> negatively correlated# c) Polychaeta and Mollusca: almost perpendicular -> uncorrelated
10.4.2 Clustering
10.4.2.1 K-means
We want to use clustering methods to define 3 environmental types. Use the function kmeans() with centers = 3 (number of clusters to be generated). Remember to set a seed; the choice of the first k centers is random.
set.seed(467456)cl =#...
Compare the three clusters with the result of the PCA using the following code that uses another plotting framework called ggplot2:
The colors of the points represent the three clusters. Answer the following question on elearning-extern (Q 4-5):
How is it possible that four observations in the middle (in red - if you have used the same seed) belong to a different cluster than the observations around them (in black)?
Which environmental variables are on average high within the black cluster (cluster 1)?
Solution
# Create clusters using all environmental variablesset.seed(467456)cl <-kmeans(dat[,7:16], centers =3) # 3 clusters to be generated# Plotlibrary(ggfortify)autoplot(pca_environment, colour = cl$cluster, size =2, loadings =TRUE, loadings.colour ='black',loadings.label =TRUE, loadings.label.size =3, loadings.label.colour ="black")
# To understand why the four observations in the middle belong to a different cluster than the observations around them, we must take into account, that the biplot is only a 2-dimensional representation of a more-than-2 dimensional space. Therefore, the four points in the middle will be dissimilar to the points around them with respect to a variable that is not well represented by the first two PCA axes.# Environmental variables that are on average high within cluster 1:# -> looking at the plot we find that exposure and grain size are high on average within cluster 1
10.4.2.2 Hierarchical clustering
Now we want to hierarchically cluster the samples with respect to their species richness, as shown in the following plot:
Create the same plot using the functions:
hclust()
plot()
as.phylo()
Load the package ape. Have a look at the help for hclust() to read what the function does and look at the examples for further help on how to use the function. Then have a look at what the function as.phylo() does. Now, color the labels using the variable week. You can do this using the argument “tip.color =”.
Choose the correct statement(s) about the species richness and its sampling on elearning-extern(Q6). To be able to read the plot more easily, you can click Zoom in the top pane of the Plots window in RStudio.
