Zhluková analýza (klaster Analysis)

---

Úvod

Klastrová (zhluková) analýza patrí medzi najpoužívanejšie metódy exploratívnej štatistiky. V praxi sa využíva všade tam, kde je potrebné rozdeliť pozorovania do homogénnych celkov - napríklad pri segmentácii zákazníkov v marketingu, identifikácii podobných krajín v makroekonomických ukazovateľoch, hodnotenízdravotných rizík, klasifikácii biologických vzoriek či v geoinformatike pri zoskupovaní priestorových sobjektov. Jej výhodou je, že pracuje s viacerými premennými naraz a dokáže odhaliť vzory, ktoré by pri samotnom hodnotení jednotlivých ukazovateľov zostali skryté. Správne zvolená metrika vzdialenosti a metóda zhlukovania umožňujú odhaliť skryté vzťahy v dátach, čím poskytujú cenný podklad pre rozhodovanie v rôznych oblastiach aplikovaného výskumu.

My predstavíme zhlukovú analýzu pri analýze údajov pochádzajúcich z meraní meteorologických premenných, konkrétne teploty vzduchu, atmosférického tlaku, rýchlosti vetra a úhrnu zrážok. Budeme využívať údaje za rok 2020. V Tab. 1. uvádzame celú nami používanú databázu.

library(knitr)
library(kableExtra)

rm(list=ls())
udaje <- read.csv2("Temperature_2020.csv",dec = ",", sep = ";",fileEncoding = "Windows-1250", stringsAsFactors = FALSE)

udaje2020 <- subset(udaje, Year == 2020)

udaje2020 <- udaje2020[, c("Year", "Temperature...C.", "Air.Pressure..hPa.", "Wind.Speed..m.s.", "Precipitation..mm.")]
rownames(udaje2020) <- seq_len(nrow(udaje2020))
udaje2020 <- subset(udaje2020, select = -Year)

Table 1.

udaje2020

Hierarchická zhluková analýza pracuje s mierami vzdialenosti medzi pozorovaniami. Aby boli tieto vzdialenosti porovnateľné, je potrebné, aby všetky premenné boli definované na rovnakejškále. Používame pritom tzv. z-škálovanie, pričom transformované \(z\) hodnoty (skóre) vypočítame nasledovne

\[z = \frac{x-\mu}{\sigma}\]

kde \(\mu\) je stredná hodnota a \(\sigma\) je štandardná odchýlka pozorovaní \(x\). Predpokladáme pritom, že súbor údajov už neobsahuje NA hodnoty, ktoré boli ošetrené v predchádzajúcich krokoch.

Touto operáciou získame škálované pozorovania, pričom ich rozloženie je znázornené nasledovne:

# =======================================================
## 1) Príprava údajov a data.frame so šlálovanými údajmi
## ======================================================

udaje_complete <- na.omit(udaje2020)
udaje_scaled <- scale(udaje_complete)

Obr. 1.

num_vars <- as.data.frame(udaje_scaled)
num_plots <- ncol(num_vars)

par(mfrow = c(ceiling(sqrt(num_plots)), ceiling(num_plots / ceiling(sqrt(num_plots)))))
par(mar = c(4, 4, 2, 1))

for (col in names(num_vars)) {
  boxplot(num_vars[[col]],
          main = col,
          col = "lightblue",
          horizontal = TRUE)
}

mtext("Boxploty numerických premenných (rok 2020)", outer = TRUE, cex = 1.3, font = 2)

Tentokrát odľahlé hodnoty nevylúčime, nakoľko predstavujú skutočné namerané pozorovania a sú súčasťou reálnych podmienok merania.

1. Teplota vzduchu (Temperature…C.): Rozdelenie teplôt je pomerne symetrické, bez extrémnych odľahlých bodov. Medián sa nachádza približne uprostred intervalu a teploty sú rozptýlené rovnomerne. Hodnoty nevykazujú výrazné anomálie.

2. Atmosférický tlak (Air.Pressure..hPa.): Atmosférický tlak má mierne asymetrické rozdelenie, no bez odľahlých pozorovaní. Väčšina hodnôt sa sústreďuje okolo mediánu a variabilita je primeraná. Údaje pôsobia stabilne a bez extrémnych výkyvov.

3. Rýchlosť vetra (Wind.Speed..m.s.): Pri rýchlosti vetra sa vyskytujú dva zreteľné odľahlé body – jeden veľmi nízky a jeden veľmi vysoký. Tieto hodnoty však predstavujú reálne namerané situácie (napr. náhle zosilnenie vetra alebo bezvetrie), preto ich ponechávame v analýze. Ostatné hodnoty sú sústredené v úzkom intervale okolo mediánu, čo naznačuje nízku variabilitu väčšiny meraní.

4. Úhrn zrážok (Precipitation..mm.): Zrážky vykazujú širšie rozptyl než ostatné premenné, čo je typické pre túto veličinu, keďže zrážky môžu byť veľmi variabilné. Rozdelenie je bez odľahlých bodov a medián sa nachádza v spodnej časti boxu, čo naznačuje, že väčšina meraní mala relatívne nízke úhrny zrážok.

Pri zhlukovej analýze je dôležitá korelačná matica premenných. Vysoká korelácia zvýhodňuje pri zhlukovej analýze korelované premenné. Preto pri korelácii nad 0,8 alebo 0.9 vylúčime jednu z korelovaných premenných. V Tab. 2. sa však takáto vysoká korelácia nenachádza, preto sa nemusíme ďalej s problémom zaoberať. > V prípade, ak máme väčší počet významne korelovaných premenných, sa odporúča i transformácia pomocou Analýzy hlavných komponentov (Principal Component Analysis)

Tab. 2

cor_mat <- cor(udaje_scaled, use="pairwise.complete.obs")
cor_mat <- round(cor_mat,2)
print(cor_mat)

                   Temperature...C. Air.Pressure..hPa. Wind.Speed..m.s.
Temperature...C.               1.00              -0.17            -0.23
Air.Pressure..hPa.            -0.17               1.00             0.07
Wind.Speed..m.s.              -0.23               0.07             1.00
Precipitation..mm.             0.54              -0.48            -0.42
                   Precipitation..mm.
Temperature...C.                 0.54
Air.Pressure..hPa.              -0.48
Wind.Speed..m.s.                -0.42
Precipitation..mm.               1.00

Každému pozorovaniu zodpovedá jeden riadok meraných hodnôt. Vzdialenosť medzi pozorovaniami \(i\) a \(j\) je:

\[ d^{ij} = \sqrt{\sum_k (x^i_k - x^j_k)^2} \] kde \(x^i_k\) je hodnota \(k\)tej premennej (Temperature…C., Air.Pressure..hPa., Wind.Speed..m.s., Precipitation..mm.) pre pozorovanie \(i\). Tento typ vzdialenosti nazývame aj Euklidovská vzdialenosť. Vzdialenosti medzi jednotlivými pozorovaniami sa súhrnne vyjadrujú aj v matici vzdialeností, čo je v našom prípade uvedené v Tab.3.. Analýzou tejto tabuľky zistíme, že najväčšia vzdialenosť je medzi mesiacmi 8 a 11, čo znamená, že tieto mesiace sa podľa meteorologických ukazovateľov najviac líšia – majú odlišnú teplotu, tlak, zrážky aj rýchlosť vetra. Naopak, najmenšia vzdialenosť je medzi mesiacmi 4 a 5, ktoré sú si veľmi podobné a vykazujú takmer rovnaké klimatické podmienky. Tieto rozdiely môžeme vysvetliť sezónnosťou – letné mesiace sa prirodzene výrazne odlišujú od jesenných či zimných mesiacov, zatiaľ čo susediace jarné mesiace sú si klimaticky veľmi blízke.

Tab. 3

## ============================
## 3) Distance matrix
## ============================
rownames(udaje_scaled) <- rownames(udaje_complete)
dist_mat <- round(dist(udaje_scaled, method = "euclidean"), 2)
dist_mat

      1    2    3    4    5    6    7    8    9   10   11
2  3.92                                                  
3  1.93 4.42                                             
4  1.61 3.34 1.96                                        
5  1.81 2.91 2.44 1.06                                   
6  3.83 2.89 3.36 2.76 2.27                              
7  2.96 3.54 2.97 1.47 1.52 2.16                         
8  4.11 3.09 3.81 2.98 2.40 0.66 2.10                    
9  3.49 2.25 3.67 2.69 1.84 1.26 2.31 1.23               
10 3.17 3.33 2.09 2.34 2.32 1.57 2.48 2.17 2.36          
11 0.77 3.90 2.50 1.54 1.66 3.89 2.71 4.06 3.42 3.51     
12 3.20 3.46 2.43 2.42 3.08 3.32 3.15 3.83 3.77 2.26 3.54

Princíp hierarchického zhlukovania (Wardova metóda)

Zhlukovanie v prípade Wardovej metódy prebieha zdola smerom nahor, t.j. začíname s jednočlennými klastrami, ktoré postupne zlučujeme. Táto metóda patrí teda medzi aglomeratívne hierarchické metódy. Minimalizuje nárast vnútornej variability pri spojení dvoch klastrov, pričom využíva nasledovné výpočty:

Wardová metóda minimalizuje sumu štvorcov chýb (Error sum of Squares - ESS)

\[ESS(C) = \sum_{i \in C} \lVert x_i - \bar{x}_C \rVert^2\] kde \(C\) je zvažovaný klaster (zhluk). V každom kroku zlučovania dvoch klasterov, Wardova metóda hľadá minimálny prírastok sumy štvorcov chýb (\(\Delta ESS\)), pričom

\[\Delta ESS = ESS(A \cup B) - ESS(A) - ESS(B)\] Dvojica zhlukov, ktoré tejto podmienke o minimalizácii vyhovuje, je následne zlúčená a prechádza sa k ďalšiemu kkroku. To spravidla vedie k vytváraniu homogénnych zhlukov, pričom nedochádza k odtrhávaniu odľahlých hodnôt tak, ako pri iných zhlukovacích metódach.

Obr. 2. Hierarchické zhlukovanie - dendogram. Červená čiara určuje rez definujúci tri klastre.

## ============================
## 4) Hierarchical klastering
## ============================

hc <- hclust(dist_mat, method = "ward.D2")

plot(hc, labels = rownames(udaje_scaled),
     main = "Hierarchical klastering of observations (Ward.D2)",
     xlab = "", sub = "")

k <- 3
h_cut <- hc$height[length(hc$height) - (k - 1)]
abline(h = h_cut, col = "red", lwd = 2, lty = 2)

klaster_membership <- cutree(hc, k = k)

udaje_klasters <- data.frame(
  Observation = rownames(udaje_complete),
  udaje_complete,
  klaster = factor(klaster_membership)
)

Tab.4. Príslušnosť pozorovaní do klastrov.

data_prac <- data.frame( Observation = rownames(udaje_klasters),
  klaster     = udaje_klasters$klaster)
colnames(data_prac) <- c("Observation","klaster")
data_prac

Na základe vykonanej zhlukovej analýzy boli mesiace rozdelené do troch klastrov. Z výslednej tabuľky klastrov vidíme toto rozdelenie: klaster 1: mesiace 1, 4, 5, 7, 11 klaster 2: mesiace 2, 6, 8, 9 klaster 3: mesiace 3, 10, 12. Tieto klastre reprezentujú mesiace, ktoré majú podobné hodnoty teploty, tlaku, rýchlosti vetra a zrážok. Klaster 1 združuje mesiace, ktoré majú podobné a mierne stabilné meteorologické podmienky — ide najmä o mesiace s podobnými teplotami a rýchlosťou vetra. Klaster 2 zahŕňa mesiace s vyššou podobnosťou najmä v úrovni zrážok a tlaku vzduchu, pričom tieto mesiace sa od klastrov 1 a 3 odlišujú práve kombináciou týchto ukazovateľov. Klaster 3 obsahuje mesiace, ktoré majú odlišný profil — či už z hľadiska teploty alebo zrážok — a preto sa od zvyšných mesiacov jasne oddeľujú. Výsledkom je teda logické rozdelenie mesiacov do troch skupín, kde každý klaster spája mesiace s podobnými meteorologickými charakteristikami.

Deskriptívne štatistiky výsledkov

Na základe Tab. 5 možno vyhodnotiť separačnú silu jednotlivých premenných vo vytvorených klastroch nasledovne: Premenná Precipitation..mm. má najvyšší podiel medzi-klastrového rozptylu (Prop_Between ≈ 0.80), čo znamená, že veľmi dobre odlišuje jednotlivé klastre – zrážky sú teda jedným z najsilnejších faktorov, podľa ktorých sa mesiace prirodzene zoskupujú. Podobne aj Wind.Speed..m.s. vykazuje vysokú hodnotu Prop_Between (≈ 0.65), čo naznačuje, že rýchlosť vetra sa medzi klastrami významne líši a predstavuje dobrý separátor. Premenná Air.Pressure..hPa. má strednú separačnú schopnosť (≈ 0.58), čo znamená, že tlak vzduchu prispieva k odlišovaniu klastrov, ale už menej výrazne ako predchádzajúce dve premenné. Najslabším separátorom je Temperature…C., ktorej podiel medzi-klastrového rozptylu je najnižší (≈ 0.29). Znamená to, že teplota sa síce medzi mesiacmi líši, ale nie dostatočne na to, aby zásadne prispievala k vytváraniu odlišných klastrov. Najviac k rozlíšeniu klastrov prispievajú zrážky a rýchlosť vetra, zatiaľ čo teplota má v tomto type analýzy iba slabý separačný efekt.

Tab. 5. Vysvetlenie vnútroklastrovej variability z hľadiska jednotlivých premenných

## ============================
## 5) Variability measures
## ============================

ssq <- function(x, m) sum((x - m)^2)

var_names <- colnames(udaje_scaled)

TSS <- sapply(var_names, function(v) ssq(udaje_scaled[, v], mean(udaje_scaled[, v])))

WSS <- sapply(var_names, function(v) {
  x <- udaje_scaled[, v]
  tapply(x, klaster_membership, function(z) ssq(z, mean(z))) |> sum()
})

BSS <- TSS - WSS

ss_table <- data.frame(
  Variable = var_names,
  TSS = TSS,
  WSS = WSS,
  BSS = BSS,
  Prop_Between = BSS / TSS
)

ss_table

attach(udaje2020)
udaje2020 <- data.frame(cbind(udaje2020,udaje_klasters$klaster))
colnames(udaje2020) <- c("Temperature_C","AirPressure_hPa","WindSpeed_ms","Precipitation_mm","klaster")

Tab. 6. Centroidy - priemerné hodnoty sledovaných premenných

library(dplyr)

descriptives <- udaje_klasters %>%
  group_by(klaster) %>%
  summarise(
    across(
      .cols = where(is.numeric),
      .fns = list(
        mean = ~mean(.x, na.rm = TRUE)
      ),
      .names = "{.col}_{.fn}"
    )
  )
descriptives

NA

Prvý klaster má mierne nadpriemernú teplotu, najvyšší priemerný atmosférický tlak (≈ 993.8 hPa) a strednú úroveň rýchlosti vetra. Druhý klaster je charakterizovaný najvyššou priemernou teplotou (≈ 14.8 °C), nízkym tlakom a najnižšou rýchlosťou vetra spomedzi všetkých klastrov. Tretí klaster má najnižšie priemerné teploty, najnižší atmosférický tlak a zároveň najvyššiu priemernú rýchlosť vetra, čo naznačuje odlišné poveternostné podmienky oproti zvyšným dvom klastrom. Z toho vyplýva, že najdôležitejším rozlišovacím faktorom medzi klastrami je teplota, zatiaľ čo tlak vzduchu a rýchlosť vetra dopĺňajú charakteristiku jednotlivých skupín.

Záver

Predložená analýza sa zaoberá zoskupovaním jednotlivých mesiacov roka na základe ich priemerných meteorologických charakteristík – konkrétne teploty vzduchu, atmosférického tlaku, rýchlosti vetra a množstva zrážok. Na základe týchto premenných boli mesiace rozdelené do troch klastrov, pričom každý klaster predstavuje odlišný typ poveternostných podmienok. Takto získaná klasifikácia umožňuje lepšie pochopiť sezónne rozdiely v počasí, identifikovať mesiace s podobnými klimatickými profilmi a môže slúžiť ako podklad pri plánovaní aktivít závislých od počasia – napríklad pri energetickom plánovaní, poľnohospodárstve či hodnotení klimatických trendov.