NMST539 | Cvičenie 6

Metóda hlavných komponent | Aplikačná časť

LS 2020/2021 | 13/04/21 | (online výuka)

zdrojový Rmd súbor (kódovanie UTF8)

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?



3. Miera informácie obsiahnutej v hlavných komponentách

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).

  • Pre tieto účely je možné využiť akýkoľvek jpg súbor – napríklad túto fotografiu:
    http://msekce.karlin.mff.cuni.cz/~maciak/NMST539/foto.jpg

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

    library("jpeg")
foto <- readJPEG("foto.jpg")

  • 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 <- 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)

  • 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):

    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 = ""))
}

  • 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:




  • 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”.
  • 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.


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.