NMST539 | Cvičenie 9 Matice vzdialenosti a mnohorozmerné škálovanie LS 2020/2021 | 27/04/21 | (online výuka) zdrojový Rmd súbor (kódovanie UTF8) Outline deviateho cvičenia: teoretické základy metódy mnohorozmerného škálovania; rôzne typy vzdialenosti pre výpočet matice vzdialenosti; metrická a nemetrická metóda mnohorozmerného škálovania; Š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 mnohorozmerného škálovania Mnohorozmerné škálovanie je štatistická metóda určená primárne pre vizualizáciu mnohorozmerných dat. Na rozdiel od vizualizácie pomocou PCA, alebo FA, ktoré sa zameriavali na zobrazenie dat z hľadiska maximálnej variability, resp. celkovej korelačnej štruktúry, mnohorozmerné škálovanie využíva pri zobrazovaní princíp vzájomnej podobnosti/nepodobnosti jednotlivých pozorovaní. Opäť sa ale jedná o určitú redukciu dimenzionality pôvodných dat aj keď vo výsledku vizualizujeme vpodstate iné data, než boli tie pôvodne. Idea mnohorozmerného škálovania spočíva v identifikácii/určení/nájdení menšieho množstva dimenzii (teda redukcia dimenzionality), tak aby bolo možné pôvodné data čo najvierohodnejšie zobraziť (reprezentovať) v menšom počte dimenzii, pričom dôraz je kladený na zachovanie vzájomnej podobnosti medzi pozorovaniami (objektami). Miera podobnosti/nepodobnosti je kvantifikovaná pomocou tzv. matice vzdialenosti. Existujú samozrejme rôzne postupy, ako mieru podobnosti definovať, resp. akú vzdialenosť pri určovaní nepodobnosti využiť. Špecifické voľby miery podobnosti dávajú špecifické varianty mnohorozmerného škálovania. Z určitého pohľadu sa dá povedať, že ak je podobnosť reprezentována pomocou klasickej korelačnej matice, tak metóda mnohorozmerného škálovania je vlastne faktorovou analýzou. Na druhej strane, ak sa miera nepodobnosti kvantifikuje vzhľadom k Euklidovej vzdialenosti, tak získame metódu hlavných komponent. Samozrejme existujú mnohé ďalšie voľby vzdialenosti. Základným vstupom pre algoritmus mnohorozmerného škálovania je tzv. matica vzdialenosti (resp. matica podobnosti/nepodobnosti) medzi všetkými možnými dvojicami pozorovaní, pričom sa predpokladá, že pôvodná dimenzia pozorovaní je neznáma. Úloha škálovacieho algoritmu je nájsť množinu bodov v priestore určitej (ideálne čo najnižšej) dimenzie, tak aby vzdialenosti medzi bodmi množiny čo najviac korespondovali s konfiguráciou v pôvodných datach. V určitých prípadoch je dokonca možné metódu mnohorozmerného škálovania aplikovať na základe matice podobnosti/nepodobnosti, ktorá nie je postavená na vzdialenosti (definícia vzdialenosti ako metriky), ale napr. na základe úrčitých skórov, alebo na základe poradí (resp. nejakej vhodnej monotónnej funkcie poradí). V takomto prípade hovoríme o tzv. nemetrickom mnohorozmernom škálovaní (nonmetric MDS). V programe R je v štandardnej inštalácii k dispozícii funkcia dist(), ktorá pre vstupné data počíta maticu vzdialenosti (resp. maticu podobnosti), pričom k dispozícii sú rôzne vzdialenosti (viď help k danej funkcii pre ďalšie podrobnosti). Matica vzdialenosti následne slúži ako vstup pre algoritmus mnohorozmerného škálovania – v programe R je k dispozícii funkcia cmdscale(). Pre ilustráciu opäť využijeme data o ekologickej a biologickej kvalite vodných tokov v Českej republike. Biologická a ekologická diverzita jednotlivých lokalít (celkovo 65 rôzných lokalít v ČR v blízkosti významných vodných tokov) je posudzovaná na základe 17 rôznych biologických indexov (definovaných podľa certifikovanej metodiky MŽP ČR). Pre lepšie vizuálne výsledky je premenná ‘Ntaxon’ faktorizovaná do troch kategórii. rm(list = ls()) bioData <- read.csv("http://msekce.karlin.mff.cuni.cz/~maciak/NMST539/bioData.csv", header = T) bioData$Ntaxonfac <- cut(bioData$Ntaxon, breaks = c(0, 12, 18, 40), labels = c("low", "medium", "high")) Dmatrix <- dist(bioData[,2:18]) DmatrixScaled <- dist(scale(bioData[,2:18])) dim(as.matrix(Dmatrix)) ## [1] 65 65 Matica vzdialenosti samozrejme závisí na měřítku pôvodných dat. Metóda mnohorozmerného škálovania preto nie je invariantná voči škálovaniu. Samotný škálovací algoritmus je ale možno aplikovať aj na maticu vzdialenosti pôvodných pozorovaní, aj na maticu vzdialenosti škálovaných (napr. štandardizovaných) dat. Pomocou škálovacieho algoritmu môžeme následne data reprezentovať graficky. Z formálneho hľadiska sa ale v grafe dívame na množinu nových bodov, ktoré určil škálovaci algoritmus tak, aby bola táto množina čo najviac podobná pôvodným datam vzhľadom k vstupnej matici vzdialenosti. mds <- data.frame(cmdscale(DmatrixScaled)) library(ggplot2) ggplot(mds, aes(x = X1, y = X2)) + geom_label(aes(fill = bioData$Ntaxonfac, label = bioData[,1])) Za základe grafu, ktorý zobrazuje v dvoch rozmeroch škálované pôvodné data (resp. nové data, ktoré sa snažia zachovať mieru podobnosti a nepodobnosti v pôvodných datach) dokážeme získať celkom dobrú predstavu o štruktúre dat z hľadiska podobnosti jednotlivých pozorovaní. V určitom zmysle by sa dalo povedať, že žiaden iný dvojrozmerný scatterplot (dve rôzne premenné z pôvodného datového súboru) neposkytne tak komplexnú (a čo najmenej zkreslenú) informáciu ohľadom celkovej podobnosti, ako je tomu v prípade mnohorozmerného škálovania. Uvedomte si, že informácia o pôvodnom rozmere (dimenzionalite) dat nie je vstupnou informáciou do algoritmu mnohorozmerného škálovania. Matica \(\mathbb{D}\) má rozmery \(n \times n\) a teda implicitne obsahuje informáciu o počte pozorovaní, avšak nie o počte premenných. Samostatne Ako by ste interpretovali jednotlivé hodnoty v matici \(\mathbb{D}\) (Dmatrix)? Štandardne používa funkcia dist() Euklidovú vzdialenosť medzi jednotlivými pozorovaniami. Dokážete jednotlivé hodnoty zrekonštruovať – t.j. získať manuálnym výpočtom? Napríklad: sqrt(sum((bioData[1,2:18] - bioData[2,2:18])^2)) ## [1] 20.55428 Príslušná hodnota z matice \(\mathbb{D}\) (Dmatrix): as.matrix(Dmatrix)[1:2,1:2] ## 1 2 ## 1 0.00000 20.55428 ## 2 20.55428 0.00000 Premenná Ntaxon v určitom zmysle súhrnne popisuje kvalitu jednotlivých lokalít. Nízke hodnoty môžu byť interpretované ako slabá ekologická kvalita, výsoké hodnoty reprezentujú lokality s veľmi dobrou ekologickou kvalitou a biologickou diverzitou. Koresponduje táto interpretácia s intuitívnym záverom získanym z grafu? Viď napr. riečné toky v okolí Ostravy, Mohelníc, prípadne Bečva – a porovnajte s vodnými tokmi napr. na Vysočine, alebo na Šumave. Premennú Ntaxon je samozrejme možné kategorizovať aj iným spôsobom (napr. detailnejšie) a následne opäť aplikovať MDS. Mnohorozmerné škálovanie na základe Euklidovej vzdialenosti / teoretické aspekty Matica vzdialenosti \(\mathbb{D} = (d_{ij})_{i, j = 1}^{n}\) je definovaná na základe Euklidovskej metriky \[ d_{ij}^2 = \sum_{k = 1}^{p} (x_{i k} - x_{j k})^2, \] kde predpokládame \(p\)-rozmerné body \(\boldsymbol{x}_i = (x_{i1}, \dots, x_{ip})^\top \in \mathbb{R}^p\) a \(\boldsymbol{x}_j = (x_{j1}, \dots, x_{jp})^\top \in \mathbb{R}^p\). Pre každé \(i, j = 1, \dots, n\) definujme dodatočné kvantity \[ a_{ij} = - \frac{d_{ij}^2}{2}, \quad\quad a_{i \bullet} = \frac{1}{n}\sum_{j = 1}^{n} a_{ij}, \quad\quad a_{\bullet j} = \frac{1}{n}\sum_{i = 1}^{n} a_{ij}, \quad\quad a_{\bullet \bullet} = \frac{1}{n^2} \sum_{i = 1}^n \sum_{j = 1}^n a_{ij}. \] S využitím značenia definovaného vyššie sa da ukázať, že platí následovné: Pre maticu \(\mathbb{B} = (b_{ij})_{i,j = 1}^n\), kde \(b_{ij} = a_{ij} - a_{i \bullet} -a_{\bullet j} + a_{\bullet \bullet}\), platí, že \(\mathbb{B} = \mathcal{X}\mathcal{X}^\top\), pričom \(\mathcal{X}\) je pôvodná matica dat s rozmermi \(n \times p\), t.j. \(\mathcal{X} = (\boldsymbol{x}_1, \dots, \boldsymbol{x}_n)^\top\), kde \(\boldsymbol{x}_i = (x_{i1}, \dots, x_{ip})^\top \in \mathbb{R}^p\). Ukážte, že platí \(b_{ii}= a_{\bullet \bullet} - 2a_{i \bullet}\), pre \(i = 1, \dots, n\). Ukážte, že platí \(\mathbb{B} = -\frac{1}{2} \mathcal{H} \mathbb{D}^2 \mathcal{H}\), kde \(\mathcal{H}\) je tzv. centrovacia matica definovaná ako \(\mathcal{H} = (\mathbb{I}_n - \frac{1}{n} \boldsymbol{1}_n\boldsymbol{1}_n^\top)\). Vzdialenosti a miery podobnosti v R funkcii dist() Pre podrobnosti viď help k funkcii dist(). Euklidova vzdialenosť - defaultná voľba; Maximová vzdialenosť - maximová, resp. supremová norma; Manhattan vzdialenosť - tzv. \(L_1\) norma; Canberra vzdialenosť - vážena verzia Manhattan vzdialenosti; Binárna vzdialenosť - proporcia nenulových elementov vo vektore \(\boldsymbol{x}−\boldsymbol{y}\); Minkowského vzdialenosť - štandardná \(L_p\)-norma (defaultná voľba je \(p=2\), čo dáva Euklidovú vzdialenosť); Naďalej využijeme defaultnú voľbu Euklidovej vzdialenosti, ale pozrieme sa na rozdiely, ktoré vzniknú štandardizáciou dat. MDS1 <- cmdscale(Dmatrix, k = 2) dim(MDS1) ## [1] 65 2 MDS2 <- cmdscale(DmatrixScaled, k = 2) dim(MDS2) ## [1] 65 2 Čo reprezentuju jednotlivé výstupy? Ako by ste hodnoty interpretovali? Je nutné pamätať na to, že data sú mnohorozmerné. Z pohľadu na scatterplot pôvodných data síce môžeme získať určitú informáciu, ale táto informácia je výrazne limitovaná mnohorozmerným charakterom dat (často môže byť dokonca hodne zavádzajúca). Ilustrujeme to na príklade troch konkrétnych lokalít, pričom data vizualizujeme štandardným scatterplotom, ale zakaždým zvolime iné dve premenné. par(mfrow = c(2,2)) cols <- rep("blue", dim(bioData)[1]) cols[which(bioData[,1] == "Olse")] <- "red" cols[which(bioData[,1] == "DesnSudk")] <- "green" cols[which(bioData[,1] == "DyjeVran")] <- "green" plot(bioData[,2], bioData[,3], xlab="SaprInd", ylab="Lital", type="n", main = "Czech River Localities") text(bioData[,2], bioData[,3], labels = bioData[,1], cex=.5, col = cols) plot(bioData[,2], bioData[,11], xlab="SaprInd", ylab="Ntaxon", type="n", main = "Czech River Localities") text(bioData[,2], bioData[,11], labels = bioData[,1], cex=.5, col = cols) plot(bioData[,2], bioData[,17], xlab="SaprInd", ylab="MDS_1", type="n", main = "Czech River Localities") text(bioData[,2], bioData[,17], labels = bioData[,1], cex=.5, col = cols) plot(bioData[,17], bioData[,18], xlab="MDS_1", ylab="MDS_2", type="n", main = "Czech River Localities") text(bioData[,17], bioData[,18], labels = bioData[,1], cex=.5, col = cols) Je zrejme, že na základe scatterplotov vyššie, nie je ľahké odpovedať na otázku, či vodné lokality Olše, Desna, a Dyje sú podobné – t.j. blízko seba vzhľadom k datam, nie vzhľadom k ich skutočnej geografickej polohe. Tieto lokality sa zdajú byť dosť podobné vzhľadom k niektorým premenným (napr. prvý graf), ale na druhej strane sú dosť rozdielné vzhľadom k iným premenným – napr. posledný graf. S využitím dvojrozmerných scatterplotov je informácia, ktorú získame o štruktúre dat (čo sa týka aj vzájomnej podobnosti/nepodobnosti jednotlivých pozorovaní) výrazne limitovaná. Ako alternatívu ale môžeme využiť metódu mnohorozmerného škálovania, ktorá vizualizuje (síce iné) data v nižšom rozmere, ale zároveň sa snaží o čo najdôveryhodnejšie zobrazenie miery podobnosti a nepodobnosti z pôvodných dat. Grafické zobrazenie pomocou mnohorozmerného škálovania (škálovanie do dvoch rozmerov): plot(MDS1[,1], MDS1[,2], xlab="Coordinate 1", ylab="Coordinate 2", type="n", main = "Czech River Localities") text(MDS1[,1], MDS1[,2], labels = bioData[,1], cex=.7, col = cols) A verzia na základe štandardizovaných dat: plot(MDS2[,1], MDS2[,2], xlab="Coordinate 1", ylab="Coordinate 2", type="n", main = "Czech River Localities") text(MDS2[,1], MDS2[,2], labels = bioData[,1], cex=.7, col = cols) Výsledky sa vizuálne hodne líšia v závislosti na tom, či sú pôvodné data štandardizované, alebo nie. Na druhej strane interpretácia (vzhľadom k celkovej podobnosti a nepodobnosti jednotlivých lokalít) je viacmenej totožná. Graf reprezentuje škálovane vzdialenosti medzi jednotlivými lokalitami (vzialenosti vhľadom k podobnosti jednotlivých lokalít, nie geografické vzdialenosti). Napr. z pohľadu Olše sú jej najpodobnejšia lokality označené ako Moravka, Tyra a VidnVidn – a to nezávisle na tom, či používame pôvodné data, alebo štandardizované. Celkový charakter troch lokalít (Olše, Desna, and Dyje) nie je možné posúdiť na základe dvojrozmerného scatterplotu. Na druhej strane, najviac relevantná informácia ohľadom podobnosti a nepodobnosti týchto troch lokalít je k dispozícii v scatterplotoch získaných na báze mnohorozmerného škálovania vyššie. Z tohto pohľadu sa teda zdá, že dané lokality sú skôr výrazne rozdielné – určite výraznejšie, než sme získali pocit napr. na základe porovnania vzhľadom k indexom Saprind a Lital – t.j. v prvom zo štvorice grafov. Samostatne Aplikujte metódu mnohorozmerného škálovania na dané data, ale využijte inú vzdialenosť pri aplikovaní funkcie dist(). Výsledky zobrazte a porovnajte. Stejné lokality lze porovnať na základe rôznych chemických koncentrácii, ktoré boli na lokalitách namerané. Data s chemickými koncentráciami je možné načítať do programu R pomocou príkazu: chemData <- read.csv("http://msekce.karlin.mff.cuni.cz/~maciak/NMST539/chemData.csv", header = T) Chemické a biologické premenné je možné zkombinovať dohromady a vizualizovať mieru podobnosti a nepodobnosti jednotlivých lokalít na základe komplexnejších bio-chemických premenných. V následujúcom aplikujeme metódu mnohorozmerného škálovania samostatne na biologické premenné a samostatne na chemické premenné. Výsledok ale vykreslíme v jednom spoločnom grafe. Každá lokalita sa v grafe objavuje dvakrát – čiernou farbou pre naškálovane lokality vzhľadom k podobnosti pôvodných biologických premenných a červenou farbou pre naškálovane lokality vhľadom k podobnosti nových, chemických premenných. chemData <- read.csv("http://msekce.karlin.mff.cuni.cz/~maciak/NMST539/chemData.csv", header = T) ind <- match(chemData[,1], bioData[,1]) bioData0 <- bioData[ind, ] MDS3 <- cmdscale(dist(bioData0[,2:18]), k = 2) MDS4 <- cmdscale(dist(chemData[,2:8]), k = 2) plot(c(MDS3[,1], MDS4[,1]), c(MDS3[,2], MDS4[,2]), xlab="Coordinate 1", ylab="Coordinate 2", type="n", main = "Czech River Localities") text(MDS3[,1], MDS3[,2], labels = bioData[,1], cex=.7) text(MDS4[,1], MDS4[,2], labels = bioData[,1], cex=.7, col = "red") legend(50, -20, legend = c("Bio Metric Measurements", "Chemical Concentrations"), col = c("black", "red"), pch = 15) Otázky na zamyslenie Nakoľko jednoznačné je zobrazenie prezentované v obrázku vyššie? Aké závery by bolo možné získať – jednak z komplexného hľadiska (bio/chemické premenné dohromady) a tiež z pohľadu buď na chemické, resp. iba na biologické premenné? 2. Metrické škálovanie vs. nemetrické škálovanie Základný rozdiel medzi oboma postupmi je v tom, že prvý využíva maticu vzdialenosti – pričom vzdialenosť je definovaná v zmysle štandardnej metriky – zatial čo druhý postup využíva pouze maticu poradí (resp. maticu nejakej vhodnej monotónnej funkcie poradí) – pričom poradia nesplňuju definíciu metriky. Pripomeňte si základné definitorické vlastnosti metriky/vzdialenosti. Opäť existujú rôzne môžnosti pre voľbu vhodnej monotónnej funkcie poradí. V následujúcom si stručne ilustrujeme nemetrické mnohorozmerné škálovanie – tzv. Kruskalov algoritmus (Kruskal’s Non-metric Multidimensional Scaling). Z dôvodu lepšieho porovnania ale použijeme štandardnú maticu vzdialenosti. library(MASS) MDS2 <- isoMDS(Dmatrix, k=2) ## initial value 7.466727 ## iter 5 value 6.387632 ## iter 10 value 6.118217 ## iter 15 value 5.959384 ## final value 5.817192 ## converged Grafická vizualizácia: plot(MDS2$points[,1], MDS2$points[,2], xlab="Coordinate 1", ylab="Coordinate 2", type="n", main = "Czech River Localities") text(MDS2$points[,1], MDS2$points[,2], labels = bioData[,1], cex=.7) Otázky V akých situáciach je vhodné použiť nemetrické škálovanie – non-metric MDS? V akých situáciach by bola preferované štandardná verzia škálovania (na základe vzdialenosti/metriky)? V programe R je množstvo ďalších možností, ako aplikovať metódu mnohorozmerného škálovania na data. Okrem funkcii cmdscale() a isoMDS(), ktoré sme ilustrovali vyššie, sú k dispozícii napr. funkcia smacofSym() v R knižnici ‘smacof’, alebo funkcia wcmdscale() v knižnici ‘vegan’, alebo pco() v R knižnici ‘ecodist’. Domáca (samostatná) úloha (Deadline: Cvičenie č.10 / 04.05.2021) Vyberte si vhodný datový súbor (napr. data, ktoré ste uvažovali v prvej, štvrtej, alebo predchádzajúcich dvoch úlohach); Aplikujte metódu mnohorozmerného škálovania a výsledky vizualizujte. Pokúste sa ilustrovať podobný prípad, kde z pohľadu niektorých dvoch premenných v pôvodných datach budú niektoré pozorovania vyzerať hodne podobne, zatiaľ čo mnohorozmerný charakter dat – zrejmý po mnohorozmernom škálovaní – bude opačný. Riešenie umiestnite na svoju webovú stránku, najneskôr v utorok, 04.05.2021, do 14:00.