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

## <b>Principal Component Analysis</b>
###  (application in R)

### <font color = "#880000">LS 2017 | Tuesday 04/04/2017 </font>
###### <a href = "cvicenie2017_5.Rmd">Rmd file (UTF8 encoding)</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>  (<a href = "NMFM301/Rmanual2.pdf">PDF súbor</a>)</li>
<li>Komárek, A.: <i>Základy práce s R.</i>  (<a href = "http://msekce.karlin.mff.cuni.cz/~komarek/vyuka/dataRko/Rmanual1.pdf">PDF súbor</a>)</li>
<li>Kulich, M.: <i>Velmi stručný úvod do R.</i> (<a href = "NMFM301/uvodr.pdf" >PDF súbor</a>)</li>
<li>De Vries, A. a Meys, J.: <i>R for Dummies.</i> (ISBN-13: 978-1119055808) </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 \leq p$ principal components instead of the whole $p \in \mathbb{N}$ dimensional random vector  (dimensionality reduction). The principel 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>

<ul>
<li>For a graphical illustration of the eigen values and eigen vectors and their role in the matrix calculations see the following example:<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>
<br><br>

### <font color = "#880000"><b>Biometrics on some river localities (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 recordid (the status of the biological life within each locality expressed in terms of some well defined and internationally recognized metrics - indexes).

The data file can be obtained from the following address: 

```{r, fig.width = 10}
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}
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)`). However, it does not work with the most recent R version (v.3.2.4). 
 

<br><br>

### <font color = "#880000"><b>Two R functions: prcomp() vs. princomp()</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)
```
<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 immage 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>
<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 = "foto.jpg", width = "900"></tr>
</table>

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


### <font color = "#880000">Various graphical representations (visualization)</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}
biplot(pc, scale = 0.5, choices = 1:2, cex = 0.8)
```

Or one can bring it all to a more detailed version by 

```{r, fig.width = 10, warning=F, message=F}
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")
ggbiplot(princomp(PCdata), elipse = T, groups = Ntaxon, circle = T)
```


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 = 'Ntaxon', loadings = T)
```


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

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

#### <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. 
You can use google to find some others not mentioned in this text. 
<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: sixth lab session / 11.04.2007)

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>