Celem projektu jest identyfikacja segmentów klientów na podstawie
cech demograficznych i zakupowych z wykorzystaniem metod klasteryzacji
niesuperwizowanej. W analizie porównuję rozwiązania uzyskane
metodami:
k-means (rozwiązanie bazowe),
PAM (odporniejsze na obserwacje odstające – medoidy),
klasteryzacja hierarchiczna (complete oraz Ward.D2),oraz wykonuję
elementy walidacji i praktycznego zastosowania: dobór liczby
klastrów,
porównanie metod, predykcja przynależności nowych obserwacji.
Wykorzystano zbiór danych Mall_Customers.csv, w którym jednostką
analizy jest #pojedynczy klient. Do klastrowania wybrano zmienne:
Age – wiek,
Annual Income (k$) – dochód roczny,
Spending Score (1–100) – wskaźnik intensywności wydatków.
Wykorzystane biblioteki
library(dplyr)
library(ggplot2)
library(factoextra)
library(cluster)
library(flexclust)
library(psych)
library(dendextend)
library(hopkins)
library(mclust)
ustawiam seed w celu powtarzalności wyników
set.seed(123)
Uporządkowanie nazw kolumn (bez spacji i znaków specjalnych)
mall <- read.csv("Mall_Customers.csv", stringsAsFactors = FALSE)
names(mall) <- make.names(names(mall))
glimpse(mall)
## Rows: 200
## Columns: 5
## $ CustomerID <int> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, …
## $ Gender <chr> "Male", "Male", "Female", "Female", "Female", "…
## $ Age <int> 19, 21, 20, 23, 31, 22, 35, 23, 64, 30, 67, 35,…
## $ Annual.Income..k.. <int> 15, 15, 16, 16, 17, 17, 18, 18, 19, 19, 19, 19,…
## $ Spending.Score..1.100. <int> 39, 81, 6, 77, 40, 76, 6, 94, 3, 72, 14, 99, 15…
summary(mall)
## CustomerID Gender Age Annual.Income..k..
## Min. : 1.00 Length:200 Min. :18.00 Min. : 15.00
## 1st Qu.: 50.75 Class :character 1st Qu.:28.75 1st Qu.: 41.50
## Median :100.50 Mode :character Median :36.00 Median : 61.50
## Mean :100.50 Mean :38.85 Mean : 60.56
## 3rd Qu.:150.25 3rd Qu.:49.00 3rd Qu.: 78.00
## Max. :200.00 Max. :70.00 Max. :137.00
## Spending.Score..1.100.
## Min. : 1.00
## 1st Qu.:34.75
## Median :50.00
## Mean :50.20
## 3rd Qu.:73.00
## Max. :99.00
W klasteryzacji opartej o odległości i centroidy konieczne jest
ujednolicenie skali zmiennych. Dlatego dane są standaryzowane (scale()),
aby żadna zmienna (np. dochód) nie dominowała wyniku.
X <- mall %>% select(Age, Annual.Income..k.., Spending.Score..1.100.)
colSums(is.na(X))
## Age Annual.Income..k.. Spending.Score..1.100.
## 0 0 0
str(X)
## 'data.frame': 200 obs. of 3 variables:
## $ Age : int 19 21 20 23 31 22 35 23 64 30 ...
## $ Annual.Income..k.. : int 15 15 16 16 17 17 18 18 19 19 ...
## $ Spending.Score..1.100.: int 39 81 6 77 40 76 6 94 3 72 ...
X_s <- scale(X)
Zanim zastosuję algorytmy klasteryzacji, sprawdzam czy dane wykazują
tendencję do tworzenia klastrów. W tym celu obliczam statystykę
Hopkinsa.
X_hopkins <- as.data.frame(X_s)
stopifnot(sum(!is.finite(as.matrix(X_hopkins))) == 0)
m <- min(50, max(10, floor(nrow(X_hopkins) / 10)))
H <- hopkins(X_hopkins, m = m)
H
## [1] 0.967832
Interpretacja: Uzyskano wartość Hopkinsa 0.967832, co wskazuje na
bardzo silną tendencję klastrową danych.
Dobór liczby klastrów
Liczbę klastrów dobieram na podstawie dwóch kryteriów:
silhouette – ocenia jakość rozdzielenia klastrów,
gap statistic – porównuje wynik do struktury losowej (bardziej
konserwatywne kryterium).
fviz_nbclust(X_s, FUNcluster = kmeans, method = "silhouette") + ggtitle("Silhouette (k-means)")
fviz_nbclust(X_s, FUNcluster = kmeans, method = "gap_stat", nboot = 50) + ggtitle("Gap statistic (k-means)")
Wyniki: silhouette sugeruje k = 8, gap statistic sugeruje k = 6. W
dalszej analizie k = 6 traktuję jako rozwiązanie główne (konserwatywne i
interpretowalne),k = 8 jako wariant porównawczy.
K-means: rozwiązanie główne (k = 6) i wariant porównawczy (k =
8)
K-means dla k = 6
km6 <- eclust(X_s, "kmeans", k = 6, hc_metric = "euclidean")
fviz_cluster( km6, geom = "point", ggtheme = theme_minimal(base_size = 12), main = "K-means clustering (k = 6)" )
fviz_silhouette(km6)
## cluster size ave.sil.width
## 1 1 13 0.22
## 2 2 57 0.37
## 3 3 22 0.56
## 4 4 22 0.31
## 5 5 40 0.54
## 6 6 46 0.26
Profilowanie klastrów (średnie wartości zmiennych w
segmentach):
mall_km6 <- mall %>% mutate(cluster_k6 = factor(km6$cluster))
mall_km6 %>% group_by(cluster_k6) %>% summarise( Age = mean(Age), Income = mean(Annual.Income..k..), Spending = mean(Spending.Score..1.100.), n = n() )
## # A tibble: 6 × 5
## cluster_k6 Age Income Spending n
## <fct> <dbl> <dbl> <dbl> <int>
## 1 1 30.5 89.5 13.4 13
## 2 2 55.5 47.9 41.9 57
## 3 3 25.3 25.7 79.4 22
## 4 4 48.3 87.5 19.6 22
## 5 5 32.9 86.1 81.5 40
## 6 6 27.7 49.6 44.3 46
K-means dla k = 8 (porównanie)
km8 <- eclust(X_s, "kmeans", k = 8, hc_metric = "euclidean")
fviz_cluster( km8, geom = "point", ggtheme = theme_minimal(base_size = 12), main = "K-means clustering (k = 8)" )
fviz_silhouette(km8)
## cluster size ave.sil.width
## 1 1 11 0.29
## 2 2 18 0.35
## 3 3 23 0.49
## 4 4 8 0.30
## 5 5 39 0.48
## 6 6 36 0.39
## 7 7 44 0.46
## 8 8 21 0.29
Profilowanie klastrów (średnie wartości zmiennych w
segmentach):
mall_km8 <- mall %>% mutate(cluster_k8 = factor(km8$cluster))
mall_km8 %>% group_by(cluster_k8) %>% summarise( Age = mean(Age), Income = mean(Annual.Income..k..), Spending = mean(Spending.Score..1.100.), n = n() )
## # A tibble: 8 × 5
## cluster_k8 Age Income Spending n
## <fct> <dbl> <dbl> <dbl> <int>
## 1 1 28 77.7 14.2 11
## 2 2 49.4 82.1 18.4 18
## 3 3 25 25.3 77.6 23
## 4 4 38.8 112. 22.6 8
## 5 5 32.7 86.5 82.1 39
## 6 6 27.2 55.9 50.1 36
## 7 7 56.3 53.7 49.4 44
## 8 8 45.5 26.3 19.4 21
Porównanie jakości k=6 vs k=8 (średnia silhouette)
mean_sil_k6 <- mean(silhouette(km6$cluster, dist(X_s))[, 3])
mean_sil_k8 <- mean(silhouette(km8$cluster, dist(X_s))[, 3])
mean_sil_k6
## [1] 0.3838699
mean_sil_k8
## [1] 0.4105837
Wniosek: k=8 ma nieco wyższą średnią silhouette (0.4105837) w
stosunku do k=6 (0.3838699), jednak ze względu na kryterium gap oraz
interpretowalność segmentów, k = 6 pozostaje rozwiązaniem głównym, a k =
8 pełni rolę porównawczą
PAM (k = 6) i porównanie z k-means
pam6 <- eclust(X_s, "pam", k = 6)
fviz_cluster( pam6, geom = "point", ggtheme = theme_minimal(base_size = 12), main = "PAM clustering (k = 6)" )
fviz_silhouette(pam6)
## cluster size ave.sil.width
## 1 1 25 0.45
## 2 2 20 0.32
## 3 3 45 0.49
## 4 4 38 0.39
## 5 5 39 0.50
## 6 6 33 0.33
mall_pam6 <- mall %>% mutate(cluster_pam6 = factor(pam6$cluster))
mall_pam6 %>% group_by(cluster_pam6) %>% summarise( Age = mean(Age), Income = mean(Annual.Income..k..), Spending = mean(Spending.Score..1.100.), n = n() )
## # A tibble: 6 × 5
## cluster_pam6 Age Income Spending n
## <fct> <dbl> <dbl> <dbl> <int>
## 1 1 25.5 25.5 75.4 25
## 2 2 46.2 26.8 18.4 20
## 3 3 56.3 54.3 49.1 45
## 4 4 26.7 57.6 47.8 38
## 5 5 32.7 86.5 82.1 39
## 6 6 41.9 88.9 17.0 33
Porównanie przypisań k-means vs PAM wykonuję przy pomocy Adjusted
Rand Index (ARI) – im bliżej 1, tym bardziej zgodne segmentacje.
comp_km_pam <- comPart(km6$cluster, pam6$cluster)
comp_km_pam
## ARI RI J FM
## 0.7428825 0.9220603 0.6535627 0.7922923
Wynik: Adjusted Rand Index wyniósł 0.74, a miary Jaccarda i
Fowlkesa–Mallowsa odpowiednio 0.65 i 0.79, co wskazuje, że obie metody
identyfikują bardzo podobne segmenty klientów.
Predykcja przynależności do klastrów (train/test)
Aby pokazać praktyczną użyteczność segmentacji, wykonuję podział na
zbiór treningowy i testowy. Następnie uczę k-means na danych
treningowych i przypisuję klaster dla nowych obserwacji.
n <- nrow(X_s)
idx_train <- sample.int(n, size = floor(0.8 * n))
X_train <- X_s[idx_train, , drop = FALSE]
X_test <- X_s[-idx_train, , drop = FALSE]
mall_train <- mall[idx_train, ]
mall_test <- mall[-idx_train, ]
k-means na zbiorze treningowym
km_train <- kmeans(X_train, centers = 6, nstart = 25)
predykcja: przypisanie nowych obserwacji do istniejących
klastrów
km_train_kcca <- as.kcca(km_train, X_train)
cluster_pred <- predict(km_train_kcca, X_test)
mall_test_pred <- mall_test %>% mutate(cluster_pred = factor(cluster_pred))
table(cluster_pred)
## cluster_pred
## 1 2 3 4 5 6
## 11 8 4 3 3 11
head(mall_test_pred)
## CustomerID Gender Age Annual.Income..k.. Spending.Score..1.100. cluster_pred
## 10 10 Female 30 19 72 5
## 15 15 Male 37 20 13 3
## 18 18 Male 20 21 66 5
## 19 19 Male 52 23 29 3
## 28 28 Male 35 28 61 5
## 33 33 Male 53 33 4 3
mall_test_pred %>% group_by(cluster_pred) %>% summarise( Age = mean(Age), Income = mean(Annual.Income..k..), Spending = mean(Spending.Score..1.100.), n = n() )
## # A tibble: 6 × 5
## cluster_pred Age Income Spending n
## <fct> <dbl> <dbl> <dbl> <int>
## 1 1 26.5 57.9 46.9 11
## 2 2 32.5 77.9 86.9 8
## 3 3 47.8 28.8 18.5 4
## 4 4 42.3 101 16 3
## 5 5 28.3 22.7 66.3 3
## 6 6 58.5 56.4 48.8 11
dodaje etykiety do train
mall_train$cluster_train <- factor(km_train$cluster)
ggplot() + geom_point( data = mall_train, aes(x = Annual.Income..k.., y = Spending.Score..1.100., color = cluster_train), alpha = 0.7 ) + geom_point( data = mall_test_pred, aes(x = Annual.Income..k.., y = Spending.Score..1.100.), shape = 1, size = 2 ) + theme_minimal(base_size = 12) + labs( title = "Predykcja przynależności nowych klientów do klastrów", subtitle = "Wypełnione: dane treningowe, puste kółka: dane testowe", color = "Cluster" )
Klasteryzacja hierarchiczna i porównanie z k-means
Klasteryzacja hierarchiczna umożliwia analizę struktury danych w
postaci dendrogramu. Porównuję dwa standardowe warianty łączenia:
complete linkage, Ward.D2 (minimalizacja wariancji
wewnątrzklastrowej)
dm <- dist(X_s, method = "euclidean")
hc_complete <- hclust(dm, method = "complete")
hc_ward <- hclust(dm, method = "ward.D2")
plot(hc_complete, main = "Klasteryzacja hierarchiczna: kompletne połączenie", cex = 0.7)
plot(hc_ward, main = "Klasteryzacja hierarchiczna: Ward.D2", cex = 0.7)
plot(hc_ward, labels = FALSE, main = "Klasteryzacja hierarchiczna: Ward.D2 (ukryte etykiety)")
Cięcie drzewa do k = 6 i silhouette
k <- 6
cl_hc_complete <- cutree(hc_complete, k = k)
cl_hc_ward <- cutree(hc_ward, k = k)
sil_complete <- silhouette(cl_hc_complete, dm)
sil_ward <- silhouette(cl_hc_ward, dm)
mean_sil_complete <- mean(sil_complete[, 3])
mean_sil_ward <- mean(sil_ward[, 3])
mean_sil_complete
## [1] 0.3745611
mean_sil_ward
## [1] 0.420117
Wynik: Ward.D2 (0.420117) ma wyższą silhouette niż complete
(0.3745611) oznacza to lepszą jakość segmentacji hierarchicznej
Wizualizacja klastrów (Ward.D2)
fviz_cluster(list(data = X_s, cluster = cl_hc_ward), ggtheme = theme_minimal(base_size = 12), main = "HC Ward.D2 clusters (k = 6)")
Porównanie HC vs k-means (ARI)
ARI_hc_ward_vs_km6 <- mclust::adjustedRandIndex(cl_hc_ward, km6$cluster)
ARI_hc_complete_vs_km6 <- mclust::adjustedRandIndex(cl_hc_complete, km6$cluster)
ARI_hc_ward_vs_km6
## [1] 0.69221
ARI_hc_complete_vs_km6
## [1] 0.8080992
Wynik: ARI (Ward.D2 vs k-means) = 0.69 oznacza to wysoką zgodność
pomiędzy klasteryzacją hierarchiczną metodą Ward.D2, a rozwiązaniem
k-means dla k = 6. ARI (complete vs k-means) = 0.81 co wskazuje na
bardzo wysoką zgodność struktur klastrów uzyskanych metodą complete
linkage oraz k-means. Metoda complete w tym przypadku lepiej odtwarza
podział uzyskany przez k-means, mimo że wcześniej charakteryzowała się
niższą wartością silhouette.
WNIOSKI
1 Dane wykazują silną tendencję klastrową (Hopkins 0.967832)
2 Dobór liczby klastrów dał konflikt kryteriów: silhouette sugerował
k=8, gap k=6. Przyjęto k=6 jako rozwiązanie główne i k=8 jako
porównawcze.
3 Segmentacja k-means (k=6) jest spójna z PAM (wysoki ARI).
4 W klasteryzacji hierarchicznej Ward.D2 osiąga lepszą silhouette,
natomiast complete ma wyższą zgodność z k-means (ARI)
5 Pokazano możliwość przypisywania nowych klientów do istniejących
segmentów (podejście train/test).