Spis treści

Wprowadzenie
Opis zbioru danych
Przygotowanie danych
Eksploracyjna analiza danych (EDA)
- 4.1 Statystyki opisowe
- 4.2 Histogramy rozkładów zmiennych
- 4.3 Analiza korelacji między zmiennymi
Klastrowanie danych w oryginalnej przestrzeni zmiennych
- 5.1 Standaryzacja danych
- 5.2 Wybór optymalnej liczby klastrów
- 5.3 Klastrowanie metodą k-means
- 5.4 Klastrowanie metodą PAM (Partitioning Around Medoids)
- 5.5 Klastrowanie hierarchiczne
Redukcja wymiarów
- 6.1 Principal Component Analysis (PCA) — teoria
- 6.2 Wyniki PCA i interpretacja składowych
- 6.3 Wizualizacja komponentów głównych
- 6.4 Multidimensional Scaling (MDS)
Klastrowanie po redukcji wymiarów
- 7.1 Klastrowanie na komponentach PCA
- 7.2 Wizualizacja klastrów w przestrzeni PCA
Porównanie wyników klastrowania
- 8.1 Indeksy jakości klastrowania (silhouette)
- 8.2 Zestawienie wyników
Wnioski

Wprowadzenie

Uczenie maszynowe nienadzorowane (ang. unsupervised learning) obejmuje metody analizy danych, które pozwalają na odkrywanie ukrytych struktur w zbiorach danych bez wykorzystania wcześniej zdefiniowanych etykiet klas. W przeciwieństwie do uczenia nadzorowanego, gdzie model uczy się na podstawie znanych odpowiedzi, metody nienadzorowane koncentrują się na identyfikacji naturalnych wzorców i zależności występujących w danych.

Do najważniejszych metod uczenia nienadzorowanego należą:

Techniki klastrowania — grupowanie obserwacji tak, aby obiekty w tej samej grupie były jak najbardziej do siebie podobne, a różniły się od obiektów z innych grup. Wyróżniamy metody podziałowe (k-means, PAM) oraz hierarchiczne (aglomeracyjne, podziałowe).
Metody redukcji wymiarów — umożliwiają reprezentację danych wielowymiarowych w przestrzeni o mniejszej liczbie wymiarów przy zachowaniu jak największej ilości informacji. Należą do nich m.in. PCA, MDS, t-SNE czy UMAP.

Celem niniejszego projektu jest kompleksowa analiza zbioru danych opisujących właściwości chemiczne win przy użyciu metod nienadzorowanych. Projekt obejmuje:

Eksploracyjną analizę danych (EDA).
Klastrowanie w oryginalnej przestrzeni zmiennych metodami k-means, PAM i hierarchiczną.
Redukcję wymiarów metodami PCA i MDS.
Klastrowanie po redukcji wymiarów i porównanie uzyskanych wyników.

# =============================================================================
# BLOK 2: Ładowanie bibliotek
# Oddzielony od instalacji — gwarantuje poprawne załadowanie po instalacji.
# =============================================================================


library(tidyverse)   # meta-pakiet: ggplot2, dplyr, tidyr, readr i inne
library(cluster)     # PAM, silhouette
library(ggrepel)
library(factoextra)  # fviz_pca_*(), fviz_cluster(), fviz_nbclust()
library(corrplot)    # corrplot()
library(NbClust)     # NbClust()
library(gridExtra)
# Ustawienie globalnego ziarna losowości dla reprodukowalności wyników
set.seed(123)

Opis zbioru danych

W analizie wykorzystano zbiór danych WineQT zawierający informacje o właściwościach fizykochemicznych czerwonych win. Każda obserwacja reprezentuje jedno wino, a zmienne opisują jego parametry chemiczne mierzone laboratoryjnie.

Zmienne w zbiorze danych:

Zmienna	Opis
`fixed.acidity`	Kwasowość stała (g/dm³)
`volatile.acidity`	Kwasowość lotna (g/dm³) — wysoka wartość = ocet
`citric.acid`	Kwas cytrynowy (g/dm³)
`residual.sugar`	Cukier resztkowy (g/dm³)
`chlorides`	Chlorki — zawartość soli (g/dm³)
`free.sulfur.dioxide`	Wolny dwutlenek siarki (mg/dm³)
`total.sulfur.dioxide`	Całkowity dwutlenek siarki (mg/dm³)
`density`	Gęstość (g/cm³)
`pH`	Wartość pH
`sulphates`	Siarczany (g/dm³)
`alcohol`	Zawartość alkoholu (% obj.)
`quality`	Ocena jakości wina (skala 0–10)
`Id`	Identyfikator obserwacji (usuwany)

# Wczytanie zbioru danych z pliku CSV
wine <- read.csv("WineQT.csv")

# Podgląd pierwszych 6 wierszy
head(wine)

##   fixed.acidity volatile.acidity citric.acid residual.sugar chlorides
## 1           7.4             0.70        0.00            1.9     0.076
## 2           7.8             0.88        0.00            2.6     0.098
## 3           7.8             0.76        0.04            2.3     0.092
## 4          11.2             0.28        0.56            1.9     0.075
## 5           7.4             0.70        0.00            1.9     0.076
## 6           7.4             0.66        0.00            1.8     0.075
##   free.sulfur.dioxide total.sulfur.dioxide density   pH sulphates alcohol
## 1                  11                   34  0.9978 3.51      0.56     9.4
## 2                  25                   67  0.9968 3.20      0.68     9.8
## 3                  15                   54  0.9970 3.26      0.65     9.8
## 4                  17                   60  0.9980 3.16      0.58     9.8
## 5                  11                   34  0.9978 3.51      0.56     9.4
## 6                  13                   40  0.9978 3.51      0.56     9.4
##   quality Id
## 1       5  0
## 2       5  1
## 3       5  2
## 4       6  3
## 5       5  4
## 6       5  5

# Liczba wierszy (obserwacji) i kolumn (zmiennych)
cat("Wymiary zbioru danych:", dim(wine)[1], "obserwacji,", dim(wine)[2], "zmiennych\n")

## Wymiary zbioru danych: 1143 obserwacji, 13 zmiennych

# Lista nazw zmiennych
cat("Nazwy zmiennych:\n")

## Nazwy zmiennych:

print(names(wine))

##  [1] "fixed.acidity"        "volatile.acidity"     "citric.acid"         
##  [4] "residual.sugar"       "chlorides"            "free.sulfur.dioxide" 
##  [7] "total.sulfur.dioxide" "density"              "pH"                  
## [10] "sulphates"            "alcohol"              "quality"             
## [13] "Id"

Przygotowanie danych

Przed przystąpieniem do analizy należy oczyścić dane: usunąć zmienną identyfikacyjną (Id), która nie niesie informacji merytorycznej, oraz sprawdzić kompletność danych.

# Usunięcie zmiennej Id — pełni tylko rolę identyfikatora, nie jest cechą chemiczną
wine <- wine[, !names(wine) %in% "Id"]

cat("Liczba brakujących wartości (NA) w zbiorze:", sum(is.na(wine)), "\n")

## Liczba brakujących wartości (NA) w zbiorze: 0

# Wymiary po usunięciu kolumny
cat("Wymiary po czyszczeniu:", dim(wine)[1], "x", dim(wine)[2], "\n")

## Wymiary po czyszczeniu: 1143 x 12

Brak wartości brakujących potwierdza kompletność zbioru — nie jest wymagana imputacja.

Eksploracyjna analiza danych (EDA)

Statystyki opisowe

Statystyki opisowe pozwalają ocenić zakres, rozkład oraz potencjalne wartości odstające dla każdej zmiennej. Widoczne są znaczące różnice w skalach zmiennych (np. total.sulfur.dioxide osiąga wartości rzędu setek, podczas gdy citric.acid mieści się w przedziale 0–1), co uzasadnia standaryzację przed klastrowaniem.

# Podstawowe statystyki opisowe dla wszystkich zmiennych
summary(wine)

##  fixed.acidity    volatile.acidity  citric.acid     residual.sugar  
##  Min.   : 4.600   Min.   :0.1200   Min.   :0.0000   Min.   : 0.900  
##  1st Qu.: 7.100   1st Qu.:0.3925   1st Qu.:0.0900   1st Qu.: 1.900  
##  Median : 7.900   Median :0.5200   Median :0.2500   Median : 2.200  
##  Mean   : 8.311   Mean   :0.5313   Mean   :0.2684   Mean   : 2.532  
##  3rd Qu.: 9.100   3rd Qu.:0.6400   3rd Qu.:0.4200   3rd Qu.: 2.600  
##  Max.   :15.900   Max.   :1.5800   Max.   :1.0000   Max.   :15.500  
##    chlorides       free.sulfur.dioxide total.sulfur.dioxide    density      
##  Min.   :0.01200   Min.   : 1.00       Min.   :  6.00       Min.   :0.9901  
##  1st Qu.:0.07000   1st Qu.: 7.00       1st Qu.: 21.00       1st Qu.:0.9956  
##  Median :0.07900   Median :13.00       Median : 37.00       Median :0.9967  
##  Mean   :0.08693   Mean   :15.62       Mean   : 45.91       Mean   :0.9967  
##  3rd Qu.:0.09000   3rd Qu.:21.00       3rd Qu.: 61.00       3rd Qu.:0.9978  
##  Max.   :0.61100   Max.   :68.00       Max.   :289.00       Max.   :1.0037  
##        pH          sulphates         alcohol         quality     
##  Min.   :2.740   Min.   :0.3300   Min.   : 8.40   Min.   :3.000  
##  1st Qu.:3.205   1st Qu.:0.5500   1st Qu.: 9.50   1st Qu.:5.000  
##  Median :3.310   Median :0.6200   Median :10.20   Median :6.000  
##  Mean   :3.311   Mean   :0.6577   Mean   :10.44   Mean   :5.657  
##  3rd Qu.:3.400   3rd Qu.:0.7300   3rd Qu.:11.10   3rd Qu.:6.000  
##  Max.   :4.010   Max.   :2.0000   Max.   :14.90   Max.   :8.000

Histogramy rozkładów zmiennych

Histogramy pokazują kształt rozkładu każdej zmiennej. Większość cech chemicznych wykazuje rozkład prawostronnie skośny (długi ogon w prawo), co jest typowe dla danych fizykochemicznych. Zmienna quality przyjmuje wartości dyskretne w wąskim przedziale.

# Przekształcenie do formatu długiego (tidy) i rysowanie histogramów
wine %>%
  pivot_longer(cols = everything(), names_to = "zmienna", values_to = "wartosc") %>%
  ggplot(aes(x = wartosc)) +
  geom_histogram(bins = 30, fill = "steelblue", color = "white", alpha = 0.85) +
  facet_wrap(~zmienna, scales = "free") +    # osobna skala osi X dla każdej zmiennej
  labs(
    title = "Rozkłady zmiennych fizykochemicznych win",
    x = "Wartość", y = "Liczba obserwacji"
  ) +
  theme_minimal(base_size = 11) +
  theme(strip.text = element_text(face = "bold"))

Rozkłady wszystkich zmiennych w zbiorze danych o winach

Analiza korelacji między zmiennymi

Macierz korelacji pozwala zidentyfikować liniowe zależności między zmiennymi. Silna korelacja między parą zmiennych może świadczyć o redundancji informacji i przemawia za zastosowaniem PCA.

Zauważalne korelacje: - free.sulfur.dioxide ~ total.sulfur.dioxide — silna korelacja dodatnia (oczekiwana) - fixed.acidity ~ density — korelacja dodatnia - fixed.acidity ~ pH — korelacja ujemna (wyższe kwasowość → niższe pH) - alcohol ~ density — korelacja ujemna

# Obliczenie macierzy korelacji Pearsona
cor_matrix <- cor(wine)

# Wizualizacja macierzy — kolor niebieski = korelacja dodatnia, czerwony = ujemna
corrplot(
  cor_matrix,
  method  = "color",          # kolorowe prostokąty
  type    = "upper",          # tylko górna trójkąt (symetria)
  order   = "hclust",         # grupowanie zmiennych wg podobieństwa
  addCoef.col = "black",      # wyświetlaj wartości współczynników
  number.cex  = 0.7,          # rozmiar czcionki współczynników
  tl.cex      = 0.9,          # rozmiar etykiet zmiennych
  tl.col      = "black",
  title = "Macierz korelacji Pearsona",
  mar   = c(0, 0, 2, 0)
)

Macierz korelacji Pearsona między zmiennymi

Obecność silnych korelacji potwierdza, że dane zawierają redundantną informację. PCA pozwoli skondensować tę informację do mniejszej liczby niezależnych składowych.

Klastrowanie danych w oryginalnej przestrzeni zmiennych

Standaryzacja danych

Przed klastrowaniem niezbędna jest standaryzacja zmiennych do średniej 0 i odchylenia standardowego 1 (transformacja Z-score). Bez standaryzacji zmienne o większych zakresach wartości (np. total.sulfur.dioxide) dominowałyby w obliczeniach odległości euklidesowych, zniekształcając wyniki klastrowania.

# Standaryzacja: każda zmienna -> srednia=0, sd=1
# Używamy tylko zmiennych chemicznych (bez quality, która jest zmienną docelową)
wine_features <- wine %>% select(-quality)
wine_scaled   <- scale(wine_features)

# Weryfikacja standaryzacji
cat("Średnie po standaryzacji (powinny być ~0):\n")

## Średnie po standaryzacji (powinny być ~0):

round(colMeans(wine_scaled), 3)

##        fixed.acidity     volatile.acidity          citric.acid 
##                    0                    0                    0 
##       residual.sugar            chlorides  free.sulfur.dioxide 
##                    0                    0                    0 
## total.sulfur.dioxide              density                   pH 
##                    0                    0                    0 
##            sulphates              alcohol 
##                    0                    0

cat("\nOdchylenia standardowe (powinny być ~1):\n")

## 
## Odchylenia standardowe (powinny być ~1):

round(apply(wine_scaled, 2, sd), 3)

##        fixed.acidity     volatile.acidity          citric.acid 
##                    1                    1                    1 
##       residual.sugar            chlorides  free.sulfur.dioxide 
##                    1                    1                    1 
## total.sulfur.dioxide              density                   pH 
##                    1                    1                    1 
##            sulphates              alcohol 
##                    1                    1

Wybór optymalnej liczby klastrów

Dobór liczby klastrów k to kluczowy krok, od którego zależą wyniki całej analizy. Stosujemy trzy uzupełniające się metody:

Metoda łokcia (Elbow method) — analizuje sumę kwadratów wewnątrz klastrów (WSS) w funkcji k. Szukamy punktu, w którym dalsze zwiększanie k nie przynosi istotnej poprawy (charakterystyczny “łokieć”).

# Metoda łokcia: suma WSS dla k = 1..10
fviz_nbclust(
  wine_scaled, 
  kmeans, 
  method = "wss",
  k.max  = 10
) +
  labs(
    title    = "Metoda łokcia (Elbow Method)",
    subtitle = "Szukamy punktu przegięcia na krzywej WSS",
    x = "Liczba klastrów k",
    y = "Suma kwadratów wewnątrz klastrów (WSS)"
  ) +
  theme_minimal()

Metoda łokcia — wybór liczby klastrów

Metoda sylwetkowa (Silhouette) — mierzy, jak dobrze każda obserwacja pasuje do swojego klastra w porównaniu z sąsiednimi klastrami. Wartości bliskie 1 oznaczają dobrze przypisane obserwacje, bliskie 0 — obserwacje na granicy klastrów.

# Metoda sylwetkowa: średnia szerokość sylwetki dla k = 2..10
fviz_nbclust(
  wine_scaled, 
  kmeans, 
  method = "silhouette",
  k.max  = 10
) +
  labs(
    title    = "Metoda sylwetkowa (Average Silhouette Width)",
    subtitle = "Wyższa wartość = lepiej rozdzielone klastry",
    x = "Liczba klastrów k",
    y = "Średnia szerokość sylwetki"
  ) +
  theme_minimal()

Metoda sylwetki — wybór liczby klastrów

# NbClust analizuje 30 różnych indeksów i rekomenduje k przez głosowanie większościowe
# Parametr: distance="euclidean", method="kmeans"
nb_result <- NbClust(
  data      = wine_scaled,
  distance  = "euclidean",
  min.nc    = 2,
  max.nc    = 6,
  method    = "kmeans",
  index     = "all"
)

# Wyświetlenie rekomendowanej liczby klastrów wg NbClust
cat("Rekomendowana liczba klastrów wg NbClust:", 
    nb_result$Best.nc["Number_clusters"], "\n")

## Rekomendowana liczba klastrów wg NbClust: NA

porównanie różnej ilości klastrów

# Testujemy k = 3, 5, 6 i porównujemy indeks sylwetki oraz WSS
k_values <- c(3, 5, 6)

# Funkcja pomocnicza — klastrowanie i obliczenie metryk dla danego k
test_k <- function(k) {
  set.seed(123)
  km <- kmeans(wine_scaled, centers = k, nstart = 25)
  sil <- silhouette(km$cluster, dist(wine_scaled))
  data.frame(
    k                 = k,
    avg_silhouette    = round(mean(sil[, 3]), 4),
    wss               = round(km$tot.withinss, 1),
    between_ratio     = round(km$betweenss / km$totss * 100, 1)
  )
}

# Zestawienie wyników dla wszystkich testowanych k
comparison_df <- do.call(rbind, lapply(k_values, test_k))

knitr::kable(
  comparison_df,
  caption = "Porównanie jakości klastrowania dla różnych k",
  col.names = c("k (liczba klastrów)", "Śr. sylwetka", "WSS", "Wyj. wariancja (%)")
)

Porównanie jakości klastrowania dla różnych k
k (liczba klastrów)	Śr. sylwetka	WSS	Wyj. wariancja (%)
3	0.1846	9051.1	27.9
5	0.1885	7249.6	42.3
6	0.1918	6672.3	46.9

# Generowanie wykresów klastrów dla każdego k obok siebie
library(gridExtra)

plots <- lapply(k_values, function(k) {
  set.seed(123)
  km <- kmeans(wine_scaled, centers = k, nstart = 25)
  fviz_cluster(
    km,
    data         = wine_scaled,
    geom         = "point",
    ellipse.type = "convex",
    palette      = "jco",
    ggtheme      = theme_minimal(),
    main         = paste0("k = ", k)
  )
})

grid.arrange(grobs = plots, ncol = 3)

Wizualizacja klastrów dla k=3, 5 i 6

# Wykresy sylwetki dla każdego k — pozwalają zobaczyć rozkład jakości przypisań
sil_plots <- lapply(k_values, function(k) {
  set.seed(123)
  km  <- kmeans(wine_scaled, centers = k, nstart = 25)
  sil <- silhouette(km$cluster, dist(wine_scaled))
  fviz_silhouette(sil, palette = "jco", ggtheme = theme_minimal(), print.summary = FALSE) +
    labs(title = paste0("Sylwetka — k = ", k)) +
    theme(legend.position = "none")
})

grid.arrange(grobs = sil_plots, ncol = 3)

Wykresy sylwetki dla k=3, 5 i 6


Na podstawie przeprowadzonych analiz wybór optymalnej liczby klastrów nie jest 
jednoznaczny. Metoda łokcia wskazuje na punkt przegięcia przy k=3, natomiast 
metoda sylwetkowa i NbClust również rekomendują k=3 jako punkt wyjścia.

Testując dodatkowo k=5 i k=6 zaobserwowano, że średnia sylwetka rośnie wraz 
z liczbą klastrów (k=3: 0.184, k=5: 0.188, k=6: 0.191), jednak przyrost jest 
marginalny — różnica między k=3 a k=6 wynosi zaledwie 0.007. Warto podkreślić, 
że wszystkie uzyskane wartości sylwetki mieszczą się w przedziale 0.18–0.20, 
co świadczy o tym, że dane nie posiadają wyraźnie rozdzielonych naturalnych 
grup — klastry częściowo nachodzą na siebie, co jest typowe dla ciągłych 
danych fizykochemicznych.

Zwiększanie liczby klastrów do k=5 lub k=6 poprawia dopasowanie statystyczne 
(wyjaśniona wariancja wzrasta z 27.9% do 46.9%), jednak kosztem interpretowalności 
— przy k=6 klastry stają się bardzo małe i trudne do opisania merytorycznie.

Z powyższych względów przyjmujemy **k = 3** jako optymalną liczbę klastrów, 
kierując się zasadą parsymonii: przy braku istotnej poprawy jakości klastrowania 
preferujemy prostsze rozwiązanie, które umożliwia sensowną interpretację 
wyodrębnionych grup win

## Klastrowanie metodą k-means

Algorytm k-means działa iteracyjnie: losuje k centrów, przypisuje każdą obserwację 
do najbliższego centrum, aktualizuje centra jako środki ciężkości klastrów i powtarza 
proces do osiągnięcia zbieżności. Parametr `nstart = 25` oznacza 25 losowych startów — 
wybierany jest najlepszy wynik, co ogranicza wpływ inicjalizacji.

``` r
# Klastrowanie k-means: k=3, 25 losowych startów dla stabilności wyników
set.seed(123)
kmeans_result <- kmeans(wine_scaled, centers = 3, nstart = 25)

# Liczebność klastrów
cat("Liczebność klastrów k-means:\n")

## Liczebność klastrów k-means:

print(table(kmeans_result$cluster))

## 
##   1   2   3 
## 530 264 349

# Odsetek wyjaśnionej wariancji (between_SS / total_SS)
cat("\nOdsetek wyjaśnionej wariancji:",
    round(kmeans_result$betweenss / kmeans_result$totss * 100, 1), "%\n")

## 
## Odsetek wyjaśnionej wariancji: 27.9 %

# Wizualizacja klastrów w przestrzeni 2D (factoextra używa PCA do redukcji wymiarów)
fviz_cluster(
  kmeans_result,
  data         = wine_scaled,
  geom         = "point",
  ellipse.type = "convex",        # otoczka wypukła wokół klastra
  palette      = "jco",
  ggtheme      = theme_minimal(),
  main         = "Klastry k-means (k=3) — wizualizacja PCA"
)

Klastry k-means zwizualizowane na dwóch pierwszych składowych PCA

Klastrowanie metodą PAM

PAM (Partitioning Around Medoids) jest wariantem k-means bardziej odpornym na wartości odstające. W odróżnieniu od k-means, centra klastrów (medoidy) muszą być rzeczywistymi obserwacjami ze zbioru danych, co czyni wyniki bardziej interpretowalnym.

# Klastrowanie PAM: k=3
set.seed(123)
pam_result <- pam(wine_scaled, k = 3)

# Liczebność klastrów
cat("Liczebność klastrów PAM:\n")

## Liczebność klastrów PAM:

print(table(pam_result$clustering))

## 
##   1   2   3 
## 519 260 364

# Indeks sylwetki dla PAM
cat("\nŚrednia szerokość sylwetki (PAM):",
    round(pam_result$silinfo$avg.width, 3), "\n")

## 
## Średnia szerokość sylwetki (PAM): 0.188

fviz_cluster(
  pam_result,
  data         = wine_scaled,
  geom         = "point",
  ellipse.type = "convex",
  palette      = "jco",
  ggtheme      = theme_minimal(),
  main         = "Klastry PAM (k=3) — wizualizacja PCA"
)

Klastry PAM zwizualizowane na dwóch pierwszych składowych PCA

Klastrowanie hierarchiczne

Klastrowanie hierarchiczne (aglomeracyjne) nie wymaga z góry podanej liczby klastrów — buduje drzewo połączeń (dendrogram), które następnie “tniemy” na pożądanej wysokości. Używamy metody Warda, która minimalizuje wariancję wewnątrz klastrów.

# Obliczenie macierzy odległości euklidesowych
dist_matrix <- dist(wine_scaled, method = "euclidean")

# Klastrowanie hierarchiczne metodą Warda (minimalizacja wariancji)
hclust_result <- hclust(dist_matrix, method = "ward.D2")

# Rysowanie dendrogramu z zaznaczeniem 3 klastrów
plot(
  hclust_result,
  labels = FALSE,    # ukrywamy etykiety (1142 obserwacji = nieczytelne)
  hang   = -1,
  main   = "Dendrogram — klastrowanie hierarchiczne (Ward.D2)",
  xlab   = "Obserwacje", ylab = "Odległość (wysokość połączenia)"
)
rect.hclust(hclust_result, k = 3, border = c("red", "blue", "green"))

Dendrogram klastrowania hierarchicznego (metoda Warda)

# Wycinamy 3 klastry z dendrogramu
hclust_clusters <- cutree(hclust_result, k = 3)

cat("Liczebność klastrów hierarchicznych:\n")

## Liczebność klastrów hierarchicznych:

print(table(hclust_clusters))

## hclust_clusters
##   1   2   3 
## 588 320 235

Redukcja wymiarów

Principal Component Analysis (PCA) — teoria i zastosowanie

PCA (Analiza Składowych Głównych) to liniowa technika redukcji wymiarów. Polega na znalezieniu nowego układu ortogonalnych osi (składowych głównych, PC), które maksymalizują wyjaśnianą wariancję. Pierwsza składowa PC1 wyjaśnia najwięcej wariancji, kolejne — coraz mniej.

Matematycznie: PCA rozwiązuje zagadnienie własne macierzy kowariancji standaryzowanych danych, gdzie wektory własne wyznaczają kierunki składowych, a wartości własne — ilość wyjaśnianej wariancji.

# PCA na standaryzowanych danych (scale.=FALSE bo już standaryzowaliśmy)
# Używamy wine_scaled bez zmiennej quality
pca_result <- prcomp(wine_scaled, scale. = FALSE)

# Podsumowanie: odchylenia standardowe i proporcje wyjaśnionej wariancji
summary(pca_result)

## Importance of components:
##                           PC1    PC2    PC3    PC4     PC5     PC6     PC7
## Standard deviation     1.7766 1.3705 1.2545 1.1007 0.97888 0.81570 0.74672
## Proportion of Variance 0.2869 0.1708 0.1431 0.1101 0.08711 0.06049 0.05069
## Cumulative Proportion  0.2869 0.4577 0.6007 0.7109 0.79798 0.85847 0.90916
##                           PC8     PC9    PC10    PC11
## Standard deviation     0.6473 0.58694 0.42099 0.24178
## Proportion of Variance 0.0381 0.03132 0.01611 0.00531
## Cumulative Proportion  0.9473 0.97857 0.99469 1.00000

Wyniki PCA i interpretacja składowych

# Scree plot: procent wyjaśnionej wariancji dla każdej składowej
fviz_eig(
  pca_result,
  addlabels = TRUE,       # dodaj procenty na słupkach
  ylim      = c(0, 50),
  barfill   = "steelblue",
  barcolor  = "white",
  linecolor = "darkred",
  ggtheme   = theme_minimal()
) +
  labs(
    title    = "Scree Plot — Wariancja wyjaśniona przez składowe PCA",
    subtitle = "Szukamy 'łokcia' — punktu, po którym przyrost wariancji jest nieistotny",
    x = "Składowa główna",
    y = "Procent wyjaśnionej wariancji (%)"
  )

Wykres osypiska (Scree Plot) — udział wariancji wyjaśnionej przez składowe

# Obliczenie skumulowanej wyjaśnionej wariancji
variance_explained <- summary(pca_result)$importance[2, ]   # proporcja
cumulative_var     <- summary(pca_result)$importance[3, ]   # skumulowana

cat("Wariancja wyjaśniona przez kolejne składowe:\n")

## Wariancja wyjaśniona przez kolejne składowe:

print(round(variance_explained * 100, 1))

##  PC1  PC2  PC3  PC4  PC5  PC6  PC7  PC8  PC9 PC10 PC11 
## 28.7 17.1 14.3 11.0  8.7  6.0  5.1  3.8  3.1  1.6  0.5

cat("\nSkumulowana wariancja wyjaśniona:\n")

## 
## Skumulowana wariancja wyjaśniona:

print(round(cumulative_var * 100, 1))

##   PC1   PC2   PC3   PC4   PC5   PC6   PC7   PC8   PC9  PC10  PC11 
##  28.7  45.8  60.1  71.1  79.8  85.8  90.9  94.7  97.9  99.5 100.0

# Ile składowych potrzeba do wyjaśnienia >= 80% wariancji?
n_components_80 <- which(cumulative_var >= 0.80)[1]
cat("\nLiczba składowych do wyjaśnienia >= 80% wariancji:", n_components_80, "\n")

## 
## Liczba składowych do wyjaśnienia >= 80% wariancji: 6

Wizualizacja komponentów PCA

# Biplot: jednocześnie pokazuje obserwacje (punkty) i zmienne (strzałki)
# Długość strzałki = siła wpływu zmiennej na daną składową
# Kierunek strzałki = interpretacja składowej
fviz_pca_biplot(
  pca_result,
  repel      = TRUE,          # ggrepel: etykiety bez nakładania
  geom.ind   = "point",
  col.ind    = "steelblue",
  col.var    = "darkred",
  alpha.ind  = 0.4,
  label      = "var",         # tylko etykiety zmiennych
  ggtheme    = theme_minimal()
) +
  labs(
    title    = "Biplot PCA (PC1 vs PC2)",
    subtitle = "Strzałki = kierunki zmiennych; punkty = obserwacje"
  )

Biplot PCA — obserwacje i ładunki zmiennych

# Ładunki (rotacje) — wkład każdej zmiennej w składowe PC1 i PC2
loadings_df <- as.data.frame(pca_result$rotation[, 1:2])
loadings_df$variable <- rownames(loadings_df)

# Wizualizacja ładunków jako wykres słupkowy
loadings_df %>%
  pivot_longer(cols = c(PC1, PC2), names_to = "Component", values_to = "Loading") %>%
  ggplot(aes(x = reorder(variable, Loading), y = Loading, fill = Loading > 0)) +
  geom_col(show.legend = FALSE) +
  facet_wrap(~Component) +
  coord_flip() +
  scale_fill_manual(values = c("TRUE" = "steelblue", "FALSE" = "tomato")) +
  labs(
    title    = "Ładunki zmiennych dla PC1 i PC2",
    subtitle = "Niebieski = ładunek dodatni, czerwony = ujemny",
    x = "Zmienna", y = "Ładunek (loading)"
  ) +
  theme_minimal()

Ładunki zmiennych dla PC1 i PC2

Multidimensional Scaling (MDS)

MDS (Wielowymiarowe Skalowanie) to alternatywna metoda redukcji wymiarów. W odróżnieniu od PCA, MDS bazuje bezpośrednio na macierzy odległości między obserwacjami i stara się zachować te odległości w przestrzeni 2D. Przy odległościach euklidesowych na standaryzowanych danych wyniki są zbliżone do PCA.

# Klasyczne MDS (cmdscale) na macierzy odległości euklidesowych
# Używamy próbki 300 obserwacji dla wydajności (dist() na pełnym zbiorze = duża macierz)
set.seed(123)
sample_idx <- sample(nrow(wine_scaled), 300)
dist_sample <- dist(wine_scaled[sample_idx, ], method = "euclidean")

mds_result <- cmdscale(dist_sample, k = 2, eig = TRUE)

# Procent wyjaśnionej wariancji przez MDS
mds_var <- round(mds_result$eig[1:2] / sum(abs(mds_result$eig)) * 100, 1)
cat("Wariancja wyjaśniona przez MDS: Dim1 =", mds_var[1], "%, Dim2 =", mds_var[2], "%\n")

## Wariancja wyjaśniona przez MDS: Dim1 = 27.3 %, Dim2 = 18.8 %

# Wizualizacja MDS
mds_df <- as.data.frame(mds_result$points)
colnames(mds_df) <- c("Dim1", "Dim2")

ggplot(mds_df, aes(x = Dim1, y = Dim2)) +
  geom_point(alpha = 0.5, color = "steelblue", size = 1.5) +
  labs(
    title    = "Multidimensional Scaling (MDS) — próbka 300 obserwacji",
    subtitle = paste0("Dim1: ", mds_var[1], "% | Dim2: ", mds_var[2], "% wyjaśnionej wariancji"),
    x = "Wymiar 1", y = "Wymiar 2"
  ) +
  theme_minimal()

Wizualizacja MDS — obserwacje w przestrzeni 2D

Klastrowanie po redukcji wymiarów

Klastrowanie na komponentach PCA

Po redukcji wymiarów do pierwszych k składowych PCA przeprowadzamy ponowne klastrowanie. Używamy składowych wyjaśniających łącznie >= 80% wariancji. Klastrowanie w przestrzeni PCA jest często skuteczniejsze, gdyż:

eliminuje szum i kolinearność obecną w oryginalnych zmiennych,
redukuje efekt “klątwy wymiarowości”,
składowe PCA są ortogonalne (niezależne od siebie).

# Wyodrębnienie pierwszych n_components_80 składowych PCA
pca_scores <- pca_result$x[, 1:n_components_80]

cat("Klastrowanie w przestrzeni", n_components_80, "składowych PCA\n")

## Klastrowanie w przestrzeni 6 składowych PCA

cat("(wyjaśniają łącznie", round(cumulative_var[n_components_80] * 100, 1), "% wariancji)\n\n")

## (wyjaśniają łącznie 85.8 % wariancji)

# k-means na składowych PCA
set.seed(123)
kmeans_pca <- kmeans(pca_scores, centers = 3, nstart = 25)

cat("Liczebność klastrów k-means na PCA:\n")

## Liczebność klastrów k-means na PCA:

print(table(kmeans_pca$cluster))

## 
##   1   2   3 
## 348 288 507

cat("\nOdsetek wyjaśnionej wariancji:",
    round(kmeans_pca$betweenss / kmeans_pca$totss * 100, 1), "%\n")

## 
## Odsetek wyjaśnionej wariancji: 32.4 %

Wizualizacja klastrów w przestrzeni PCA

# Wizualizacja klastrów PCA z użyciem fviz_cluster
fviz_cluster(
  kmeans_pca,
  data         = pca_scores,
  geom         = "point",
  ellipse.type = "norm",        # elipsy normalne (rozkład Gaussa)
  palette      = "jco",
  ggtheme      = theme_minimal(),
  main         = "Klastry k-means w przestrzeni PCA"
) +
  labs(subtitle = paste0("PC1–PC2 (", 
                          round(cumulative_var[2] * 100, 1), 
                          "% wariancji)"))

Klastry k-means w przestrzeni pierwszych dwóch składowych PCA

# Dodanie etykiet klastrów do oryginalnych (niestandaryzowanych) danych
wine_clustered <- wine_features %>%
  mutate(Klaster = factor(kmeans_pca$cluster))

# Obliczenie średnich dla każdego klastra i zmiennej
cluster_profiles <- wine_clustered %>%
  group_by(Klaster) %>%
  summarise(across(everything(), mean), .groups = "drop") %>%
  pivot_longer(-Klaster, names_to = "Zmienna", values_to = "Srednia")

# Normalizacja do [0,1] dla czytelnego porównania na jednym wykresie
cluster_profiles <- cluster_profiles %>%
  group_by(Zmienna) %>%
  mutate(Srednia_norm = (Srednia - min(Srednia)) / (max(Srednia) - min(Srednia) + 1e-10)) %>%
  ungroup()

ggplot(cluster_profiles, aes(x = Zmienna, y = Srednia_norm, fill = Klaster)) +
  geom_col(position = "dodge") +
  scale_fill_brewer(palette = "Set1") +
  labs(
    title    = "Profile klastrów — znormalizowane średnie wartości zmiennych",
    subtitle = "Porównanie 3 klastrów dla każdej cechy chemicznej",
    x = "Zmienna chemiczna", y = "Znormalizowana średnia [0–1]"
  ) +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

Profile chemiczne klastrów — średnie wartości zmiennych w każdym klastrze

Porównanie wyników klastrowania

Indeksy jakości klastrowania (Silhouette)

Indeks sylwetki (silhouette) jest jedną z najczęściej stosowanych miar jakości klastrowania. Dla każdej obserwacji i obliczany jest wskaźnik:

\[s(i) = \frac{b(i) - a(i)}{\max(a(i), b(i))}\]

gdzie: a(i) = średnia odległość do obserwacji w tym samym klastrze (spójność), b(i) = minimalna średnia odległość do obserwacji z sąsiedniego klastra (separacja).

Wartości s(i) bliskie 1 oznaczają dobre przypisanie, bliskie 0 — obserwacje na granicy.

# Sylwetka dla k-means na oryginalnych danych
sil_kmeans <- silhouette(kmeans_result$cluster, dist(wine_scaled))

fviz_silhouette(sil_kmeans, palette = "jco", ggtheme = theme_minimal()) +
  labs(title = "Wykres sylwetki — k-means (dane oryginalne)")

##   cluster size ave.sil.width
## 1       1  530          0.25
## 2       2  264          0.12
## 3       3  349          0.13

Wykres sylwetki dla k-means

# Sylwetka dla PAM (wbudowana w obiekt pam_result)
fviz_silhouette(pam_result, palette = "jco", ggtheme = theme_minimal()) +
  labs(title = "Wykres sylwetki — PAM")

##   cluster size ave.sil.width
## 1       1  519          0.26
## 2       2  260          0.08
## 3       3  364          0.16

Wykres sylwetki dla PAM

# Sylwetka dla k-means na danych PCA
sil_pca <- silhouette(kmeans_pca$cluster, dist(pca_scores))

fviz_silhouette(sil_pca, palette = "jco", ggtheme = theme_minimal()) +
  labs(title = "Wykres sylwetki — k-means po redukcji PCA")

##   cluster size ave.sil.width
## 1       1  348          0.15
## 2       2  288          0.18
## 3       3  507          0.27

Wykres sylwetki dla k-means po PCA

Zestawienie wyników

# Podsumowanie indeksów jakości dla wszystkich metod
avg_sil_kmeans <- mean(sil_kmeans[, 3])
avg_sil_pam    <- pam_result$silinfo$avg.width
avg_sil_pca    <- mean(sil_pca[, 3])

results_table <- data.frame(
  Metoda = c(
    "k-means (oryg.)",
    "PAM (oryg.)",
    "k-means po PCA"
  ),
  `Liczba_klastrow` = c(3, 3, 3),
  `Srednia_sylwetka` = round(c(avg_sil_kmeans, avg_sil_pam, avg_sil_pca), 4),
  `Wyjasniona_wariancja` = c(
    paste0(round(kmeans_result$betweenss / kmeans_result$totss * 100, 1), "%"),
    "—",
    paste0(round(kmeans_pca$betweenss / kmeans_pca$totss * 100, 1), "%")
  )
)

knitr::kable(
  results_table,
  caption = "Porównanie metod klastrowania — indeksy jakości",
  col.names = c("Metoda", "Liczba klastrów", "Śr. sylwetka", "Wyj. wariancja")
)

Porównanie metod klastrowania — indeksy jakości
Metoda	Liczba klastrów	Śr. sylwetka	Wyj. wariancja
k-means (oryg.)	3	0.1846	27.9%
PAM (oryg.)	3	0.1883	—
k-means po PCA	3	0.2103	32.4%

# Tabela krzyżowa k-means (oryg.) vs k-means (PCA) — sprawdzamy zgodność przypisań
cat("Zgodność przypisań — k-means oryginalne vs k-means po PCA:\n")

## Zgodność przypisań — k-means oryginalne vs k-means po PCA:

print(table(
  "k-means oryg."  = kmeans_result$cluster,
  "k-means po PCA" = kmeans_pca$cluster
))

##              k-means po PCA
## k-means oryg.   1   2   3
##             1   1  22 507
##             2   2 262   0
##             3 345   4   0

Wnioski

Na podstawie przeprowadzonej analizy można sformułować następujące wnioski:

1. Eksploracja danych:
Zbiór danych zawiera 1142 win opisanych 11 cechami chemicznymi. Zmienne wykazują znaczne zróżnicowanie skal oraz wyraźne korelacje (szczególnie free.sulfur.dioxide i total.sulfur.dioxide), co uzasadniało standaryzację i zastosowanie PCA.

2. Klastrowanie w oryginalnej przestrzeni:
Wszystkie trzy metody (k-means, PAM, hierarchiczna) wskazały na k=3 jako punkt wyjścia, potwierdzony przez metodę łokcia i NbClust. Przeprowadzono również analizę porównawczą dla k=5 i k=6 — choć wyższe k nieznacznie poprawiało indeks sylwetki (wzrost o 0.007 między k=3 a k=6), różnica ta jest statystycznie nieistotna, a wszystkie wartości sylwetki poniżej 0.20 wskazują na brak wyraźnej naturalnej struktury klastrowej w danych. Ostatecznie wybrano k=3 zgodnie z zasadą parsymonii — mniejsza liczba klastrów zapewnia lepszą interpretowalność przy porównywalnej jakości statystycznej. Wyodrębnione trzy grupy win różnią się przede wszystkim zawartością alkoholu, kwasowością i poziomem siarki, co może odpowiadać różnym stylom i przedziałom jakości win.

3. Redukcja wymiarów:
PCA wykazała, że pierwszych 6 składowych wyjaśnia łącznie ponad 80% wariancji danych. Scree Plot potwierdził, że większość informacji koncentruje się w kilku pierwszych składowych, a pozostałe zawierają głównie szum. Interpretacja ładunków pokazuje, że PC1 związana jest z kwasowością i gęstością, zaś PC2 — z zawartością siarki.

4. Klastrowanie po PCA:
Klastrowanie w przestrzeni składowych PCA dało wyniki porównywalne z klastrowaniem na surowych danych, przy jednoczesnym zmniejszeniu wymiarowości i eliminacji kolinearności. Średnia sylwetka wskazuje na umiarkowanie dobrą separację klastrów (typową dla danych ciągłych bez wyraźnych granic).

5. Podsumowanie:
Połączenie PCA z klastrowaniem stanowi efektywne podejście do analizy wielowymiarowych danych chemicznych: PCA redukuje szum i kolinearność, a klastrowanie identyfikuje naturalne grupy win o podobnych właściwościach. Wyniki mogą stanowić punkt wyjścia do dalszej analizy, np. oceny czy wyodrębnione klastry korelują z oceną jakości (quality).

Redukcja wymiarów i klastrowanie danych o winach

Lukasz Gawor

2026-03-15