12  Metódy neasistovaného učenia

12.1 Úvod

  • V neasistovanom (štatistickom/strojovom) učení odozva Y nie je definovaná, všetky vlastnosti pozorovaného objektu sú rovnocenné náhodné premenné X_1,\ldots,X_p.
  • Cieľom analýzy tak už nie je predpovedať hodnotu (alebo rozdelenie premennej), ale skôr objavovať zaujímavé veci, napr. ako vhodne zobraziť údaje, alebo nájdenie skupín medzi pozorovanými objektami (riadky datasetu) či ich vlastnosťami (stĺpce datasetu).
  • Preto sa neasistované učenie často využíva v rámci prieskumnej analýzy údajov (exploratory data analysis, EDA).
  • Má tendenciu byť viac subjektívna, a pretože nemá jeden cieľ (predpoveď), je ťažšie vyhodnotiť úspešnosť.
  • Základným prístupom ku hľadaniu skupín medzi pozorovaniami je zhluková analýza (clustering). Naopak, hľadaniu skupín užitočných a menej užitočných vlastností sa venuje oblasť metód redukcie rozmernosti (dimensionality reduction, ordination), ktorá je užitočné aj pre optimálne zobrazenie údajov.
  • Redukcia rozmernosti v širšom kontexte zahŕňa
    • výber vlastností (feature selection, častejšie v asistovanom učení – napr. regularizačná technika LASSO) a
    • extrahovanie vlastností (feature extraction).
  • Základnou metódou používanou na extrakciu užitočných vlastností v oblasti neasistovného učenia je analýza hlavných komponentov (principal component analysis, PCA).1

12.2 Analýza hlavných komponentov

12.2.1 Princíp

  • Vlastnosti X_1,\ldots,X_p tvoria osi p-rozmerného priestoru a skúmané objekty (štatistické jednotky) sa v ňom prostredníctvom pozorovaných hodnôt ich vlastností x_{ij} (kde i=1,\ldots,n, a j=1,\ldots,p) zobrazia ako body.
  • Ak by sme aj dokázali vidieť vo viac ako 3-rozmernom priestore tvar mraku bodov (charakter skúmaného javu), potrebujeme nájsť spôsob, ako túto informáciu sprostredkovať ďalším ľuďom (či už obrázkom alebo modelom).

  • Hľadáme podpriestor, do ktorého projekcia “najlepšie” reprezentuje mrak bodov.
  • V prípade PCA, os Z “najlepšieho” (1D) podpriestoru vedie v smere najväčšieho rozptylu bodov.
  • Ekvivalentne, kolmá vzdialenosť bodov od tohto podpriestoru je najmenšia.

  • Os Z sa nazýva hlavný (principal) komponent a je lineárnou kombináciou osí X_j pôvodneho priestoru, Z = v_1 X_1 + \ldots + v_pX_p pričom pre vektor \mathbf{v} tzv. záťaží (loadings) platí ||\mathbf{v}||=1. Fixovaním dĺžky vektora sa zafixuje aj rozptyl hlavného komponentu, čo je nevyhnutné pre nájdenie maxima.

12.2.2 Optimalizácia

  • Nech stĺpce matice \mathbf{X} obsahujú (centrované) súradnice bodov v pôvodnom priestore a vektor \mathbf{z} = \mathbf{X}\mathbf{v} obsahuje polohu bodov na osi Z. V žargóne sa tieto nové súradnice nazývajú scores (zásahy?).
  • Vyjadrime rozptyl novej premennej Z, ktorý chceme maximalizovať: \begin{split} var(Z) &= E[Z^2] \\ & = \frac1n\mathbf{z}^T\mathbf{z} \\ & = \frac1n\mathbf{v}^T\mathbf{X}^T\mathbf{X}\mathbf{v} \\ & = \mathbf{v}^T\left(\frac1n\mathbf{X}^T\mathbf{X}\right)\mathbf{v} \\ & = \mathbf{v}^T\mathbf{\Sigma}\mathbf{v} \end{split} kde \mathbf{\Sigma} je kovariančná matica náhodného vektora (X_1,...,X_p).
  • Táto nová premenná pritom nemusí byť jediná, vo všeobecnosti sa ich dá nájsť rovnaký počet ako je vlastností. Označíme ich (Z_1,...,Z_p) alebo \{PC_j\}_{j=1}^p (kde j=1,\ldots,p), každému zodpovedá smerový vektor \mathbf{v}_j a sú na seba kolmé (rovnako osi X_j), pretože z geometrického hľadiska ide o obyčajnú rotáciu súradnicového systému.

  • Rozmernosť podpriestoru (počet hlavných komponentov) je subjektívna voľba, hlavné komponenty sú zoradené podľa variability, ktorú “vysvetľujú”.

  • Vektor záťaží \mathbf{v}_1 prvého hlavného komponentu maximalizuje celkový rozptyl s uplatnením jedinej podmienky, \mathbf{v}_1 = \arg\max\limits_{\mathbf{v}}\left(\mathbf{v}^T\mathbf{\Sigma}\mathbf{v}\right) \quad\text{pre}\quad \mathbf{v}^T\mathbf{v}=1 Uplatnením metódy Lagrangeových multiplikátorov \begin{split} L(\mathbf{v}) = \mathbf{v}^T\mathbf{\Sigma}\mathbf{v}-\lambda(\mathbf{v}^T\mathbf{v}-1) & \\ \frac{\partial L}{\mathbf{v}} = 2\mathbf{\Sigma}\mathbf{v}-2\lambda\mathbf{v} \qquad\qquad & \\ \mathbf{\Sigma}\mathbf{v}_1-\lambda_1\mathbf{v}_1 &= 0 \\ \mathbf{\Sigma}\mathbf{v}_1 &= \lambda_1\mathbf{v}_1 \end{split} dostaneme riešenie, ktorým je (normovaný) vlastný vektor zodpovedajúci najväčšiemu vlastnému číslu (\lambda_1) matice \mathbf{\Sigma}.
  • Vektor záťaží \mathbf{v}_2 druhého hlavného komponentu musí pri maximalizácii spĺňať okrem jednotkovej dĺžky aj podmienku kolmosti na \mathbf{v}_1, čiže \mathbf{v}_2 = \arg\max\limits_{\mathbf{v}}\left(\mathbf{v}^T\mathbf{\Sigma}\mathbf{v}\right) \quad\text{pre}\quad \mathbf{v}^T\mathbf{v}=1 \quad\text{a}\quad \mathbf{v}_1^T\mathbf{v}=0 Riešenie vo forme (normovaného) vlastného vektora matice \mathbf{\Sigma}, ktorý zodpovedá druhému najväčšiemu vlastnému číslu (\lambda_2), spĺňa aj druhú podmienku, pretože vlastné vektory tvoria ortonormálnu bázu.
  • Takto môžeme hlavné komponenty nájsť všetky naraz, a to diagonalizáciou kovariančnej matice \mathbf{\Sigma}, čiže rozkladom (eigenvalue decomposition) \mathbf{\Sigma} = \mathbf{V}\mathbf{\Lambda}\mathbf{V}^T kde
    • \mathbf{\Lambda} je diagonálna matica s vlastnými číslami (\lambda_1,\ldots,\lambda_p) na diagonále, predstavuje kovariančnú maticu hlavných komponentov,
    • \mathbf{V}=(\mathbf{v}_1,\ldots,\mathbf{v}_p) je matica s vlastnými vektormi v stĺpcoch, pričom tie sú ortonormálne, takže platí \mathbf{V}^T\mathbf{V}=1.
  • V dvojrozmernom prípade je \mathbf{V} jednoduchá matica rotácie o uhol \theta,

\mathbf{V} = \begin{pmatrix} \cos\theta & -\sin\theta \\ \sin\theta & \cos\theta \end{pmatrix}

  • Symetrické matice majú p reálnych vlastných čísel (vrátane duplicít), ktorým zodpovedá p navzájom lineárne nezávislých vlastných vektorov. Množina vlastných čísel sa nazýva aj spektrum matice.
  • Smer vektorov záťaží \mathbf{v}_1,\ldots,\mathbf{v}_p je unikátny, ale ich orientácia nie, dve implementácie metódy teda môžu dať hodnoty s opačnými znamienkami.
  • Ak sú vlastnosti X_1,\ldots,X_p pozorované v rôznych jednotkách, je okrem centrovania potrebné aj ich preškálovanie a to podelením štandardnou odchýlkou. V opačnom prípade by vlastnosti s najväčším rozptylom zásadne ovplyvnili transformáciu.

12.2.3 Aplikácie

  • Transformáciu medzi vlastnosťami (features) v matici \mathbf{X} a zásahmi (scores) v matici \mathbf{Z}=(\mathbf{z}_1,\ldots,\mathbf{z}_p) realizujeme pomocou záťaží (loadings) \mathbf{V} pomocou vzťahu \begin{split} \underset{n\times p}{\mathbf{Z}} &= \underset{n\times p}{\mathbf{X}}\cdot\underset{p\times p}{\mathbf{V}} \\ \mathbf{X} &= \mathbf{Z}\cdot\mathbf{V}^T \end{split} takže pri použití všetkých hlavných komponentov dokážeme obnoviť pôvodné pozorovania vlastností.
  • Zmysel PCA je však v redukcii rozmernosti, preto sa v aplikáciách nepoužívajú všetky PC, ale iba tie “najvýznamnejšie”.
  • Ak je treba zobraziť mnohorozmerný priestor, zvyčajne siahneme po prvých dvoch, nanajvýš prvých troch PC.
  • PCA sa dá prirodzene použiť na filtrovanie dát (odstránenie šumu), \underset{n\times p}{\mathbf{X}} = \underset{n\times q}{\mathbf{Z}}\cdot\underset{q\times p}{\mathbf{V}^T} + \underset{n\times p}{\mathbf{R}} \qquad\text{pre}\quad q<p kde matica \mathbf{R} obsahuje rezíduá.
  • Hlavné komponenty môžu poslúžiť ako vysvetľujúce premenné v klasickej lineárnej regresii (principal component regression, PCR) Y = \beta_0 + \beta_1Z_1 + \ldots + \beta_qZ_q + \varepsilon za predpokladu, že smer najväčšej variability vlastností dokáže vysvetliť variabilitu odozvy. Aj keď to nebýva splnené úplne, predpovedná schopnosť takého modelu je často dobrá a zároveň rieši problém multikolinearity. Počet hlavných komponentov zaradených do regresie sa určí napr. metódou krížovej validácie.
  • Pozn.: Smery vysvetľujúcich premenných v PCR sú určené neasistovaným prístupom. Asistovanú alternatívu poskytuje metóda parciálnych najmenších štvorcov (partial least squares, PLS).
  • Iná aplikácia sa núka v oblasti dopočítania chýbajúcich hodnôt (missing value imputation). Tie však
    • musia byť náhodne chýbajúce (napr. vybité baterky v elektrickej váhe),
    • nie systematicky (neschopnosť pacienta zúčastniť sa merania v dôsledku extrémnej nadváhy).
    Štandardný postup, v ktorom sa chýbajúce hodnoty nahradia priemerom v danom stĺpci (vlastnosti), zanedbáva korelácie medzi stĺpcami. Môže však poslúžiť ako počiatočný stav \mathbf{X}^{(0)} v nasledujúcom algoritme (s indexom iterácie i):
    1. pomocou \mathbf{X}^{(i-1)} sa vypočítajú matice \mathbf{Z} a \mathbf{V} (pre prvých q komponentov),
    2. z nich pomocná matica \mathbf{X} a
    3. z tej sa doplnia hodnoty na chýbajúcich pozíciách do \mathbf{X}^{(i-1)}, takže vznikne aktualizovaná matica \mathbf{X}^{(i)}.
  • Týmto spôsobom sa dajú “dopočítať” napr. aj preferenčné matice, v ktorých riadky predstavujú zákazníkov, stĺpce zodpovedajú produktom (napr. film) a samotné prvky obsahujú hodnotenie produktu zákazníkom. Také matice sú veľmi riedke (nutne chýbajúce hodnoty), ale ak je dostatočný prekryt riadkov (každý film videlo viac zákazníkov), odporúčací algoritmus (recommender system) zákazníkov zaradí do “skupiniek” a filmy do “žánrov”. Jedna skupinka obľubuje jeden žáner (spolu je ich toľko ako hlavných komponentov), potom zásahy (scores) reprezentujú silu príslušnosti zákaznikov do skupiniek, a záťaže (loadings) reprezentujú intenzitu, s akou film patrí do žánru.

12.2.4 Určenie počtu PC

  • Čím menej hlavných komponentov by sme mohli použiť, tým ľahšie by sa dali interpretovať (a my by sme ľahšie porozumeli dátam). Koľko “málo” je optimálne?
  • Naneštastie neexistuje jediná alebo jednoduchá odpoveď.
  • V rozhodnutí však pomôže tzv. scree plot (suťový graf?). Na horizontálnej osi má index (poradie) hlavného komponentu, na vertikálnej osi y môže mať buď priamo varianciu \lambda_k vysvetlenú k-tym komponentom, alebo jej pomer ku celkovej variancii (proportion of variance explained, PVE) PVE(k) = \frac{\lambda_k}{\sum_{j=1}^p\lambda_j}
  • Za q zvolíme také najväčšie k,
    • pre ktoré \lambda_k'>1 (Kaiser criterion, \lambda_k' je vl.číslo korelačnej matice), alebo
    • od ktorého je pokles PVE výrazne pomalší (“lakeť” klesajúcej funkcie).

Príklad (zločiny)
  • Štatistický súbor údajov datasets::USArrests o počte zatknutí na 100 tis. obyvateľov v roku 1973 v 50 štátoch USA podľa typu zločinu (vražda, napadnutie, znásilnenie). Dataset navyše obsahuje podiel obyvateľov žijúcich v mestách (UrbanPop).
  • Zobrazenie datasetu je už aj pre 4 premenné pomerne náročné.
Kód 
plotly::plot_ly(USArrests, 
                x = ~Murder, y = ~Assault, z = ~Rape,
                color = ~UrbanPop,
                size = 2) |> 
  plotly::add_markers()
Kód 
c("mean", "var") |> 
  sapply(function(x) apply(USArrests, 2, FUN = x)) |> 
  t() |> round(1)
     Murder Assault UrbanPop Rape
mean    7.8   170.8     65.5 21.2
var    19.0  6945.2    209.5 87.7
  • Rozptyl počtu napadnutí (Assault) je oveľa väčší ako u ostatných. Hoci ide o bezrozmerné premenné, škálovanie na rovnaký rozptyl je nutné, inak by bol prvý hlavný komponent takmer výhradne v smere napadnutia a nič zaujímavé by sme sa nedozvedeli.
  • V Rku na analýzu PC slúži napr. funkcia stats::prcomp. Jej výstup obsahuje
    • (sdev) štandardnú odchýlku vysvetlenú jednotlivými PC,
    • (rotation) maticu záťaží \mathbf{V}, ktorou sa rotuje súradnicový systém vlastností do súr.systému hlavných komponentov,
    • (center) centroid datasetu,
    • (scale) štandardnú odchýlku jednotlivých vlastností, ktorou sa škálujú na rovnakú mierku a
    • (x) maticu zásahov \mathbf{Z}.
    Zobrazíme maticu záťaží a časť matice zásahov:
Kód 
# manuálny výpočet:
# USArrests |> cor() |> eigen() |> getElement("values") |> sqrt()  # std.dev.
# V <- USArrests |> cor() |> eigen() |> getElement("vectors")  # loadings
# Z <- USArrests |> as.matrix() |> apply(2, scale) |> (`%*%`)(V)  # scores

# automatika:
fit <- prcomp(USArrests, scale = TRUE)

fit$rotation |> round(3)
            PC1    PC2    PC3    PC4
Murder   -0.536 -0.418  0.341  0.649
Assault  -0.583 -0.188  0.268 -0.743
UrbanPop -0.278  0.873  0.378  0.134
Rape     -0.543  0.167 -0.818  0.089
Kód 
fit$x |> head() |> round(2)
             PC1   PC2   PC3   PC4
Alabama    -0.98 -1.12  0.44  0.15
Alaska     -1.93 -1.06 -2.02 -0.43
Arizona    -1.75  0.74 -0.05 -0.83
Arkansas    0.14 -1.11 -0.11 -0.18
California -2.50  1.53 -0.59 -0.34
Colorado   -1.50  0.98 -1.08  0.00
  • Dáta s osami oboch systémov znázorníme pomocou biplotu.
Kód 
biplot(fit, scale = 0, cex = 0.7)  # scale=0 aby šípky zodpovedali záťažiam

  • Vidíme, že všetky druhy zločinu majú v prvom PC približne rovnaké váhy, zatiaľčo podiel mešťanov iba polovičnú. To znamená, že PC1 približne zodpovedá celkovej miere závažných zločinov.
  • Druhý PC naložil väčšinu záťaže na UrbanPop, takže predstavuje stupeň urbanizácie jednotlivých štátov.
  • Premenné zodpovedajúce zločinom sú sústredené blízko pri sebe, to znamená, že sú navzájom značne korelované: štáty s vysokou kriminalitou jedného druhu majú problém aj s inými druhmi kriminality. S urbanizovanosťou súvisia menej.
  • z grafu zásahov je vidieť, že najvyššiu kriminalitu mala Kalifornia, Nevada a Florida. Kalifornia má tiež vysoký podiel mestského obyvateľstva.
  • Podľa suťového grafu by nám stačil prvý, nanajvýš aj s druhým hlavným komponentom.
Kód 
screeplot(fit, type = "lines", main = "")
abline(h = 1, lty = 2)
plot(fit$sdev^2 / sum(fit$sdev^2), type = "b", ylab = "PVE", xlab = "PC")

  • Prvý PC vysvetľuje asi 62% variability, druhý už iba 25%, spolu teda 87% celkovej variablity údajov.

12.2.5 Použitá literatúra

12.3 Zhluková analýza

12.3.1 Úvod

12.3.1.1 Princíp

  • Zhlukovanie (angl. clustering) označuje širokú množinu metód neasistovaného učenia určených na hľadanie podskupín (zhlukov) v súbore údajov.
  • Zvyčajne sa hľadajú zhluky pozorovaní (riadky tabuľky údajov), no môže sa použiť aj na hľadanie skupín medzi vlastnosťami (stĺpce).
  • Objekty (štatistické jednotky - riadky, alebo vlastnosti - stĺpce) v skupine sú si navzájom “celkom podobné” (within-cluster homogeneity) zatiaľčo medzi skupinami sa “dosť líšia” (between-cluster heterogeneity).
  • Napríklad
    • v marketingu je cieľom objaviť skupiny zákazníkov na cielenie konkrétnej reklamy,
    • v medicíne sa hľadajú skupiny pacientov vhodných pre špecifické liečebné postupy, resp. zhlukovaním sú predbežne diagnostikované varianty ochorenia.
    • v astronómii je záujem zoskupiť podobné hviezdy a galaxie,
    • v genetike zas objaviť skupiny génov s podobným procesom expresie (prenosu informácie na RNA alebo proteíny).

12.3.1.2 Podobnosť a odlišnosť

  • Aby sa zhluky dali od seba odlíšiť, je potrebné definovať mieru “podobnosti”, resp. “odlišnosti”.
  • Niektoré miery vzdialenosti už poznáme (euklidovskú, manhattanskú, Mahalanobisovu,…), z nich napr. euklidovská je síce často používaná, ale je aj citlivá na odľahlé hodnoty.
  • Pri numerických vlastnostiach daných v rôznej mierke sa hodí vzdialenosť založená na korelácii.

  • V text mining-u sa bežne používa kosínusová vzdialenosť.
  • Gowerova miera je vhodná pre zmiešané údaje (mixed data):
    • pre kvantitatívnu (intervalovú) premennú X je definovaná ako manhattanská miera normovaná rozsahom, d_{i,j} = \frac{|x_i-x_j|}{\max_l(x_l)-min_l(x_l)}
    • ordinálna je najprv prekódovaná na poradia (s korekciou pre zhody),
    • nominálna premenná je prekódovaná na pomocné binárne premenné a z nich vypočítaný kockový (Sorensen-dice) koeficient podobnosti S = \frac{2a}{2a+b}, kde a je počet zhôd na hodnote 1, b je počet nezhôd (a zhody na 0 sa nerátajú). Miera vzdialenosti (nepodobnosti) je potom doplnok do 1. Celková vzdialenosť je priemerom čiastkových vzdialeností zo všetkých uvažovaných vlastností (kvantitatívných a kvalitatívnych). Môže byť vážený pre korekciu neproporcionality medzi počtom jedničiek kvalitatívnych voči kvantitatívnym.

12.3.1.3 Metódy

  • Zhlukovacie metódy sa delia podľa toho, či sú objekty do zhlukov zaradené výlučne (jednoznačná príslušnosť, hard methods) alebo s určitou mierou príslušnosti (soft, fuzzy methods).
  • Bežnejšie sú metódy z prvej skupiny. Tam sa ďalej líšia podľa spôsobu separácie zhlukov na
    • metódy s priamym delením, napr. k-means alebo k-medoids clustering,
    • metódy s hierarchickým delením.

  • Hard metódy používajú heuristické algoritmy, ktoré sa musia nejak vysporiadať s komplexnosťou problému zaradenia n objektov do K (výlučných) skupín.
    • Počet možných delení do K tried je tzv. Stirlingovo číslo druhého druhu {n \brace K} = \frac1{K!}\sum_{k=0}^K (-1)^k{K \choose k}(K-k)^n,
    • a počet všetkých delení n objektov zas tzv. Bellovo číslo B_n=\sum_{k=0}^n{n \brace k}=\frac1e\sum_i^\infty\frac{i^n}{i!}. Napríklad B_{10}=115975, no už B_{100}=4.8\times10^{115} je viac, než počet atómov v pozorovateľnom vesmíre (cca 10^{80}).
  • Zo soft metód sú obľúbené tie, ktoré modelujú zmes (normálnych) rozdelení, angl. Gaussian mixture models.

12.3.2 K-means

12.3.2.1 Princíp

  • Jednoduchá a elegantná metóda.

  • Počet tried K treba definovať vopred.

  • Matematicky je problém formulovaný nasledovne: ak C_1,\ldots,C_K označujú množiny indexov objektov (pozorovaní) v každom zhluku, pričom

    • C_1\cup C_2\cup\ldots\cup C_K=\{1,\ldots,n\},
    • C_k\cap C_{k'}=\emptyset pre všetky k\neq k',

    čiže každý objekt patrí práve do jedného zhluku, potom cieľom je minimalizovať vnútro-zhlukový rozptyl W, \underset{C_1,\ldots,C_K}{minimize}\left\{\sum_{k=1}^K W(C_k)\right\},\qquad\text{kde}\quad W(C_k)=\frac1{n_k}\sum_{i\in C_k}\sum_{j=1}^p(x_{ij}-g_{kj})^2, g_{kj}=\frac1{n_k}\sum_{i\in C_k}x_{ij} a použitá bola euklidovská vzdialenosť.

  • Nájsť presné riešenie je veľmi náročné.

  • Našťastie existuje jednoduchý algoritmus, ktorým sa dá nájsť celkom dobré riešenie, aj keď predstavuje iba lokálne optimum.

  • S danými počiatočnými hodnotami centroidov \mathbf{g}_1,\ldots,\mathbf{g}_K sa opakujú nasledujúce kroky:

    1. zaradenie objektov do zhluku podľa najmenšej vzdialenosti od centroidu,
    2. prepočítanie centroidov,

    až kým sa zaradenie v bode 1 nezmení.

  • Počiatočné hodnoty sa určia buď

    • (Forgy-ho metódou) náhodným stotožnením centroidov s K objektami, alebo
    • náhodným zaradením objektov do K tried a výpočtom ich centroidov.
Ilustrácia (Forgy)

  • Algoritmus vo všeobecnosti nájde skôr lokálne, než globálne minimum. Výsledok závisí od počiatočných (náhodných) podmienok.
  • Preto je dôležité nechať odhad zbehnúť viackrát, vždy s inými počiatočnými hodnotami centroidov. Nakoniec sa vyberie najlepšie riešenie.

12.3.2.2 Alternatívy

  • Ak dáta obsahujú odľahlé hodnoty, alternatívou ku k-priemerom môže byť metóda k-medoidov (partitioning around medians, PAM).
  • Robustnosť je však kompenzovaná vyššou výpočtovou náročnosťou. V Rku sa použije funkcia cluster::pam
  • Súradnice medoidov sú mediány jednotlivých vlastností, navyše sú vždy totožné s niektorými z pozorovaní.)
  • S nárastom počtu objektov sa zhlukovanie spomaľuje. Preto sa pre veľké datasety oplatí použiť algoritmus (clustering large applications, CLARA), ktorý rozdelí dataset na viacero menších (rovnako veľkých) vzoriek. Na každej aplikuje PAM, získa medoidy ale zhluky urobí nad celým datasetom a získa celkovú mieru vnútrozhlukovej vzdialenosti. Vyhrávajú medoidy (tej vzorky), pre ktoré je táto miera najmenšia.
  • Algoritmus je implementovaný vo funkcii cluster::clara.
  • Klasické zhlukovanie predpokladá konvexné hranice zhlukov. Komplikovanejšie štruktúry vyžadujú komplikovanejšie metódy, napr. spektrálne zhlukovanie. V Rku ich implementuje napr. balík kernlab.

Príklad (gejzír)
  • Dataset MASS::geyser s dĺžkou trvania gejzírov (duration) a čakacou dobou na erupciu.
Kód 
# skrátenie výpisu z objektu vráteného funkciou kmeans
print.kmeans <- function (x, ...) 
{
    cat("K-means clustering with ", length(x$size), " clusters of sizes ", 
        paste(x$size, collapse = ", "), "\n", sep = "")
    cat("\nCluster means:\n")
    print(x$centers, ...)
    cat("\nClustering vector:\n")
    print(head(x$cluster,10), ...)
    cat("...\n")
    cat("\nWithin cluster sum of squares by cluster:\n")
    print(x$withinss, ...)
    ratio <- sprintf(" (between_SS / total_SS = %5.1f %%)\n", 
        100 * x$betweenss/x$totss)
    cat(sub(".", getOption("OutDec"), ratio, fixed = TRUE), "Available components:\n", 
        sep = "\n")
    print(names(x))
    if (!is.null(x$ifault) && x$ifault == 2L) 
        cat("Warning: did *not* converge in specified number of iterations\n")
    invisible(x)
}
Kód 
dat <- MASS::geyser
fit <- kmeans(dat, centers = 3, nstart =  10)
print(fit)  # modifikovaná metóda pre skrátenie niektorých dlhších výpisov
K-means clustering with 3 clusters of sizes 123, 77, 99

Cluster means:
   waiting duration
1 76.54472 3.218835
2 88.02597 2.589827
3 54.83838 4.438889

Clustering vector:
 1  2  3  4  5  6  7  8  9 10 
 1  1  3  1  1  1  3  2  1  3 
...

Within cluster sum of squares by cluster:
[1] 2020.012 1444.276 2819.925
 (between_SS / total_SS =  89.1 %)

Available components:

[1] "cluster"      "centers"      "totss"        "withinss"     "tot.withinss"
[6] "betweenss"    "size"         "iter"         "ifault"
Kód 
dat |> 
  transform(cluster = as.factor(fit$cluster)) |> 
  ggplot() + aes(x = waiting, y = duration) +
  geom_point(aes(color = cluster)) + 
  geom_point(data = as.data.frame(fit$center), shape = 2)

  • Zhluková analýza je zjavne ovplyvnená rozdielnym rozsahom časov, preto ju zopakujeme pre normované dáta.
Kód 
# inverse to function scale()
# - x is matrix or data frame to be unsclaed,
# - scaled is the output of scale()
unscale <- function(x, scaled = x) {
  par <- attributes(scaled)[c("scaled:center", "scaled:scale")] |>
    setNames(c("mean","sd")) 
  n <- nrow(x)
  x * rep(par$sd, each=n) + rep(par$mean, each=n)
}

# clustering on scaled data
sdat <- MASS::geyser |> 
  scale()
fit <- kmeans(sdat, centers = 3, nstart =  10)
# plot unscaled data
centroids <- fit$centers |>  
  unscale(sdat) |> 
  as.data.frame() |> 
  tibble::rownames_to_column("cluster")
dat |> 
  transform(cluster = as.factor(fit$cluster)) |>
  ggplot() + aes(x = waiting, y = duration) +
  geom_point(aes(color = cluster)) + 
  geom_point(data = centroids, shape = 13, size = 3, show.legend = FALSE)

12.3.3 Hierarchické zhlukovanie

12.3.3.1 Dendrogram

  • Nevýhodou metódy k-centroidov (k-means) je nutnosť špecifikovať počet zhlukov K vopred.
  • Hierarchické zhlukovanie toto nevyžaduje, navyše výstupom je prehľadný stromový diagram, ktorý sa volá dendrogram.
  • Každý list predstavuje jeden objekt, navzájom sú spájané vetvami v rôznej výške. Vertikálna súradnica spojov predstavuje mieru vzdialenosti (odlišnosti) skupín objektov.
Kód 
dat1 <- data.frame(
  X1 = c(-1.5, -1.4, -1.0, -0.6, -0.1, 0.1, 0.7, 1.3, 1.5),
  X2 = c(-0.4, -1.5, -1.1, -1.0, 0.9, -0.8, -0.2, -0.3, 0.0)
  )

plot(X2 ~ X1, dat1, asp = 1, type = "n")
text(X2 ~ X1, dat1, labels = row.names(dat1))
fit <- dat1 |> 
  dist() |>  # distance matrix calculation
  hclust()  # hierarchical clustering
plot(fit, ann = FALSE)

  • Odlišnosť dvoch objektov tak môžeme určiť podľa výšky, v ktorej sa ich vetvy spájajú. Pozor, nemôžeme súdiť na základe horizontálnej vzdialenosti!
  • Identifikácia zhlukov v dendrograme prebieha pomocou rezov v konkrétnej výške.
Kód 
# generate data
set.seed(1)
dat <- list(A = c(-5,0), B = c(0,0), C = c(1,3)) |> 
  lapply(mvtnorm::rmvnorm, n = 15, sigma = diag(2)) |> 
  lapply(as.data.frame) |> 
  lapply(setNames, nm = paste0("X",1:2)) |> 
  dplyr::bind_rows(.id = "class") |> 
  dplyr::mutate(class = as.factor(class))

# fit model  
fit <- dat |> 
  dplyr::select(-class) |> 
  dist() |> 
  hclust()

# --- visualize result

# basic dendrogram
#fit |> 
#  plot(hang = -1, cex = 0.6, main="", ylab="")

# cut to get 2 clusters
fit |> 
  as.dendrogram() |> 
  dendextend::set("labels_col", 2:3, k=2) |> 
  dendextend::set("labels_cex", 0.6) |> 
  plot()
abline(h = 8, lty = 2)
# plot with dots or labels
plot(X2 ~ X1, dat,
     col = cutree(fit, k = 2) + 1, 
     pch = 16, cex=1.5)
#plot(X2 ~ X1, dat, type = "n")
#text(X2 ~ X1, dat, 
#     col = cutree(fit, k = 3) + 1, 
#     labels = rownames(dat))


# cut to get 3 clusters
fit |> 
  as.dendrogram() |> 
#  dendextend::set("branches_k_col", 2:4, k=3) |> 
  dendextend::set("labels_col", 2:4, k=3) |> 
  dendextend::set("labels_cex", 0.6) |> 
  plot()
abline(h = 5, lty = 2)
# plot with dots or labels
plot(X2 ~ X1, dat,
     col = cutree(fit, k = 3) + 1, 
     pch = 16, cex=1.5)
#plot(X2 ~ X1, dat, type = "n")
#text(X2 ~ X1, dat, 
#     col = cutree(fit, k = 3) + 1, 
#     labels = rownames(dat))


# use package dendextend to color leaves and possibly branches
# https://talgalili.github.io/dendextend/articles/dendextend.html
# http://www.sthda.com/english/wiki/beautiful-dendrogram-visualizations-in-r-5-must-known-methods-unsupervised-machine-learning

  • Výška rezu tak slúži rovnakému účelu ako voľba počtu zhlukov v metóde k-means.
  • Hierarchické zhlukovanie nemusí vždy zodpovedať realite. Napr. ak by naše pozorovania tvorili dve skupiny podľa pohlavia a tri skupiny podľa národnosti, potom nedáva zmysel, aby tretia skupina vznikla rozdelením jednej z tých dvoch. Vtedy môže hierarchické zhlukovanie dávať menej presné výsledky ako k-means.

12.3.3.2 Metódy

  • Cieľom hierarchického zhlukovania je nájsť vnorené delenia objektov.
  • Na to existujú dva prístupy:
    • zoskupujúci (agglomerative)
    • rozkladný (divisive)
  • Aglomeratívna metóda AGNES (AGglomerative NESting) začína zdola:
    1. Každý objekt je na začiatku považovaný za samostatný zhluk (list).
    2. V každom kroku iteračného algoritmu sa spoja dva najpodobnejšie zhluky.
    3. Táto procedúra skončí, keď sú všetky objekty súčasťou jedného zhluku (koreň).
  • Divizívna metóda DIANA (DIvisive ANAlysis clustering) začína z vrchu hierarchie (od koreňa). V každom kroku iteračného algoritmu je zhluk rozdelený na také dva menšie, ktoré sú zo všetkých možných čo najviac heterogénne. Cyklus končí vtedy, keď každý zhluk obsahuje jediný objekt.

  • Hierarchické zhlukovanie je implementované v dvoch predinštalovaných balíkoch Rka.
    • Aglomeratívne vo funkciách stats::hclust a cluster::agnes.
    • Divizívne v cluster::diana

12.3.3.3 Vzdialenosť zhlukov

  • Odlišnosť jednotlivých objektov už vieme merať, napr. euklidovskou alebo manhattanskou vzdialenosťou.
  • Ale ako vyjadriť rozdielnosť dvoch zhlukov?
  • Na to vzniklo niekoľko metód (linkage):
    • maximum (complete linkage) - Najväčšia vzdialenosť medzi objektami z jedného a z druhého zhluku. Má tendenciu tvoriť zhluky rovnakého priemeru (diameter, najv.vzdialenosť), je citlivá na odľahlé hodnoty. Takže objekty môžu byť oveľa bližšie ku cudzím než ku vlastným.
    • minimum (single linkage) - Najmenšia vzdialenosť medzi objektami z jedného a z druhého zhluku. Má tendenciu tvoriť reťaze (pridávať prvky po jednom alebo asociovať také zhluky, ktoré sú prepojené reťazou).
    • priemer (average linkage) - Priemerná vzdialenosť medzi objektami z jedného a z druhého zhluku. Kompromis medzi predošlými. Avšak mení sa s monotónnou transformáciou vzdialenosti objektov, zatiaľčo max a min závisia len od poradia.
    • centroid - Vzdialenosť medzi centroidmi. Môže viesť k inverziám, kedy dva zhluky sú spojené v menšej výške, než sú jednotlivo (problém s vizualizáciou a interpretáciou).
    • Ward - Minimalizuje vnútro-zhlukový rozptyl a tak produkuje kompaktné zhluky.
Ilustrácia (linkage)
Kód 
plot(X2 ~ X1, dat1, asp = 1, type = "n")
text(X2 ~ X1, dat1, labels = row.names(dat1))
dat1 |> 
  dist() |>  # distance matrix calculation
  hclust(method = "complete")  |> # hierarchical clustering
  plot(xlab = "", ylab = "", main = "max (complete)", sub = "")
dat1 |> 
  dist() |>  # distance matrix calculation
  hclust(method = "single")  |> # hierarchical clustering
  plot(xlab = "", ylab = "", main = "min (single)", sub = "")
dat1 |> 
  dist() |>  # distance matrix calculation
  hclust(method = "average")  |> # hierarchical clustering
  plot(xlab = "", ylab = "", main = "average", sub = "")
dat1 |> 
  dist() |>  # distance matrix calculation
  hclust(method = "centroid")  |> # hierarchical clustering
  plot(xlab = "", ylab = "", main = "centroid", sub = "")
dat1 |> 
  dist() |>  # distance matrix calculation
  hclust(method = "ward.D")  |> # hierarchical clustering
  plot(xlab = "", ylab = "", main = "Ward", sub = "")

  • Funkcia cluster::agnes vypočíta aj hodnotu tzv. aglomeračného koeficientu (AC), ktorý v intervale [0,1] vyjadruje silu zhlukovej štruktúry. Vyvážené zhluky majú hodnotu blízke 1. Nevýhodou je, že rastie s n, takže sa nedajú porovnať dáta príliš rozdielnej dĺžky.
Kód 
c("complete", "single", "average", "ward") |> 
  sapply(function(x) cluster::agnes(dat1, method = x)$ac) |> 
  round(2) |> 
  rbind(AC = _)
   complete single average ward
AC     0.73   0.54    0.65  0.8
Príklad (Agnes)
  • Dataset MASS::geyser s dĺžkou trvania gejzírov (duration) a čakacou dobou na erupciu.
Kód 
dat <- MASS::geyser
fit <- dat |> 
  cluster::agnes(stand = TRUE, method = "ward")
plot(fit, which = 2, labels = FALSE)
abline(h=15, lty=2)
#fit |> as.dendrogram() |> plot()

dat |> 
  transform(
    cluster = fit |> cutree(k = 3) |> as.factor()
    ) |> 
  ggplot() + aes(x = waiting, y = duration) +
  geom_point(aes(color = cluster))

12.3.4 Zmes normálnych rozdelení

12.3.4.1 Princíp

  • Tradičné zhlukovacie metódy ako k-means a hierarchické sú heuristické. Zhluky tvoria priamo z dát, teda bez použitia nejakej pravdepodobnostnej miery.
  • Metódy založené na modeli pravdepodobnosti poskytujú neostrú hranicu medzi zhlukmi (soft assignment), každý objekt má priradenú pravdepodobnosť príslušnosti ku každému zhluku.
  • Tento prístup navyše rieši problém identifikácie optimálneho počtu zhlukov automaticky (implicitne).
  • Najpopulárnejšia je trieda modelov so zmiešaným gaussovským rozdelením (Gaussian mixture model, GMM).
  • GMM predpokladá, že každý objekt (jeho pozorované vlastnosti) pochádza z jedného z K (p-rozmerných) normálnych rozdelení s hustotou f_k(\mathbf{x})=\frac1{(2\pi)^{p/2}|\mathbf{\Sigma}_k|^{1/2}}\exp\left(-\frac12(\mathbf{x}-\mathbf{\mu}_k)^T\mathbf{\Sigma}_k^{-1}(\mathbf{x}-\mathbf{\mu}_k)\right) kde \mathbf{\mu}_k je centroid a \mathbf{\Sigma}_k kovariančná matica k-teho zhluku.
  • Potom GMM je hustotou zmesi všetkých K normálnych rozdelení, f(\mathbf{x}) = \sum_{k=1}^K\pi_k f_k(\mathbf{x}) kde \{\pi_k\} sú váhy.
  • Váhy spolu s centroidmi a kovariančnými maticami tvoria parametre GMM.
Príklad (kyslosť)
  • Príklad: dataset mclust::acidity, rozdelenie pravdepodobnosti indexu kyslosti vody v 155 jazerách v USA.
Kód 
# 1D
fit <- mclust::acidity |> 
  mclust::densityMclust(plot = FALSE)
fit |> 
  plot(what = "density", data = mclust::acidity, breaks = 15, xlab = "acidity")
fit |> 
  plot(what = "diagnostic", type = "cdf")

Príklad (gejzír)
  • Príklad: datasets::faithful s dĺžkou trvania gejzírov (duration) a čakacou dobou na erupciu, gejzír Old Faithful v Yellowstonskom národnom parku.
Kód 
# 2D
fit <- faithful |> 
  mclust::densityMclust(plot = FALSE)
fit |> 
  plot(what = "density", type = "hdr", 
       data = faithful, points.cex = 0.5)
fit |> 
  plot(what = "density", type = "persp")

  • Zdá sa, že aj tieto gejzírové dáta vykazujú trimodalitu, jeden menší zhluk je však veľmi blízko druhého, väčšieho.

12.3.4.2 Typy kovariancií

  • Rozkladom kovariančnej matice \mathbf{\Sigma}_k = \lambda_k\mathbf{D}_k\mathbf{A}_k\mathbf{D}_k^\intercal získame kontrolu nad geometriou zhlukov, konkrétne
    • objemom elipsoidu (volume), prostredníctvom skalárnej hodnoty \lambda_k,
    • tvarom (shape), cez diagonálnu maticu \mathbf{A}_k, pričom |\mathbf{A}_k|=1,
    • orientáciou, pomocou ortogonálnej matice \mathbf{D}_k.
  • Taxonómia kóduje parametrizáciu modelu trojpísmenovou skratkou v poradí podľa jednotlivých komponentov kovariančnej matice. Jednotlivé písmená potom značia:
    • E = equal (rovnaké),
    • I = identity matrix (jednotková matica),
    • V = variable (rôzne).
  • V jednorozmernom prípade sú iba dva modely: E pre rovnaký rozptyl a V pre variabilný rozptyl.´
  • Prehľad taxonómie modelov dáva nasledujúca tabuľka a po nej ilustrácie zodpovedajúcich elíps.
Model \mathbf{\Sigma}_k Distribution Volume Shape Orientation
EII \lambda\mathbf{I} Spherical Equal Equal —
VII \lambda_k\mathbf{I} Spherical Variable Equal —
EEI \lambda\mathbf{A} Diagonal Equal Equal Coordinate axes
VEI \lambda_k\mathbf{A} Diagonal Variable Equal Coordinate axes
EVI \lambda\mathbf{A}_k Diagonal Equal Variable Coordinate axes
VVI \lambda_k\mathbf{A}_k Diagonal Variable Variable Coordinate axes
EEE \lambda\mathbf{D}\mathbf{A}\mathbf{D}^\intercal Ellipsoidal Equal Equal Equal
EVE \lambda\mathbf{D}\mathbf{A}_k\mathbf{D}^\intercal Ellipsoidal Equal Variable Equal
VEE \lambda_k\mathbf{D}\mathbf{A}\mathbf{D}^\intercal Ellipsoidal Variable Equal Equal
VVE \lambda_k\mathbf{D}\mathbf{A}_k\mathbf{D}^\intercal Ellipsoidal Variable Variable Equal
EEV \lambda\mathbf{D}_k\mathbf{A}\mathbf{D}_k^\intercal Ellipsoidal Equal Equal Variable
VEV \lambda_k\mathbf{D}_k\mathbf{A}\mathbf{D}_k^\intercal Ellipsoidal Variable Equal Variable
VVV \lambda_k\mathbf{D}_k\mathbf{A}_k\mathbf{D}_k^\intercal Ellipsoidal Variable Variable Variable

Príklad (gejzír)
  • Hustota zhlukov v datasete MASS::geyser má (v reze) eliptický tvar s približne rovnakým objemom, varianciou aj orientáciou (s osami).
  • Model umožňuje zobraziť neistotu v zaradení jednotlivých objektov do zhlukov.
Kód 
library(mclust)
fit <- MASS::geyser |> 
  mclust::Mclust(G = 3)

plot(fit, what = "density")
points(MASS::geyser, pch=19, cex=0.5)
plot(fit, what = "uncertainty")
fit |> summary()

---------------------------------------------------- 
Gaussian finite mixture model fitted by EM algorithm 
---------------------------------------------------- 

Mclust EEI (diagonal, equal volume and shape) model with 3 components: 

 log-likelihood   n df      BIC       ICL
      -1371.823 299 10 -2800.65 -2814.577

Clustering table:
  1   2   3 
 91 107 101

12.3.5 Použitá literatúra

  1. Ďalšími ordinačnými metódami sú napr. principal coordinate analysis (PCoA) resp. multidimensional scaling (MDS), uniform manifold approximation and projection for dimension reduction (UMAP), t-distributed stochastic neighbor embedding (t-SNE)…↩︎