NMST539 | Cvičenie 6
Metóda hlavných komponent | Aplikačná časť
LS 2020/2021 | 13/04/21 | (online výuka)
Outline siedmeho cvičenia:
-
metóda hlavných komponent v programe R;
-
dôležté aplikačné aspekty PCA a príklady;
Štatisticky program R je k dispozícii (GNU public licence) na adrese https://www.r-project.org
RStudio – “user-friendly” interface (jeden z mnohých, ktoré na internete nájdete): RStudio.
Užitočné návody a manuály pre prácu s programom R:
-
Bína, V., Komárek, A. a Komárková, L.: Jak na jazyk R. (PDF súbor)
-
Komárek, A.: Základy práce s R. (PDF súbor)
-
Kulich, M.: Velmi stručný úvod do R. (PDF súbor)
-
De Vries, A. a Meys, J.: R for Dummies. (ISBN-13: 978-1119055808)
Odborná literatúra:
-
Hardle, W. and Simar.L.: Applied Multivariate Statistical Analysis. Springer, 2015
-
Mardia, K., Kent, J., and Bibby, J.:Multivariate Analysis, Academic Press, 1979.
1. Metóda hlavných komponent v programe R (opakovanie)
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}).\]
2. R funkcie ‘princomp()’ a ‘prcomp()’
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.
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).
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() :
summary(prcomp(PCdata))
## Importance of components:
## PC1 PC2 PC3 PC4 PC5 PC6 PC7
## Standard deviation 23.0661 8.1033 4.50909 3.42042 2.17783 1.82627 1.65944
## Proportion of Variance 0.8273 0.1021 0.03162 0.01819 0.00738 0.00519 0.00428
## Cumulative Proportion 0.8273 0.9294 0.96104 0.97923 0.98661 0.99180 0.99608
## PC8 PC9 PC10 PC11 PC12 PC13 PC14
## Standard deviation 1.37629 0.60761 0.42677 0.22385 0.13114 0.09806 0.02037
## Proportion of Variance 0.00295 0.00057 0.00028 0.00008 0.00003 0.00001 0.00000
## Cumulative Proportion 0.99902 0.99960 0.99988 0.99996 0.99998 1.00000 1.00000
summary(princomp(PCdata))
## Importance of components:
## Comp.1 Comp.2 Comp.3 Comp.4 Comp.5
## Standard deviation 22.8879807 8.0407225 4.4742688 3.39400594 2.161016463
## Proportion of Variance 0.8273203 0.1021054 0.0316157 0.01819215 0.007375218
## Cumulative Proportion 0.8273203 0.9294257 0.9610414 0.97923357 0.986608787
## Comp.6 Comp.7 Comp.8 Comp.9
## Standard deviation 1.812162414 1.646620645 1.36566514 0.602917613
## Proportion of Variance 0.005186244 0.004281992 0.00294542 0.000574083
## Cumulative Proportion 0.991795031 0.996077023 0.99902244 0.999596525
## Comp.10 Comp.11 Comp.12 Comp.13
## Standard deviation 0.4234770324 2.221225e-01 1.301262e-01 9.730192e-02
## Proportion of Variance 0.0002832164 7.791911e-05 2.674164e-05 1.495208e-05
## Cumulative Proportion 0.9998797419 9.999577e-01 9.999844e-01 9.999994e-01
## Comp.14
## Standard deviation 2.021332e-02
## Proportion of Variance 6.452592e-07
## Cumulative Proportion 1.000000e+00
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.
PCdata <- scale(PCdata, center = TRUE, scale = TRUE)
summary(prcomp(PCdata))
## Importance of components:
## PC1 PC2 PC3 PC4 PC5 PC6 PC7
## Standard deviation 2.9922 1.3774 1.2203 0.83936 0.57503 0.38570 0.37509
## Proportion of Variance 0.6395 0.1355 0.1064 0.05032 0.02362 0.01063 0.01005
## Cumulative Proportion 0.6395 0.7750 0.8814 0.93170 0.95532 0.96594 0.97599
## PC8 PC9 PC10 PC11 PC12 PC13 PC14
## Standard deviation 0.31946 0.25953 0.24952 0.21451 0.19375 0.11730 0.08427
## Proportion of Variance 0.00729 0.00481 0.00445 0.00329 0.00268 0.00098 0.00051
## Cumulative Proportion 0.98328 0.98809 0.99254 0.99583 0.99851 0.99949 1.00000
summary(princomp(PCdata))
## Importance of components:
## Comp.1 Comp.2 Comp.3 Comp.4 Comp.5
## Standard deviation 2.9690582 1.3667152 1.2108763 0.83287359 0.57059435
## Proportion of Variance 0.6395033 0.1355069 0.1063665 0.05032265 0.02361893
## Cumulative Proportion 0.6395033 0.7750102 0.8813767 0.93169932 0.95531825
## Comp.6 Comp.7 Comp.8 Comp.9 Comp.10
## Standard deviation 0.38272054 0.37219820 0.316989404 0.257526895 0.24759625
## Proportion of Variance 0.01062598 0.01004972 0.007289451 0.004811168 0.00444727
## Cumulative Proportion 0.96594423 0.97599395 0.983283402 0.988094570 0.99254184
## Comp.11 Comp.12 Comp.13 Comp.14
## Standard deviation 0.212856261 0.192254476 0.1163912359 0.0836145071
## Proportion of Variance 0.003286837 0.002681379 0.0009827565 0.0005071876
## Cumulative Proportion 0.995828677 0.998510056 0.9994928124 1.0000000000
Porovnajte směrodatné chyby a odmocniny vlastných čísel spočítaných z výberovej variančnej-kovariančnej matice:
sqrt(eigen(var(PCdata))$values)
## [1] 2.99216408 1.37735123 1.22029957 0.83935520 0.57503484 0.38569896
## [7] 0.37509473 0.31945628 0.25953103 0.24952310 0.21451275 0.19375064
## [13] 0.11729702 0.08426521
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?
4.Grafická vizualizácia PCA
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:
pc <- princomp(PCdata)
par(mfrow = c(1,2))
plot(pc, type = "l", col = "red", main = "")
plot(pc, col = "red", main = "")

summary(pc)
## Importance of components:
## Comp.1 Comp.2 Comp.3 Comp.4 Comp.5
## Standard deviation 2.9690582 1.3667152 1.2108763 0.83287359 0.57059435
## Proportion of Variance 0.6395033 0.1355069 0.1063665 0.05032265 0.02361893
## Cumulative Proportion 0.6395033 0.7750102 0.8813767 0.93169932 0.95531825
## Comp.6 Comp.7 Comp.8 Comp.9 Comp.10
## Standard deviation 0.38272054 0.37219820 0.316989404 0.257526895 0.24759625
## Proportion of Variance 0.01062598 0.01004972 0.007289451 0.004811168 0.00444727
## Cumulative Proportion 0.96594423 0.97599395 0.983283402 0.988094570 0.99254184
## Comp.11 Comp.12 Comp.13 Comp.14
## Standard deviation 0.212856261 0.192254476 0.1163912359 0.0836145071
## Proportion of Variance 0.003286837 0.002681379 0.0009827565 0.0005071876
## Cumulative Proportion 0.995828677 0.998510056 0.9994928124 1.0000000000
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() :
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)

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:
Domáca (samostatná) úloha
(Deadline: Cvičenie č.8 / 20.04.2021)
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.
-
Vyberte si vhodný datový súbor (napr. data, ktoré ste uvažovali v prvej, alebo štvrtej domácej úlohe);
-
Z podkladového datového súboru uvažujte vhodné (numerické) premenné. Dané premenné štandardizujte a aplikujte metódu hlavných komponent.
-
Hlavné komponenty a pôvodné premenné pomocou nástrojov v programe R vhodne vizualizujte.
-
Riešenie umiestnite na svoju webovú stránku, najneskôr v utorok, 20.04.2021, do 14:00.
|