= read.table(file = "http://evol.bio.lmu.de/_statgen/Multivariate/11SS/RIKZGroups.txt", header = T)
dat head(dat)
Exercise - Multivariate statistics
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.
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
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:
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!
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
andpca_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 theimportance
, here the second row contains the proportion of variance of all components. So what we want to plot is:summary(pca_species)$importance[2,]
. - Don’t forget to give your plots a title.
Now, use your results to answer the questions on elearning-extern (Q 1-3) (“07_Test for Exercise in R”).
- 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
Clustering
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:
library(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")
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)?
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.