Na dzisiejszym laboratorium zajmiemy się technikami LASSO oraz regresji grzbietowej.
Ponieważ, nie licząc kilku gotowych funkcji, zadania wykorzystują wyłącznie to, co już poznaliśmy, i są dość szczegółowo opisane, to nie zamieszczam w tym skrypcie zadania z przykładowym rozwiązaniem.
Zaczniemy natomiast od czytanki dotyczącej typów zmiennych - gdzie tym razem przez “zmienną” rozumiem rodzaj danych, z jakim mamy do czynienia w tabelce, a nie pojęcie programistyczne.

Rodzaje zmiennych w statystyce

Na poprzednich laboratoriach spotkaliśmy się z pojęciem “zmienne ilościowe” oraz “zmienne jakościowe”. Jaka jest właściwie pomiędzy nimi różnica?

Zmienne ilościowe oznaczają na ogół wynik jakiegoś pomiaru lub zliczenia, podczas gdy zmienne jakościowe (kategoryczne) oznaczają pewną, na ogół niemierzalną, cechę jakiegoś obiektu. W przypadku tych pierwszych mają sens operacje arytmetyczne, a w przypadku tych drugich - nie. Można mieć dwa razy więcej podgrzybków, ale nie można być dwa razy bardziej podgrzybkiem. Zmienne jakościowe, w przeciwieństwie do ilościowych, nie mają również określonego porządku. Napis “2 kg podgrzybków > 1 kg podgrzybków” ma sens, ale napis “podgrzybek > maślak” nie (w każdym razie nie obiektywny). Ogólnie rzecz biorąc, poszczególne wartości zmiennej jakościowej nie mają ze sobą żadnego związku.

Operacje arytmetyczne oraz istnienie porządku to główne sposoby na określenie jakiego typu jest zmienna. Nie ma znaczenia natomiast to, w jaki sposób jest kodowana - jeśli w bazie danych gatunki grzybów zostaną oznaczone liczbami naturalnymi to wcale nie oznacza, że gatunek nagle stanie się zmienną ilościową. Dlatego wczytując dane w R należy zawsze sprawdzić, czy zmienne jakościowe są kodowane jako factory, a zmienne ilościowe jako wektory numeryczne.

Pomylenie typu zmiennej ma, oczywiście, wiele potencjalnych negatywnych konsekwencji.
Potraktowanie zmiennej jakościowej jako zmiennej ilościowej sprawia, że błędnie przypisujemy strukturę geometryczną poszczególnym wartościom (czyli arbitralnie ustawiamy je na prostej rzeczywistej). W praktyce zdarza się, że to nawet działa, ale równie dobrze może sprawić, że model straci jakikolwiek sens.
Potraktowanie zmiennej ilościowej jako zmiennej jakościowej sprawia, że tracimy informację wynikającą z jej geometrii (czyli tego, że wiemy, co to znaczy “dwa razy więcej”). Sprawia to również, że nasz model ma dużo więcej parametrów (ponieważ w przypadku zmiennych jakościowych każda wartość jest kodowana jako osobna zmienna), więc ich estymacja jest dużo mniej dokładna, a do tego ryzykujemy przeuczenie. A w końcu, kiedy mamy wykonać predykcję, bardzo możliwe że nowa wartość zmiennej będzie inna niż te, które widzieliśmy w danych. Warto zdawać sobie sprawę z tego, że jeśli mamy do czynienia z rzeczywistą aparaturą pomiarową, to praktycznie nigdy nie otrzymamy dwóch dokładnie takich samych wyników.

Zmienne kategoryczne i ilościowe nie wyczerpują wszystkich możliwości. Na ostatnich laboratoriach spotkaliśmy się na przykład ze zmiennymi rangowymi - posiadającymi określony porządek, ale brak struktury geometrycznej (okazało się natomiast, że w tym przypadku mogliśmy je potraktować jako zmienne ilościowe bez znacznej utraty skuteczności modelu). Poniższa tabela podsumowuje kilka podstawowych rodzajów zmiennych zgodnie z rosnącą ilością informacji na temat związku pomiędzy dwiema wartościami. Wymieniam ponadto przykłady technik statystycznych które można (lub nie można) stosować do danego typu zmiennej - z niektórymi nie spotkamy się na tym przedmiocie, ale część z Was prawdopodobnie przerabiała je na innych przedmiotach.

Typ zmiennej Przykład Co działa Co nie działa
Kategoryczna
Categorical
Płeć
Gatunek
Test niezależności \(\chi^2\)
Test dwumianowy
ANOVA (w pewnym sensie)
Mediana
Średnia
Test Wilcoxona
Rangowa
Ordinal
Poziom wykształcenia
Skala Likerta
Mediana
Test Wilcoxona
Regresja rangowa
Średnia
Test t Studenta
Regresja liniowa
Interwałowa
(właśnie dlatego daję odnośniki do angielskiej Wikipedii…)
Interval
Stopnie Celsjusza
Data
Średnia
Mediana
Odchylenie standardowe
Cokolwiek co wymaga mnożenia lub dzielenia
Regresja liniowa
(teoretycznie - bo w praktyce działa doskonale)
Ilorazowa
Ratio
Liczba (np. podgrzybków)
Masa
Czas od wyróżnionej chwili 0
Wszystko

Istnieją również inne typy zmiennych, a nawet zupełnie inne sposoby ich klasyfikacji. Nieco więcej można o tym przeczytać tutaj oraz tutaj.
Na tych zajęciach przez zmienne jakościowe rozumiemy zmienne kategoryczne oraz rangowe, czyli wszystkie te, na których nie są określone operacje arytmetyczne. Przez zmienne ilościowe rozumiemy natomiast zarówno zmienne ilorazowe, jak i interwałowe - czyli wszystko, na czym jest określona jakakolwiek operacja arytmetyczna.

Wybór modelu metodą LASSO

Zadanie 1. W tym zadaniu porównamy techniki LASSO oraz regresji grzbietowej z kryteriami informacyjnymi (AIC, BIC) w przypadku regresji liniowej.
Naszym celem jest zbadanie, jak dobrze każda z tych technik oddaje prawdziwy model oraz jak dobrze nadaje się do predykcji. Metoda LASSO została omówiona na Wykładzie 10. Więcej informacji na jej temat można przeczytać m.in. tutaj, a obejrzeć tutaj (ten tutorial dotyczy regresji logistycznej, ale myślę, że na tym etapie domyślicie się już, jak to się ma do regresji liniowej). Jedną z największych zalet metody LASSO, poza tym, że umożliwia całkiem skuteczny wybór modelu, jest to, że można ją stosować do danych w których mamy więcej cech niż obserwacji.

W pierwszym kroku stworzymy zbiór danych, których użyjemy do porównania procedur. Jako zmienne objaśniające wykorzystamy zmienne ilościowe z danych Boston z pakietu MASS.
Żeby wiedzieć, jak dobrze działa dana procedura, musimy wiedzieć jak wygląda prawdziwy model. Dlatego wektor parametrów \(\beta\) i zmienną objaśnianą \(Y\) musimy stworzyć sami. Można to zrobić na wiele różnych sposobów, na przykład wybrać wartości parametrów samodzielnie. Ja proponuję następujący wektor, który wylosowałem z rozkładu \(Unif(-3, 3)\), a następnie wyzerowałem co trzecią obserwację:

B <- c(crim=1.52274943515658, zn=-2.79081122064963, indus=0, 
       nox=1.4332089680247, rm=1.03212527744472, age=0, 
       dis=2.75923800934106, tax=2.43884606240317, ptratio=0, 
       black=-0.371096117421985, lstat=1.09251574007794, medv=0)

Następnie trzeba obliczyć \(Y = X_i \beta + \varepsilon_i\), gdzie \(X_i\) to i-ty wiersz tabeli Boston, a \(\varepsilon\) to zmienna o standardowym rozkładzie normalnym. Najłatwiej jest obliczyć cały wektor \(X\beta\) korzystając z mnożenia macierzowego: as.matrix(X) %*% B, a następnie dodać do niego wektor niezależnych zmiennych gaussowskich otrzymany funkcją rnorm().

Szczegółowy opis proponowanej procedury utworzenia zmiennej objaśnianej wygląda następująco:

W następnych punktach sprawdzimy, jak poszczególne metody wyboru modelu radzą sobie z przewidywaniem \(Y_i\) oraz \(\beta_j\) na zbiorze testowym. W tym przypadku wybór modelu oznacza utworzenie modelu liniowego wykorzystującego wyłącznie te parametry, które są niezerowe - przyjmujemy, że one stanowią prawdziwy model, a parametry zerowe do tego prawdziwego modelu nie należą.

Żeby ocenić poszczególne metody, w pierwszym kroku podziel dane na zbiór treningowy rozmiaru 60 oraz zbiór testowy złożony z pozostałych obserwacji. Jeśli masz ochotę zrobić to porządniej, to możesz przeprowadzić walidację krzyżową, ale w tym zadaniu nie jest to obowiązkowe.

Pierwsze dwie metody pojawiły się na wcześniejszych zajęciach:

1. Zwykły model liniowy. Stwórz model liniowy za pomocą lm. Czy model poprawnie wykrył, od których zmiennych zależy \(Y\)? Porównaj wartości estymatorów \(\hat{\beta}_j\) z prawdziwymi wartościami \(\beta_j\). Zbadaj błąd średniokwadratowy predykcji na zbiorze testowym (pamiętaj, że znasz prawdziwe wartości zmiennej \(Y\) na zbiorze testowym). Porównaj go z minimalnym błędem jaki możemy otrzymać, czyli \(Y_i - X_i \beta\).

2. Kryteria informacyjne. Wybierz model za pomocą funkcji stepAIC. Czy kryterium informacyjne poprawiło model i predykcję? Czy jesteśmy blisko minimalnego błędu? Na ile uzyskany model jest poprawny? Poprawność możemy ocenić na przykład poprzez sprawdzenie, ile spośród wyzerowanych parametrów jest w modelu, a ile niezerowych do niego nie trafiło. Porównaj kryterium AIC z kryterium BIC. To ostatnie, zgodnie z teorią, prowadzi do wyboru poprawnego modelu, o ile możemy porównać ze sobą wszystkie potencjalne modele (zamiast wykorzystywać strategię krokową) oraz posiadamy nieskończenie wiele obserwacji.

Następne metoda bazuje na wprowadzeniu kary za rozmiar modelu bezpośrednio w optymalizowaną funkcję. W regresji LASSO szukamy estymatora \(\hat{\beta}_{LASSO} = \text{arg min}_{\beta} ||Y- X\beta ||^2 + \lambda |\beta|\). Optymalną wartość \(\lambda\) wybiera się za pomocą walidacji krzyżowej - tworzy się tak zwaną siatkę parametrów, czyli wektor potencjalnych wartości parametru \(\lambda\), a następnie dla każdej wartości ocenia się błąd testowy. Oczywiście nie będziemy implementować tej procedury samodzielnie, tylko wykorzystamy gotowe funkcje.

3. Regresja LASSO. Dopasuj model LASSO do zbioru treningowego za pomocą funkcji cv.glmnet z biblioteki glmnet. Funkcja ta sama wybierze pewną liczbę wartości \(\lambda\) do przetestowania, i dla każdej z nich przeprowadzi walidację krzyżową aby ocenić błąd predykcji. Ponieważ nasze dane są już wystandaryzowane, ustaw w funkcji cv.glmnet atrybut standardize=F; ponadto ustaw thresh=1e-16 aby otrzymać dokładniejszą estymację. Porównaj błąd na zbiorze testowym wybranym na początku tego zadania. Może być on większy niż w przypadku innych technik, ponieważ metoda LASSO może wprowadzić spore obciążenie do estymowanych parametrów - w następnym zadaniu zobaczymy, jak to poprawić. Uwaga: funkcja cv.glmnet i inne funkcje z tego pakietu wymagają na wejściu macierzy, a nie ramki danych - należy wykorzystać funkcję as.matrix().

Kilka wskazówek:

  1. Wynik funkcji cv.glmnet(), zapisany w zmiennej simple.cv, możemy przedstawić na następującym wykresie. Przedstawia on błąd średniokwadratowy estymowany za pomocą 10-krotnej walidacji krzyżowej w zależności od logarytmu parametru \(\lambda\). Pionowa kreska po lewej stronie zaznacza wartość dla której błąd jest minimalny. Możemy również wpisać w konsolę simple.cv, aby uzyskać informacje na temat optymalnej wartości \(\lambda\) oraz optymalnego błędu testowego.
plot(simple.cv)

  1. Wartości parametrów dla optymalnego modelu otrzymamy funkcją coef. W tym przypadku kropka oznacza, że estymator parametru wynosi zero.
coef(simple.cv)
## 13 x 1 sparse Matrix of class "dgCMatrix"
##                        1
## (Intercept) -0.043513055
## crim         1.267752539
## zn          -2.619723970
## indus        0.075594732
## nox          1.307328703
## rm           0.857750048
## age          .          
## dis          2.444809390
## tax          2.358589551
## ptratio      0.004343864
## black       -0.393510182
## lstat        1.104596738
## medv        -0.050041554
  1. Predykcję wykonujemy, jak zwykle, korzystając z funkcji predict(), która automatycznie obliczy wektor \(X \hat{\beta}\) dla parametrów \(\hat{\beta}\) z powyższej tabeli. Różnica w stosunku do predykcji które jak dotąd wykonywaliśmy jest taka, że w tym przypadku funkcji predict() trzeba podać nowe obserwacje w macierzy, a nie ramce danych.
Y_lasso <- predict(simple.cv, X.test)

Istnieje jeszcze jedna metoda wizualizacji modelu LASSO - jest to tak zwana ścieżka LASSO, którą możemy uzyskać w następujący sposób:

# Uwaga, tutaj X.train to macierz, a Y.train to wektor numeryczny! 
# Biblioteka glmnet nie lubi zmiennych typu data.frame!
simple.lasso <- glmnet(X.train, Y.train)
plot(simple.lasso)

Wykres przedstawia wartości poszczególnych współczynnikow \(\beta\) w zależności od rozmiaru modelu mierzonego normą L1 (sumą wartości bezwzględnych współczynników). Jest ona tym większa, im mniejszy jest parametr \(\lambda\).
Warto zwrócić uwagę, czy w środku ścieżki pojawiają się współczynniki, które następnie ponownie wypadają z modelu.
Jest naturalne zjawisko w przypadku LASSO.

Zadanie 2. Pełna procedura LASSO. W poprzednim zadaniu zastosowaliśmy metodę LASSO w bardzo prosty sposób. Ma to kilka wad. Po pierwsze, wiadomo, że LASSO dość szybko (tzn. dla stosunkowo małych wartości kar \(\lambda\)) wprowadza do modelu zmienne które tak naprawdę do niego nie należą. Po drugie, estymator \(\hat{\beta}\) w metodzie LASSO jest obciążony (co jest oczywiste, jeśli się pomyśli, że regresja liniowa daje estymator nieobciążony, a my ją zaburzamy dodając karę za rozmiar modelu). To oznacza, że estymatory są zaburzone, i to tym bardziej, im wyższa jest wartość kary \(\lambda\). To z kolei wpływa na predykcję.

Nieco bardziej skomplikowana procedura pozwala w wielu przypadkach na znaczne poprawienie wyników. Składa się z trzech etapów:

Przeprowadzimy je w następujący, nieco uproszczony sposób:

  1. Utwórz zbiór walidacyjny, do którego przeniesiesz 20 obserwacji ze zbioru treningowego (w zbiorze treningowym powinno pozostać 40 obserwacji).
  2. Wybierz optymalny parametr \(\lambda\) na zbiorze treningowym za pomocą walidacji leave-one-out cv.glmnet.
  3. Utwórz wektor wartości progowych, poniżej których będziesz zerować parametry \(\beta\), na przykład za pomocą seq(0., 1, length.out=1000), oraz pusty wektor takiej samej długości co wektor progów, w którym zapiszesz błędy walidacyjne (funkcja numeric()).
  4. Dla każdej wartości progowej stwórz wektor progowanych parametrów. Następnie korzystając z nich oblicz predykcję na zbiorze walidacyjnym i porównaj z prawdziwymi wartościami zmiennej \(Y\) dla obserwacji z tego zbioru. Zapisz błąd średniokwadratowy w wektorze stworzonym w punkcie 3. Rozwiązując ten punkt możesz wykorzystać pętlę for.
  5. Znajdź próg minimalizujący błąd średniokwadratowy na zbiorze walidacyjnym. Następnie wyproguj parametry i wybierz te, które są niezerowe.
  6. Na pełnym zbiorze treningowym (rozmiaru 60) utwórz model liniowy, który wyjaśni zmienną \(Y\) za pomocą tych predyktorów, które odpowiadają niezerowym parametrom z punktu 5. Porównaj estymatory otrzymane w ten sposób z prawdziwymi wartościami \(\beta\).

W praktyce zarówno optymalną wartość \(\lambda\), jak i optymalny próg do wyzerowania współczynników, powinno wybierać się metodą walidacji krzyżowej na trzech zbiorach: treningowym, walidacyjnym i testowym. W tej metodzie dzielimy dane na 10 kawałków, tak jak w przypadku zwykłej walidacji krzyżowej. Następnie dla każdej pary parametrów (\(\lambda\) oraz progu) bierzemy jeden kawałek jako zbiór testowy, inny kawałek jako zbiór walidacyjny, a pozostałe obserwacje jako zbiór treningowy. Oceniamy błąd testowy i powtarzamy dla każdego możliwego wyboru zbioru walidacyjnego i testowego, a następnie dla każdej pary parametrów.
Ma to oczywiście bardzo duży koszt obliczeniowy, dlatego w tym zadaniu zbiór walidacyjny wybraliśmy tylko raz.

Skorelowane zmienne

Jak wiadomo, regresja liniowa nie radzi sobie zbyt dobrze z silnie skorelowanymi zmiennymi. Można to wyjaśnić na dwa sposoby. Formalny dowód jest taki, że odwracanie macierzy \(X^TX\) jest wówczas źle uwarunkowane numerycznie. Nieformalny powód jest taki, że współczynniki regresji liniowej mierzą wpływ danej zmiennej objaśniającej na zmienną objaśnianą, gdy wartości pozostałych zmiennych są ustalone. Wysoka korelacja oznacza natomiast, że dwie zmienne przyjmują duże wartości jednocześnie - nie możemy zwiększyć jednej, nie zwiększając drugiej. Im większa korelacja, tym mniejszy sens pojedynczych współczynników. Jeszcze inne wyjaśnienie można przeczytać tutaj.

Korelacja między zmiennymi utrudnia estymację parametrów oraz interpretację modelu. Z tego powodu opracowano wiele metod radzenia sobie z nią, między innymi regresję grzbietową (ang. ridge regression). W regresji grzbietowej szukamy \(\hat{\beta}_{RIDGE} = \text{arg min}_{\beta} || Y - X\beta ||^2 + \lambda || \beta ||^2\); Podobnie jak w przypadku LASSO, parametr \(\lambda\) wybieramy poprzez sprawdzenie zestawu potencjalnych parametrów za pomocą walidacji krzyżowej. Pomimo łudząco podobnej postaci, regresja grzbietowa ma zupełnie inne zastosowanie niz LASSO:

Skorelowane zmienne to nie jedyne zastosowanie regresji grzbietowej. Służy ona również do poprawy wyników w sytuacji, gdy głównym problemem jest wariancja estymatorów - pozwala na ograniczenie tej wariancji kosztem wprowadzenia pewnego obciążenia.
Więcej na temat różnych aspektów tej techniki, wraz z przykładami w R, można przeczytać tutaj, a nieco bardziej matematyczne podejście do tematu jest przedstawione tutaj.
Jeszcze więcej informacji, wraz z przykładami w języku Python, znajduje się tutaj.
Inne podejścia do skorelowanych zmiennych są opisane m.in. tutaj.

Zadanie 3. W tym zadaniu zilustrujemy działanie regresji grzbietowej na wysoce skorelowanych danych.

Regresja grzbietowa jest bardzo popularnym narzędziem do radzenia sobie ze skorelowanymi zmiennymi objaśniającymi i ma wielu zwolenników. Podobnie jednak jak w przypadku krokowego wyboru modelu, niektórzy statystycy krytykują tę technikę. Jest do tego kilka powodów, między innymi:

Alternatywnym podejściem radzenia sobie z bardzo silnie skorelowanymi zmiennymi (\(\rho>0.99\)) jest potraktowanie ich jako jednej zmiennej. Podejście to można zrealizować i porównać z regresją grzbietową następująco: