Na dzisiejszych zajęciach zajmiemy się bardziej zaawansowanymi metodami estymacji błędu testowego w zadaniach klasyfikacji i regresji. Przed przystąpieniem do zadań zainstaluj biblioteki caret. Jest to bardzo rozbudowana biblioteka służąca do budowy i weryfikacji różnych modeli. Więcej informacji na jej temat można przeczytać tutaj.

Walidacja krzyżowa

Walidacja krzyżowa (ang. k-fold cross-validation) polega na podziale danych na \(k\) rozłącznych zbiorów testowych. Następnie dla każdego zbioru testowego trenujemy model korzystając z pozostałych danych i oceniamy jego błąd. Otrzymujemy w ten sposób \(k\) estymatorów błędu testowego.
Ważne jest to, żeby obserwacje przypisywać do zbiorów testowych w sposób losowy. Wynika to z faktu, że w praktyce dane są często w jakiś sposób uporządkowane. Biorąc do zbioru testowego kolejne “kawałki” uporządkowanych danych zaburzamy estymację błędu testowego. Trzeba jednocześnie pamiętać że każda obserwacja ma trafić co najwyżej do jednego zbioru testowego.

Ilustracja: K-krotna walidacja krzyżowa. Źródło: Wikipedia, Cross-validation (statistics)

Ilustracja: K-krotna walidacja krzyżowa. Źródło: Wikipedia, Cross-validation (statistics)

Alternatywnym podejściem jest metoda Monte Carlo cross-validation, w której tworzymy wiele (np. 1000) zbiorów testowych, zawierających np. po 1/10 obserwacji, za każdym razem wybierając zbiór testowy losowo (bez zwracania). Tak jak w poprzedniej metodzie, każdemu zbiorowi testowemu odpowiada zbiór treningowy złożony z pozostałych obserwacji.

W praktyce zamiast trzymać zbiory testowe i treningowe w osobnych zmiennych, warto pamiętać wyłącznie numery wierszy odpowiadające danemu zbiorowi testowemu. Pozwala to łatwo otrzymywać odpowiednie zbiory testowe i treningowe, a jednocześnie znacznie zaoszczędzić pamięć komputera. Przydaje się tu funkcja cut z atrybutem labels=FALSE.

Którą metodę wybrać? To zależy. Zaletą metody k-fold jest to, że każda obserwacja jest wykorzystana zarówno do treningu jak i testowania. W metodzie Monte Carlo pewne obserwacje mogą w ogóle nie wystąpić w żadnym zbiorze testowym, a inne mogą wystąpić w wielu. Otrzymujemy za to więcej estymatorów błędu, co pozwala lepiej przedstawić jego rozkład.

W rzeczywistości nie zawsze jest tak łatwo. Na dzisiejszych zajęciach będziemy pracować na prostym zbiorze danych, w którym możemy po prostu wybierać obserwacje losowo. W wielu praktycznych zastosowaniach metoda walidacji krzyżowej jest dużo bardziej złożona. Ma to miejsce na przykład wtedy, gdy dane mają pewną specyficzną strukturę - na przykład, kiedy mamy dane pochodzące od zbioru pacjentów, a od każdego pacjenta pobierane jest kilka próbek. Taką strukturę należy uwzględnić w walidacji krzyżowej.
Niestety, jest to dużo bardziej złożony problem, i na ogół rozwiązanie trzeba wymyślić dla każdego zastosowania z osobna. Przykłady prac naukowych dotyczących optymalnej walidacji krzyżowej w poszczególnych zastosowaniach można zobaczyć tutaj i tutaj.

Kolejnym praktycznym problemem w walidacji krzyżowej są mocno niezbalansowane klasy - na przykład, kiedy budujemy klasyfikator binarny na danych, w których mamy 1000 obserwacji typu 1 i 10 obserwacji typu 2. Zwykła walidacja krzyżowa nie działa w takim przypadku z dwóch powodów: po pierwsze, przypisanie wszystkim obserwacjom klasy 1 daje bardzo dobry klasyfikator, a po drugie, ponieważ w niektórych zbiorach testowych (a nawet treningowych) mogą się w ogóle nie pojawić obserwacje typu 2.
Proste wprowadzenie do walidacji krzyżowej na niezbalansowanych danych, z przykładami w języku Python, można przeczytać tutaj. Warto jednak zdawać sobie sprawę z tego, że tematyka niezbalansowanych danych jest na tyle złożona, że istnieje niejedna książka na ten temat.

Zadanie 1. Zadanie przykładowe. W tym zadaniu zbudujemy model logistyczny przewidujący złośliwość nowotworu na podstawie innych cech. Następnie przeprowadzimy 10-krotną walidację krzyżową w celu oszacowania błędu predykcji w tym modelu. Wykorzystamy w tym celu pętlę for (warto spróbować zrobić to zadanie za pomocą arytmetyki wektorowej!).

Wczytaj i obejrzyj dane biopsy z biblioteki MASS oraz przeczytaj ich dokumentację. Następnie wyrzuć obserwacje z brakującymi danymi (NA). Wykorzystaj funkcję na.omit. Sprawdź, czy w danych są cechy, których nie powinniśmy brać pod uwagę budujac model. Na podstawie pozostałych cech zbuduj model regresji logistycznej przewidujący złośliwość nowotworu.

Oceń błąd modelu za pomocą \(k\)-krotnej walidacji krzyżowej, z samodzielnie wybraną wartością \(k\).
Porównaj wyniki swojej metody z funkcją train z pakietu caret. Przykłady jej użycia dla regresji liniowej znajdziesz tutaj; wykorzystanie do regresji logistycznej jest bardzo zbliżone.

Uwaga. W tym zadaniu dla uproszczenia wykozystamy regresję logistyczną, czyli potraktujemy zmienne objaśniające jako ilościowe. Tak naprawdę nie jest to poprawne podejście, ponieważ cechy nie są mierzone, ocenione w uporządkowanej skali od 1 do 10. To oznacza, że są to zmienne rangowe, pośrednie pomiędzy zmiennymi jakościowymi a ilościowymi. Poprawne podejście do tych danych wymagało by wykorzystania innych technik, takich jak regresja rangowa, które jednak znacznie wykraczają poza program tych zajęć. Jak jednak zobaczymy, model logistyczny, mimo, że teoretycznie nie nadaje się do tych danych, daje bardzo dobre rezultaty. Ot, różnica między teorią a praktyką.

Rozwiązanie. Wczytujemy i podglądamy dane:

library(MASS)
?biopsy
dim(biopsy)
## [1] 699  11
head(biopsy)
##        ID V1 V2 V3 V4 V5 V6 V7 V8 V9     class
## 1 1000025  5  1  1  1  2  1  3  1  1    benign
## 2 1002945  5  4  4  5  7 10  3  2  1    benign
## 3 1015425  3  1  1  1  2  2  3  1  1    benign
## 4 1016277  6  8  8  1  3  4  3  7  1    benign
## 5 1017023  4  1  1  3  2  1  3  1  1    benign
## 6 1017122  8 10 10  8  7 10  9  7  1 malignant

Mamy 11 cech (kolumn) i 699 obserwacji (wierszy). Czyscimy dane:

biopsy <- na.omit(biopsy)
dim(biopsy)
## [1] 683  11

Po usunięciu wierszy zawierających brakujące wartości zostały nam 683 obserwacje. Sprawdźmy poziom bazowy zmiennej class:

levels(biopsy$class)
## [1] "benign"    "malignant"

Zmienna class ma dwa poziomy: ‘benign’ (nowotwór łagodny) i ‘malignant’ (nowotwór złośliwy). Poziomem bazowym jest ‘benign’.
Sprawdźmy zbalansowanie klas:

table(biopsy$class)
## 
##    benign malignant 
##       444       239

Prywatna opinia autora skryptu jest taka, że mamy w miarę dobre zbalansowanie. Na jeden nowotwór złośliwy przypadają dwa niezłośliwe. To oznacza, że klasyfikator przypisujący każdemu nowotworowi typ niezłośliwy może uzyskać dokładność na poziomie 65%, co nie jest dużo lepsze niż klasyfikator losowy. Możemy również oszacować jak zbalansowane zbiory testowe otrzymamy.
Załóżmy, że wybieramy losowo 60 obserwacji do zbioru testowego. Niech \(X\) będzie liczbą wylosowanych nowotworów złośliwych. Dla uproszczenia przyjmijmy na chwilę, że obserwacje do zbioru testowego wybieramy ze zwracaniem. Wowczas możemy przybliżyć rozkład \(X\) jako \(X \sim \text{Bin}(60, 239/(239+444)) = \text{Bin}(60, 0.35)\).
Wylosowana liczba nowotworów złożliwych na 95% mieści się w przedziale od qbinom(0.025, 60, 0.35) == 14 do qbinom(0.975, 60, 0.35)) == 28, co daje stosunkowo niewielkie odchylenie od oczekiwanej wartości równej 0.35*60 = 21: długość otrzymanego przedziału względem wartości oczekiwanej wynosi (28 - 14)/21 == 0.66.
Dla porównania, w przypadku idealnie zbalansowanych klas względna długość takiego przedziału wyniosłaby 0.533. Natomiast gdybyśmy mieli jedynie 10 nowotworów złośliwych i 673 nowotwory łagodne, to w zbiorze testowym otrzymalibyśmy od qbinom(0.025, 60, 0.015) == 0 do qbinom(0.975, 60, 0.015) == 3, i względna różnica względem wartości oczekiwanej wyniosłaby 3/0.9 = 3.33.

Ja osobiście zaczynam się martwić w momencie, gdy na jedną obserwację typu 1 przypada 10 obserwacji typu 2.
Ale to jest ocena subiektywna, i ktoś inny mógłby uznać że przy naszych danych już należy brać poprawkę na niezbalansowane klasy.

Budujemy wstępny model, żeby zobaczyć jak poszczególne zmienne mogą wpływać na typ nowotworu (krok opcjonalny). Oczywiście nie chcemy w modelu wykorzystać ID pacjenta, więc usuwamy je z formuły:

model <- glm(class ~ . - ID, biopsy, family=binomial)
summary(model)
## 
## Call:
## glm(formula = class ~ . - ID, family = binomial, data = biopsy)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -3.4841  -0.1153  -0.0619   0.0222   2.4698  
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -10.10394    1.17488  -8.600  < 2e-16 ***
## V1            0.53501    0.14202   3.767 0.000165 ***
## V2           -0.00628    0.20908  -0.030 0.976039    
## V3            0.32271    0.23060   1.399 0.161688    
## V4            0.33064    0.12345   2.678 0.007400 ** 
## V5            0.09663    0.15659   0.617 0.537159    
## V6            0.38303    0.09384   4.082 4.47e-05 ***
## V7            0.44719    0.17138   2.609 0.009073 ** 
## V8            0.21303    0.11287   1.887 0.059115 .  
## V9            0.53484    0.32877   1.627 0.103788    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 884.35  on 682  degrees of freedom
## Residual deviance: 102.89  on 673  degrees of freedom
## AIC: 122.89
## 
## Number of Fisher Scoring iterations: 8

W celu wstępnej oceny modelu sprawdźmy błąd treningowy. Jeśli model ma wysoki błąd treningowy, to od razu będziemy wiedzieć, że nie ma sensu sprawdzać go dalej - błąd testowy będzie na pewno wyższy. W pierwszym kroku robimy predykcję prawdopodobieństwa tego, że nowotwór jest złośliwy. Następnie przypisujemy klasę ‘malignant’ tam, gdzie to prawdopodobieństwo jest wyższe niż 1/2. Na poprzednich laboratoriach wykorzystaliśmy w tym celu funkcję ifelse, a tym razem zrobimy to w inny sposób:

malignant_probability <- predict(model, biopsy, 'response')
predicted_class <- rep('benign', nrow(biopsy))  
predicted_class[malignant_probability > 0.5] <- 'malignant'
mean(predicted_class != biopsy$class)
## [1] 0.03074671

Błąd treningowy jest bardzo niski, czyli nasz model ma potencjał. Należy jednak sprawdzić go dokładniej, i ocenić, jak dobrze radzi sobie na danych które nie były wykorzystane do treningu. W tym celu przeprowadzamy 10-krotną walidację krzyżową. Na początku, aby ułatwić sobie pisanie formuł, usuniemy zmienną ID z danych:

biopsy$ID <- NULL  # najprostszy sposób na usunięcie kolumny 

Aby podzielić nasze dane wykorzytamy funkcję cut, która zwróci wektor zawierający indeksy zbiorów testowych do których należą obserwacje. Następnie permutujemy otrzymany wektor aby zapewnić losowość przypisania obserwacji do zbioru testowego. Wykorzystamy w tym celu funkcję sample, która domyślnie zwraca permutację.

podzial <- cut(1:nrow(biopsy), 10, labels = F)
podzial <- sample(podzial) 

Dla każdego zbioru testowego trenujemy nowy model. W tym celu wykorzystamy pętlę for. Jej składnia w R to for(x in X) {...}, gdzie X to dowolny wektor, x to zmienna przymująca kolejne wartości z tego wektora, a w nawiasach klamrowych znajdują się komendy wykorzystujące zmienną x. Więcej na temat tej pętli i przykłady jej zastosowania można przeczytać w rozdziale 9 manualu R.

bledy.kfold <- numeric(10)  # tu zapiszemy bledy testowe
bledy.treningowe <- numeric(10)  # a tu bledy treningowe
for (i in 1:10) {
  treningowe <- which(podzial != i)  # wektor zawierający indeksy obserwacji ze zbioru treningowego
  model <- glm(class ~ ., biopsy, family=binomial, subset = treningowe)  # model zbudowany na danych treningowych
  
  predykcja <- predict(model, biopsy, 'response')  # robimy predykcje na całych danych
  predykcja <- ifelse(predykcja > 0.5, 'malignant', 'benign')  # zamieniamy prawdopodobieństwa na klasy
  
  predykcja_treningowa <- predykcja[treningowe]  # wybieramy predykcje odpowiadające zbiorowi treningowemu
  predykcja_testowa <- predykcja[-treningowe]
  
  # poniżej łączymy dwa typy indeksowania tabeli biopsy: indeksowanie wektorem całkowitoliczbowym oraz przez nazwę kolumny
  bledy.kfold[i] <- mean(predykcja_testowa != biopsy[-treningowe, 'class'])  
  bledy.treningowe[i] <- mean(predykcja_treningowa != biopsy[treningowe, 'class'])  
}

Przedstawmy błędy na wykresie:

library(ggplot2)
## Registered S3 methods overwritten by 'ggplot2':
##   method         from 
##   [.quosures     rlang
##   c.quosures     rlang
##   print.quosures rlang
dane_do_wykresu <- data.frame('Fold' = as.integer(rep(1:10, 2)),
                              'Blad' = c(bledy.kfold, bledy.treningowe),
                              'Typ' = rep(c('Testowy', 'Treningowy'), each=10))
ggplot(dane_do_wykresu,
       aes(x=Fold, y=Blad, color=Typ)) + 
  ggtitle('Porównanie błędu treningowego i testowego', 
          subtitle = paste0('Średni błąd testowy = ', round(mean(bledy.kfold), 3), ' treningowy = ', round(mean(bledy.treningowe), 3))) + 
  geom_line() + geom_point() + scale_x_continuous(breaks = c(1,3,5,7,9)) +  theme_minimal() 

Średni błąd testowy, równy 0.03363, jest jedynie nieznacznie wyższy od średniego błędu treningowego równego 0.02944. Nasz klasyfikator ma bardzo wysoką skuteczność.

Na wykresie widoczne jest dość nietypowe zjawisko: w niektórych powtórzeniach błąd testowy był nawet niższy niż błąd treningowy. Jest to spowodowane losowym przypisaniem obserwacji do zbioru testowego, oraz tym, że zbiory testowe są mało liczne. W niektórych powtórzeniach nasz klasyfikator myli się tylko raz - natomiast w tych, w których myli się trzy razy, na ogół błąd testowy przekracza już błąd treningowy. Takie zjawisko jest bardzo rzadkie, i świadczy o tym, że 5-krotna walidacja krzyżowa mogłaby być w tym przypadku lepszym wyborem - większe zbiory testowe zapewniłyby mniejszą wariancję estymatora błędu testowego, ponieważ pomylenie się w klasyfikacji pojedynczej obserwacji miałoby mniejszy wpływ na średni błąd.

Porównajmy teraz nasze wyniki z implementacją walidacji krzyżowej z biblioteki caret. W pierwszym kroku tworzymy obiekt typu trainControl, który określa, w jaki sposób chcemy ocenić nasz model.

library(caret)
## Loading required package: lattice
train_control <- trainControl(method='cv', number=10)

Teraz wywołujemy funkcję train, która wytrenuje nasz model i przetestuje go za pomocą walidacji krzyżowej.

kfold_train <- train(class~., data=biopsy, 
                     method='glm', family=binomial, trControl=train_control)

Podsumowanie procedury trenowania wyświetlamy następująco:

print(kfold_train)
## Generalized Linear Model 
## 
## 683 samples
##   9 predictors
##   2 classes: 'benign', 'malignant' 
## 
## No pre-processing
## Resampling: Cross-Validated (10 fold) 
## Summary of sample sizes: 615, 615, 615, 614, 615, 614, ... 
## Resampling results:
## 
##   Accuracy  Kappa    
##   0.966283  0.9256307

Dokładność klasyfikatora wyniosła 0.966283, a więc jest bardzo bliska tej, którą otrzymaliśmy przeprowadzając walidację krzyżową ręcznie. Możemy ją również wyświetlić wywołując komendę kfold_train$results$Accuracy. Zmienna kfold_train$results zawiera również odchylenie standardowe estymatora dokładności klasyfikatora. Wyniki poszczególnych iteracji otrzymamy pisząc kfold_train$resample, a wytrenowany model pisząc kfold_train$finalModel. Powinien on być taki sam, jak model który stworzyliśmy na samym początku na pełnych danych.

Wybór modelu

Model, który opracowaliśmy w poprzednim zadaniu, sprawdza się bardzo dobrze.
Wykorzystuje jednak wiele zmiennych nieistotnych statystycznie, które mogą nie tylko być niepotrzebne, ale nawet zwiększać błąd testowy ze względu na przeuczenie (ponieważ model może wykorzystywać fałszywe informacje).
Warto to poprawić z dwóch powodów:

Oczywiście powyższe rozważania są czysto teoretyczne, ponieważ nasz klasyfikator nie będzie wykorzystywany w diagnostyce. Niedługo jednak część z Was być może będzie pracować nad zastosowaniami uczenia maszynowego w medycynie. Warto wiedzieć, jak to robić tak, by nie zabijać ludzi.

Do wyboru modelu (czyli wyboru w pewnym sensie optymalnego podzbioru zmiennych objaśniających) wykorzystamy następujące metody, gdzie \(k\) to liczba wybranych predyktorów, \(n\) to liczba obserwacji, a \(L\) to wiarygodność modelu:

Mimo, że te techniki są niemal identyczne, kryteria AIC i BIC mają inne własności. BIC sprawdza się lepiej w analizie eksploracyjnej na (bardzo) dużych zbiorach danych, gdy szukamy “prawdziwego modelu” - czyli chcemy odkryć rzeczywiste prawa rządzące zmiennością wyjaśnianej zmiennej. AIC natomiast sprawdza się lepiej w zadaniu predykcji, gdy chcemy znaleźć model który da mały błąd średniokwadratowy na nowych danych. Więcej informacji na ten temat znajduje się w artykułach podlinkowanych w powyższych punktach.

Do znalezienia optymalnego modelu wykorzystamy strategię krokową, czyli stopniowe dodawanie i usuwanie zmiennych. Dla kryteriów AIC i BIC jest ona zaimplementowana w funkcji stepAIC z pakietu MASS. Pozwala ona na wyszukiwanie w przód (czyli stopniowe dodawanie predyktorów), wstecz (czyli stopniowe usuwanie predyktorów w pełnym modelu) oraz w obie strony.
Zauważmy, że jeśli mamy więcej cech niż obserwacji, to nie możemy stosować przeszukiwania wstecz, ponieważ nie da się w takiej sytuacji zbudować pełnego modelu. Możemy natomiast stosować stratgię w przód i strategię naprzemienną, zaczynającą od pustego modelu.

Krokowy wybór modelu, choć bardzo popularny, jest nieco kontrowersyjną techniką. Pozwala na znaczne przyspieszenie obliczeń w porównaniu ze sprawdzeniem wszystkich możliwych podzbiorów zmiennych objaśniających, ale nie gwarantuje znalezienia optymalnego modelu.
Żeby dowiedzieć się więcej na ten temat, warto zapoznać się z dyskusją w tym wątku oraz z uwagami w tym artykule. Na następnych zajęciach zapoznamy się z alternatywnym sposobem wyboru modelu w regresji liniowej. Niestety w przypadku innych modeli często nie mamy dobrej alternatywy.

Innym sposobem wyboru modelu jest minimalizacja błędu testowego estymowanego przez walidację krzyżową. W przypadku regresji liniowej można ją zrealizować korzystając z funkcji train oraz trainControl z pakietu caret. Przykłady ich użycia do wyboru modelu można znaleźć tutaj. O ile wiem, pakiet caret niestety nie pozwala na wykorzystanie tej metody dla regresji logistycznej.

Zadanie 2. Utwórz klasyfikator typu nowotworu w oparciu o regresję logistyczną i dane biopsy. Przeprowadź wybór modelu poprzez minimalizację AIC korzystając z funkcji stepAIC. Następnie, korzystając z tej samej funkcji, przeprowadź wybór modelu poprzez minimalizację BIC. Porównaj te trzy modele poprzez k-krotną walidację krzyżową, z samodzielnie wybraną wartością parametru k, korzystając z pakietu caret.

Metoda Bootstrap

Zadanie 3. Włodzimierz Bielski zajmuje się produkcją niebieskich landrynek, które sprzedaje Meksykaninowi Tuco. Tuco zauważył, że zawartość cukru w dostarczanych landrynkach waha się. Chce oszacować średnią zawartość cukru w landrynkach uzyskiwanych metodą Włodzimierza.

Zawartość cukru w landrynkach znajduje się w pliku walter.csv. Załaduj plik i wyestymuj średnią zawartość cukru w landrynkach.

Tuco zna się na statystyce i martwi się, czy otrzymany estymator mu wystarcza. Chciałby zbadać jego wariancję, żeby dowiedzieć się, czy na innych próbach dostanie podobne wyniki. Pomóż mu, korzystając z metody bootstrap:

Wskazówka: Funkcję z pierwszego kroku możesz utworzyć jako srednia_podzbioru <- function(dane, indeksy) mean(dane[indeksy]).

Zadanie 4. Włodzimierz przysięga, że średnia zawartość cukru w jego landrynkach jest wyższa niż 75%. Jeśli jest ona mniejsza, Tuco będzie się na niego bardzo gniewał. Co prawda Włodzimierz nie rozumie, co to jest rozkład prawdopodobieństwa (a tym bardziej test statystyczny), więc nie wie jak sprawdzić czy jego hipoteza jest prawdziwa, ale zarzeka się, że “to widać”. Pomóż Włodzimierzowi, korzystając z metody bootstrap:

Więcej ciekawych zastosowań metody bootstrap można znaleźć w Google, szukając np. “bootstrap hypothesis testing”, “bootstrap confidence intervals” itp.

Zadania dodatkowe (nieobowiązkowe).

Zadanie 5. W tym zadaniu zautomatyzujemy procedurę walidacji krzyżowej w modelu liniowym, umożliwiając wykonanie jej dla innych danych i parametrów. Napisz funkcję kfold(formula, dane, y, k), która przyjmuje jako argumenty:

Funkcja powinna zwracać średni testowy błąd kwadratowy modelu liniowego. Domyślnie funkcja powinna wykonywać walidację leave-one-out, gdzie parametr k jest równy liczbie obserwacji.

Wektor wartości zmiennej predykowanej podajemy poza ramką danych dla ułatwienia - nie trzeba będzie go wyciągać do liczenia błędu testowego.

Na danych biopsy zbuduj model liniowy przewidujący zmienną V2 w zależności od pozostałych zmiennych od V1 do V9. Porównaj działanie Twojej funkcji i funkcji train z pakietu caret.