Na dzisiejszych laboratoriach poznamy nowe podejście do regresji i klasyfikacji, czyli metody drzewiaste takie jak drzewa regresji i klasyfikacji i lasy losowe.
Będziemy korzystać z metod zaimplementowanych w pakietach randomForest i rpart.

Drzewa CART

W pierwszym zadaniu zbadamy, jak drzewo klasyfikacji radzi sobie z przewidywaniem Drzewa klasyfikacji i regresji (ang. CART) zostały omówione na Wykładzie 11. Proste wprowadzenie do tej tematyki można również obejrzeć tutaj.

Zadanie 1. Zadanie przykładowe. Załaduj dane biopsy z pakietu MASS. Stwórz drzewo klasyfikacji które pomoże ocenić złośliwość nowotworu. Zbadaj błąd testowy za pomocą k-krotnej walidacji krzyżowej. Przytnij drzewo i sprawdź, czy poprawiło to klasyfikator. Wykorzystaj funkcje rpart, plotcp oraz prune z pakietu rpart.

Rozwiązanie. Dane biopsy oglądaliśmy już na laboratorium 10, więc odpuścimy sobie ten etap. Ładujemy pakiet rpart (oczywiście należy go wcześniej zainstalować) i budujemy drzewo. Oczywiście musimy w tym celu usunąć z modelu zmienną ID, ponieważ budowanie klasyfikatora wykorzystującego ID pacjenta nie ma żadnego sensu! Usuniemy również brakujące obserwacje (choć akurat drzewa decyzyjne dość dobrze sobie z nimi radzą, o czym można przeczytać m.in. w manualu do pakietu rpart).
Ponieważ kolumna class ma typ factor, to funkcja rpart domyśli się, że należy stworzyć drzewo klasyfikacji a nie regresji. Mimo to wyspecyfikujemy to za pomocą argumentu method. Domyślnie do stworzenia drzewa wykorzystany jest indeks Giniego. Liście drzewa są kolejno dzielone dopóki zapewnia to odpowiednio wysoką poprawę błędu klasyfikacji i dopóki liczba obserwacji odpowiadających liściowi nie spadnie poniżej określonego progu. Próg można zmodyfikować przekazując funkcji rpart odpowiednie argumenty w polu control, co opisuje dokładniej jej dokumentacja.

library(MASS)
library(rpart)
biopsy$ID <- NULL  
biopsy <- na.omit(biopsy)
drzewo <- rpart(class~., data=biopsy, method='class')
plot(drzewo, main="Klasyfikacja złośliwości nowotworu")
text(drzewo, use.n=TRUE, all=T, cex=.8)

Widzimy, że w pierwszym kroku nowotwory są rozdzielane za pomocą zmiennej V2, czyli jednorodności rozmiaru komórek, co jest zgodne z wiedzą dotyczącą histopatologii nowotworów. Próbki tkanki, dla których ta zmienna jest mniejsza niż 2.5 (a więc komórki mają stosunkowo jednorodne rozmiary), są następnie oceniane na podstawie zmiennej V6, czyli obecność w pewnym sensie nietypowych jąder komórkowych. Komórki dla których ta zmienna jest mniejsza niż 5.5 są oceniane jako nowotwór łagodny. W tym liściu mamy 405 obserwacji odpowiadających nowotworom łagodnym i 5 odpowiadających nowotworom złośliwym, a więc - gdyby drzewo było idealnie poprawne - skuteczność klasyfikacji w tym węźle wynosiłaby około 98.8%. Komórki, dla których zmienna V6 jest większa niż 5.5, oceniamy jako nowotwór złośliwy, ze skutecznością 7/(1+7) = 87.5%.

Przycięcie drzewa, czyli (w tym przypadku) złączenie niektórych liści, może poprawić błąd testowy, ponieważ zmniejsza wariancję predykcji - liście z większą liczbą obserwacji są mniej podatne na losowe błędy klasyfikacji. Funkcja plotcp umożliwia sprawdzenie błędu testowego w zależności od parametru złożoności drzewa (cp, ang. complexity parameter). Im mniejszy parametr cp, tym bardziej złożone drzewo - dla nieskończonej wartości otrzymujemy jedynie korzeń, czyli klasyfikator, który każdej obserwacji przypisuje klasę która najczęściej występuje w danych.

plotcp(drzewo)

Wykres przedstawia błąd walidacji krzyżowej w odniesieniu do błędu drzewa z jednym liściem (dlatego wartości Inf odpowiada błąd równy 1 - jest to kwestia skalowania osi Y, a nie tego, że klasyfikator za każdym razem się myli). Dla parametru cp=0.011 otrzymujemy pełne drzewo. Przerywana linia pozioma reprezentuje błąd na pełnym drzewie plus jedno odchylenie standardowe tego błędu.
Wybieramy karę cp=0.037, ponieważ zapewnia ona mały błąd przy jak najprostszym modelu. Wykorzystujemy tę karę aby przyciąć drzewo, które następnie wizualizuemy ponownie:

przyciete_drzewo <- prune(drzewo, cp=0.037)
plot(przyciete_drzewo, main="Klasyfikacja złośliwości nowotworu")
text(przyciete_drzewo, use.n=TRUE, all=T, cex=.8)

Jak widać, otrzymane drzewo jest dużo prostsze. Między innymi omawiany wcześniej węzeł klasyfikujący nowotwory w zależności od zmiennej V6 został połączony w jeden liść, w którym dokładność klasyfikacji wynosi teraz 406/(406+12) = 97%.

Powyższą analizę można również przeprowadzić za pomocą pakietu caret. Dodatkowo pozwoli to na uzyskanie bezwzględnego błędu walidacji krzyżowej, zamiast błędu względem drzewa z jednym liściem. Taki błąd pozwoli nam porównać nasz model z modelem logistycznym, który badaliśmy na laboratorium 10.
Aby otrzymać błąd bezwzględny możemy również przeprowadzić walidację krzyżową ręcznie, podobnie jak w laboratorium 10. Niestety nie udało mi się znaleźć funkcji która zrobiła by to za nas, czyli przeprowadziła walidację krzyżową na gotowym drzewie (pakiet caret liczy drzewo ponownie) - jeśli komuś uda się taką funkcję znaleźć, będę wdzięczny za informację.

library(caret)
## Loading required package: lattice
## Loading required package: ggplot2
walidacja_caret <- train(class ~ ., 
                          data=biopsy, 
                          method="rpart", 
                          trControl = trainControl(method = "cv", number=10),
                         tuneLength=10)  # ten parametr kontroluje liczbe wartosci parametru cp do przetestowania
walidacja_caret
## CART 
## 
## 683 samples
##   9 predictors
##   2 classes: 'benign', 'malignant' 
## 
## No pre-processing
## Resampling: Cross-Validated (10 fold) 
## Summary of sample sizes: 615, 614, 615, 615, 615, 615, ... 
## Resampling results across tuning parameters:
## 
##   cp          Accuracy   Kappa    
##   0.00000000  0.9501705  0.8911613
##   0.08786611  0.9194160  0.8241398
##   0.17573222  0.9194160  0.8241398
##   0.26359833  0.9194160  0.8241398
##   0.35146444  0.9194160  0.8241398
##   0.43933054  0.9194160  0.8241398
##   0.52719665  0.9194160  0.8241398
##   0.61506276  0.9194160  0.8241398
##   0.70292887  0.9194160  0.8241398
##   0.79079498  0.7482097  0.3090500
## 
## Accuracy was used to select the optimal model using the largest value.
## The final value used for the model was cp = 0.

Najwyższa dokładność, odpowiadająca parametrowi cp = 0 (czyli pełnemu drzewu), wyniosła ~95%. Jest zatem nieco niższa niż w przypadku modelu logistycznego, dla którego wniosła 96.6%. Na podstawie wyników możemy również oszacować, że dla przyjętej wcześniej wartości cp=0.037 dokładność wynosi około 92%. Zauważmy, że kosztem nieco niższej dokładności otrzymujemy zdecydowanie prostszy model.

Zadanie 2. Drzewo klasyfikacyjnie domyślnie przypisuje nowej obserwacji tę klasę, która występuje najczęściej w odpowiednim liściu. Jakkolwiek rozwiązanie to jest bardzo intuicyjne, wiąże się z nim pewien problem. Jeśli test diagnostyczny mówi, że pacjent na 51% jest zdrowy, a na 49% ma złośliwy nowotwór, to niekoniecznie chcemy uznać że jest fajnie i odesłać go do domu.
Żeby uzyskać trochę bezpieczniejszą klasyfikację wykorzystujemy ważenie danych. Klasie \(k\) przypisujemy wagę \(w_k\), a nową obserwację przypisujemy do “cięższej” klasy. Na przykład, w liściu w którym mamy 1 nowotwór złośliwy i 5 łagodnych, nowej obserwacji przypiszemy klasę ‘złośliwy’ jeśli waga tej klasy będzie co najmniej 5 razy wyższa niż waga klasy ‘łagodny’. Jest to analogiczne podejście do progu prawdopodobieństwa, którego używaliśmy w regresji logistycznej. Jeśli wagi sumują się do 1, mówimy o prawdopodobieństwach klas a priori - interpretujemy je jako ogólne prawdopodobieństwo, że jakaś obserwacja należy do danej klasy (gdy nie znamy odpowiadających jej wartości predyktorów).
Aby przypisać nową obserwację do konkretnej klasy, mnożymy prawdopodobieństwa a priori przez prawdopodobieństwa empiryczne, czyli te, które otrzymujemy obliczając proporcje klas w danym liściu - jest to to samo, co unormowana “waga” danej klasy. Obserwacji przypisujemy tę klasę, dla której otrzymujemy najwyższą wartość takiego iloczynu.
W przypadku dwóch klas jest to równoznaczne z ustaleniem progu prawdopodobieństwa, tak samo jak w przypadku regresji logistycznej. Warto jednak wiedzieć, że drzewa decyzyjne pozwalają łatwo uogólnić tę procedurę na większą liczbę klas.

Podziel dane biopsy na zbiór testowy i treningowy. Zbuduj drzewo decyzyjne na zbiorze treningowym i przeprowadź predykcję na zbiorze testowym za pomocą funkcji predict. Funkcja ta zwróci prawdopodobieństwa empiryczne poszczególnych klas.
Utwórz wektor progów od 0 do 1. Dla każdego progu zaklasyfikuj obserwacje, przypisując klasę ‘malignant’ jeśli jej prawdopodobieństwo przekroczy próg, i oblicz True Positive Rate (procent nowotworów złośliwych zaklasyfikowanych jako łagodne) oraz False Positive Rate (procent nowotworów łagodnych zaklasyfikowanych jako złośliwe). W tym celu przyda się pętla for. W celu obliczenia TPR i FPR można stworzyć macierz błędu korzystając z funkcji table.
Utwórz krzywą ROC. Sprawdź, jaki procent nowotworów złośliwych zaklasyfikujemy jako łagodne w sytuacji, gdy (na przykład ze względów na koszty związane z dalszą terapią) możemy pozwolić sobie na zaklasyfikowanie jedynie 10% nowotworów łagodnych jako złośliwe.

Opcja dla ambitnych: Wykorzystując \(k\)-krotną walidację krzyżową utwórz \(k\) krzywych ROC i przedstaw je na jednym wykresie. W tym przypadku trzeba przeprowadzić walidację krzyżową ręcznie, na przykład korzystając z funkcji cut oraz which.

Uwaga. W odróżnieniu od regresji logistycznej pakiet rpart może wykorzystać wagi i macierz kosztu nie tylko podczas klasyfikowania nowej obserwacji, ale również podczas budowania drzewa, tak aby zoptymalizować pod tym kątem predykcję. Można to zrealizować komendą rpart(class~., data=biopsy, method='class', parms=list('prior'=wektor_unormowanych_wag)).
Prowadzi to jednak do sytuacji, w której dla każdego zestawu prawdopodobieństw otrzymujemy inny klasyfikator. Jedną z konsekwencji tego faktu jest to, że nie możemy zbudować sensownej krzywej ROC.
Bardziej ogólnym sposobem ważenia danych jest podanie macierzy kosztu \(L\), w której wyraz \(L_{ij}\) oznacza koszt predykcji klasy \(j\), gdy prawdziwą klasą jest \(i\). W typowych zastosowaniach \(L_{kk} = 0\), a wszystkie pozostałe wyrazy są dodatnie. Funkcja rpart umożliwia również korzystanie z macierzy kosztu podczas budowania drzewa.

library(ggplot2)
ggplot(roc_data, aes(x=FPR, y=TPR, group=CV)) + geom_line(alpha=0.4, color='orange') + theme_minimal()

Lasy losowe

Znajomy chemik kwantowy, gdy pierwszy raz zobaczył wzór na regresję LASSO, skonstatował: “Aha. Czyli chodzi o to, żeby dopasować prostą do punktów, tylko tak niezbyt dobrze. Wy naprawdę nie macie nic lepszego do roboty na tej statystyce?”.

Lasy losowe to kolejna metoda w której poprawiamy poprzez psucie. W pierwszym kroku na podstawie naszych danych losujemy pewną liczbę \(B\) prób bootstrapowych (czyli \(B\)-krotnie losujemy z naszych danych \(n\) obserwacji ze zwracaniem, gdzie \(n\) zawsze jest równe wyściowej liczbie obserwacji). W następnym kroku na każdej próbie bootstrapowej budujemy drzewo losowe, ale żeby przypadkiem nie zbudować go zbyt dobrze, za każdym razem dzieląc węzeł losujemy \(m\) predyktorów, i dzielimy węzeł korzystając jedynie z nich. Okazuje się, że to poprawia błąd testowy. Formalny dowód tego faktu można znaleźć m.in. w oryginalnej pracy Leo Breimana z 2001 r. Z kolei proste wprowadzenie do tematyki lasów losowych można obejrzeć tutaj.

Do budowania lasów losowych w R służy funkcja randomForest z biblioteki o tej samej nazwie.
Cała procedura wygląda tak samo jak w przypadku pozostałych metod - wpisanie w konsolę nazwy utworzonego modelu wyświetli podstawowe informacje, działają również funkcje summary, plot, predict itd.

Walidację krzyżową lasu losowego można przeprowadzić ręcznie (patrz Lab10) lub korzystając z pakietu caret. Sposób jego użycia jest praktycznie taki sam jak w przypadku pozostałych modeli - szczegóły można znaleźć wpisując w Google “random forests in caret”.
Korzystając z funkcji randomForest warto przekazać argument importance = TRUE, aby funkcja oceniła istotność zmiennych - czyli intuicyjnie mówiąc to, jak duży wpływ ma każda z nich na predykcję, lub jak mocno zmieniła by się predykcja gdyby tę zmienną usunąć z danych. Jest to co innego niż istotność statystyczna - zbieżność nazw w języku polskim jest przypadkowa. Dokładniejsze wyjaśnienie istotności cech w lasach losowych można znaleźć m.in. tutaj lub w książce The Elements of Statistical Learning.

Zadanie 3. Zbuduj las losowy przewidujący złośliwość nowotworu na danych biopsy. Sprawdź, czy jego dokładność w walidacji krzyżowej przebije 96%, czyli tę którą osiągnęliśmy korzystając z modelu logistycznego. Estymowany błąd testowy znajdziesz w podstawowych informacjach dotyczących utworzonego modelu, które uzyskasz wpisując jego nazwę w konsolę. Jest on opisany jako ‘OOB estimate of error rate’ - jest to tak zwany błąd out-of-bag, który przybliża błąd walidacji krzyżowej. Przedstaw istotność zmiennych na wykresie, wykorzystując funkcję varImpPlot lub bibliotekę ggplot2. Możesz również przetestować różne liczby zmiennych losowanych przy podziale węzłów modyfikując atrybut mtry. Możesz też zmienić liczbę drzew za pomocą argumentu ntrees. Może być tu pomocna funkcja plot, która wywołana na utworzonym modelu przedstawi zależność błędu testowego od liczby użytych drzew. W praktyce bierze się po prostu “dostatecznie dużo” drzew (na przykład 500) i nie przeprowadza dodatkowych analiz dotyczących ich optymalnej liczby, ponieważ ma to bardzo nieznaczny wpływ na błąd. Więcej informacji na ten temat można przeczytać tutaj.

Dla chętnych: Przeprowadź ręcznie walidację krzyżową i porównaj z estymatorem OOB.

Dla ułatwienia, stworzony model powinien wyglądać następująco:

las
## 
## Call:
##  randomForest(formula = class ~ ., data = biopsy, importance = TRUE,      mtry = 4, ntree = 400) 
##                Type of random forest: classification
##                      Number of trees: 400
## No. of variables tried at each split: 4
## 
##         OOB estimate of  error rate: 2.93%
## Confusion matrix:
##           benign malignant class.error
## benign       432        12  0.02702703
## malignant      8       231  0.03347280

Poniżej zamieszczam również wartość średnią oraz wykres pudełkowy błędu testowego obliczonego na podstawie ręcznie przeprowadzonej 10-krotnej walidacji krzyżowej:

mean(bledy_walidacji)
## [1] 0.0445735
ggplot(data.frame('Blad'=bledy_walidacji)) + geom_boxplot(aes(y=Blad)) + theme_minimal()

Tutaj natomiast mamy ręcznie utworzony wykres słupkowy przedstawiający istotność zmiennych (atrybut importance utworzonego modelu, zapisanego w zmiennej las, jest nieco nietypową zmienną, dlatego trzeba się trochę nakombinować żeby przekształcić go do ramki danych):

waznosc <- las$importance
ggplot(data.frame('Zmienna' = rownames(waznosc), 'Waznosc' = as.numeric(waznosc))) + geom_col(aes(x=Zmienna, y=Waznosc)) + coord_flip() + theme_minimal() + ggtitle('Ważność zmiennych w lesie losowym')

Zadania dodatkowe.

Zadanie 4. Indeks Giniego. Załóżmy, że mamy 20 obserwacji - 10 zielonych, 10 fioletowych. Rozpatrzmy drzewo klasyfikacyjne, które ma dwa liście.

  1. Załóżmy, że do każdego liścia przypisujemy 10 obserwacji. Stwórz wykres przedstawiający indeks Giniego dla pierwszego liścia w zależności od liczby przypisanych mu fioletowych obserwacji.
  2. Przedstaw zależność indeksu Giniego od liczby obserwacji fioletowych przypisanych do pierwszego liścia w sytuacji, gdy do liscia przypisujemy dowolną liczbę obserwacji (od 1 do 19). Możesz zrobić to na dwa sposoby:

Zadanie 5. Boosting. Obejrzyj dane Khan z pakietu ISLR. Będziemy zajmować się przewidywaniem typu nowotworu na podstawie ekspresji genów. Zauważ, że dane są listą podzieloną już na zbiór treningowy i testowy.

Wykonaj procedurę boosting do predykowania typu nowotworu (zmienna y), korzystając z funkcji gbm z biblioteki o tej samej nazwie. Użyj zbioru train do trenowania, a test do testowania. Wygeneruj 500 drzew (domyślna wartość 100 jest trochę mała). Przetestuj różne wartości parametru shrinkage, najlepiej wykorzystując siatkę parametrów. Przykładową siatkę parametrów długości n możesz stworzyć komendą seq(0.001, 0.9, length.out = n). Narysuj wykres zależności błędu testowego w zależności od parametru shrinkage.

Dla chętnych: za pomocą walidacji krzyżowej oszacuj, najlepszą wartość parametru shrinkage oraz liczby generowanych drzew.

Uwaga: w przypadku lasów losowych i baggingu nie da się przetrenować modelu zwiększając liczbę drzew; w przypadku boostingu owszem. Liczba drzew zatem również powinna zostać ostrożnie dobrana, np. za pomocą walidacji krzyżowej.

Zadanie 6. Porównaj drzewa klasyfikacji i lasy losowe na danych Khan z pakietu ISLR:

  1. Wykonaj klasyfikację za pomocą drzewa klasyfikującego korzystając z funkcji rpart z pakietu o tej samej nazwie. Zwizualizuj drzewo za pomocą funkcji plot. Które predyktory zostały wykorzystane? Porównaj ich średnie i wariancje między klasami.
  2. Przytnij uzyskane drzewo tak, żeby miało trzy liście. Wykorzystaj funkcję prune.
  3. Zbuduj las losowy, wykorzystując w każdej iteracji maksymalnie 20 predyktorów. Policz średni błąd testowy. Wykorzystaj funkcję randomForest z biblioteki o tej samej nazwie; najlepiej uruchomić z argumentem importance = TRUE.
  4. Powtórz punkt 3 od kilku do kilkunastu razy dla różnej liczby rozważanych predyktorów. Narysuj wykres zależności błędu testowego od liczby predyktorów. Czy podobne predyktory zostały uznane za najistotniejsze przez las losowy i zwykłe pojedyncze drzewo? Czy różna liczba uwzględnianych predyktorów przy budowaniu lasu wpłynęła na wybór najistotniejszych predyktorów?
  5. Narysuj wykres istotności predyktorów za pomocą funkcji varImpPlot, najlepiej dla kilku różnych lasów losowych.

Po wykonaniu tego zadania zastanów się, jak przy użyciu funkcji randomForest wykonać procedurę bagging.