Grupowanie hierarchiczne – praktyka

analiza skupień, grupowanie, klasteryzacja, klastry, klastrowanie, grupowanie

W pierwszym wpisie na temat grupowania hierarchicznego skupiłem się na omówieniu podstawowych pojęć i teoretycznym wprowadzeniu. Dziś kolej na część praktyczną.

Z tego artykułu dowiesz się: 1. Jak wykonać grupowanie hierarchiczne? 2. Jakie biblioteki wspierają grupowanie hierarchiczne? 3. O czym należy pamiętać, przeprowadzając proces grupowania?

Mam nadzieję, że część teoretyczna była dla Ciebie w miarę przystępna i jasna. Jeśli udało Ci się przez nią przebrnąć i interesują Cię algorytmy hierarchiczne, to w tym wpisie będzie już tylko ciekawiej, bo skupia się ona tylko i wyłącznie na praktyce. 😉 Nie przedłużając, wczytajmy zatem biblioteki, zbiór i zaczynajmy! 🙂

1. Wczytanie podstawowych bibliotek.

In [1]:
# biblioteki do manipulacji zbiorem
import pandas as pd
import numpy as np
# wizualizacja danych
import seaborn as sns
from matplotlib import pyplot as plt
# redukcja wymiarów zbioru - na potrzeby wizualizacji grup
from sklearn.decomposition import PCA
# skalowanie zmiennych - wskazane przy metodach korzystających z miar odległości
from sklearn.preprocessing import StandardScaler
# biblioteki do grupowania hierarchicznego: sklearn i scipy
from sklearn.cluster import AgglomerativeClustering
from scipy.cluster.hierarchy import dendrogram, linkage, fcluster

Będziemy operować na dosyć "szerokim" zbiorze, dlatego zmieniam liczbę wyświetlanych kolumn przez bibliotekę Pandas na 30.

In [2]:
pd.options.display.max_columns = 30

2. Wczytanie zbioru.

Na potrzeby tego wpisu posłużę się zbiorem, który już pojawiał się na blogu przy okazji projektu dotyczącego przewidywanie defaultu wśród posiadaczy kart kredytowych.

In [3]:
ef = pd.ExcelFile('data/data.xls')
df = ef.parse('Data', skiprows=1, names = ['id', 'lim_kredytu', 'plec', 'wyksztalcenie', 'stan_cywilny', 'wiek', 'opozn_plat_wrz', 'opozn_plat_sie', 'opozn_plat_lip', 'opozn_plat_cze', 'opozn_plat_maj', 'opozn_plat_kwi', 'kwota_wyciagu_wrz', 'kwota_wyciagu_sie', 'kwota_wyciagu_lip', 'kwota_wyciagu_cze', 'kwota_wyciagu_maj', 'kwota_wyciagu_kwi', 'platnosc_wrz', 'platnosc_sie', 'platnosc_lip', 'platnosc_cze', 'platnosc_maj', 'platnosc_kwi', 'y'])
df.drop('id', axis = 1, inplace = True) # usuwam zbędną kolumnę - id
df.drop('y', axis = 1, inplace = True) # usuwam zbędną kolumnę - y
In [4]:
df.shape
Out[4]:
(30000, 23)

Na wstępie mam 30 tysięcy obserwacji i 23 zmienne. Szczególnie ta pierwsza wartość jest dosyć wysoka jak na możliwości mojej maszyny (w pierwszym wpisie wspominałem o wysokiej złożoności obliczeniowej algorytmów hierarchicznych). W kolejnych krokach pozbędę się części obserwacji i wykonam kilka podstawowych transformacji.

3. Wstępne przygotowanie zbioru.

In [5]:
df.plec.replace([1,2], [0, 1], inplace = True) # zamieniam typ zmiennej na binarny
df.wyksztalcenie.replace([0, 1, 2, 3, 4, 5, 6], ['nieznane', 'wyzsze_pelne', 'wyzsze', 'srednie', 'inne', 'nieznane', 'nieznane'], inplace = True)
df.stan_cywilny.replace([0, 1, 2, 3], ['nieznany', 'w_zwiazku', 'kawaler_panna', 'inny'], inplace = True)
df = pd.get_dummies(df) # wykonuję transformację zmiennych kategorycznych do kodowania dummy

Podejrzyjmy, jak wygląda zbiór po wstępnych transformacjach.

In [6]:
df.head()
Out[6]:
lim_kredytu plec wiek opozn_plat_wrz opozn_plat_sie opozn_plat_lip opozn_plat_cze opozn_plat_maj opozn_plat_kwi kwota_wyciagu_wrz kwota_wyciagu_sie kwota_wyciagu_lip kwota_wyciagu_cze kwota_wyciagu_maj kwota_wyciagu_kwi platnosc_wrz platnosc_sie platnosc_lip platnosc_cze platnosc_maj platnosc_kwi wyksztalcenie_inne wyksztalcenie_nieznane wyksztalcenie_srednie wyksztalcenie_wyzsze wyksztalcenie_wyzsze_pelne stan_cywilny_inny stan_cywilny_kawaler_panna stan_cywilny_nieznany stan_cywilny_w_zwiazku
0 20000 1 24 2 2 -1 -1 -2 -2 3913 3102 689 0 0 0 0 689 0 0 0 0 0 0 0 1 0 0 0 0 1
1 120000 1 26 -1 2 0 0 0 2 2682 1725 2682 3272 3455 3261 0 1000 1000 1000 0 2000 0 0 0 1 0 0 1 0 0
2 90000 1 34 0 0 0 0 0 0 29239 14027 13559 14331 14948 15549 1518 1500 1000 1000 1000 5000 0 0 0 1 0 0 1 0 0
3 50000 1 37 0 0 0 0 0 0 46990 48233 49291 28314 28959 29547 2000 2019 1200 1100 1069 1000 0 0 0 1 0 0 0 0 1
4 50000 0 57 -1 0 -1 0 0 0 8617 5670 35835 20940 19146 19131 2000 36681 10000 9000 689 679 0 0 0 1 0 0 0 0 1
In [7]:
df.shape
Out[7]:
(30000, 30)

Po transformacjach liczba zmiennych urosła do 30. By zobrazować proces grupowania, będę wizualizować obserwacje, wraz z informacją o grupie, do której zostały przypisane. Mając na uwadze liczbę zmiennych, musiałbym użyć trzydziestowymiarowej przestrzeni, co może być dosyć problematyczne. 😉 Na potrzeby tego wpisu zmniejszę liczbę wymiarów do dwóch. Użyję do tego odpowiedniego algorytmu - PCA. 🙂

4. Zmniejszenie liczby wymiarów.

In [8]:
pca = PCA(n_components=2)
df_pca = pca.fit_transform(df)
df = pd.DataFrame(df_pca, columns=['c1', 'c2'], index=df.index)
In [9]:
df.head()
Out[9]:
c1 c2
0 -166488.191086 -75538.153664
1 -114226.976817 9780.671250
2 -98432.362181 -33471.387331
3 -71230.675903 -95224.273104
4 -114834.618365 -68729.185234

5. Przygotowanie zbioru do grupowania.

Przed rozpoczęciem procesu grupowania wskazane jest wykonanie kilku operacji na zbiorze:

  1. Usunięcie ujemnych wartości (przed transformacją logarytmiczną).
  2. Usunięcie skośności zmiennych (transformacja logarytmiczna).
  3. Centrowanie i skalowanie zmiennych.

Kolejność powyższych operacji ma oczywiście znaczenie. 🙂 Wykonanie powyższych transformacji (szczególnie skalowanie zmiennych - usuwanie skośności jest istotniejsze w przypadku np. algorytmu k-średnich, który będzie omawiany z użyciem tego samego zbioru, w jednym z kolejnych wpisów) jest istotne ze względu na sposób, w jaki działają algorytmy hierarchiczne - mierzą odległości pomiędzy obserwacjami. Odległości te mogą nieco przekłamane, jeśli wartości na poszczególnych wymiarach są znacząco różne (np. jedna zmienna została wyrażona w złotówkach, a druga w milionach złotych).

Sprawdzam podstawowe statystyki zbioru.

In [10]:
df.agg(['mean', 'median', 'std', 'min', 'max']).round(2)
Out[10]:
c1 c2
mean 0.00 -0.00
median -52210.70 -44280.93
std 166517.94 115827.56
min -179141.27 -457588.61
max 2196733.98 599288.02

Po rzucie oka na podstawowe statystyki zbioru mam kilka zastrzeżeń:

  • zbiór nie powinien zawierać ujemnych wartości (ze względu na transformacje mające na celu usunięcie skośności),
  • obie zmienne mają prawoskośny rozkład,
  • odchylenie standardowe powinno wynosić 0.
In [11]:
sns.distplot(df.c1, color = '#eb6c6a').set(title = 'Wykres gęstości - zmienna c1')
plt.show()
sns.distplot(df.c2, color = '#eb6c6a').set(title = 'Wykres gęstości - zmienna c2')
plt.show()

Usuwam wartości ujemne.

By pozbyć się wszystkich wartości ujemnych, zwiększam wartość poszczególnych zmiennych o wartość bezwzględną z najmniejszej wartości występującą w danej zmiennej. By uniknąć zer, dodaję jeszcze 1. 🙂

In [12]:
for col in df:
    if df[col].min() <= 0:
        df[col] = df[col] + np.abs(df[col].min()) + 1

Usuwam skośność zmiennych.

In [13]:
df = np.log(df)

Usuwam wartości odstające.

Rozkład powinien już przypominać rozkład normalny, dlatego zastosuję metodę 1.5 odstępu międzykwartylowego.

In [14]:
q1 = df.quantile(0.25)
q3 = df.quantile(0.75)
iqr = q3 - q1

low_boundary = (q1 - 1.5 * iqr)
upp_boundary = (q3 + 1.5 * iqr)
num_of_outliers_L = (df[iqr.index] < low_boundary).sum()
num_of_outliers_U = (df[iqr.index] > upp_boundary).sum()
outliers = pd.DataFrame({'lower_boundary':low_boundary, 'upper_boundary':upp_boundary,'num_of_outliers__lower_boundary':num_of_outliers_L, 'num_of_outliers__upper_boundary':num_of_outliers_U})
In [15]:
outliers
Out[15]:
lower_boundary num_of_outliers__lower_boundary num_of_outliers__upper_boundary upper_boundary
c1 9.544414 355 47 13.998835
c2 12.266477 34 18 13.731749

W pętli usuwam wszystkie obserwacje, które nie spełniają dwóch warunków.

In [16]:
for row in outliers.iterrows():
    df = df[(df[row[0]] >= row[1]['lower_boundary']) & (df[row[0]] <= row[1]['upper_boundary'])]

Wykonuję standaryzację zmiennych.

In [17]:
scaler = StandardScaler()
scaler.fit(df)
df_std = scaler.transform(df)
df = pd.DataFrame(data=df_std, index=df.index, columns=df.columns)

Po wykonaniu kilku transformacji ponownie sprawdzam rozkłady zmiennych i podstawowe statystyki zbioru.

In [18]:
sns.distplot(df.c1, color = '#eb6c6a').set(title = 'Wykres gęstości - zmienna c1')
plt.show()
sns.distplot(df.c2, color = '#eb6c6a').set(title = 'Wykres gęstości - zmienna c2')
plt.show()
df.agg(['mean', 'std', 'max', 'min']).round(2)
Out[18]:
c1 c2
mean 0.00 -0.00
std 1.00 1.00
max 2.70 3.07
min -2.68 -3.17

Jest zdecydowanie lepiej. 🙂 Przy okazji ograniczę nieco liczebność zbioru. Jak wcześniej wspomniałem, jedną z wad grupowania hierarchicznego jest dosyć wysoka złożoność obliczeniowa. Ponadto, zbyt dużą liczba obserwacji nieco utrudni obserwację klastrów na wykresie. Pozostawię zatem 1000 obserwacji.

In [19]:
df = df.sample(1000, random_state=20190314)

Sprawdźmy jeszcze, jak prezentują się wylosowane obserwacje na płaszczyźnie dwuwymiarowej.

In [20]:
sns.lmplot('c1', 'c2', data = df, scatter_kws={"color": "#eb6c6a"}, fit_reg=False).set(title = 'Wykres punktowy zbioru')
plt.show()

Wszystko wygląda dobrze. Przechodzę zatem do wykonania grupowania.

6. Grupowanie obserwacji.

6.1. Grupowanie z użyciem biblioteki sklearn.

Wykonam grupowanie, w którym miarą odległości będzie odległość euklidesowa, natomiast jako metodę łączenia grup obiorę kryterium Warda. Sklearn daje możliwość wyboru liczby grup. Ustalam zatem ten parametr na 3.

In [21]:
model_1 = AgglomerativeClustering(linkage='ward', affinity='euclidean', n_clusters=3).fit(df)

Model zbudowany. By łatwiej zobrazować grupy na wykresie punktowym, dodam nową zmienną do zbioru. Będzie ona oznaczać grupę, do jakiej przypisano daną obserwację.

In [22]:
df['label_1'] = model_1.labels_

Mając tak przygotowany zbiór, wizualizuję go w przestrzeni dwuwymiarowej. 🙂

In [23]:
sns.lmplot('c1', 'c2', data = df, hue = 'label_1', fit_reg=False, palette = ['#eb6c6a', '#6aeb6c', '#6c6aeb']).set(title='Wizualizacja grup')
plt.show()
6.2. Grupowanie z użyciem biblioteki SciPy.

Proces grupowania można również przeprowadzić, używając biblioteki SciPy. Ma ona względem sklearn kilka plusów:

  1. Podziału na grupy można dokonać, korzystając z innego kryterium niż maksymalna liczba grup.
  2. Umożliwia ona wizualizację drzewa budowanego w procesie grupowania. Dzięki temu możemy zobaczyć, jak przebiegało grupowanie (możliwa jest ocena jakości grupowania po kształcie drzewa) oraz zdecydować, po jakim dystansie dzielącym obserwację od grupy będziemy budować nową grupę, zamiast łączyć je w jedną.

Buduję zatem drugi model z użyciem SciPy.

In [24]:
model_2 = linkage(df.iloc[:,0:2], method = 'ward', metric = 'euclidean')

Buduję dendrogram, by zobrazować działanie algorytmu.

In [25]:
fig = plt.figure(figsize=(25, 10))
dn = dendrogram(model_2)
plt.show()

Dzielę obserwacje na grupy według wybranego kryterium (maksymalnej liczby grup). Podobnie jak w przypadku poprzedniego modelu, przypisuję nazwy grup do obserwacji.

In [26]:
clusters = fcluster(model_2, 3, criterion='maxclust')
df['label_2'] = clusters
In [27]:
sns.lmplot('c1', 'c2', data = df, fit_reg=False, hue = 'label_2', palette = ['#eb6c6a', '#6aeb6c', '#6c6aeb'])
plt.show()

Podsumowanie

Udało mi się otrzymać identyczne wyniki w przypadku obu modeli, co może potwierdzać, że w obu bibliotekach algorytm został zaimplementowany poprawnie. 😉 Mam nadzieję, że dzięki temu wpisowi nieco przybliżyłem Ci ideę grupowania hierarchicznego. 🙂

Jeśli masz jakieś pytania, to proszę podziel się nimi w komentarzu pod wpisem - zapraszam do dyskusji. Jeśli artykuł przypadł Ci do gustu, to proszę, podziel się nim w mediach społecznościowych ze swoimi znajomymi. Będę bardzo wdzięczny. 🙂

Podobał Ci się ten artykuł?

Jeśli tak, to zarejestruj się, by otrzymywać informacje o nowych wpisach. Dodatkowo w prezencie wyślę Ci bezpłatny poradnik 🙂

Bądź pierwszy, który skomentuje ten wpis!

Dodaj komentarz

Twój adres email nie zostanie opublikowany.


*