NMST539 | Cvičenie 10 Zhluková analýza LS 2020/2021 | 04/05/21 | (online výuka) zdrojový Rmd súbor (kódovanie UTF8) Outline desiateho cvičenia: základy zhlukovej analýzy v programe R; princípy hierarchického a nehierarchického zhlukovania; K-means algoritmus; Š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. Základné princípy zhlukovej analýzy Zhluková analýza (Cluster Analysis, CA) je štatistická metóda určená k rozlíšeniu jednotlivých objektov (pozorováni) pomocou vytvorenia určitého počtu konkrétnych a disjunktných zhlukov (tzv. clustrov), do ktorých sú pozorovania jednoznačne priradené na základe určitej miery podobnosti/nepodobnosti. Táto miera vzájomnej podobnosti/nepodobnosti je väčšinou meraná pomocou konkrétnej metriky (tzv. proximity measure) – t.j. napr. nejakej vhodnej vzdialenosti – v programe R k tomuto účelu slúži funkcia dist(). Je ale nutné pripomenúť, že zhluková analýza nepredstavuje jednú konkrétnu metódu klasifikácie subjektov (pozorovaní) do skupín (clustrov), ale jedná sa o súhrnné označenie množstva rôznych postupov a algoritmov, ktoré sú určené k tomuto účelu. Zaroveň je nutné zdôrazniť základný a principiálny rozdiel medzi zhlukovaním a medzi diskrimináciou/klasifikáciou (resp. zhlukovou analýzou a diskriminačnou a klasifikačnou analýzou). Zhlukovanie funguje bez apriornej informácie o skutočnom zaradení subjektov do jednotlivých skupín. Skupiny pozorovaní (zhluky) sú neznáme a cielom zhlukovacieho algoritmu je využiť mieru podobnosti a nepodobnosti medzi objektami a tieto zhluky definovať (a k ním priradiť príslušné pozorovania). Klasifikácia funguje na princípe datového súboru, u ktorého je priradenie pozorovaní do jednotlivých skupín predom známe (tzv. group label) a cieľom je pochopiť princíp priradenia jednotlivých pozorovaní do jednotlivých skupín, aby bolo možné následne nové pozorovanie bez ‘group label’ zaradiť do správnej kategórie. Z teoretického hľadiska rozlišujeme dva základne postupy (algoritmické princípy) v zhlukovaní: Hierarchické algoritmy - priradenie subjektov do jednotlivých skupín je stabilné a v priebehu algoritmu sa už nemení. Nehierarchické (tzv. partitioning) algoritmy - priradenie subjektov do jednotlivých skupín sa môže v priebehu algoritmu meniť. Finálný výsledok zhlukovacieho algoritmu je ovplyvnený konkrétnou voľbou vzdialenosti (proximity measure) a a konrkétnou voľbou zhlukovacieho algoritmu. K dispozícii je množstvo rôznych zhlukovacích algoritmov, ktoré pracuju na rôznych principoch. Výstupom z algoritmu je ale zakaždym priradenie jednotlivých pozorovaní do určitého počtu disjunktných zhlukov (clustrov). Zhlukovacie algoritmy môžeme rozdeliť napr. následovne: Hierarchické algoritmy Partitioning algoritmy Algoritmy založené na strede zhlukov (tzv. centroid based) Algoritmy založene na distribúcii (tzv. distribution based) Algoritmy založené na hustote (tzv. density based) alebo iné… Poznámka Zhluková analýza patrí prevažne do exploratívnej časti štatistických metód. V určitom zmysle je ale možné zhlukovanie chápať aj ako metódu redukcie dimenzionality. Z interpretačného hľadiska má zmysel hľadať počet jednotlivých zhlukov tak, aby zhluky umožňovali zmysluplnú intepretáciu. Po určení konkrétnych zhlukov je následne možné jednotlivé pozorovania reprezentovať pomocou ich príslušnosti do zhlukov (redukcia dimenzionality – t.j. namiesto pôvodných premenných využiť novú kategorickú premennú – tzv. group label). 2. Hierarchické zhlukovanie Hierarchické aglomeratívne zhlukovanie funguje na princípe ‘zdola-hore’, kedy sa najprv vytvárajú malé zhluky (jednoprvkové zhluky na začiatku algoritmu), ktoré sa následne v postupnosti krokov zlučujú do väčších zhlukov, až je vytvorený jeden celkový zhluk obsahujúci všetky pozorovania. Pri zlučovaní menších zhlukov v jednotlivých krokoch sa využíva princíp vzájomnej podobnosti (medzi pozorovaniami navzájom, medzi pozorovaniami a zhlukmi, resp. medzi zhlukmi navzájom), definovaný na základe použitej vzdialenosti, resp. proximity measure. Aglomeratívne hierarchické algoritmy sú výpočetne veľmi jednoduché najme vďaka faktu, že podobnosť (resp. vzdialenosť) medzi zhlukmi v jednotlivých krokoch je možné napočítať z matice vzdialenosti v predchádzajích krokoch (vstupom pre algoritmus je teda matica vzdialenosti medzi pôvodnýcmi pozorovaniami, z ktorých každé tvorí vlastný, samostatný zhluk, cluster). Hierarchické deliace algoritmy pracujú na princípe ‘zhora-dole’, kedy je najprv vytvorený jeden komplexný zhluk obsahujúci všetky pozorovania dohromady a následne v jednotlivých krokoch dochádza k deleniu existujúcich zhlukov, opäť na princípe vzájomnej nepodobnosti medzi zhlukmi (resp. medzi zhlukmi a pozorovaniami). Deliace (divisive) aglomeratívne algoritmy bývajú výpočetne náročné, nakoľko celková hierarchia rôznych delení musí byť vybudovaná za účelom dosiahnutia finálneho delenia. Výstup zhlukovacieho algoritmu je možné graficky vizualizovať, čo výrazne pomôže pri voľbe finálneho počtu zhlukov. Na druhu stranu, vizualizácia je praktický neinformatívna, ak sa jedná o zhlukovanie väčšieho počtu pozorovaní. Klučovým problémom teda zostáva voľba finálneho počtu zhlukov. . Pre ilustráciu zhlukovania v programe R využijeme datový súbor mtcars. data <- mtcars attach(data) head(data) ## mpg cyl disp hp drat wt qsec vs am gear carb ## Mazda RX4 21.0 6 160 110 3.90 2.620 16.46 0 1 4 4 ## Mazda RX4 Wag 21.0 6 160 110 3.90 2.875 17.02 0 1 4 4 ## Datsun 710 22.8 4 108 93 3.85 2.320 18.61 1 1 4 1 ## Hornet 4 Drive 21.4 6 258 110 3.08 3.215 19.44 1 0 3 1 ## Hornet Sportabout 18.7 8 360 175 3.15 3.440 17.02 0 0 3 2 ## Valiant 18.1 6 225 105 2.76 3.460 20.22 1 0 3 1 Využijeme všetky premenné, ktoré sú v datovom súbore mtcars k dispozícii a najrpv spočítame maticu vzdialenosti medzi jednotlivými pozorovaniami pomocou funkcie mtcars. Použijeme defaultne nastavenie funkcie – teda Euklidovskú vzdialenosť. D <- dist(mtcars) ### for the euclidian distance by default Následne použijeme maticu vzdialenosti \(\mathbb{D}\) ako vstup pre hierarchický aglomeratívny algoritmus – funkcia hclust() v programe R (funkcia je k dispozícii pod štandardnou inštaláciou programu R). Podrobnosti v helpe funkcie> ?hclust(). HC1 <- hclust(D) plot(HC1, xlab = "Observations", ylab = "Proximity measure") Defaultnou zhlukovacou metódou implementovanou vo funkcii hclust() je tzv. ‘complete linkage’ prístup, ktorý využíva pri vytváraní skupín princíp najvzdialenejišieho suseda. Ďalšie metódy sú k dispozícii pod voľbou dodatočného parametru ‘method’: ‘method = ’complete’ (najvzdialenejší sused | default) ‘method = ’single’ (najbližší sused) ‘method = ’average’ (priemer zhluku, resp. kompromis medzi najvzdialenejším a najbližším susedom) ‘method = ’median’ (robustná verzia metódy s priemerom) ‘method = ’centroid’ (postup na základe geometrickej vzdialenosti) ‘method = ’mcquitty’ (vážena verzia metódy priemerov) ‘method = ’ward.D’ (postup založený na najmenšom náraste tzv. measure of heterogeneity) Rôzne zhlukovacie algoritmy sú k dispozícii aj v dodatočných rozšíreniach a R balíčkoch. Napr. knižnica ‘sparcl’ (inštalácia pomocou príkazu install.packages('sparcl')), ktoráá ponúka funkciu ColorDendrogram(). Ak je k dispozícii informácia o konečnom počete zhlukov, je možné počet zhlukov definovať nasledujúcim spôsobom: plot(HC1, xlab = "Observations", ylab = "Proximity measure") groups <- cutree(HC1, k=3) rect.hclust(HC1, k=3, border="red") Ekvivalentná verzia, pouze vylepšená z vizuálneho hľadiska: plot(HC1, xlab = "Observations", ylab = "Proximity measure", hang = -1) groups <- cutree(HC1, k=3) rect.hclust(HC1, k=3, border="red") alebo: cutHC1 = cutree(HC1, 3) labelColors = c("red", "green", "blue") # function to get color labels colLab <- function(n) { if (is.leaf(n)) { a <- attributes(n) labCol <- labelColors[cutHC1[which(names(cutHC1) == a$label)]] attr(n, "nodePar") <- c(a$nodePar, lab.col = labCol) } n } # using dendrapply clusDendro = dendrapply(as.dendrogram(HC1), colLab) plot(clusDendro, main = "Cool Dendrogram", type = "triangle") Pomocou trochu sofistikovanejšieho kódu je možné v obrázku znázorniť aj dodatočné informácie (napr. niektoré premenné) z pôvodného datového súboru (napr. počet prevodových stupňov – gear). local({ colLab <<- function(n) { if(is.leaf(n)) { a <- attributes(n) i <<- i+1 attr(n, "edgePar") <- c(a$nodePar, list(col = mycols[i], lab.font= i%%3)) } n } mycols <- mtcars$gear i <- 0 }) x.clust.dend1 <- dendrapply(as.dendrogram(HC1), colLab) plot(x.clust.dend1) Samostatne Využijte defaultný zhlukovací algoritmus vo funkcii hclust(), ale pre maticu \(\mathbb{D}\) použijte rôzne varianty pre definíciu vzdialenosti. Porovnajte výsledky. Zvoľte vhodnu vzdialenosti pri aplikovaní funkcie dist() a využijte rôzne postupy zhlukovania ktoré sú k dispozícii vo funkcii hclust() pod voľbou parametru ‘method = …’. Opäť porovnajte výsledky. Výsledky ilustrujte a porovnajte pomocou dendogramov. 2. Nehierarchické zhlukovanie Princíp fugovania nehierarchického zhlukovania sme už uviedli v úvode. Najznámejším a asi aj načastejšie používaným nehierarchickým zhlukovacím algoritmom je tzv. K means algoritmus. Myšlienka tohto algoritmu spočíva vo vytvorení predom definovaného počtu \(K\) zhlukov (clustrov), tak aby bol celkový súčet štvorcov vzdialenosti jednotlivých pozorovaní od stredu im príslušiacého zhluku minimálny. Predpokládajme, že pozorovania \(\boldsymbol{x}_i\) pre \(i = 1,\dots, n\) formujú \(K\) rôznych zhlukov v zmysle im príslušných \(K\) stredných hodnôt, čo môžeme teoreticky zapísať ako \[ E[\boldsymbol{x}_i | k_i] = \boldsymbol{\mu}_{k_{i}}, \] kde \(k_{i} \in \{1, \dots, K\}\). Celkový počet zhlukov je teda \(K \in \mathbb{N}\) a tento počet sa v priebehu algoritmu nemení. Napríklad, pre \(k_i = 1\), dostaneme, že pozorovanie \(\boldsymbol{x}_i\) náleži do zhluku č.1, teda jeho stredná hodnota je \(\boldsymbol{\mu}_{k_{i}} \equiv \boldsymbol{\mu}_1\). Jednotlivé parametre stredných hodnôt \(\boldsymbol{\mu}_{k}\), pre \(k = 1, \dots, K\) sú neznáme parametre, ktorých odhad sa v priebehu algoritmu môže výrazne meniť (odhad totiž závísí od jednotlivých pozorovaní, ktoré sú zaradené do daného zhluku). Výstupom algoritmu je konfigurácia, ktorá pri daných pozorovaniach a pri danom počete zhlukov minimalizuje celkový súčet štvorcov chýb (resp. vzdialenosti od stredu), t.j. výraz \[ \sum_{k =1}^{K} \sum_{i \in \pi(k)} (\boldsymbol{x}_i - \boldsymbol{\mu}_k)^2 \] je najmenší možný, kde \(\pi(k) = \{i \in \{1, \dots, n\};~E[\boldsymbol{x}_i | k] = \boldsymbol{\mu}_k\}\) je množina indexov pozorovaní, ktoré patria do zhluku \(k \in \{1, \dots, K\}\). Jednotlivé stredné hodnoty majú určitú ‘’váhu’’ (určité zastúpenie) v celkovom súbore pozorovaní (v tzv. celkovej ‘’mixture’’) a pre nové pozorovanie \(\boldsymbol{x}\), bez znalosti jeho priradenia do konkrétneho zhluku, možeme pre strednú hodnotu napísať \[ E[x] = P(k = 1) \mu_1 + \dots + P(k = K) \mu_{K}. \] V programe R je pre zhlukovací K-means algoritmus určená funkcia kmeans(), opäť k dispozícii pod štandardnou inštaláciou programu R. Cieľom algoritmu je nájsť, resp. zostrojiť \(K \in \mathbb{N}\) disjunktných clustrov (hodnota \(K\) je zvolená predom a je vstupom pre K-means algoritmus), tak aby celková variabilia vrámci clustrov (within-cluster variability) bola minimálna. Stredy jednotlivých clustrov, teda parametre \(\boldsymbol{\mu}_k\) sú neznáme. Pripomína to klasický problém sliepka-vajce a otázku, čo bolo skôr: ak by sme poznali správne priradenie pozorovaní \(\boldsymbol{x}_i\) do jednotlivých clustrov \(k_i\), problém nájdenia (odhadnutia) neznámych priemerov \(\boldsymbol{\mu}_k\) by bol triviálny. Na druhú stranu, ak by sme poznali skutočné hodnoty neznámych parametrov \(\boldsymbol{\mu}_k\), úloha priradenia pozorovaní medzi jednotlivé clustre by bola opäť triviálna. Riešenie v K-means algoritme sa hľada postupne vo forme iteratívnych krokov. V niektorých prípadoch môže byť algoritmus pomerne časovo náročný. Pre ilustráciu opôť využijeme datový súbor mtcars a budeme hľadať tri neznáme clustre: Kmeans <- kmeans(mtcars, centers = 3) colVector <- as.numeric(Kmeans$cluster) plot(mtcars$mpg ~ mtcars$hp, bg = colVector, xlab = "Horse Power", ylab = "Miles per Gallon", pch = 21, col = "black") points(Kmeans$centers[,1] ~ Kmeans$centers[,4], col = 1:3, pch = 8, cex = 2) text((mtcars$mpg + 1) ~ mtcars$hp,labels=rownames(mtcars), col=colVector, cex = 0.5) Je asi nutné pripomenúť, že sa opäť dívame na mnohorozmerné data zobrazené prostredníctvom dvojrozmerného \(xy\) scatterplotu. Interpretácia zhlukov na základe predchádzajúceho obrázku asi nie je jednoduchá. Situácia sa ale zmení, ak vykreslíme tie isté data a tie isté zhluky prostredníctvom iných dvoch premenných. Dodatočne zahrnieme aj informáciu o počte válcov. plot(mtcars$mpg ~ mtcars$disp, bg = colVector, xlab = "Displacement", ylab = "Miles per Gallon", pch = 19 + cyl/2, col = "black") points(Kmeans$centers[,1] ~ Kmeans$centers[,3], col = "black", pch = 21, cex = 2, bg = 1:3) text((mtcars$mpg + 1) ~ mtcars$disp,labels=rownames(mtcars), col=colVector, cex = 0.5) legend(420, 34, legend = c("4 Cylinders", "6 Cylinders", "8 Cylinders"), pch = c(21, 22, 23), bg = "grey") Poznámka Z popisu algoritmov je zrejme, že škálovanie hraje podstatnú rolu: ak nahradíme niektorú premennú, napr. \(X_{j}\) novou premennou \(2 X_{j}\) (t.j. preškálujeme a použijeme polovičné měřítko), tak potom príslušná dimenzia sa započíta dvojnásobne pri určovaní vzdialenosti medzi pozorovaniami. Následne to bude mať aj vplyv na priradenie pozorovaní do zhlukov. Doporučuje sa používať štandardizované data. Pre centroidy jednotlivých zhlukov \(\mu_{k j}\) pre \(k = 1, \dots, K\) a pôvodné premenné \(j = 1, \dots, p\) môžeme clustre interpretovať v zmysle štandardných odchýliek od marginálných priemerov. Aplikujte K-means algoritmus analogicky ako v predchádzajúcom prípade, ale najprv štandardizujte data. Výsledky porovnajte. Zhluková analýza pripomína v určitom zmysle zhlukové modely – tzv. mixture models – na ktoré sa pozrieme za dva týždne. Využijeme ešte jeden ilustračný príklad. Data reprezentujú konzumáciu proteinu v Európe. protein <- read.csv("http://msekce.karlin.mff.cuni.cz/~maciak/NMST539/protein.csv", header = T, row.names=1) sprotein <- scale(protein) head(protein) ## RedMeat WhiteMeat Eggs Milk Fish Cereals Starch Nuts Fr.Veg ## Albania 10.1 1.4 0.5 8.9 0.2 42.3 0.6 5.5 1.7 ## Austria 8.9 14.0 4.3 19.9 2.1 28.0 3.6 1.3 4.3 ## Belgium 13.5 9.3 4.1 17.5 4.5 26.6 5.7 2.1 4.0 ## Bulgaria 7.8 6.0 1.6 8.3 1.2 56.7 1.1 3.7 4.2 ## Czechoslovakia 9.7 11.4 2.8 12.5 2.0 34.3 5.0 1.1 4.0 ## Denmark 10.6 10.8 3.7 25.0 9.9 21.9 4.8 0.7 2.4 Pre K-means algoritmus použijeme len informáciu o konzumácii bieleho a červeného mäsa – t.j. prvé dve premenné z pôvodných dat (ktoré pre použitie zhlukovacieho algoritmu štandardizujeme). subset<-sprotein[,(1:2)] plot(subset,pch=19,col=4,cex=0.2) text(subset,labels=rownames(protein)) Definujeme tri zhluky a pomocou K-means algoritmu sa pokusíme roztriediť Európske krajiny podľa miery konzumácie bieleho a červeného mäsa. grpMeat<- kmeans(subset, center=3, nstart=10) plot(subset,pch=19,col=4,cex=0.2) text(subset,labels=rownames(protein), col=grpMeat$cluster+1) means<-grpMeat$centers points(means,col=c(2,3,4),cex=2,pch=4) Celkový počet clustrov bol zvolený predom. V niektorých prípadoch ale chýba táto informácia a nie je jednoduché počet clustrov zvoliť predom. Pre voľbu počtu clustrov sa doporučujú rôzne postupy: Zhluková analýza sa často využíva ako exploratívnýy nástroj. Z tohto pohľadu je užitočné voliť hodnotu \(K\) tak, aby celkový počet clustrov dával zmysel z hľadiska interpretácie. Iný postup pre voľbu \(K \in \mathbb{N}\) je založený na rôznych modeloch, kedy sa zhlukovací algoritmus aplikuje pre postupnosť hodnôt \(K_1 < K_2 < \dots, K_M\) a výsledne sa zvolí model, ktorý dáva najlepšie výsledky (napr. z pohľadu vhodne zvoleného tzv. ‘’model selection’’ kritéria). Využiť lze napr. BIC alebo AIC kritérium (alebo im podobné), prípadne aplikovať tzv. krížove overovanie (cross-validation). Cekový počet stupňov voľnosti (resp. parametrov) je \(K \times p\), kde \(K\) je počet zhlukov a \(p \in \mathbb{N}\) je počet premenných v pôvodných datach (t.j. jeden neznámy parameter pro strednú hodnotu každého zhluku v každej dimenzii). Celkový súčet štvorcov pre data o konzumácii bieleho a červeného mäsa a tri clustr (K-means algoritmus): mean_overall<-apply(subset,2,mean) means_between<-(t(means)-mean_overall)^2 sum(t(means_between)*grpMeat$size) ## [1] 35.96535 Variabilitu vrámci jednotlivých clustrov môžeme vizualizovať: grpMeat<- kmeans(subset, center=3, nstart=10) plot(subset,pch=19,col=4,cex=0.2) text(subset,labels=rownames(protein), col=grpMeat$cluster+1) means<-grpMeat$centers points(means,col=c(2,3,4),cex=2,pch=4) for(i in (1:nrow(protein))){ if (grpMeat$cluster[i]==1){ lines(c(subset[i,1],means[1,1]), c(subset[i,2],means[1,2]),col=2) } if (grpMeat$cluster[i]==2){ lines(c(subset[i,1],means[2,1]), c(subset[i,2],means[2,2]),col=3) } if (grpMeat$cluster[i]==3){ lines(c(subset[i,1],means[3,1]), c(subset[i,2],means[3,2]),col=4) }} Analogicky, variabilita medzi jednotlivými clustrami: plot(subset,pch=19,col=grpMeat$cluster+1,cex=0.2) means<-grpMeat$centers points(means,col=c(2,3,4),cex=2,pch=4) points(mean_overall[1],mean_overall[2],pch=4,lwd=3) lines(c(mean_overall[1],means[3,1]), c(mean_overall[2],means[3,2]), col=4) lines(c(mean_overall[1],means[2,1]), c(mean_overall[2],means[2,2]), col=3) lines(c(mean_overall[1],means[1,1]), c(mean_overall[2],means[1,2]), col=2) text(c(10,12,9),c(10,8,6),grpMeat$size) Celkový súčet štvorcov je (podobne ako v prípade lineárnej regresie) rozložený do troch hlavných častí: variabilita vrámci jednotlivých clustrov (within clusters sum of squares) grpMeat$tot.withinss ## [1] 12.03465 variabilita medzi jednotlivými zhlukmi (between clusters sum of squares) grpMeat$betweenss ## [1] 35.96535 a na záver, celková variabilita dat (total sum of squares) definovaná ako súčet predchádzajúcich dvoch kvantít grpMeat$tot.withinss+grpMeat$betweenss ## [1] 48 plot(subset,pch=19,col=4,cex=0.2) text(subset,labels=rownames(protein), col=grpMeat$cluster+1) means<-apply(subset, 2, mean) points(means[1], means[2],bg="black",cex=2,pch=21) for(i in (1:nrow(subset))){ lines(c(subset[i,1],means[1]), c(subset[i,2],means[2]),col="black", lty = 3) } Samostatne Aplikujte K-means algoritmus s inými počiatočnými hodnotami. Aký to ma vplyv na výsledné zhluky? Pomocou internetu preskúmajte ďalšie možností zhlukovej analýzy, ktoré sú v rôznych R knižniciach a rozšíreniach. Domáca (samostatná) úloha (Deadline: Cvičenie č.11 / 11.05.2021) Vyberte si vhodný datový súbor (napr. data, ktoré ste uvažovali v prvej, štvrtej, alebo predchádzajúcich troch úlohach); Použijte K-means algoritmus a pomocou vhodného rozhodovacieho kritéria rozhodnite o počte zhlukov. Následne využijte hierarchický aglomeratívny zlukovací algoritmus a zvoľte stejný počet zhlukov, ako pri K-means algoritme. Porovnajte výsledky a jednotlivé zhluky. Riešenie umiestnite na svoju webovú stránku, najneskôr v utorok, 11.05.2021, do 14:00.