W pierwszym wpisie na temat grupowania hierarchicznego skupiłem się na omówieniu podstawowych pojęć i teoretycznym wprowadzeniu. Dziś kolej na część praktyczną.
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.¶
# 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.
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.
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
df.shape
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.¶
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.
df.head()
df.shape
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.¶
pca = PCA(n_components=2)
df_pca = pca.fit_transform(df)
df = pd.DataFrame(df_pca, columns=['c1', 'c2'], index=df.index)
df.head()
5. Przygotowanie zbioru do grupowania.¶
Przed rozpoczęciem procesu grupowania wskazane jest wykonanie kilku operacji na zbiorze:
- Usunięcie ujemnych wartości (przed transformacją logarytmiczną).
- Usunięcie skośności zmiennych (transformacja logarytmiczna).
- 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.
df.agg(['mean', 'median', 'std', 'min', 'max']).round(2)
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.
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. 🙂
for col in df:
if df[col].min() <= 0:
df[col] = df[col] + np.abs(df[col].min()) + 1
Usuwam skośność zmiennych.
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.
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})
outliers
W pętli usuwam wszystkie obserwacje, które nie spełniają dwóch warunków.
for row in outliers.iterrows():
df = df[(df[row[0]] >= row[1]['lower_boundary']) & (df[row[0]] <= row[1]['upper_boundary'])]
Wykonuję standaryzację zmiennych.
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.
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)
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.
df = df.sample(1000, random_state=20190314)
Sprawdźmy jeszcze, jak prezentują się wylosowane obserwacje na płaszczyźnie dwuwymiarowej.
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.
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ę.
df['label_1'] = model_1.labels_
Mając tak przygotowany zbiór, wizualizuję go w przestrzeni dwuwymiarowej. 🙂
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:
- Podziału na grupy można dokonać, korzystając z innego kryterium niż maksymalna liczba grup.
- 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.
model_2 = linkage(df.iloc[:,0:2], method = 'ward', metric = 'euclidean')
Buduję dendrogram, by zobrazować działanie algorytmu.
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.
clusters = fcluster(model_2, 3, criterion='maxclust')
df['label_2'] = clusters
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 :-)
Cześć Mateusz
Z dużym zainteresowaniem śledzę Twoją stronę. Świetne artykuły. Powiedz mi proszę czym było podyktowane to iż nie zastosowałeś standaryzacji danych (np. klasyczny sklearn.cluster.StandardScaler()) przed redukcją ilości zmiennych niezależnych (PCA). Brak takiej standaryzacji faworyzuje dane w PCA które maja wartości dużo większe i większy wpływ na wynik analizy PCA. Dobrego dnia 🙂