Exercises

This page contains tutorials/exercises to deepen your understanding of PCA. I would recommend starting with the challenge in Exercise 1 (can be discussed within your group) but feel free to choose the exercise which is the most interesting for you (you will probably not have time during the course to do all exercises). A description is provided as introduction of each exercise.

Exercise 1: challenge

Exercise 1 is based on a real question I once got from a collaborator (RNAseq data) where simply performing a PCA was sufficient to answer their question.

Let's set the exercise as follows:

A collaborator comes to you (because you just took a course in biostats) and says their data analysis results do not make sense
They have done an experiment on drosophila with 3 conditions, 3 replicates per condition, 2 types of RNA. The 3 conditions correspond to a control condition, a condition where gene A has been suppressed, a condition where gene B has been suppressed. The 2 types of RNA (which do not matter for this exercise) are cytosolic RNA and Ribosome-Protected-Fragments RNA.
You ran a PCA and see the following representation of the samples (different conditions are represented as different colors, the RNA types are represented as different shapes and the numbers next to each sample dot is the replicate identifier (1, 2, 3).

Can you guess what can be the reason why their results do not make sense? Can you interpret PC1, PC2?
If needed, check this hint.

Exercise 2

This exercise/tutorial is there in case you want to go through the example from the PCA lecture (2004 Olympic games decathlon) and run the R code at your own pace, testing different parameters etc. The provided R file contain the entire lecture with the actual R code used and where slide content are provided as comments (Slide titles are shown after ##). In order for the code from the lecture to run, you will first need to install the following R packages: summarytools, dplyr, corrplot, ggbiplot, tidyr, tibble, forcats, GGally, RColorBrewer, knitr.

If you prefer to use base R which is sufficient for performing PCA (and avoid having to install the packages mentioned above), you can use the code provided below. To get the data, you can either download it from Canvas (option 1; data) or install the R package in which it is saved and pre-process it yourself (option 2). The code below shows both options to import the data and shows how to compute the PCA.


## To get the data, you can either download it from Canvas (option 1) or install
## the R package in which it is saved and pre-process it yourself (option 2)

### --- OPTION 1 --- ###

### Download the data from Canvas

### Read the data into R
decathlon_OG_new <- read.csv("decathlon_OG_new.csv", 
                             row.names = 1) # to specify that the first column of the csv file contain rownames

### Check the data (you might notice that variables with a name starting by a 
### number have been automatically modified but this will not have any impact on
### the analysis. You can continue with running the PCA)
decathlon_OG_new

### run PCA with scale.=TRUE
pca <- prcomp(decathlon_OG_new, scale. = TRUE)

### plot the results
biplot(pca)




### --- OPTION 2 --- ###

install.packages("ClustOfVar")
library(ClustOfVar)

### Load data
data("decathlon")

### Keep Olympic games results only (dataset also contains another competition)
decathlon_OG <- subset(decathlon, Competition == "OlympicG")

### Keep only the 10 sport variables
decathlon_OG <- decathlon_OG[, !colnames(decathlon_OG) %in% c("Rank", "Points", "Competition")]

### Reformatting the data so that the higher the score in each sport, the better
vec <- c(-1, 1, 1, 1, -1, -1, 1, 1, 1, -1)
decathlon_OG_new <- as.matrix(decathlon_OG) %*% diag(vec)
colnames(decathlon_OG_new) <- colnames(decathlon_OG)

### PCA
pca <- prcomp(decathlon_OG_new, scale. = TRUE)
biplot(pca)

Exercise 3

Part of this exercise is a guide through some manual mathematical steps of PCA. These steps (e., f., g., i., j., k.) are optional so if you prefer to focus on running and interpreting the PCA, feel free to skip them. Note that step h. is not optional.

Let’s start principal component analysis (PCA) with Edgar Anderson’s iris data, which is already available as iris in R.

Examine the contents and data types of the data set iris. R functions str, summary and pairs can be useful for this.
As Setosa is somewhat different from the other two species, we will limit our analysis to these two. Leave the data of Versicolor and Virginica only, and examine the contents of the data.
The first 4 columns contains 4 different variables and the last column has species. Split the data into a numeric matrix m and a character vector species.
Transform each variable to have the same mean 0 and the same standard deviation 1 using scale.
(optional) Compute eigenvectors and eigenvalues of the variance-covariance matrix of the transformed (or scaled) m using cov and eigen.
(optional) Compute the scores of the first principal component. Matrix dot product in R is %*% as A %*% B for \(\mathbf{A} \mathbf{B}\). Hint: The scores of the principal components can be computed from taking the matrix dot product of the scaled data by the eigenvectors of the covariance matrix.
(optional) Compute the proportion of variance explained by each PC and find the rotation matrix and loading matrix.
Using prcomp, perform PCA. Don’t forget to add scale. = TRUE in the arguments of the function. The default is scale. = FALSE though the help of the function states scaling is recommended.
(optional) Compare the scores of PC1 from our manual computation (f.) and from prcomp.
(optional) Compute the proportion of explained variance using the output from prcomp. Compare the results with that of g).
(optional) Find the rotation matrix or loading matrix from the output of prcomp and compare it with that of g).
Draw a PCA score plot for PC1 and PC2. Add colors by species. What does the plot show?
Draw a PCA score plot for PC2 and PC3.
Draw a scree plot.
Draw a bi-plot. Can you tell
1. which variables are the most important for PC1?
2. which variables are highly correlated?

Exercise 4

This exercise provides an additional example of how to perform PCA, this time on gene expression data from cell lines.

NCI60 data set consists of gene expression values for 6830 genes for 64 cell lines. Using this data set, repeat the analysis using prcomp of Exercise 1. You can download the data using the command below.

m <- t(as.matrix(read.table(
  url("https://web.stanford.edu/~hastie/ElemStatLearn/datasets/nci.data.csv"),
  sep = ",", row.names = 1, header = TRUE
)))

Examine the contents and data types of the data set nci. Check the dimension of the data first. This is much larger than iris.
Using prcomp, perform PCA.
Compute the proportion of variance explained by each PC.
Find the rotation matrix. Show the first 5 columns and rows only.
Draw a PCA score plot for PC1 and PC2.
Draw a PCA score plot for PC2 and PC3.
Draw a scree plot.
Draw a bi-plot. Note it can take long time.

You may find the loadings displayed in the biplot are in different scale from those you computed. It is primarily due to the scaling by function biplot. Please refer to https://stats.stackexchange.com/questions/276645/arrows-of-underlying-variables-in-pca-biplot-in-r for the question.

Solutions of exercises 3 and 4