<table width = 900 border = 0><tr><td>
# <b>NMST539 | Cvičenie 6</b>

## <b>Metóda hlavných komponent | Aplikačná časť</b>
### <font color = "#880000"><b>LS 2020/2021 |  13/04/21 | (online výuka)</b></font>
###### <a href = "cvicenie2021_7.Rmd">zdrojový Rmd súbor (kódovanie UTF8)</a>

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

<b>Outline siedmeho cvičenia:</b> 

<ul>
<li>metóda hlavných komponent v programe R;</li>
<li>dôležté aplikačné aspekty PCA a príklady;</li>
</ul><br>


<b>Štatisticky program R</b> je k dispozícii (GNU public licence) na adrese <a href = "https://www.r-project.org">https://www.r-project.org</a><br><br>

<b>RStudio</b> -- "user-friendly" interface (jeden z mnohých, ktoré na internete nájdete): <a href = "https://www.rstudio.com">RStudio</a>.<br><br>

<b>Užitočné návody a manuály pre prácu s programom R:</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><br>



<b>Odborná literatúra:</b><br>
<ul>
<li>Hardle, W. and Simar.L.: Applied Multivariate Statistical Analysis. Springer, 2015</li>
<li>Mardia, K., Kent, J., and Bibby, J.:Multivariate Analysis, Academic Press, 1979.</li>
</ul>

<br><br>

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


### <font color = "#880000"><b>1. Metóda hlavných komponent v programe R (opakovanie)</b></font>

Metóda hlavných komponent je v programe R implementovaná pomocou príkazu `princomp()` 
ktorý využíva spektrálný rozklad (EIV) teoretickej (resp. empirickej) variančnej-kovariančnej matice. 
Pre konkrétne data $\mathcal{X} = (\boldsymbol{X}_1, \dots, \boldsymbol{X}_p)$ s odhadnutou variančnou-kovariančnou maticou 
$\mathcal{S}$ získame spektrálny rozklad -- t.j. $EIV(\mathcal{S})$ (analogicky, ak je známa teoretická matica, získame $EIV(\mathcal{\Sigma})$). 

Alternatívne a pri praktických (t.j. empirických) úlohach aj ekvivalentne je možné v programe R využiť aj príkaz `prcomp()`, 
ktorý funguje výhradne na empirickom princípe a využíva singulárný rozklad (SVD)
datovej matice $\mathcal{X}$ -- teda $SVD(\mathcal{X})$). 

Označme oba rozklady ako
$$
EIV(\mathcal{S}) = \Gamma \Lambda \Gamma^\top \quad \quad \textrm{a analogicky} \quad \quad SVD(\mathcal{X}) = UDV^\top.
$$

Ekvivalentnosť oboch rozkladov je viac-menej okamžitá: 
$$
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}),
$$ 
pre $\Gamma \equiv V$ a $\widetilde{\Lambda} \equiv D^2$. Taktiež platí, že
$$EIV(n\mathcal{S}) = \Gamma \widetilde{\Lambda} \Gamma^\top = \Gamma (n \Lambda) \Gamma^\top = n \Gamma \Lambda \Gamma^\top = n EIV(\mathcal{S}).$$

<br><br>




### <font color = "#880000"><b>2. R funkcie 'princomp()' a 'prcomp()'</b></font>

V programe R je k dispozícii množstvo rôznych knižníc (libraries) a funkcii určených pre aplikáciu metódy hlavných komponent. 
Napríklad funkcia `PCA()` v knižnici 'FactoMineR', alebo `pca()` v knižnici 'PCAtools', alebo funkcia `acp()` v knižnici 'amap'. 
V štandardnej inštalácii programu R (knižnica 'stats') sú defaultne k dispozícii dve základne funkcie:  `prcomp()` a `princomp()`. 

Funkcia `prcomp()` využíva singulárny rozklad (SVD rozklad) matice dat (viď teoretické odvodenie vyššie). Druhá funkcia funguje 
na princípe spektrálnej dekompozície teoretickej, resp. empirickej variančnej kovriančnej matice (EIV rozklad). 

Konkrétna implementácia jednotlivých príkazov je podrobnejšie špecifikované v helpe oboch funcií -- t.j. `?prcomp()`, alebo `?princomp()`.

V nasledujúcom R kóde využijeme ilustračné data a obe defaultne funkcie porovnáme. <br><br>

Ilustračné data reprezentuju 65 rôznych lokalít (t.j. počet pozorovaní, resp. riadkov v datovej matici) v Českej republike (jedná sa o geografické lokality blízko významných vodných tokov) a pre každú lokalitu je zaznamenaných 17 hodnôt rôznych biologických metrík a indexov 
hodnotiacich ekologickú a biologickú kvalitu danej lokality (t.j. 17 premenných, resp. stĺpcov v datovej matici). 


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

### v datach ponechame pouze numericke premenne
PCdata <- data[,-c(1,11,12,14)]

library(corrplot)
corrplot(cor(PCdata), method="ellipse")
```


Aplikácia PCA na data a porovnanie R funkcii `prcomp()` a `princomp()`:

```{r}
summary(prcomp(PCdata))
summary(princomp(PCdata))
```

Výsledky PCA nie sú invariantné voči základným transformáciam ako napr. štandardizácia. Porovnajte predchádzajúce výsledky 
s výsledkami po štandardizácii nižšie. Z tohto dôvodu sa metóda hlavných komponent doporučuje aplikovať už na štandardizované data. 

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

Porovnajte směrodatné chyby a odmocniny vlastných čísel spočítaných z výberovej variančnej-kovariančnej matice: 

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

Dokážete vysvetliť, prečo směrodatné chýby, ktoré sú výstupom oboch funkcií `prcomp()` a `princomp()` sú rozdielné? Ktoré z nich 
korespondujú s Vaším intuitívnym očakávaním? 
Prečo jednotlivé proporcie vysvetlenej variability vzájomne korespondujú? Využijte help k obom funkciám (napr. `help(princomp)`) 
a zamerajte sa na konkrétnu formu štandardizície, ktorá je aplikovaná pri počítaní výberovej variančnej-kovariančnej matice.
Aky faktor je použitý pri delení v oboch prípadoch? 

<br><br>

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

### <font color = "#880000"><b>3. Miera informácie obsiahnutej v hlavných komponentách</b></font>

Metóda hlavných komponent sa najčastejšie využíva k dvom podstatným úloham v štatistike: Buď sa jedná o exploratívnu 
metódu, pomocou ktorej sa snažime data vyzualizovať a pochopiť ich základnú štruktúru, alebo ide o metódu, pomocou 
ktorej chceme dosiahnuť podstatnú redukciu dimensionality (tzv. 'curse of dimensionality'). V druhom prípade sú dokonca na mieste 
rôzne štatistické postupy (napr. štatistické testy), ktorých cieľom je získanú redukciu dimenzie štatisticky vyhodnotiť. 

V prípade druhej varianty aplikácie PCA sa snažíme dosiahnúť čo najväčšiu redukciu dimenzie, ale zároveň zachovať čo najväčšiu mieru  celkovej informácie z pôvodných dat. Celkovú mieru informácie obsiahnutej v jednotlivých hlavných komponentách si budeme ilustrovať na nasledujúcom príklade s maticou 
dat, ktorá reprezentuje obrázok/fotografiu (t.j. vlastne tri matice dat - jedna matica pre každý farebných RGB kanál). 

<ul>
<li>Pre tieto účely je možné využiť akýkoľvek <i>jpg</i> súbor -- napríklad túto fotografiu:<br>
<a href = "http://msekce.karlin.mff.cuni.cz/~maciak/NMST539/foto.jpg">http://msekce.karlin.mff.cuni.cz/~maciak/NMST539/foto.jpg</a> </li>

<li>Následne je možné pomocou R knižnice `library("jpeg")` obrázok načítať do programu R: 

```{r}
library("jpeg")
foto <- readJPEG("foto.jpg")
```
</li>


<li>Obrázok (data) rozdelíme na tri matice pre každý farebný kanál samostatne a na každú datovú 
maticu (t.j. farebný kanál) aplikujeme metódu hlavných komponent: 

```{r, eval = F}
r <- foto[,,1]
g <- foto[,,2]
b <- foto[,,3]

dim(r)
dim(g)
dim(b)

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>Následne môžeme pomocou prvých niekoľko hlavných komponent obrázok spätne zrekonštruovať. Postupne použijeme 
3, 10, 20, 50, 100 a 500 prvých hlavných komponent  (z celkového počtu 682 premenných -- t.j. stĺpcov v jednotlivých datových maticiach `r`, `g` a `b`): 

```{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>Hlavné komponentty je pomocou rotácie možné vrátiť späť do priestoru pôvodných obrázkov 
a obrázky vzniknuté po (výraznej) redukcii dimenzionality uložiť a vizuálne porovnať s pôvodným obrázkom. 
Získame tak aspoň približný obraz jednak o celkovej miere redukcie dimenzionality, aj o celkovej miere zachovanej informáciez pôvodného obrázku: <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>Podobný postup kompresie obrázkov a porovnávania s originálom (napr. na základe metódy hlavných komponent) sa 
uplatňuje pri tzv. "image processingu". </li>

<li>V predchádzajúcom príklade sme automaticky využili predpoklad, že všetky tri farebné kanály sú rovnako dôležité. 
Toto ale nemusí byť pravda a niekedy je možno výhodnejšie aplikovať výraznejšiu redukciu dimenzionality na jeden kanál 
a ostatné kanály zatažiť kompresiou menej.  </li>
</ul>
<hr size="3" width="100%" noshade style="color:#440000" align="left" />
<br>




### <font color = "#880000"><b>4.Grafická vizualizácia PCA</b></font>

Základná informácia v súvislosti s metódou hlavných komponent je celková miera vysvetlenej 
variability, ktorú ponúka prvých $k \leq p \in \mathbb{N}$ hlavných komponent (kde $p$ značí celkový počet 
premenných v pôvodných datach). Grafickú predstavu môžeme získať napr. pomocou takýchto obrázkov: 

```{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)
```


Druhá informácia sa týka samotných hlavných komponent a ich usporiadania vzhľadom k pôvodným premenným. 
Takúto informáciu získame napr. pomocou jednoduchého biplotu -- t.j. funkcie `biplot()`: 

```{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)
```

<br><br>

Analogické grafy/obrázky (aj vo výrazne sofistikovanejších a komplexnejších prevedeniach) 
je možné získať pomocou mnohých iných príkazov a dodatočných R knižníc. Pre stručnosť aspoň niekoľko 
konkrétnnych odkazov:

<ul>
<li> <a href = "https://cran.r-project.org/web/packages/ggfortify/vignettes/plot_pca.html">pomocou knižníc 'ggfortify' a 'ggplot2'</a>;</li>

<li> <a href = "https://www.datacamp.com/community/tutorials/pca-analysis-r">pomocou knižníce 'ggbiplot'</a>;</li>

<li><a href = "https://cran.r-project.org/web/packages/geomorph/vignettes/geomorph.PCA.html">pomocou ďalších balíčkov a knižníc v programe R</a>;</li>

<li>pomocou ďalších nástrojov, ktoré určite nájdete na <a href = "www.google.com">google.com;</a></li>
</ul>




<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>Domáca (samostatná) úloha </b></font>
### (Deadline: Cvičenie č.8 / 20.04.2021)

<hr size="3" width="100%" noshade style="color:#440000" align="left" />
Vyberte si niektorú z vyššie uvedených možnosti pre vizualizáciu hlavných komponent a aplikujte 
metódu hlavných komponent na konkrétne data. <br><br>

<ul>
<li>Vyberte si vhodný datový súbor (napr. data, ktoré ste uvažovali v prvej, alebo štvrtej domácej úlohe); </li>

<li>Z podkladového datového súboru uvažujte vhodné (numerické) premenné. Dané premenné štandardizujte a aplikujte 
metódu hlavných komponent. </li>

<li>Hlavné komponenty a pôvodné premenné pomocou nástrojov v programe R vhodne vizualizujte. </li>

<li> Riešenie umiestnite na svoju webovú stránku, <font color = "red">najneskôr v utorok, 20.04.2021, do 14:00</font>.</li>
</ul>
<hr size="3" width="100%" noshade style="color:#440000" align="left" />
<br><br>



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