Lab 9 - klasyfikacja

Regresja logistyczna

Zadanie 1. Zadanie przykładowe. W tym zadaniu spróbujemy przewidzieć, czy wynik giełdy w danym dniu będzie dodatni czy ujemny na podstawie wyników z pięciu poprzednich dni. W tym celu wykorzystamy regresję logistyczną. Regresja logistyczna została omówiona na Wykładzie 8. Elementarne wprowadzenie do regresji logistycznej można również usłyszeć tutaj. Więcej na temat tej techniki można przeczytać tutaj oraz tutaj.

Zainstaluj i załadu pakiet ISLR i dane Smarket, zawierające wartości indeksu giełdowego S&P 500. Przeczytaj informacje na temat danych korzystając z funkcji help(Smarket) oraz obejrzyj dane korzystając albo z funkcji head, albo View. Sprawdź, ile jest obserwacji oraz zmiennych, które zmienne są jakościowe, a które ilościowe. Co oznacza i ile wartości przyjmuje zmienna Direction? Co oznaczają zmienne o nazwie Lag?

Zbuduj model regresji logistycznej, który na podstawie zmiennych Lag1, Lag2, …, Lag5 wykona predykcję dla zmiennej Direction. Wykorzystaj funkcję glm z argumentem family = binomial.

Sprawdź podstawowe informacje o stworzonym modelu, korzystając z funkcji summary. Możesz również sprawdzić jakie informacje zawiera model za pomocą funkcji attributes. W jaki sposób wynik giełdy w danym dniu zależy od wyniku sprzed dwóch dni?

Oblicz prawdopodobieństwa przynależności każdej z obserwacji do klasy ”Up” za pomocą funkcji predict. Za pomocą wykresów pudełkowych przedstaw zależność pomiędzy estymowanym prawdopodobieństwem a prawdziwą wartością zmiennej Direction. Czy wyniki zgadzają się z intuicją?

Wygeneruj wektor predykowanych wartości ”Up” i ”Down”, przyjmując jako poziom odcięcia wartość prawdopodobieństwa 0.5. Oblicz średni błąd klasyfikatora (w tym zadaniu zbadamy jedynie błąd treningowy naszego klasyfikatora). Utwórz macierz błędu i oblicz czułość oraz specyficzność wykrywania wzrostu na giełdzie. Zinterpretuj wyniki.

Powtórz klasyfikację dla 10 poziomów odcięcia od 0.4 do 0.6. Dla każdego poziomu oblicz False Positive Rate i True Positive Rate (uwaga, nie mylić z False Discovery Rate!). Utwórz krzywą ROC, czyli zależność TPR od FPR.

Załóżmy, że chcemy zainwestować jeśli model wskaże wzrost na giełdzie. Na podstawie krzywej ROC oceń, jakiego typu wartość progową należy przyjąć, jeśli chcemy grać bezpiecznie i nie lubimy ryzyka. A jaką, jeśli chcemy grać agresywnie i maksymalizować szanse na zysk?

Uwaga. Przewidywanie wartości procesu losowego (na przykład wyniku giełdy) w danym dniu, korzystając z obserwacji w kilku poprzednich dniach, to tak zwane one-step forecasting. W przypadku, gdy chcemy przewidzieć wyniki na kilka dni do przodu, mówi się o multi-step forecasting. Więcej na temat tego typu predykcji można dowiedzieć się na kursie Szeregi Czasowe.

Rozwiązanie. W pierwszym kroku ładujemy i poglądamy dane.

library(ISLR)
data("Smarket")
head(Smarket)

##   Year   Lag1   Lag2   Lag3   Lag4   Lag5 Volume  Today Direction
## 1 2001  0.381 -0.192 -2.624 -1.055  5.010 1.1913  0.959        Up
## 2 2001  0.959  0.381 -0.192 -2.624 -1.055 1.2965  1.032        Up
## 3 2001  1.032  0.959  0.381 -0.192 -2.624 1.4112 -0.623      Down
## 4 2001 -0.623  1.032  0.959  0.381 -0.192 1.2760  0.614        Up
## 5 2001  0.614 -0.623  1.032  0.959  0.381 1.2057  0.213        Up
## 6 2001  0.213  0.614 -0.623  1.032  0.959 1.3491  1.392        Up

Mamy do czynienia z 9 zmiennymi, z których ostatnia jest jakościowa, a pozostałe są ilościowe. Sprawdzamy poziomy zmiennej jakościowej Direction:

levels(Smarket$Direction)

## [1] "Down" "Up"

Możemy również sprawdzić wymiary naszej ramki danych komendą dim:

dim(Smarket)

## [1] 1250    9

Mamy 1250 obserwacji i 9 cech.

W dokumentacji zbioru danych Smarket, wyświetlonej komendą ?Smarket, możemy przeczytać, że zmienne Lag to wartość indeksu w poprzedzających dniach. Przechodzimy do budownia modelu. Zwróćmy uwagę, że poziomem bazowym zmiennej Direction jest Down. To oznacza, że wynikiem modelowania będzie prawdopodobieństwo tego, ze zmienna przyjmie wartość Up.

Do zbudowania modelu wykorzystujemy funkcję glm. Jest to funkcja służąca do budowania tzw. uogólnionych modeli liniowych, do których należy m.in. zwykła regresja liniowa oraz regresja logistyczna. Domyślnie funkcja ta zbuduje zwykły model liniowy. Atrybut family = binomial sprawi, że zbudujemy model logistyczny.

model <- glm(Direction ~ Lag1 + Lag2 + Lag3 + Lag4 + Lag5, Smarket, family = binomial)
summary(model)

## 
## Call:
## glm(formula = Direction ~ Lag1 + Lag2 + Lag3 + Lag4 + Lag5, family = binomial, 
##     data = Smarket)
## 
## Deviance Residuals: 
##    Min      1Q  Median      3Q     Max  
## -1.376  -1.204   1.071   1.145   1.347  
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)
## (Intercept)  0.074163   0.056674   1.309    0.191
## Lag1        -0.071325   0.050104  -1.424    0.155
## Lag2        -0.044136   0.050025  -0.882    0.378
## Lag3         0.009229   0.049879   0.185    0.853
## Lag4         0.007211   0.049898   0.145    0.885
## Lag5         0.009311   0.049490   0.188    0.851
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 1731.2  on 1249  degrees of freedom
## Residual deviance: 1728.3  on 1244  degrees of freedom
## AIC: 1740.3
## 
## Number of Fisher Scoring iterations: 3

Przypomnienie z wykładu. W regresji logistycznej rozpatrujemy binarną zmienną objaśnianą $Y \in \{0, 1\}$ i modelujemy prawdopodobieństwo zdarzenia $Y=1$ jako $$\mathbb{P}(Y = 1) = \frac{1}{1 + \exp(-\beta_0 - \beta_1 X_1 - ... - \beta_k X_k)}.$$ To oznacza, że jeśli $\beta_i > 0$, to prawdopodobieństwo rośnie do 1 wraz ze wzrostem $X_i$, i spada do 0 wraz ze spadkiem $X_i$. Wartość parametru $\beta_i$ kontroluje tempo wzrostu prawdopodobieństwa. Parametr $\beta_0$ kontroluje prawdopodobieństwo zdarzenia $Y=1$ w sytuacji, gdy wszystkie zmienne $X_i$ mają wartość 0. Im większy ten parametr, tym większe jest to prawdopodobieństwo.

Na podstawie wyniku funkcji summary widzimy, że prawdopodobieństwo dodatniego wyniku giełdy w danym dniu rośnie tym bardziej, im gorszy był wynik w dniu poprzednim. Wartości parametrów $\beta$ są jednak nieistotne statystycznie. To oznacza, że nie możemy wykluczyć tego, że w rzeczywistości wartość indeksu w danym dniu nie zależy od wartości w dniu poprzednim, a nawet, że prawdziwa zależność jest odwrotna niż to sugeruje model - a więc, że dodatni wynik w danym dniu zwiększa szanse na dodatni wynik następnego dnia. Możemy uzyskać więcej informacj konstruując przedziały ufności dla parametrów $\beta$:

confint(model)

## Waiting for profiling to be done...

##                   2.5 %     97.5 %
## (Intercept) -0.03686767 0.18534962
## Lag1        -0.17004940 0.02662310
## Lag2        -0.14254775 0.05381678
## Lag3        -0.08864473 0.10714767
## Lag4        -0.09071165 0.10515270
## Lag5        -0.08777508 0.10649382

Jak widać, wszystkie przedziały ufności zawierają zero. Mimo tego, że wyestymowana wartość parametru odpowiadającego zmiennej Lag1 wynosi -0.07, to nie możemy wykluczyć, że rzeczywista wartość to 0, -0.1, lub nawet 0.02. Musimy zatem traktować predykcje takiego modelu z dużą dozą nieufności.

Obliczmy prawdopodobieństwa tego, że zmienna Direction przyjmie wartość Up, dla każdego wiersza naszych danych. Wykorzystamy do tego funkcję predict. Domyślnie funkcja ta zwraca wartość $\log(p/(1-p)) = \beta_0 + \beta_1 X_1 + ... + \beta_k X_k$, gdzie $p = \mathbb{P}(Y=1)$. Argument type = 'response' sprawia, że funkcja wraca wartość $\mathbb{P}(Y=1)$.

prawdopodobienstwa <- predict(model, type = "response")

Zwizualizujemy rozkłady przewidywanych prawdopodobieństw za pomocą wykresów kubełkowych:

library(ggplot2)

## Registered S3 methods overwritten by 'ggplot2':
##   method         from 
##   [.quosures     rlang
##   c.quosures     rlang
##   print.quosures rlang

ggplot(data = data.frame('P' = prawdopodobienstwa, 'Direction' = Smarket$Direction)) + geom_boxplot(aes(x = Direction, y = P)) + theme_minimal()

Wykres kubełkowy (boxplot) czytamy następująco: Pozioma gruba kreska zaznacza medianę, pudełko zaznacza przedział międzykwartylowy (czyli przedział od kwantylu na poziomie 0.25 do kwantylu na poziomie 0.75), “wąsy” mają długość równą półtorej długości przedziału międzykwartylowego (chyba że miałyby wykroczyć poza najmniejszą lub największą wartość w danych, wtedy je się skraca), a kropki oznaczają pozostałe obserwacje. W przypadku rozkładu normalnego zakres wąsów odpowiada 99% prawdopodobieństwa. Wykres kubełkowy to bardzo powszechny sposób przedstawiania danych, ale może być bardzo zdradliwy. Więcej na ten temat można przeczytać tutaj.

Z powyższych wykresów widać, że model działa raczej średnio - nawet gdy zmienna Direction przyjmuje wartość Down, to mediana prawdopodobieństwa ‘Up’ jest większa niż 1/2. Jedyne co nas pociesza, to że mediana tego prawdopodobienstwa jest przynajmniej nieco wyższa gdy prawdziwą wartością rzeczywiście jest ‘Up’.

Zrobimy teraz predykcję zmiennej Direction przyjmując poziom odcięcia 0.5. To oznacza, że przyjmujemy, że $Y=1$ wtedy, gdy $\mathbb{P}(Y=1) \geq 1/2$. Wykorzystamy w tym celu funkcję ifelse(A, B, C), która przyjmuje wektor logiczny A, i dla każdej współrzędnej tego wektora zwraca B gdy współrzędna ta ma wartość TRUE, a C, gdy współrzędna ma wartość FALSE.

predykcja <- ifelse(prawdopodobienstwa  > 0.5, "Up", "Down")

Tworzymy macierz błędu (a.k.a. macierz konfuzji). Zwróćmy uwagę, że jeśli uznamy że nasz klasyfikator ma za zadanie wykryć wzrost na giełdzie, to wówczas wzrost interpretujemy jako wynik pozytywny.

conf.matrix <- table('Pred'=predykcja, 'True'=Smarket$Direction)
conf.matrix

##       True
## Pred   Down  Up
##   Down  116  98
##   Up    486 550

Obliczamy dokładność (a.k.a. accuracy), czułość (a.k.a. sensitivity, recall, True Positive Rate), specyficzność (a.k.a. specificity, True Negative Rate) oraz precyzję (a.k.a. precision, Positive Predictive Value):

accuracy <- sum(diag(conf.matrix))/sum(conf.matrix)  # suma wartości na przekątnej dzielona przez sumę wszystkich wartości macierzy
sensitivity <- conf.matrix[2, 2]/(conf.matrix[1,2]+conf.matrix[2,2])
specificity <- conf.matrix[1, 1]/(conf.matrix[1,1]+conf.matrix[2,1])
precision <- conf.matrix[2, 2]/(conf.matrix[2,1] + conf.matrix[2, 2])
c('Dokładność' = accuracy, 'Czułość' = sensitivity, 'Specyficzność' = specificity, "Precyzja" = precision)

##    Dokładność       Czułość Specyficzność      Precyzja 
##     0.5328000     0.8487654     0.1926910     0.5308880

Jak widać, ogółem nasz klasyfikator jest właściwie beznadziejny - jego dokładność wynosi 0.53, czyli jest niewiele lepsza niż klasyfikatora losowego. Ma jednak jedną mocną stronę: wysoka czułość oznacza, że skutecznie wykrywa możliwość wzrostu na giełdzie. Niestety, nie oznacza to, że możemy w pełni ufać jego wskazaniom: nawet jeśli nasz klasyfikator sygnalizuje wzrost, to mamy jedynie 53% prawdopodobieństwa, że wzrost rzeczywiście nastąpi. Z kolei jeżeli na giełdzie wystąpi spadek, to mamy jedynie 19% prawdopodobieństwa że klasyfikator nas o tym ostrzeże.

Powyższe wyniki zależą oczywiście od przyjętego progu odcięcia, który mogliśmy wybrać bardzo niefortunnie. Żeby ocenić ogólną jakość klasyfikatora, niezależnie od progu odcięcia, stosuje się krzywą ROC (ang. receiver operator characteristic curve). Żeby ją utworzyć, potrzebujemy dwóch miar jakości klasyfikatora: spotkane już True Positive Rate, czyli stosunek wyników prawdziwie pozytywnych do wszystkich obserwacji pozytywnych w danych, oraz False Positive Rate, czyli stosunek wyników fałszywie pozytywnych do wszystkich obserwacji negatywnych w danych. Krzywa ROC to wykres opisany poniższym wzorem: $$\text{ROC} = \big \{ \big( \text{FPR}(p), \text{TPR}(p) \big) \colon p \in [0, 1] \big \}.$$ Wartość $p$ powyżej to przyjęty próg odcięcia. Dla $p=0$ wszystkim obserwacjom jest przypisywany wynik pozytywny, wobec czego zarówno TPR oraz FPR są równe 1: klasyfikator idealnie wykrywa obserwacje pozytywne i idealnie myli się w przypadku obserwacji negatywnych. Dla $p=1$ sytuacja się odwraca, i TPR oraz FPR są równe 0. Dla innych progów krzywa ROC mówi nam, z jaką ilością wyników fałszywie pozytywnych musimy się pogodzić (FPR), jeśli chcemy poprawnie wykryć daną proporcję obserwacji pozytywnych (TPR).

Stwórzmy teraz wektor zawierający 20 kolejnych progów odcięcia. W naszym przypadku prawdopodobieństwa zdarzenia $Y = 1$ wahają się od 0.40 do 0.62, więc nie musimy rozpatrywać wszystkich progów odcięcia - próg 0.1 da taką samą klasyfikację, jak próg 0.2, i oba progi będą odpowiadały temu samemu punktowi na krzywej ROC.

progi <- seq(0.4, 0.65, length.out = 20)

Najkrótszym sposobem obliczenia TPR jest wykorzystanie funkcji sapply w następujący sposób:

macierz_predykcji <- sapply(progi, function(p) prawdopodobienstwa > p) 

Powyższy kawałek kodu utworzył macierz rozmiaru 1250 x 20, w której każda kolumna odpowiada klasyfikacji przy zadanym progu odcięcia. Wartość TRUE w macierzy oznacza, że prawdopodobieństwo jest powyżej progu odcięcia, a więc że przyjmujemy klasyfikację pozytywną ($Y = 1$, Direction = Up).
Teraz korzystamy z funkcji apply, aby obliczyć FPR i TPR dla każdej kolumny:

P <- sum(Smarket$Direction == "Up")  
N <- sum(Smarket$Direction == "Down")
TPR <- apply(macierz_predykcji, 2, function(x) sum(x & Smarket$Direction=='Up'))/P
FPR <- apply(macierz_predykcji, 2, function(x) sum(x & Smarket$Direction=='Down'))/N

Oprator & to wektorowa koniunkcja: dla dwóch wektorów logicznych x oraz y napis x & y zwraca wektor, w którym i-ta współrzędna odpowiada koniunkcji i-tych współrzędnych tych dwóch wektorów. Sumowanie wektora logicznego zwraca liczbę wystąpień wartości TRUE. Altrernatywnym sposobem obliczenia TPR i FPR to napisanie funkcji, która obliczy te wartości przy zadanym progu odcięcia, i wywołanie jej za pomocą sapply() na wektorze progów.

Mamy teraz dwa wektory, TPR oraz FPR, odpowiadające kolejnym progom odcięcia. Możemy zatem stworzyć krzywą ROC:

ggplot(data.frame('TPR' = TPR, 'FPR' = FPR, 'p' = progi)) + geom_line(aes(x=FPR, y=TPR)) + geom_abline(slope=1, intercept=0, alpha=0.2) + theme_minimal() 

Klasyfikator losowy (przypisujący każdej obserwacji losową klasę) spełnia zależność $\text{FPR} \approx \text{TPR}$. Takiemu klasyfikatorowi odpowiada na powyższym wykresie szara ukośna linia. Badany klasyfikator jest tym lepszy, im bardziej “na lewo i w górę” znajduje się jego krzywa ROC. Jak widać, w naszym przypadku krzywa ROC jest jedynie nieznacznie powyżej klasyfikatora losowego. To oznacza, że tak ponure wyniki jakie otrzymaliśmy wcześniej nie były spowodowane jedynie niefortunnym wyborem progu odcięcia.

Poniższy wykres przedstawia krzywą ROC, w której oznaczono progi odcięcia odpowiadające kolejnym punktom.

ggplot(data.frame('TPR' = TPR, 'FPR' = FPR, 'p' = progi)) + geom_line(aes(x=FPR, y=TPR)) + geom_abline(slope=1, intercept=0, alpha=0.2) + theme_minimal() + geom_text(aes(x=FPR, y=TPR + 0.05, label=round(p, 2)))

Praktyczne wnioski z powyższej krzywej są następujące. Załóżmy, że wygrywamy, jeśli poprawnie założymy że wystąpi wzrost na giełdzie. Jeśli chcemy grać bezpiecznie, przyjmiemy próg który maksymalizuje specyficzność, czyli minimalizuje FPR przy zachowaniu rozsądnego poziomu czułości. To nam pozwoli ograniczyć ryzyko, ale jednocześnie wykryć jakieś szanse wygranej. Próg odcięcia będzie duży.

Jeśli chcemy maksymalizować potencjał zysku, maksymalizujemy TPR, przy zachowaniu rozsądnego poziomu FPR. To nam pozwoli wykryć możliwie dużo okazji do wygranej, ale jednocześnie w jakimś stopniu kontrolować ryzyko. Próg odcięcia będzie mały.

Co oznacza powyżej sformułowanie “rozsądny poziom FPR”, to już kwestia pozamatematematyczna - musimy sami wiedzieć, jakie ryzyko możemy dopuścić podczas inwestowania oraz jaki minimalny zysk chcemy uzyskać. Z wykresu ROC z zaznaczonymi progami widać, że przykładowe sensowne progi to np. 0.54 dla gdy bezpiecznej i 0.49 dla gry agresywnej.

Zadanie 2. Regresja logistyczna pozwala na przypisanie obserwacji do jednej z dwóch klas. Istnieją sposoby na rozszerzenie tej techniki na przypadek $K > 2$ klas. Jednym z nich jest wytrenowanie $K$ klasyfikatorów binarnych, gdzie $i$-ty klasyfikator decyduje czy obserwacja należy do klasy $i$, czy do jakiejś innej klasy. Każdej obserwacji przypisujemy następnie klasę o najwyższym prawdopodobieństwie. Jest to tak zwana klasyfikacja one-vs-all lub one-vs-rest. Innym sposobem na poradzenie sobie z większą liczbą klas, który wykracza jednak poza program tego przedmiotu, jest regresja wielomianowa.

Załaduj dane iris. Podziel dane na zbiór treningowy i testowy. Wykorzystaj w tym celu funkcję sample.split z pakietu caTools, która zadba o to, aby w obu zbiorach znalazły się podobne proporcje różnych gatunków. Jest to o tyle ważne, że jeśli do obu zbiorów przypisujemy obserwacje całkowicie losowo, to może się zdarzyć, że w zbiorze testowym dostaniemy wyłącznie jeden gatunek. To z kolei sprawi, że nasz błąd testowy nie będzie odawał rzeczywistej jakości klasyfikatora.

Dodaj do zbioru treningowego trzy kolumny binarne (czyli factory o dwóch poziomach), gdzie i-ta kolumna będzie zawierała informację czy obserwacja należy do i-tego gatunku, czy nie. Utwórz trzy modele logistyczne objaśniające te nowe kolumny za pomocą zmiennej Sepal.Length. Alternatywnie, zamiast dodawać trzy kolumny do tabeli iris, możesz utworzyć trzy nowe ramki danych. Możesz też wybrać inny zestaw zmiennych objaśniających.

Dla każdej obserwacji ze zbioru testowego oblicz, korzystając z funkcji predict(), wektory prawdopodobieństw przynależności do poszczególnych gatunków. Przypisz każdej obserwacji gatunek o najwyższym prawdopodobieństwie. Może się tu przydać funkcja max.col. Utwórz macierz błędu. Jaka jest ogólna dokładność klasyfikatora? Który gatunek jest najprecyzyjniej klasyfikowany? Możesz również sprawdzić, co się stanie, jeśli dodamy do modelu więcej zmiennych objaśniających.

Przykładowy wynik:

table('Predicted' = class, 'True' = iris[-train, 'Species'])

##             True
## Predicted    setosa versicolor virginica
##   setosa         46         11         1
##   virginica       0         20        41
##   versicolor      3         19         8

Linear Discriminant Analysis, LDA

Zadania dodatkowe

Zadanie 4. Podobnie jak dla LDA, wytrenuj model na losowej połowie danych i wykonaj predykcję dla drugiej połowy, tym razem używając metody k najbliższych sąsiadów. Wykorzystaj funkcję knn z biblioteki class; zwróć uwagę, że ta funkcja ma kompletnie inną składnię niż w dotychczasowych modelach. Narysuj wykres błędu w zależności od wartości k od 1 do 20.

Zadanie 5. Do przewidywania wyników giełdy bardzo często stosuje się jeszcze jedną technikę, z którą się już spotkaliśmy: regresję liniową. W tak zwanym modelu autoregresyjnym zakłada się, że wynik danego dnia zależy liniowo od wyników z kilku poprzednich dni.

Przeprowadź regresję liniową zmiennej Today na podstawie zmiennych Lag. Następnie na podstawie tego modelu wykonaj predykcję zmiennej Today. Na jej podstawie utwórz predykcję zmiennej Direction: przyjmij wartość Up tam, gdzie predykcja zmiennej Today jest dodatnia. Oceń takie podejście korzystając z poznanych miar jakości klasyfikatorów. Opcjonalnie, możesz dodać interakcje pomiędzy zmiennymi.

Zadanie 6. Przeprowadź predykcję wartości zmiennej Direction (lub Today) w następujący sposób:

• Przyjmij, że dzisiaj będzie to samo, co wczoraj.

Zbadaj jak dobrze działa takie podejście korzystając z poznanych już miar jakości klasyfikatorów. Możesz wybrać inną liczbę dni do klasyfikacji i przyjąć, że danego dnia będzie to, co występowało najczęściej w ostatnich dniach. Jeśli przewidujesz zmienną Today, możesz również porównać to podejście z modelem autoregresyjnym za pomocą błędu średniokwadratowego.