---
title: "Ćwiczenia 3 — Uczenie nienadzorowane: PCA, mapy próbek, skupienia (wersja R, z rozwiązaniami)"
subtitle: "Dane demo: kohorta T2D (QinJ_2012)"
output:
  html_document:
    toc: true
    toc_float: true
    code_folding: show
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE, message = FALSE, warning = FALSE)
```

Na Ćwiczeniach 2 oglądaliśmy cechy pojedynczo. Teraz patrzymy na wszystkie cechy
jednocześnie i pytamy, czy próbki układają się w jakąś strukturę: czy tworzą grupy, jakie
kierunki zmienności dominują w danych i czy ta struktura ma związek z klasą. Metody
nienadzorowane (ang. *unsupervised learning*) nie korzystają z etykiet: klasę nanosimy na
wyniki dopiero na końcu, żeby sprawdzić, czy największa zmienność w danych wiąże się
z chorobą. Jeśli tak, model nadzorowany ma łatwe zadanie; jeśli nie, różnice między klasami
są subtelne i model musi ich szukać wśród wielu innych źródeł zmienności.

Pracujemy na całym zbiorze, bez podziału na trening i test. Przy oglądaniu danych bez
etykiet wolno tak robić, bo niczego nie dopasowujemy do klasy i niczego nie oceniamy.
Decyzje dotyczące modelu (np. wybór przekształcenia albo liczby składowych jako cech)
podejmuje się później, wyłącznie na zbiorze treningowym.

# Wczytanie danych

Wczytujemy te same pliki co na Ćwiczeniach 2 i zostawiamy cechy obecne w co najmniej 10%
próbek. Arkusz trzeba zapisać w tym samym katalogu co arkusze z Ćwiczeń 1 i 2, czyli obok
katalogu `dane/`. Arkusz najlepiej uruchamiać w nowej sesji R (w RStudio: *Session → Restart R*),
blok po bloku od początku.

Pakiety metod nienadzorowanych instaluje się jednorazowo (jeśli `library()` zgłasza błąd
`there is no package called ...`):

```{r instalacja, eval=FALSE}
install.packages(c("vegan", "Rtsne", "uwot", "cluster", "mclust"))   # jednorazowo, potem tylko library()
```

```{r dane}
library(ggplot2)                        # wykresy
library(vegan)                          # vegdist(): odległość Braya–Curtisa
library(Rtsne)                          # t-SNE
library(uwot)                           # UMAP
library(cluster)                        # silhouette(): szerokość sylwetki
library(mclust)                         # adjustedRandIndex(): zgodność dwóch podziałów

# row.names = 1: pierwsza kolumna to nazwy próbek; check.names = FALSE: nazwy taksonów bez zmian
X <- read.csv("dane/t2d_features.csv", row.names = 1, check.names = FALSE)
# factor() z levels: klasa jako zmienna kategoryczna, najpierw control, potem T2D
y <- factor(read.csv("dane/t2d_labels.csv")$grupa, levels = c("control", "T2D"))

Xf <- X[, colMeans(X > 0) >= 0.10]      # cechy obecne w co najmniej 10% próbek
dim(Xf)                                 # liczba próbek i cech
table(y)                                # liczba próbek w każdej klasie

cols <- c(control = "#2C5F7C", T2D = "#D55E00")  # kolory klas na wszystkich rycinach (granatowy, pomarańczowy)
```

# 1. Przekształcenie CLR

Na Ćwiczeniach 2 okazało się, że na surowych danych odległości między próbkami wyznacza
w dużej mierze jedna cecha (*Prevotella copri*), bo ma największy rozrzut. Metody z tego
arkusza opierają się na odległościach i wariancji, więc bez przekształcenia pokazywałyby
głównie tę jedną cechę. Do tego dane są kompozycyjne (ang. *compositional data*): wartości
w każdej próbce sumują się do 100%, więc wzrost jednej cechy oznacza spadek pozostałych.
Przekształcenie CLR (ang. *centered log-ratio*) rozwiązuje oba problemy: każdą wartość
logarytmuje się i odejmuje od niej średnią logarytmów w tej samej próbce. Wynik mówi, ile
razy cecha jest liczniejsza od „typowej” cechy w danej próbce, więc liczą się proporcje,
a nie wartości bezwzględne. Logarytm z zera nie istnieje, dlatego przed logarytmowaniem
dodajemy małą stałą `1e-6`, czyli pseudoliczność (ang. *pseudocount*), tę samą co
w całym kursie.

```{r clr}
# własna funkcja: CLR dla każdego wiersza (próbki) macierzy M
clr <- function(M) {
  L <- log(as.matrix(M) + 1e-6)         # logarytm naturalny po dodaniu pseudoliczności
  L - rowMeans(L)                       # od każdej wartości odejmujemy średnią logarytmów jej wiersza
}
Z <- clr(Xf)                            # dane po przekształceniu CLR: te same próbki i cechy
dim(Z)                                  # wymiary się nie zmieniają

variances <- sort(apply(Z, 2, var), decreasing = TRUE)   # wariancja każdej cechy po CLR
round(sum(variances[1:5]) / sum(variances), 3)           # udział 5 cech o największej wariancji
d_top1 <- dist(Xf[, "species:Prevotella copri"])         # odległości liczone tylko na P. copri
cor(dist(Z), d_top1)                                     # na ile ta cecha wyznacza odległości po CLR
```

**Na co patrzeć.** Na surowych danych pięć cech skupiało około 60% całkowitej wariancji,
a po CLR skupia tylko kilka procent. Korelacja odległości z odległością liczoną na samej
*P. copri* spadła z około 0,70 (Ćwiczenia 2) do wartości bliskiej zeru: żadna pojedyncza
cecha nie wyznacza już podobieństwa próbek.

# 2. Analiza głównych składowych

Analiza głównych składowych (ang. *principal component analysis*, PCA) szuka kierunków,
wzdłuż których próbki różnią się najbardziej. Pierwsza składowa główna (PC1) to kierunek
największej zmienności, druga (PC2) kolejny, prostopadły do pierwszego, i tak dalej. Po co:
dwie pierwsze składowe pozwalają pokazać próbki opisane 190 cechami na jednej płaskiej rycinie,
a odsetek wariancji mówi, jaką część informacji ta rycina zachowuje. PCA liczymy na danych
CLR po wycentrowaniu, bez standaryzacji, jak na wykładzie 3: po CLR wszystkie cechy mają
wspólną jednostkę (logarytm ilorazu).

```{r pca}
pca <- prcomp(Z, center = TRUE, scale. = FALSE)   # center: odjęcie średniej każdej cechy; bez standaryzacji
var_prop <- pca$sdev^2 / sum(pca$sdev^2)          # sdev^2: wariancja każdej składowej; dzielimy przez sumę
round(var_prop[1:5], 3)                           # odsetek wariancji pięciu pierwszych składowych
round(cumsum(var_prop)[c(2, 10)], 3)              # cumsum(): odsetek skumulowany dla 2 i 10 składowych
```

Wykres osypiska (ang. *scree plot*) pokazuje odsetek wariancji kolejnych składowych:

```{r osypisko, fig.width=6, fig.height=3.5}
scree <- data.frame(pc = 1:20, value = 100 * var_prop[1:20])   # 20 pierwszych składowych, w procentach
ggplot(scree, aes(x = pc, y = value)) +   # wykres: składowa na osi X, odsetek wariancji na osi Y
  geom_col(fill = "grey60") +                     # słupek dla każdej składowej
  labs(x = "Składowa główna", y = "Wariancja wyjaśniona (%)") +   # podpisy osi
  theme_bw()   # jasne tło z siatką
```

**Na co patrzeć.** PC1 i PC2 wyraźnie odstają od pozostałych, a dalej słupki maleją powoli.
Wariancja rozkłada się na wiele kierunków, więc rycina z dwiema osiami pokazuje tylko
niewielką część informacji zawartej w danych.

Wykres próbek na dwóch pierwszych składowych. PCA nie zna etykiet; kolor klasy nanosimy
po obliczeniach:

```{r pca-rozrzut, fig.width=6, fig.height=4.5}
pc_df <- data.frame(PC1 = pca$x[, 1], PC2 = pca$x[, 2], group = y)   # pca$x: położenie próbek na składowych
ggplot(pc_df, aes(x = PC1, y = PC2, colour = group)) +   # osie: PC1 i PC2; kolor: klasa
  geom_point(shape = 1, stroke = 0.8) +           # jedna próbka = jedno kółko
  scale_colour_manual(values = cols, name = NULL) +   # kolory klas z wektora cols, legenda bez tytułu
  labs(x = paste0("PC1 (", round(100 * var_prop[1], 1), "%)"),     # odsetek wariancji w podpisie osi
       y = paste0("PC2 (", round(100 * var_prop[2], 1), "%)")) +   # paste0() skleja napis bez spacji
  theme_bw()   # jasne tło z siatką
```

Ładunki (ang. *loadings*) mówią, z jaką wagą każda cecha wchodzi do składowej. Pięć cech
o największym ładunku PC1 (bez względu na znak, bo znak składowej jest umowny i może się
różnić między programami):

```{r ladunki}
load_pc1 <- pca$rotation[, 1]                                # ładunki PC1: jedna liczba na cechę
top5 <- names(sort(abs(load_pc1), decreasing = TRUE))[1:5]   # pięć największych wartości bezwzględnych
data.frame(loading = round(load_pc1[top5], 3), row.names = sub("^species:", "", top5))   # tabela ładunków z krótkimi nazwami cech
```

**Na co patrzeć.** Cechy o przeciwnych znakach działają w przeciwnych kierunkach: próbki
z jednego końca PC1 mają względnie więcej cech o ładunku dodatnim, a z drugiego końca więcej
cech o ładunku ujemnym. Ładunek opisuje udział cechy w zmienności danych, a nie jej związek
z chorobą.

🔧 **Do sprawdzenia.** W bloku `pca` zmienić `scale. = FALSE` na `scale. = TRUE`
(standaryzacja: każda cecha z tą samą wagą) i uruchomić sekcję ponownie. Jak zmienia się
odsetek wariancji PC1 i wygląd wykresu próbek?

> **✏️ Zadanie 1.** Odczytać, jaki odsetek wariancji wyjaśniają razem PC1 i PC2 i ile
> wyjaśnia 10 pierwszych składowych. Czy na wykresie PC1–PC2 klasy tworzą oddzielne chmury,
> czy się nakładają? Wzdłuż której osi różnią się bardziej? Odpowiedzieć w 2–3 zdaniach.

```{r zad-1}
round(var_prop[1] + var_prop[2], 3)          # odsetek wariancji PC1 i PC2 razem
round(cumsum(var_prop)[10], 3)               # odsetek wariancji 10 pierwszych składowych
round(tapply(pca$x[, 1], y, median), 1)      # mediana położenia na PC1 w każdej klasie
round(tapply(pca$x[, 2], y, median), 1)      # mediana położenia na PC2 w każdej klasie
```

**Rozwiązanie.** PC1 i PC2 wyjaśniają razem 18,2% wariancji, a 10 pierwszych składowych 36,3%.
Na wykresie PC1–PC2 klasy się nakładają: nie ma oddzielnych chmur. Klasy różnią się nieco
bardziej wzdłuż PC2 (mediany: control −1,9, T2D 4,8) niż wzdłuż PC1 (−7,3 i −1,1), ale
przesunięcie jest małe w porównaniu z rozrzutem próbek, który na każdej osi sięga kilkudziesięciu
jednostek. Największa zmienność w danych nie wiąże się więc wyraźnie z chorobą. W Pythonie PC2
ma przeciwny znak niż w R, więc mediany PC2 mają odwrócone znaki; znak składowej jest umowny.

> **✏️ Zadanie 2.** Sprawdzić, co opisuje PC1: policzyć korelację położenia próbek na PC1
> z liczbą cech niezerowych w próbce, `cor(pca$x[, 1], rowSums(Xf > 0))`. Co oznacza ta
> wartość? Wskazówka do interpretacji: po CLR z pseudolicznością 1e-6 każde zero staje się
> bardzo małą liczbą (logarytm z 1e-6 to około −14).

```{r zad-2}
cor(pca$x[, 1], rowSums(Xf > 0))           # PC1 a liczba cech niezerowych w próbce
cor(pca$x[, 1], prcomp(Xf > 0)$x[, 1])     # PC1 z danych CLR a PC1 z samej macierzy obecności (0/1)
```

**Rozwiązanie.** Korelacja PC1 z liczbą cech niezerowych wynosi −0,64: próbki z małą liczbą
obecnych cech leżą po jednej stronie osi, a próbki z wieloma cechami po drugiej. Druga linijka
pokazuje to jeszcze wyraźniej: PC1 z danych CLR prawie pokrywa się z PC1 policzoną na samej
informacji „cecha jest / cechy nie ma” (korelacja 0,95 co do wartości bezwzględnej; znak
składowej jest umowny). Po CLR z pseudolicznością 1e-6 zero zamienia się w bardzo małą liczbę
(około −14), więc różnica między obecnością a brakiem cechy jest największą różnicą w danych
i to ją opisuje PC1, a nie klasę. Wynik zależy od pseudoliczności: przy 0,5 ten związek prawie
znika (wykład 3).

# 3. Analiza współrzędnych głównych na odległości Braya–Curtisa

PCA opiera się na odległości euklidesowej między próbkami. W danych o składzie próbek często
stosuje się inną miarę: odległość Braya–Curtisa (ang. *Bray–Curtis dissimilarity*), liczoną
na udziałach procentowych; 0 oznacza identyczny skład, 1 brak wspólnych cech. Analiza
współrzędnych głównych (ang. *principal coordinates analysis*, PCoA) rysuje próbki tak, by
odległości na rycinie jak najlepiej odtwarzały wybraną miarę. Po co: inna miara odległości
może pokazać inną strukturę danych, więc warto sprawdzić, czy wnioski z PCA się powtarzają.

```{r pcoa, fig.width=6, fig.height=4.5}
d_bc <- vegdist(Xf, method = "bray")              # odległość Braya–Curtisa dla każdej pary próbek (na udziałach)
pcoa <- cmdscale(d_bc, k = 2, eig = TRUE)         # k = 2: dwie osie do ryciny; eig = TRUE: zwraca wartości własne
axis_prop <- pcoa$eig / sum(pcoa$eig[pcoa$eig > 0])   # odsetek osi (względem sumy dodatnich wartości własnych)
round(axis_prop[1:2], 3)                          # odsetek dwóch pierwszych osi

pcoa_df <- data.frame(axis1 = pcoa$points[, 1], axis2 = pcoa$points[, 2], group = y)   # położenie próbek
ggplot(pcoa_df, aes(x = axis1, y = axis2, colour = group)) +   # osie: dwie pierwsze osie PCoA; kolor: klasa
  geom_point(shape = 1, stroke = 0.8) +   # punkty jako puste kółka; stroke: grubość obwódki
  scale_colour_manual(values = cols, name = NULL) +   # kolory klas z wektora cols, legenda bez tytułu
  labs(x = paste0("Oś 1 (", round(100 * axis_prop[1], 1), "%)"),   # podpis osi z odsetkiem osi
       y = paste0("Oś 2 (", round(100 * axis_prop[2], 1), "%)")) +   # to samo dla osi 2
  theme_bw()   # jasne tło z siatką

# czy pierwsza oś PCoA wiąże się z udziałem najliczniejszej cechy (korelacja rang)
cor(pcoa$points[, 1], Xf[, "species:Prevotella copri"], method = "spearman")
```

**Na co patrzeć.** Odległość Braya–Curtisa liczona na udziałach procentowych zależy przede
wszystkim od cech liczebnych, dlatego pierwsza oś PCoA wyraźnie wiąże się z udziałem
*P. copri*. To inny obraz niż w PCA na danych CLR: wybór miary odległości zmienia to, co
widać na mapie próbek.

# 4. Mapy t-SNE i UMAP

t-SNE (ang. *t-distributed stochastic neighbor embedding*) i UMAP (ang. *uniform manifold
approximation and projection*) układają próbki na płaszczyźnie tak, żeby próbki bliskie
w danych leżały blisko na mapie. Po co: lepiej niż PCA pokazują lokalne grupy podobnych
próbek. Interpretuje się tylko sąsiedztwo; odległości między odległymi grupami i wielkości
grup na mapie nie odpowiadają danym. Obie metody są losowe, więc ustalamy ziarno. Wejściem
jest 10 pierwszych składowych PCA, co pomija składowe o małej wariancji i skraca obliczenia.
Parametr `perplexity` w t-SNE to w przybliżeniu liczba sąsiadów, których bierze pod uwagę
każda próbka; `n_neighbors` w UMAP pełni podobną rolę.

```{r mapy, fig.width=10, fig.height=4.5}
scores10 <- pca$x[, 1:10]                         # 10 pierwszych składowych głównych
set.seed(42)                                      # ziarno: ta sama mapa przy każdym uruchomieniu
tsne <- Rtsne(scores10, dims = 2, perplexity = 30, pca = FALSE)   # pca = FALSE: wejście to już składowe PCA
set.seed(42)   # ziarno: ten sam wynik przy każdym uruchomieniu
umap_xy <- umap(scores10, n_neighbors = 15, min_dist = 0.1, n_threads = 1)   # min_dist jak w wersji Python

# tabela „długa”: obie mapy jedna pod drugą, z nazwą metody w osobnej kolumnie
maps <- rbind(
  data.frame(d1 = tsne$Y[, 1], d2 = tsne$Y[, 2], group = y, method = "t-SNE (perplexity = 30)"),
  data.frame(d1 = umap_xy[, 1], d2 = umap_xy[, 2], group = y, method = "UMAP (n_neighbors = 15)")
)
ggplot(maps, aes(x = d1, y = d2, colour = group)) +   # osie: dwa wymiary mapy; kolor: klasa
  geom_point(shape = 1, stroke = 0.8) +   # punkty jako puste kółka; stroke: grubość obwódki
  scale_colour_manual(values = cols, name = NULL) +   # kolory klas z wektora cols, legenda bez tytułu
  facet_wrap(~ method, scales = "free") +         # osobny panel dla każdej metody
  labs(x = "Wymiar 1", y = "Wymiar 2") +   # podpisy osi (bez jednostek)
  theme_bw()   # jasne tło z siatką
```

**Na co patrzeć.** Kolory pokazują, czy najbliżsi sąsiedzi próbki należą zwykle do tej
samej klasy. Jeśli w małych fragmentach mapy kolory są wymieszane, sąsiedztwo w danych słabo
wiąże się z chorobą. Skupiska widoczne na mapie należy traktować jako hipotezę do
sprawdzenia, bo t-SNE i UMAP tworzą wyraźne skupiska także w danych bez wyraźnej struktury.

🔧 **Do sprawdzenia.** W wywołaniu `umap()` zmienić `n_neighbors` na 5, a potem na 50.

> **✏️ Zadanie 3.** W wywołaniu `Rtsne()` zmienić `perplexity` na 5, a potem na 50
> (ziarno bez zmian) i za każdym razem narysować mapę ponownie. Jak zmienia się liczba
> i rozstaw skupisk? Czy przy którejś wartości klasy rozdzielają się wyraźniej? Co z tego
> wynika dla interpretacji jednej mapy t-SNE?

```{r zad-3, fig.width=12, fig.height=4}
maps_p <- NULL                                    # pusta tabela na wyniki trzech map
for (p in c(5, 30, 50)) {                         # pętla po wartościach perplexity
  set.seed(42)                                    # to samo ziarno dla każdej wartości
  ts <- Rtsne(scores10, dims = 2, perplexity = p, pca = FALSE)   # mapa t-SNE dla danej perplexity
  maps_p <- rbind(maps_p, data.frame(d1 = ts$Y[, 1], d2 = ts$Y[, 2], group = y,
                                     perplexity = paste("perplexity =", p)))   # dopisanie mapy do tabeli
}
maps_p$perplexity <- factor(maps_p$perplexity, levels = paste("perplexity =", c(5, 30, 50)))   # kolejność paneli
ggplot(maps_p, aes(x = d1, y = d2, colour = group)) +   # osie: dwa wymiary mapy; kolor: klasa
  geom_point(shape = 1, stroke = 0.8) +                 # punkty jako puste kółka
  scale_colour_manual(values = cols, name = NULL) +     # kolory klas z wektora cols
  facet_wrap(~ perplexity, scales = "free") +           # osobny panel dla każdej wartości perplexity
  labs(x = "Wymiar 1", y = "Wymiar 2") +                # podpisy osi
  theme_bw()                                            # jasne tło z siatką
```

**Rozwiązanie.** Przy `perplexity = 5` mapa rozpada się na wiele małych, zwartych grupek
złożonych z kilku próbek; przy 30 i 50 próbki tworzą jedną, mniej więcej równomierną chmurę.
W żadnym ustawieniu klasy nie tworzą oddzielnych skupisk, a kolory w małych fragmentach mapy są
wymieszane. Liczba i rozstaw skupisk na mapie t-SNE zależą więc od parametru, a nie tylko od
danych: małe skupiska przy `perplexity = 5` nie oznaczają, że w danych istnieją takie grupy.
Jednej mapy t-SNE nie należy interpretować bez sprawdzenia kilku wartości parametru.

# 5. Analiza skupień

Analiza skupień (ang. *clustering*) dzieli próbki na grupy podobnych próbek. Po co:
sprawdzamy, czy dane same dzielą się na grupy i czy te grupy odpowiadają klasom. Metoda
k-średnich (ang. *k-means*) dzieli próbki na k skupień tak, by próbki leżały jak najbliżej
środka swojego skupienia. Liczbę skupień wybieramy dwiema miarami: metodą łokcia (ang.
*elbow method*; suma kwadratów odległości od środków skupień, WSS, maleje z każdym k i szuka
się miejsca, po którym spada już wolno) i średnią szerokością sylwetki (ang. *silhouette
width*; od −1 do 1, im wyżej, tym lepiej oddzielone skupienia). Skupienia wyznaczamy na tych
samych 10 składowych PCA.

```{r kmeans, fig.width=9, fig.height=3.5}
ks <- 2:6                                         # sprawdzane liczby skupień
d_scores <- dist(scores10)                        # odległości między próbkami, potrzebne do sylwetki
wss <- numeric(length(ks))                        # puste wektory na wyniki
sil <- numeric(length(ks))   # pusty wektor na średnie sylwetki
for (i in seq_along(ks)) {                        # pętla po kolejnych wartościach k
  set.seed(42)                                    # k-średnich startuje z losowych środków
  km_k <- kmeans(scores10, centers = ks[i], nstart = 25)   # nstart = 25: 25 startów, zostaje najlepszy podział
  wss[i] <- km_k$tot.withinss                     # WSS dla tego k
  sil[i] <- mean(silhouette(km_k$cluster, d_scores)[, "sil_width"])   # średnia sylwetka wszystkich próbek
}
data.frame(k = ks, WSS = round(wss), silhouette = round(sil, 3))   # tabela wyników dla każdego k

# tabela „długa”: obie miary jedna pod drugą
k_df <- rbind(data.frame(k = ks, value = wss, measure = "WSS (metoda łokcia)"),
              data.frame(k = ks, value = sil, measure = "Średnia sylwetka"))
k_df$measure <- factor(k_df$measure, levels = c("WSS (metoda łokcia)", "Średnia sylwetka"))   # kolejność paneli
ggplot(k_df, aes(x = k, y = value)) +   # wykres: k na osi X, wartość miary na osi Y
  geom_line() + geom_point() +   # linia łącząca punkty i same punkty
  facet_wrap(~ measure, scales = "free_y") +      # osobna skala osi Y dla każdej miary
  labs(x = "Liczba skupień k", y = NULL) +   # podpis osi X, bez podpisu osi Y
  theme_bw()   # jasne tło z siatką
```

**Na co patrzeć.** Krzywa WSS opada bez wyraźnego załamania, a średnia sylwetka jest
największa dla k = 2, ale niska. Próbki nie tworzą wyraźnie oddzielonych grup; dalej
porównujemy podział na dwa skupienia z dwiema klasami.

Grupowanie hierarchiczne (ang. *hierarchical clustering*) łączy kolejno najbliższe próbki
i skupienia w drzewo, czyli dendrogram. Metoda Warda (`ward.D2`) przy każdym łączeniu wybiera
parę skupień, której połączenie najmniej zwiększa WSS. Przecięcie drzewa daje podział na
skupienia. Dendrogram rysujemy grafiką podstawową R, bo ggplot2 nie ma gotowej funkcji do
rysowania drzew.

```{r hclust, fig.width=9, fig.height=4}
hc <- hclust(d_scores, method = "ward.D2")        # drzewo na tych samych odległościach co sylwetka
plot(hc, labels = FALSE, hang = -1, xlab = "", sub = "", ylab = "Wysokość",   # labels = FALSE: bez nazw próbek
     main = "Dendrogram (metoda Warda)")   # tytuł ryciny
rect.hclust(hc, k = 2, border = "#B2182B")        # ramki wokół dwóch skupień
```

Porównanie skupień z klasami. Skorygowany indeks Randa (ang. *adjusted Rand index*, ARI)
mierzy zgodność dwóch podziałów tych samych próbek: 1 oznacza pełną zgodność, 0 zgodność
taką jak przy podziale losowym.

```{r skupienia-klasy}
set.seed(42)   # ziarno: ten sam wynik przy każdym uruchomieniu
km <- kmeans(scores10, centers = 2, nstart = 25)  # k-średnich dla k = 2
hc_cl <- cutree(hc, k = 2)                        # cutree(): przecięcie drzewa na 2 skupienia
table(cluster = km$cluster, class = y)            # tabela krzyżowa: skupienie k-średnich × klasa
table(cluster = hc_cl, class = y)                 # skupienie metody Warda × klasa
round(adjustedRandIndex(km$cluster, y), 3)        # ARI: k-średnich a klasa
round(adjustedRandIndex(hc_cl, y), 3)             # ARI: metoda Warda a klasa
round(adjustedRandIndex(km$cluster, hc_cl), 3)    # ARI: obie metody skupień między sobą
```

**Na co patrzeć.** Numery skupień są umowne i mogą się zamienić po zmianie ziarna albo
między programami; ARI od numeracji nie zależy. Gdyby skupienia odpowiadały klasom, w każdym
wierszu tabeli przeważałaby wyraźnie jedna klasa.

> **✏️ Zadanie 4.** Na podstawie tabel i wartości ARI: czy skupienia odpowiadają klasom?
> Czy obie metody skupień znajdują podobny podział? Policzyć medianę liczby cech niezerowych
> w każdym skupieniu k-średnich, `tapply(rowSums(Xf > 0), km$cluster, median)`, i ocenić,
> czym różnią się skupienia. Odpowiedzieć w 2–3 zdaniach.

```{r zad-4}
tapply(rowSums(Xf > 0), km$cluster, median)   # mediana liczby cech niezerowych w każdym skupieniu k-średnich
```

**Rozwiązanie.** Skupienia nie odpowiadają klasom: w każdym skupieniu jest prawie tyle samo
próbek control i T2D (k-średnich: 102 i 95 w większym skupieniu, 72 i 75 w mniejszym), a ARI
z klasą jest bliski zeru (−0,002 dla k-średnich, 0,003 dla metody Warda). Obie metody znajdują
za to podobny podział (ARI między nimi 0,74), więc w danych istnieje powtarzalny podział na dwie
grupy, tylko niezwiązany z chorobą.

`rowSums(Xf > 0)` liczy, ile ze 190 cech wykryto w każdej próbce. Mediana tej liczby wynosi 89
w większym skupieniu i 69 w mniejszym: typowa próbka z pierwszego skupienia ma wykrytych
89 cech, a z drugiego 69. Metoda k-średnich podzieliła więc próbki na „bogatsze” i „uboższe”
w cechy, czyli według tego samego, co opisuje PC1 (Zadanie 2), a nie według choroby. Analiza
nienadzorowana znajduje najsilniejszą strukturę w danych, a ta nie musi mieć związku z klasą,
którą chcemy przewidywać.