<table width = 900 border = 0><tr><td>
# <b>NMST539 | Lab Session 6</b>

## <b>Principal Component Analysis</b>
### (applications and examples in R)
### <b><font color = "#880000">Summer Term  2018/2019 | Monday 01/04/19</font></b>
###### <a href = "cvicenie2019_6.Rmd">Rmd file (UTF8 coding)</a>

**************

The R-software is available for download from the website: <a href = "https://www.r-project.org">https://www.r-project.org</a>

A user-friendly interface (one of many): <a href = "https://www.rstudio.com">RStudio</a>.<br><br>

<b>Manuals and introduction into <i>R</i>
(in Czech or English):</b><br>
    
<ul>
<li>Bína, V., Komárek, A. a Komárková, L.: <i>Jak na jazyk R.</i> (in Czech) | (<a href = "NMFM301/Rmanual2.pdf">PDF manual</a>)</li>
<li>Komárek, A.: <i>Základy práce s R.</i> (in Czech) | (<a href = "http://msekce.karlin.mff.cuni.cz/~komarek/vyuka/dataRko/Rmanual1.pdf">PDF manual</a>)</li>
<li>Kulich, M.: <i>Velmi stručný úvod do R.</i> (in Czech) | (<a href = "NMFM301/uvodr.pdf" >PDF manual</a>)</li>
<li>Venables, W.N., Smith, D.M., and the R Core Team: <i>An Introduction to R.</i> |  (<a href = "https://cran.r-project.org/doc/manuals/R-intro.pdf">PDF manual</a>)</li>
<li>De Vries, A. a Meys, J.: <i>R for Dummies.</i> (ISBN-13: 978-1119055808) |   (<a href = "http://www.ievbras.ru/ecostat/Kiril/R/Biblio/R_eng/R%20dummies.pdf">PDF book</a>)</li>		
</ul>

**************


The principal component analysis (PCA) - is a multivariate exploratory method based on the variance-covariance matrix of some random vector $\boldsymbol{X} \in \mathbb{R}^{p}$. Of course, instead of working with the vector itself and the theoretical variance-convariance matrix  we use a random sample drawn from the same distribution and the sample verion of the varinance covariance matrix. 

The main idea is to consider $\ell <  p$ principal components instead of the whole $p \in \mathbb{N}$ dimensional random vector  (dimensionality reduction). The principal components are orthogonal linear combinations of the original covariates, and, in addtion, they are ordered with respect to a decreasing variance. In most cases $\ell \ll p$ (a significant reduction of dimensionality). <br><br>


In general, the principal components can be seen from (four) different perspectives:

<ul>
<li>orthogonal directions within the data with the highest proportional variance;</li>
<li>best linear subspace approximation to the data;</li>
<li>eigen value decomposition of the (sample) variance-covariance matrix;</li>
<li>singular value decomposition of the data frame;</li>
</ul>

Statistically, the principal components are linear combinations of the original covariates, such that these new covariates are mutually uncorrelated, with zero mean
and the variace-covariance matrix which is diagonal, with the eigen values of $\Sigma$ on its diagonal. The matrix $\Sigma$ is the variance-covariace matrix of the original data. 

<br><br>

### <font color = "#880000"><b>1. Principal Components and Regression</b></font>

Let $\boldsymbol{Y}$ be some random vector of the length $n \in \mathbb{N}$ and, in addition, we also have some data matrix $\mathbb{X}$ of the type $n \times p$, where, in general, $p \in \mathbb{N}$ is considered to be less than $n$.  A common task in regression is to use the information obtained in $\mathbb{X}$ and to approximate the $n$ dimensional random vector $\boldsymbol{Y}$ in a lower dimensional space $\mathcal{L}(\mathbb{X})$ -- a $p$ dimensional space generated by the columns of the matrix $\mathbb{X}$. 

If we denote the orthogonal projection of $\boldsymbol{Y}$ into $\mathcal{L}(\mathbb{X})$ as $P_{\mathbb{X}}(\boldsymbol{Y})$ then the linear regression model for regessing $\boldsymbol{Y}$
on $\mathbb{X}$ in the form  $P_{\mathbb{X}}(\boldsymbol{Y}) = a_1 \boldsymbol{X}_1 + \dots + a_p \boldsymbol{X}_p$ gives the minumum error
$$
\| \boldsymbol{Y} - P_{\mathbb{X}}(\boldsymbol{Y}) \|_2^2.
$$

The corresponding projection matrix (which needs to be idempotent) is 
$$
\mathbb{H} = \mathbb{X}(\mathbb{X}^\top \mathbb{X})^{-1}\mathbb{X}^\top = \mathbb{Q}\mathbb{Q}^\top,
$$
where $\mathbb{X}$ contains the basis vectors of $\mathcal{L}(\mathbb{X})$ in its columnts, or, respectively, the colums of $\mathbb{Q}$ contain the orthogonal basis of the same linear space.  Obviously, it also holds, that 
$$
P_{\mathbb{X}}(\boldsymbol{Y}) = \mathbb{H}\boldsymbol{Y} = \mathbb{X}(\mathbb{X}^\top \mathbb{X})^{-1}\mathbb{X}^\top \boldsymbol{Y}.
$$


This concept can be easily generalized for any orthogonal projection of $\boldsymbol{Y}$ into a linear span of the columns of some general $n \times p$ matrix $\mathbb{B}$, such that $\mathbb{B}^\top\mathbb{B} = \mathbb{I}$. Let this projection be denoted as $P_{\mathbb{B}}(\boldsymbol{Y})$. Then it holds that $P_{\mathbb{B}}(\boldsymbol{Y}) = \mathbb{B}\mathbb{B}^\top\boldsymbol{Y}$ and, moreover, it can be shown that 
$$
\|\boldsymbol{Y} - P_{\mathbb{\Gamma}}(\boldsymbol{Y})\|_2^2 \leq \|\boldsymbol{Y} - P_\mathbb{B}(\boldsymbol{Y})\|_2^2
$$
where $\Gamma$ is the $n\times p$ matrix with the first eigen vectors of the variance covariance matrix $\Sigma$ (see Eckart and Young (1936) and Mirsky (1960) approximation). 


```{r, fig.width = 10, fig.height=4.5}
set.seed(1234)
x <- seq(0,1, length = 100)
y <- 0 + 2 * x + rnorm(100)

par(mfrow = c(1,3))
### regression of Y on X
plot(y ~ x, pch = 21, bg = "lightblue", xlab = "X variable", ylab = "Y variable", main = "Projection of Y onto L(X)", ylim = c(-1.5, 4))
m <- lm(y ~ x)
abline(m, col = "red", lwd = 2)

index <- c(10, 60, 80)

for (i in 1:100){
  lines(c(x[i], x[i]), c(y[i], m$fitted[i]), col = "blue", lty =2)
}


### regression of X on Y 
plot(y ~ x, pch = 21, bg = "lightblue", xlab = "X variable", ylab = "Y variable", main = "Projection of X onto L(Y)", ylim = c(-1.5, 4))
m <- lm(x ~ y)
lines((x - m$coeff[1])/m$coeff[2] ~ x, col = "red", lwd = 2)

index <- c(10, 60, 80)

for (i in 1:100){
  lines(c(x[i], m$coeff[1] + m$coeff[2] * y[i]), c(y[i], y[i]), col = "blue", lty =2)
}


### regression onto the first PC

### data standardization
means <- apply(cbind(x,y), 2, mean)
sdErs <- apply(cbind(x,y), 2, sd)
cd <- cbind((x - means[1])/sdErs[1], (y - means[2])/sdErs[2])

### first eigen vector
firstEigenVector <- eigen(cov(cd))$vectors[,1]

### projection onto the linear space generated by the first eigen vector
project2pc1 <- t(firstEigenVector) %*% t(cd)

### rotation back to the original coordinates
rotation <- as.numeric(project2pc1) %*% t(firstEigenVector)

### rescaling to the original mean/variance
rescaled <- cbind(rotation[,1] * sdErs[1] + means[1], rotation[,2] * sdErs[2] + means[2])

plot(y ~ x, pch = 21, bg = "lightblue", xlab = "X variable", ylab = "Y variable", main = "Projection onto 1.PC", ylim = c(-1.5,4), xlim = c(-0, 1))

### first PC
lines(rescaled[,2] ~ rescaled[,1], col = "red", lwd = 2)

### orthogonal projections (but they are not orthogonal)
for (i in 1:100){
  lines(c(x[i], rescaled[i,1]), c(y[i], rescaled[i,2]), col = "blue", lty =2)
}
```

Alternatively, instead of finding the spectral decomposition of the variance-covariance matrix $\Sigma$ we can use an idea of a standard linear regressin approach 
and we can alternate between two different regression problems. For some data $\boldsymbol{Y}_1, \dots, \boldsymbol{Y_n}$ we have
$$
\sum_{i = 1}^n \|\boldsymbol{Y}_i - P_\mathbb{B}(\boldsymbol{Y}_i)\|_2^2 = \sum_{i = 1}^n \|\boldsymbol{Y}_i - \mathbb{B}\mathbb{B}^\top\boldsymbol{Y}_i\|_2^2 =  \sum_{i = 1}^n \|\boldsymbol{Y}_i - \mathbb{B} \boldsymbol{v}_i\|_2^2 = \sum_{i = 1}^{n}\sum_{j = 1}^{p} |Y_{i j} - \boldsymbol{b}_j\boldsymbol{v}_{i}|,  
$$
where $\boldsymbol{v}_i = \mathbb{B}^\top\boldsymbol{Y}_i$, for $i = 1, \dots, n$, and $\boldsymbol{b}_j$ is the corresponding row of $\mathbb{B}$. However, neither $\mathbb{B}$ nor the vectors $\boldsymbol{v}_i$ are known. Therefore, we need to minimize the term on the right hand side with respect to some orthogonal matrix $\mathbb{B}$, such that $\mathbb{B}^\top\mathbb{B} = \mathbb{I}$ and some vectors $\boldsymbol{v}_1, \dots, \boldsymbol{v}_n$. This can be easily done by an alternating regression approach. <br><br>

For a brief illustration see the next example (by <a href = "https://www.stat.ubc.ca/~matias/">Matías Salibián-Barrera</a>). 
```{r}
alter.pca.k1 <- function(x, max.it = 500, eps=1e-10) {
  n2 <- function(a) sum(a^2)
  p <- dim(x)[2]
  x <- scale(x, scale=FALSE)
  it <- 0
  old.a <- c(1, rep(0, p-1))
  err <- 10*eps
  while( ((it <- it + 1) < max.it) & (abs(err) > eps) ) {
    b <- as.vector( x %*% old.a ) / n2(old.a)
    a <- as.vector( t(x) %*% b ) / n2(b)
    a <- a / sqrt(n2(a))
    err <- sqrt(n2(a - old.a))
    old.a <- a
  }
  conv <- (it < max.it)
  return(list(a=a, b=b, conv=conv))
}


set.seed(1234)
n <- 20
p <- 5
x <- matrix(rt(n*p, df=2), n, p)
```

And we can obtain the first principal component by adopting the alternating regression approach: 

```{r}
alter.pca.k1(x)$a
```

which can be directly compared with the eigen vector itself (obtained by the singular value decomposition or, equivalently, by the eigen value decomposition instead):

```{r}
svd(cov(x))$u[,1]
eigen(cov(x))$vectors[,1]
```


<br><br>

For some practical applications it can be useful to investigate the empirical performance of both approaches. This can be effectively done, for instance, by comparing
the computational times needed for both approaches (for this purpose we consider a higher dimensional example).  

```{r}
n <- 2000
p <- 500
x <- matrix(rt(n*p, df=2), n, p)

system.time( tmp <- alter.pca.k1(x) )

system.time( e1 <- svd(cov(x))$u[,1] )
```

<br><br>
The principal component analysis performed by the R command `princomp()` uses the eigen value decomposition of the sample variance-covariance matrix $\mathcal{S}$ of 
some data $\mathcal{X} = (\boldsymbol{X}_1, \dots, \boldsymbol{X}_p)$  (denoted as $EIV(\mathcal{S})$) while the second command,  `prcomp()`, uses the singular value decomposition 
of the data matrix $\mathcal{X}$ instead (denoted as $SVD(\mathcal{X})$). 

Let the eigen value decomposition (the singular value decomposition respectively) is of the form 
$$
EIV(\mathcal{S}) = \Gamma \Lambda \Gamma^\top \quad \quad \textrm{respectively} \quad \quad SVD(\mathcal{X}) = UDV^\top.
$$

Then we easily obtain, that both methods are equivalent as we have
$$
n \mathcal{S} = \mathcal{X}^\top\mathcal{X} =  VDU^\top UD V^\top = V D^2 V^\top = \Gamma \widetilde{\Lambda} \Gamma^\top = EIV(n\mathcal{S}),
$$ 
for $\Gamma \equiv V$ and $\widetilde{\Lambda} \equiv D^2$. Moreover, it holds that 
$$EIV(n\mathcal{S}) = \Gamma \widetilde{\Lambda} \Gamma^\top = \Gamma (n \Lambda) \Gamma^\top = n \Gamma \Lambda \Gamma^\top = n EIV(\mathcal{S}).$$
Thus, we we multiply the square matrix by some constant, then the eigen vectors remain the same, just the eigen values multiply correspondingly. 
<br><br>

#### <font color = "#880000">Try by yourself</font> 
<hr size="3" width="100%" noshade style="color:#440000" align="left" />
The principal component analysis relies on the eigen vector decomposition of the sample variance-covariance matrix. Alternatively, one can use a singular value decomposition of the 
data matrix itself. However, in both situations it is crucial to understand the role of the eigen vectors and the eigen values. See the following examples to get some useful insight. 

<ul>
<li>For a graphical illustration of the eigen values and eigen vectors and their role in the matrix calculations see:<br>
<a href = "http://setosa.io/ev/eigenvectors-and-eigenvalues/">http://setosa.io/ev/eigenvectors-and-eigenvalues/</a></li>

<li>For an interactive illustration of 2D and 3D situations of of the PC analysis, see the following example:<br>
<a href = "http://setosa.io/ev/principal-component-analysis/">http://setosa.io/ev/principal-component-analysis/</a></li>
</ul>
<hr size="3" width="100%" noshade style="color:#440000" align="left" />
<br><br>

### <font color = "#880000"><b>2. Biometrics on the river localities in the Czech Republic</b></font>

We will start with a data sample which represents 65 river localities in the Czech Republic. For each locality there are 17 biological metric recorded. Each biological metric is used as a qualitative/quantitative status of some biological activity (for instance, phytobenthos, macrobenthos, fishes) within each locality expressed in terms of some well defined and internationally recognized metrics - indexes). The idea was to compare the quality of the fresh water environments across all European union member states. 

The data file can be obtained from the following address: 

```{r, fig.width = 10, message=F}
data <- read.csv("http://msekce.karlin.mff.cuni.cz/~maciak/NMST539/bioData.csv", header = T)
attach(data)
```



One can use a fancy way to visualize the correlation structure within the data (using the library 'corrplot' which needs to be installed by `install.packages("corrplot")` firstly). 

```{r, fig.width = 10, message=F}
library(corrplot)
PCdata <- data[,-c(1,11,12,14)]
corrplot(cor(PCdata), method="ellipse")
```

Another nice alternative can be obtained using the `qgraph()` command from the library 'qgraph' (again, the library needs to be installed by the command `install.packages("qgraph"")`). 

```{r, fig.width = 10}
#library("qgraph")
#qgraph(cor(PCdata))
```
 

<br><br>

### <font color = "#880000"><b>3. Two R functions:  princomp() vs. prcomp()</b></font>

There are various libraries PCA available for the R environment: standard commands are `prcomp()` and `princomp()`, which are both available under the classical R installation (library 'stats'). While the first one is using the singular value decomposition of the data matrix (SVD) the second one uses the spectral decomposition of the variance-covariance (correlation) matrix  instead. For additional details see the help sessions for both functions (type `?prcomp()` and `?princomp()`).

We can compare two results obtained by the two commands mentioned above: 

```{r}
PCdata <- data[,-c(1,11,12,14)]
summary(prcomp(PCdata))
summary(princomp(PCdata))
```

The results are, however, not invariant with respect to transformations. 

```{r}
PCdata <- scale(PCdata, center = TRUE, scale = TRUE)
summary(prcomp(PCdata))
summary(princomp(PCdata))
```

Compare the standard deviation values with the square roots of the eigen values of the sample variance covariance matrix. 

```{r}
sqrt(eigen(var(PCdata))$values)
```


Can you explain why the Standard deviation values provided for two procedures (`prcomp()` and `princomp()`) are different? Why the proportions 
of variance and the cumulative proportions both corresponds? Check the help `help(princomp)` and focus on the standardization form when calculating the sample variance-covariance matrix: what is the divisor factor used in both fucntions? <br><br>


#### <font color = "#880000">Try by yourself</font> 
<hr size="3" width="100%" noshade style="color:#440000" align="left" />
A simple (visual example) on how the PC analysis actually works can be obtained from the next illustration. 
<ul>
<li>Firstly, install the library `library("jpeg")` and load it by 
```{r}
library("jpeg")
```
</li>
<li>Download a sample image from the following address: <a href = "http://msekce.karlin.mff.cuni.cz/~maciak/NMST539/foto.jpg">http://msekce.karlin.mff.cuni.cz/~maciak/NMST539/foto.jpg</a> and save it to your working environment; </li>
<li>Load the picture into the R software by using the followinng command (with the right working directory):

```{r}
foto <- readJPEG("foto.jpg")
```
</li>
<li>Decompose the picture into RGB scale and apply the PC analysis to all three:
```{r, eval = F}
r <- foto[,,1]
g <- foto[,,2]
b <- foto[,,3]

foto.r.pca <- prcomp(r, center = FALSE)
foto.g.pca <- prcomp(g, center = FALSE)
foto.b.pca <- prcomp(b, center = FALSE)

rgb.pca <- list(foto.r.pca, foto.g.pca, foto.b.pca)
```
</li>

<li>Now, we will reconstruct the original picture using 3, 10, 20, 50, 100 and 500 components (out of 682 covariates together).
```{r, eval =F}
for (i in c(3,10,20,50,100,500)) {
  pca.img <- sapply(rgb.pca, function(j) {
    compressed.img <- j$x[,1:i] %*% t(j$rotation[,1:i])
  }, simplify = 'array')
  
  writeJPEG(pca.img, paste("foto_", i, "_components.jpg", sep = ""))
}
```
</li>

<li>Compare the compressed images with the original photo and see how much dimensionality reduction (loss of information in each compressed fiture respectively) is present in each figure;<br><br>
<center>
<table width = "900" border = 0>
<tr><td width = "300"><img src = "foto_3_components.jpg", width = "290", align = "center">
    <td width = "300"><img src = "foto_10_components.jpg", width = "290", align = "center">
    <td width = "300"><img src = "foto_20_components.jpg", width = "290", align = "center"></tr>
</table><br>

<table width = "900" border = 0> 
<tr><td width = "300"><img src = "foto_50_components.jpg", width = "290", align = "center">
     <td width = "300"><img src = "foto_100_components.jpg", width = "290", align = "center">
    <td width = "300"><img src = "foto_500_components.jpg", width = "290", align = "center"></tr>
</table><br>

<table>
<tr><td><img src = "foto.jpg", width = "900"></tr>
</table><br>

</center>
</li>

<li>The image compression (a lower dimensional representation with respect to various principal components) can be now used to compare the original image with some other one in order to asses some mutual (important) similarities (image processing). 

<center>
<table>
<tr><td><img src = "foto_3_components.jpg", width = "290">
    <td><img src = "foto_10_components.jpg", width = "290">
    <td><img src = "foto_20_components.jpg", width = "290"></tr>
<tr><td><img src = "foto_50_components.jpg", width = "290">
     <td><img src = "foto_100_components.jpg", width = "290">
    <td><img src = "foto_500_components.jpg", width = "290"></tr>
</table><br>

<table>
<tr><td><img src = "foto2_3_components.jpg", width = "290">
    <td><img src = "foto2_10_components.jpg", width = "290">
    <td><img src = "foto2_20_components.jpg", width = "290"></tr>
<tr><td><img src = "foto2_50_components.jpg", width = "290">
     <td><img src = "foto2_100_components.jpg", width = "290">
    <td><img src = "foto2_500_components.jpg", width = "290"></tr>
</table><br>

</center>

</li>

<li>The PCs obtained above used the assumption that all three color channels as same important. On the other hand, this may not be true, therefore, it is also possible to combine all
three color channels into just one data matrix and to run the PC on the overal image. </li>
</ul>
<hr size="3" width="100%" noshade style="color:#440000" align="left" />
<br>








### <font color = "#880000"><b>4. Various graphical representations of PCs</b></font>

First of all, it is usefull to know the amount of variability (the proportion of the total variability)
which is explained by some first principal components. 

```{r, fig.width = 10}
pc <- princomp(PCdata)
par(mfrow = c(1,2))
plot(pc, type = "l", col = "red", main = "")
plot(pc, col = "red", main = "")
summary(pc)
```

<br><br>

#### <font color = "#880000">Try by yourself</font>
<hr size="3" width="100%" noshade style="color:#440000" align="left" />
<ul>
<li>Use the output above and plot a figure where one can clearly see the proportions of explained 
variability for each component. </li>
<li>Similarly, from the output above, create a figure with cumulative proportions of the explaied variability. </li>
<li>Apply a similar approach (PCA) to the iris dataset available in the R environment (use `data(iris)`).</li>
</ul>
<hr size="3" width="100%" noshade style="color:#440000" align="left" />
<br>


There are also other options to visualize the results. A standard graphical tool is available through the commnad `biplot()` which is available in the standard R installation. This fuction automatically projects the data on the first two PCs. Other principle components can be also used in the figure (parameter 'choices = 1:2').


```{r, fig.width = 10}
par(mfrow = c(1,2))
biplot(pc, scale = 0.5, choices = 1:2, cex = 0.8)
biplot(pc, scale = 0.5, choices = 1:2, xlim = c(-1,1), ylim = c(-1, 1), cex = 0.8)
```

A more fancy version can by obtained by installing the `ggbiplot` package from the Github repository (use  command `install_github("ggbiplot", "vqv")`). After loading in the library one can use the `ggbiplot()` command (and additional parameters).  

```{r, fig.width = 10, warning=F, message=F}
library("ggbiplot")
plot1 <- ggbiplot(princomp(PCdata), elipse = T, groups = Lital, circle = T)
plot2 <- ggbiplot(princomp(PCdata), elipse = T, groups = Bindex, circle = T)

grid.newpage()
vplayout <- function(x, y) viewport(layout.pos.row = x, layout.pos.col = y)
pushViewport(viewport(layout = grid.layout(1, 2)))
print(plot1, vp = vplayout(1, 1))
print(plot2, vp = vplayout(1, 2))
```

<br><br>


Another type of plot can be obtained using the library 'ggfortify' (use `install.packages('ggfortify', dep = T)` to install the library) and the command `autoplot()`. 


```{r, warning=F, message=F, fig.width = 10}
#library("ggfortify")
#autoplot(prcomp(PCdata), data = data,  colour = 'SaprInd', loadings = T)
```

<br><br>

Last but not least, there is also a complex command `fviz_pca()` available in the library 'factoextra' (use `install.packages("factoextra")` for installation). 

```{r, fig.width = 10, message=F}
#library("factoextra")
#par(mfrow = c(1,3))
#fviz_pca_ind(pc)
#fviz_pca_var(pc)
#fviz_pca_biplot(pc)
```


<br><br>

#### <font color = "#880000">Try by yourself </font>
<hr size="3" width="100%" noshade style="color:#440000" align="left" />
There are many more packages and libraries available for the principal component analysis in R (especially for dealing with all kinds of psychological questionaires, etc.). 
You can use google to find some others not mentioned in this text. <br><br>
<ul>
<li>Use some two dimensional artificial datasets with small amount of observations (e.g. 3 - 5)
and perform the PC analysis. Using some appropriate visualization techniques try to understand 
the algebra machinery behind.</li>
</ul>
<hr size="3" width="100%" noshade style="color:#440000" align="left" />
<br>




<br><br>
<hr size="3" width="100%" noshade style="color:#440000" align="left" />
<hr size="3" width="100%" noshade style="color:#440000" align="left" /><br><br> 




## <font color = "#880000"><b>Homework Assignment</b></font>
### (Deadline: seventh lab session / 08.04.2019)

Use the command below to instal the `SMSdata` packages with various multivariate datasets: 

```{r, eval = F}
install.packages(pkgs="http://www.karlin.mff.cuni.cz/~hlavka/sms2/SMSdata_1.0.tar.gz", repos=NULL, type="source")
```

Chose one dataset of your choise from the list of all available datasets in the package:

```{r, eval = F}
data(package = "SMSdata") 
```
 
 There are 21 different datasets and you can load each of them by typing and running 
`data(dataset_name)`, for instance:

```{r, eval = F}
data(plasma) 
```

<ul>
<li>Choose one dataset and apply the principal component analysis; </li>

<li>Discuss the right choice of the number of principal components. Provide some graphical output and interpret the results;</li>
</ul>


</td></tr></table>
<br><br><br>