“Some people use statistics as a drunk man uses lamp posts–for support rather than for illumination.”
– Prawdopodobnie Andrew Lang inspirowany tekstem napisanym w 1903 roku przez A. E. Housmana
Kolejny piątek, więc i kolejny zestaw przykładów. Dzisiaj weźmiemy na tapet (absolutnie nie na tapetę) parowe testy statystyczne.
Eksperyment
Jako, że poprzedni wpis koncentrował się na przeprowadzaniu eksperymentów z wykorzystaniem powtarzanej walidacji krzyżowej, dzisiaj nie poświęcimy temu etapowi wiele uwagi. Wyróżnimy tylko te drobne elementy, które mogą być dla nas nowe.
Przygotowanie danych
Jak zwykle zaczynamy od danych, ale tym razem nie będę już one syntetyczne. Ostatnio bliskim mi krajem stała się chwilowo Australia, skorzystamy więc ze zbioru danych australian dostępnego na repozytorium KEEL. Zestaw ten można pobrać także w postaci pliku CSV z repozytorium benchmark_datasets autorstwa dr. (tak, z kropką) Ksieniewicza oraz mojego.
Wczytajmy nasze dane:
import numpy as np
dataset = 'australian'
dataset = np.genfromtxt("%s.csv" % (dataset), delimiter=",")
X = dataset[:, :-1]
y = dataset[:, -1].astype(int)
Zdefiniowanie i inicjalizacja modeli klasyfikacji
Kiedy przeprowadzamy eksperymenty, to najczęściej chcemy coś ze sobą porównać. Z tego powodu wykorzystamy dzisiaj nie jeden, ale trzy modele klasyfikacji. Będą to (i) Naiwny klasyfikator bayesowski (GNB) (ii) klasyfikator k najbliższych sąsiadów (kNN) oraz (iii) drzewo klasyfikacyjne i regresyjne (CART). Modele inicjalizujemy z domyślnymi parametrami, pamiętając o ustawieniu random_state dla algorytmu CART, i zawieramy je w słowniku clfs:
from sklearn.naive_bayes import GaussianNB
from sklearn.neighbors import KNeighborsClassifier
from sklearn.tree import DecisionTreeClassifier
clfs = {
'GNB': GaussianNB(),
'kNN': KNeighborsClassifier(),
'CART': DecisionTreeClassifier(random_state=42),
}
Przeprowadzenie eksperymentu & zapisanie wyników
Teraz rutynowa i nudna już dla nas procedura wielokrotnej walidacji krzyżowej. Zwróćmy uwagę za zdefiniowanie scores nie jako listy, ale jako dwuwymiarowej macierzy numpy, w której liczba wierszy odpowiada liczbie testowanych przez nas modeli, a liczba kolumn równa jest liczbie wyników uzyskiwanych w procesie walidacji (tutaj wynoszącej 10).
Dodatkowo warte zauważenia jest wykorzystanie funkcji enumerate(), która pozwala nam na wprowadzenie inkrementowanego licznika bezpośrednio do pętli, a także użycie znanej z wykładu funkcji clone() do klonowania klasyfikatorów. Funkcja ta w tym wypadku nie jest konieczna, ale korzystanie z niej na pewno będzie dobrym nawykiem, który zaprocentuje w przyszłości.
n_splits = 5
n_repeats = 2
rskf = RepeatedStratifiedKFold(
n_splits=n_splits, n_repeats=n_repeats, random_state=42)
scores = np.zeros((len(clfs), n_splits * n_repeats))
for fold_id, (train, test) in enumerate(rskf.split(X, y)):
for clf_id, clf_name in enumerate(clfs):
clf = clone(clfs[clf_name])
clf.fit(X[train], y[train])
y_pred = clf.predict(X[test])
scores[clf_id, fold_id] = accuracy_score(y[test], y_pred)
Teraz przedstawmy uzyskane wyniki tak, jak było to robione w poprzednim wpisie:
mean = np.mean(scores, axis=1)
std = np.std(scores, axis=1)
for clf_id, clf_name in enumerate(clfs):
print("%s: %.3f (%.2f)" % (clf_name, mean[clf_id], std[clf_id]))
>> GNB: 0.789 (0.02)
>> kNN: 0.708 (0.03)
>> CART: 0.810 (0.03)
Jak widzimy, najwyższą wartości jakości klasyfikacji uzyskał algorytm CART, o 2% gorszy okazał się GNB, a kNN to już zupełnie najgorzej.
Jednak czy powyższe wyniki pozwalają nam na bezsprzeczne stwierdzenie, że CART jest w tym przypadku najlepszym modelem klasyfikacji? Absolutnie nie. Dlatego właśnie korzystamy z testów statystycznych.
Zanim jednak przejdziemy do właściwego tematu dzisiejszego dnia, zastosujmy starą harcerską sztuczkę z wykorzystaniem funkcji save() biblioteki numpy. Dzięki temu możemy zapisać uzyskane przez nas wyniki z poszczególnych foldów do pliku NPY i wykorzystać go w osobnym skrypcie do analizy statystycznej. Warto zapamiętać tą funkcję, ponieważ okazuje się ona bardzo przydatna podczas prowadzenia badań. Jeżeli nasze obliczenia trwają długo, to nie chcemy przecież uruchamiać ich za każdym razem od nowa, aby mieć wygodny dostęp do wyników.
np.save('results', scores)
Analiza statystyczna
W końcu nadszedł czas na clou dzisiejszego wpisu, czyli parowe testy statystyczne. Mają one na celu porównanie między sobą modeli każdy z każdym i ewentualne wskazanie w każdej z par, czy jeden z modeli jest statystycznie znacząco lepszy od drugiego.
Wczytanie wyników eksperymentu
Zaczniemy od wczytania wcześniej zapisanych wyników poszczególnych foldów walidacji krzyżowej:
import numpy as np
scores = np.load('results.npy')
print("Folds:\n", scores)
>> Folds:
>> [[0.77697842 0.76258993 0.79710145 0.80291971 0.81751825 0.75539568
>> 0.81294964 0.76086957 0.81751825 0.78832117]
>> [0.67625899 0.69064748 0.70289855 0.76642336 0.69343066 0.68345324
>> 0.67625899 0.70289855 0.72992701 0.75912409]
>> [0.79856115 0.83453237 0.81884058 0.83941606 0.75182482 0.76978417
>> 0.79856115 0.8115942 0.86131387 0.81751825]]
Wszystko jest na swoim miejscu, nic nie zginęło, możemy więc brać się do pracy.
t-statystyka oraz p-wartość
Jednymi z częściej wykorzystywanych statystycznych testów parowych, jeżeli chodzi o badania związane z uczeniem maszynowym, są test T Studenta oraz test Wilcoxona (ang. Wilcoxon signed-rank test). My skorzystamy ze znanego z wykładu testu T studenta dla prób zależnych z poziomem ufności alpha równym 0.05. Test ten implementuje funkcja ttest_rel().
Aby jak najlepiej zrozumieć to, co się za chwilę wydarzy, rozbijemy naszą analizę statystyczną na 5 kroków. Każdy z tych kroków wiąże się z utworzeniem tabeli (macierzy NxN, gdzie N to liczba porównywanych metod) zawierającej odpowiednie informacje. Wszystkie te tabele będziemy czytać wierszami.
Pierwsze dwa kroki to przedstawienie w tabelach znanych nam z wykładu wartości t-statystyki (tutaj wskazującej na różnice w jakości klasyfikacji pomiędzy metodami) oraz p-wartości (przyrównywane do wartości alpha w celu przyjęcia lub odrzucenia hipotezy zerowej) dla każdej pary klasyfikatorów. W tym celu utwórzmy dwie macierze dwuwymiarowe t_statistic oraz p_value, a następnie uzupełnijmy je wykorzystując funkcję ttest_rel(), do której będziemy podawać kolejno jakości klasyfikacji uzyskane przez każdą z par klasyfikatorów w każdym z foldów. Następnie wyświetlmy wyniki:
from scipy.stats import ttest_rel
alfa = .05
t_statistic = np.zeros((len(clfs), len(clfs)))
p_value = np.zeros((len(clfs), len(clfs)))
for i in range(len(clfs)):
for j in range(len(clfs)):
t_statistic[i, j], p_value[i, j] = ttest_rel(scores[i], scores[j])
print("t-statistic:\n", t_statistic, "\n\np-value:\n", p_value)
>> t-statistic:
>> [[ 0. 6.29911296 -1.64022221]
>> [-6.29911296 0. -6.98938694]
>> [ 1.64022221 6.98938694 0. ]]
>>
>> p-value:
>> [[1.00000000e+00 6.14348994e-06 1.18317253e-01]
>> [6.14348994e-06 1.00000000e+00 1.58542999e-06]
>> [1.18317253e-01 1.58542999e-06 1.00000000e+00]]
Wartości wydają się być w porządku, na przekątnej macierzy t-statystyki mamy zera (te same wartości jakości klasyfikacji, bo porównujemy klasyfikator sam ze sobą), a na przekątnej p-wartości mamy jedynki (z tego samego powodu).
Trudno jednak ukryć, że powyższa reprezentacja pozostawia wiele do życzenia pod względem wizualnym. Spróbujmy coś na to poradzić, korzystając z biblioteki tabulate. Biblioteka ta pozwala nam na proste sformatowanie macierzy numpy do postaci w miarę czytelnej tabeli, a osiągamy to podając do funkcji tabulate() jako argumenty naszą macierz oraz wcześniej zdefiniowane nagłówki kolumn będące akronimami wykorzystanych metod klasyfikacji (argument headers). Dodatkowo, za pomocą argumentu floatfmt możemy ustawić format wyświetlania wartości zmiennoprzecinkowych. Dla polepszenia czytelności, przed wykorzystaniem funkcji tabulate(), dodajmy z przodu kolumnę zawierającą akronimy klasyfikatorów. Na koniec zobaczmy efekty:
from tabulate import tabulate
headers = ["GNB", "kNN", "CART"]
names_column = np.array([["GNB"], ["kNN"], ["CART"]])
t_statistic_table = np.concatenate((names_column, t_statistic), axis=1)
t_statistic_table = tabulate(t_statistic_table, headers, floatfmt=".2f")
p_value_table = np.concatenate((names_column, p_value), axis=1)
p_value_table = tabulate(p_value_table, headers, floatfmt=".2f")
print("t-statistic:\n", t_statistic_table, "\n\np-value:\n", p_value_table)
>> t-statistic:
>> GNB kNN CART
>> ---- ----- ----- ------
>> GNB 0.00 6.30 -1.64
>> kNN -6.30 0.00 -6.99
>> CART 1.64 6.99 0.00
>>
>> p-value:
>> GNB kNN CART
>> ---- ----- ----- ------
>> GNB 1.00 0.00 0.12
>> kNN 0.00 1.00 0.00
>> CART 0.12 0.00 1.00
O wiele lepiej, człowiek w końcu jest w stanie na to spojrzeć i coś zrozumieć.
Przewaga
Trzecim krokiem będzie wyznaczenie tabeli przewagi, ukazującej, który z klasyfikatorów z pary osiągnął lepszą jakość klasyfikacji. W tym celu tworzymy wypełnioną zerami macierz advantage i wpisujemy do niej wartości 1 tam, gdzie w macierzy t_statistic występują wartości większe od 0:
advantage = np.zeros((len(clfs), len(clfs)))
advantage[t_statistic > 0] = 1
advantage_table = tabulate(np.concatenate(
(names_column, advantage), axis=1), headers)
print("Advantage:\n", advantage_table)
>> Advantage:
>> GNB kNN CART
>> ---- ----- ----- ------
>> GNB 0 1 0
>> kNN 0 0 0
>> CART 1 1 0
Dzięki temu wyraźnie widzimy, że GNB jest lepszy od kNN i gorszy od CART, kNN jest najgorszy, a CART osiągnął wartość jakości klasyfikacji większą niż GNB oraz kNN.
Różnice statystycznie znaczące
Przedostatnim krokiem jest utworzenie macierzy istotności ukazującej, czy różnica pomiędzy daną parą klasyfikatorów jest istotna statystycznie. Dokonujemy tego poprzez przyrównanie p-wartości w tabeli p-value do naszej wartości alpha i na tej podstawie przyjęcie bądź odrzucenie hipotezy zerowej (czyli, że pomiędzy klasyfikatorami nie ma istotnej różnicy statystycznej). Jeżeli p-wartość jest mniejsza lub równa alpha to odrzucamy hipotezę zerową, a jeżeli p-wartość jest większa od alpha, to przyjmujemy tą hipotezę.
significance = np.zeros((len(clfs), len(clfs)))
significance[p_value <= alfa] = 1
significance_table = tabulate(np.concatenate(
(names_column, significance), axis=1), headers)
print("Statistical significance (alpha = 0.05):\n", significance_table)
>> Statistical significance (alpha = 0.05):
>> GNB kNN CART
>> ---- ----- ----- ------
>> GNB 0 1 0
>> kNN 1 0 1
>> CART 0 1 0
Na podstawie tej tabeli możemy stwierdzić, że GNB i kNN oraz CART i kNN są statystycznie od siebie różne.
Wynik końcowy analizy statystycznej
Ostatnim krokiem naszej analizy będzie przedstawienie tabeli obserwacji końcowych, która ukaże nam, dla każdego z klasyfikatorów, algorytmy od których jest on statystycznie znacząco lepszy.
W tym celu wykorzystujemy dwie poprzednio utworzone tabele:
stat_better = significance * advantage
stat_better_table = tabulate(np.concatenate(
(names_column, stat_better), axis=1), headers)
print("Statistically significantly better:\n", stat_better_table)
>> Statistically significantly better:
>> GNB kNN CART
>> ---- ----- ----- ------
>> GNB 0 1 0
>> kNN 0 0 0
>> CART 0 1 0
Widzimy, że zarówno GNB jak CART są statystycznie znacząco lepsze od kNN. Nie mamy natomiast, pomimo różnicy w średniej jakości klasyfikacji, podstaw do stwierdzenia, że CART spisał się lepiej niż GNB. W ewentualnych wnioskach należałoby więc uznać różnice pomiędzy tymi dwiema metodami za nieistotne statystycznie.
Dobrnęliśmy do końca. W przyszłą środę przeprowadzimy badania na wielu zestawach danych i wykonamy globalny rankingowy test statystyczny.
