Na ostatnich dwóch laboratoriach zajmiemy się eksploracyjną analizą danych, czyli taką, w której chcemy trochę lepiej zobaczyć jak “wyglądają” skomplikowane dane.
Jedną z najważniejszych technik eksploracyjnej analizy danych (i zastosowań statystyki w ogóle) jest analiza głównych składowych (ang. principal component analysis, PCA). Technika ta zostanie (lub, w zależności od grupy, już została) omówiona szczegółowo na Wykładzie 13. W skrócie: mając wycentrowane dane \(X\), czyli takie, dla których w każdej kolumnie obserwacje sumują się do zera, szukamy unormowanej kombinacji liniowej kolumn, czyli \(X\phi = \phi_1 X[:, 1] + \phi_2 X[:, 2] + \dots + \phi_k X[:, k]\) gdzie \(|| \phi ||_2 = 1\), która jak najlepiej odwzorowuje dane. Jakość odwzorowania mierzymy następująco: dla danego wektora \(\phi\) rozpatrujemy prostą \(t\phi, t\in \mathbb{R}\), i patrzymy, jak wiele obserwacji znajduje się z grubsza na tej prostej. Liczbę obserwacji “z grubsza na prostej” mierzymy rzutując wszystkie obserwacje z \(X\) na tę prostą (czyli licząc po prostu iloczyn macierzowy \(X\phi\)) i licząc wariancję zrzutowanych obserwacji - to określa jak szeroki przekrój danych zapewnia wektor \(\phi\). Po znalezieniu optymalnego wektora \(\phi\) szukamy następnego, w dokładnie taki sam sposób, z tą różnicą, że rozpatrujemy wyłącznie wektory prostopadłe do \(\phi\). Otrzymujemy w ten sposób wektory \(\phi^{(1)}, \phi^{(2)}, \dots, \phi^{(k)}\), które kolejno wyjaśniają coraz mniej wariancji danych (tutaj liczby w nawiasach to po prostu indeksy, nie mylić z potęgowaniem).

W tym skrypcie przyjmujemy następującą terminologię: wektor \(X\phi^{(i)}\) nazywamy i-tą składową główną danych \(X\), wektor \(\phi^{(i)}\) nazywamy i-tym kierunkiem głównym, a wektor \(\phi^{(i)}\) przemnożony przez odchylenie standardowe odpowiadającej mu składowej głównej, czyli \(\text{sd}(X\phi^{(i)})\phi^{(i)}\), nazywamy ładowaniem i-tej składowej głównej (część ludzi mówi też o i-tych loadingsach…).
Niestety, we współczesnej literaturze ta terminologia nie jest do końca jednoznaczna - niektórzy autorzy rozumieją co innego przez składowe główne i ładowania, a niektórzy nazywają składowymi głównymi zarówno wektor \(\phi^{(i)}\) jak i rzut danych na ten wektor, czyli \(X\phi^{(i)}\).
W praktyce trzeba po prostu rozumieć, jaki wzór dany autor ma na myśli.

Analiza składowych głównych. Strzałki oznaczają ładowania składowych. Źródło

Analiza składowych głównych. Strzałki oznaczają ładowania składowych. Źródło

Istnieje również alternatywny, ale równoważny opis tej techniki: szukamy \(k\)-wymiarowej elipsoidy, która jak najlepiej odwzorowuje dane (a więc zamiast przeprowadzania przez chmurę punktów prostej, jak w regresji liniowej, przybliżamy te punkty elipsoidą). Ładowania pokrywają się z półosiami takiej elipsoidy.

Do oblilczenia składowych głównych w praktyce wykorzystujemy rozkład macierzy na wartości szczególne lub osobliwe (ang. singular value decomposition, SVD), który - z punktu widzenia matematyki stosowanej - jest prawdopodobnie najważniejszym wynikiem algebry liniowej. Mając macierz \(X\) wymiaru \(n \times k\), rozkład SVD polega na przedstawieniu jej jako \(X = U\Sigma V^T\), gdzie \(U\) jest macierzą ortonormalną (czyli taką, że \(U^TU = I\)) wymiaru \(n \times n\), \(V\) jest macierzą ortonormalną wymiaru \(k \times k\), a \(\Sigma\) jest macierzą \(n \times k\) z niezerowymi wyrazami wyłącznie na głównej diagonali (zwanymi wartościami szczególnymi), uporządkowanymi niemalejąco.
Przekształcając równanie na SVD otrzymujemy \(XV = U\Sigma\). Okazuje się, że macierz \(V\) zawiera w kolumnach kierunki główne, wartości na diagonali macierzy \(\Sigma\) są proporcjonalne do odchyleń standardowych rzutów \(X\) na kierunki główne, a \(U\Sigma\) zawiera składowe główne. Macierz \(V\Sigma^T\) zawiera w kolumnach ładowania.

Podobnie jak w przypadku regresji liniowej, przeprowadzając analizę składowych głównych warto zadbać o to, by wszystkie zmienne były w podobnej skali, tak by nie porównywać milimetrów z kilometrami. Najczęściej osiąga się to poprzez dzielenie kolumn przez ich odchylenie standardowe, choć istnieją również inne podejścia, o których można przeczytać np. tutaj.

Więcej na temat idei stojącej za PCA można przeczytać tutaj, a przykłady różnych zastosowań tutaj.
A jeśli ktoś by chciał wiedzieć naprawdę dużo, to może na przykład przeczytać tę oto 500-stronową książkę poświęconą wyłącznie PCA.
Więcej na temat rozkładu SVD od strony algebry liniowej można przeczytać w tym artykule oraz tej pracy.

Uwaga na boku: nowoczesne algorytmy do regresji liniowej również są oparte na SVD - w ramach ćwiczenia polecam sprawdzić, w jaki sposób rozkład SVD macierzy \(X\) pozwala uprościć równanie \(\hat{Y} = X(X^TX)^{-1}X^TY\).

W pierwszym zadaniu przeprowadzimy analizę PCA ręcznie za pomocą rozkładu SVD, a następnie porównamy wyniki z wbudowanymi funkcjami prcomp oraz princomp. Ta pierwsza wykorzystuje rozkład SVD, ta druga bardziej prymitywny algorytm oparty na macierzy korelacji. Zadanie 1. Zadanie przykładowe. Wczytaj i obejrzyj zbiór danych dotyczący sportowców ais.txt, dostępny tutaj.
Wybierz zmienne numeryczne, wycentruj je i wyskaluj. Oblicz kierunki główne otrzymanej macierzy korzystając z rozkładu SVD zaimplementowanego w funkcji svd.
Sprawdź otrzymane składowe główne oraz macierz przejścia z bazy standardowej do bazy składowych głównych (czyli macierz zawierającą kierunki główne w kolumnach). Zbadaj procent wyjaśnionej zmienności w danych w zależności od liczby składowych głównych. Czy potrzeba wielu składowych żeby dobrze odwzorować dane?
Zbadaj ładowania pierwszych dwóch składowych głównych. Czy mają one jasną interpretację?

Przedstaw dwie pierwsze składowe główne na wykresie punktowym. Kolor punktu powinien odpowiadać uprawianemu sportowi, a kształt płci. Wnioskują po wykresie, która cecha bardziej odpowiada za zmienność danych? Czy problem klasyfikacji sportowców na podstawie numerycznych pomiarów wydaje się prosty czy trudny?

Porównaj swoje wyniki z wynikami funkcji prcomp(X, retx=T) oraz princomp(X, scores=T).
Uwaga na przyszłość: domyślnie żadna z tych funkcji nie centruje ani nie skaluje danych! Można za to je do tego przekonać odpowiednimi argumentami.

Rozwiązanie. Ładujemy dane i dzielimy je na dwie ramki danych, zawierające zmienne ilościowe oraz jakościowe. Te pierwsze centrujemy i skalujemy (uwaga, bez tego etapu otrzymalibyśmy zupełnie inne wyniki!).

dane <- read.table("ais.txt", header = T)
A <- dane[, c(1, 2)]
B <- dane[,-c(1,2)]
B <- scale(B)

Robimy rozkład SVD macierzy \(B\) i dla wygody zapisujemy macierze \(U, \Sigma, V\) w osobnych zmiennych. Ponieważ macierz \(\Sigma\) jest diagonalna, funkcja svd zwraca ją w postaci wektora.

B_svd <- svd(B)
B_V <- B_svd$v
B_Sigma <- B_svd$d
B_U <- B_svd$u

Standardowo w PCA istotność danej składowej ocenia się przez procent wariancji, jaką ona wyjaśnia. W tym celu podniesiemy wartości osobliwe do kwadratu:

B_Sigma^2
##  [1] 1003.183078  514.092951  232.617259  178.730370  159.852290   87.213876
##  [7]   21.137363    8.228131    4.661558    1.065500    0.217624

Jak widać, znaczna część wariancji jest wyjaśniona przez pierwsze kilka składowych. Sprawdźmy, jaki procent wariancji danych zachowamy, jeśli weźmiemy kilka pierwszych składowych. W tym celu obliczymy unormowaną sumę kumulatywną wariancji.

cumsum(B_Sigma^2)/sum(B_Sigma^2)
##  [1] 0.4537237 0.6862397 0.7914488 0.8722857 0.9445843 0.9840298 0.9935899
##  [8] 0.9973113 0.9994197 0.9999016 1.0000000

Jak widać, wzięcie pierwszych pięciu składowych wystarcza aby wyjaśnić 94% zmienności danych.
To umożliwia zastosowanie PCA do tak zwanej redukcji wymiaru - zamiast pracować na całej macierzy \(X\), wystarczy pracować na kilku pierwszych składowych głównych. W tym celu obliczamy macierz \(XV\) i zachowujemy tylko pierwszych kolumn.

Sprawdźmy dwa pierwsze kierunki główne:

data.frame('Zmienna' = colnames(B), 'Pierwszy kierunek' = B_V[, 1], 'Drugi kierunek' = B_V[,2])
##    Zmienna Pierwszy.kierunek Drugi.kierunek
## 1      RCC        0.37449886    -0.15896540
## 2      WCC        0.07610204     0.14658880
## 3       Hc        0.38921888    -0.16930064
## 4       Hg        0.39398641    -0.14839399
## 5     Ferr        0.18096731     0.04021179
## 6      BMI        0.25689422     0.42396429
## 7      SSF       -0.17657999     0.52587090
## 8   X.Bfat       -0.23763944     0.47350830
## 9      LBM        0.39995694     0.18522542
## 10      Ht        0.29428834     0.19798604
## 11      Wt        0.33804507     0.38332531

Współrzędne pierwszego kierunku głównego mówią nam, że pierwsza składowa będzie tym większa, im większa będzie wartość zmiennych RCC, Hc i Hg (odpowiadające liczbie czerwonych krwinek i hemoglobiny we krwi), LBM (masa mięśniowa), Ht (wzorst) i BMI, a tym mniejsza, im większy procent tkanki tłuszczowej (zmienna X.Bfat) oraz grubość fałdu skóry (zmienna SSF).

Obliczmy teraz składowe główne, czyli odpowiednio przekształcone dane. Skorzystamy w tym celu z mnożenia macierzowego:

SG <- B %*% B_V
SG <- data.frame(SG)
colnames(SG) <- paste0('PC', 1:11)  # 'PC' jak 'Principal Component'
SG$Sport <- dane$Sport
SG$Sex <- dane$Sex

Przedstawmy dwie pierwsze składowe na wykresie punktowym. Osie takiego wykresu, czyli składowe główne, odpowiadają kombinacjom liniowym wyjściowych cech.

library(ggplot2)
ggplot(SG) + geom_point(aes(x=PC1, y=PC2, col=Sport, pch=Sex)) + theme_minimal()

Powyższy wykres przedstawia nasze dane w dwóch wymiarach. Jak widać, pierwsza składowa główna rozróznia sportowców różnych płci - a więc płeć odpowiada za najwięcej zmienności w naszych danych. Mężczyźni mają wysokie wartości pierwszej składowej, co wydaje się być zgodne z tym, że pierwszy kierunek główny wskazuje na większy wzrost i masę mięśniową.
Na wykresie nie widać wyraźnego rozróżnienia dyscyplin sportowych. Moze to świadczyć albo o tym, że sportowcy z różnych dyscyplin nie różnią się parametrami (co jest raczej wątpliwe), albo o tym, że żeby zauważyć różnicę pomiędzy dyscyplinami należy wziąć więcej składowych głównych.
Z drugiej strony wszystkie dyscypliny które mamy w danych wymagają podobnego typu sprawności, możliwe więc, że parametry sportowców rzeczywiście są do siebie stosunkowo podobne - gdybyśmy uwzględnili sportowców uprawiających trójbój siłowy, prawdopodobnie łatwiej byłoby zauważyć różnicę.

Możemy się jeszcze przekonać na własne oczy, że wyrazy na diagonali macierzy \(\Sigma\) rzeczywiście są proporcjonalne do odchyleń standardowych składowych głównych. Konkretnie, są one równe odchyleniu standardowemu przemnożonemu przez pierwiastek z liczby obserwacji w danych. Policzmy ręcznie wariancje składowych głównych:

apply(B %*% B_V, 2, var)
##  [1] 4.990960588 2.557676373 1.157299795 0.889205819 0.795285026 0.433899882
##  [7] 0.105161008 0.040935973 0.023191833 0.005300995 0.001082707

Łatwo teraz sprawdzić, że są one równe kwadratom wyrazów macierzy \(\Sigma\) podzielonym przez liczbę obserwacji:

B_Sigma^2/nrow(B)
##  [1] 4.966252862 2.545014609 1.151570588 0.884803810 0.791347972 0.431751863
##  [7] 0.104640409 0.040733320 0.023077022 0.005274752 0.001077347

W praktyce do przeprowadzenia analizy składowych głównych w R korzystamy z funkcji prcomp lub princomp.
Ta pierwsza wykorzystuje rozkład SVD, dzięki czemu lepiej nadaje się do dużych zbiorów danych.

pr.pca <- prcomp(B, retx=T)  
prin.pca <- princomp(B, scores=T) 

Funkcja prcomp zwraca odchylenia standardowe składowych głównych w polu sdev i kierunki główne w polu rotation. Dodatkowo, jeśli przekażemy argument retx=T, zwróci składowe główne w polu x. Domyślnie centruje dane, za to nie skaluje kolumn przez ich odchylenia standardowe - ale można ją do tego przekonać odpowiednimi argumentami.

summary(pr.pca)
## Importance of components:
##                           PC1    PC2    PC3     PC4    PC5     PC6     PC7
## Standard deviation     2.2340 1.5993 1.0758 0.94298 0.8918 0.65871 0.32429
## Proportion of Variance 0.4537 0.2325 0.1052 0.08084 0.0723 0.03945 0.00956
## Cumulative Proportion  0.4537 0.6862 0.7914 0.87229 0.9446 0.98403 0.99359
##                            PC8     PC9    PC10   PC11
## Standard deviation     0.20233 0.15229 0.07281 0.0329
## Proportion of Variance 0.00372 0.00211 0.00048 0.0001
## Cumulative Proportion  0.99731 0.99942 0.99990 1.0000

Funkcja princomp działa niemal tak samo. Odchylenia standardowe mamy w polu sdev, macierz kierunków głównych \(V\) w polu loadings (tutaj mamy właśnie do czynienia ze wspomnianą niejednoznacznością terminologii), a jeśli przekażemy jej argument scores=T to otrzymamy składowe główne w polu scores. Do obliczania odchyleń standardowych akurat ta funkcja stosuje dzielenie przez \(n\), a nie przez \(n-1\), więc otrzymamy nieco inny wynik.

summary(prin.pca)
## Importance of components:
##                           Comp.1    Comp.2    Comp.3     Comp.4     Comp.5
## Standard deviation     2.2285091 1.5953102 1.0731126 0.94064011 0.88957741
## Proportion of Variance 0.4537237 0.2325160 0.1052091 0.08083689 0.07229864
## Cumulative Proportion  0.4537237 0.6862397 0.7914488 0.87228569 0.94458433
##                            Comp.6      Comp.7      Comp.8      Comp.9
## Standard deviation     0.65707828 0.323481698 0.201824973 0.151911231
## Proportion of Variance 0.03944544 0.009560092 0.003721452 0.002108348
## Cumulative Proportion  0.98402977 0.993589863 0.997311315 0.999419664
##                             Comp.10      Comp.11
## Standard deviation     0.0726274906 3.282296e-02
## Proportion of Variance 0.0004819086 9.842787e-05
## Cumulative Proportion  0.9999015721 1.000000e+00

Warto też zwrócić uwagę na jeden ważny fakt. Porównajmy kierunki drugich składowych głównych otrzymane tymi dwiema funkcjami.

data.frame('prcomp' = pr.pca$rotation[,2], 'princomp' = prin.pca$loadings[,2])
##             prcomp    princomp
## RCC    -0.15896540  0.15896540
## WCC     0.14658880 -0.14658880
## Hc     -0.16930064  0.16930064
## Hg     -0.14839399  0.14839399
## Ferr    0.04021179 -0.04021179
## BMI     0.42396429 -0.42396429
## SSF     0.52587090 -0.52587090
## X.Bfat  0.47350830 -0.47350830
## LBM     0.18522542 -0.18522542
## Ht      0.19798604 -0.19798604
## Wt      0.38332531 -0.38332531

Kierunki są one identyczne z dokładnością do znaku. Jest to całkowicie naturalne: w tej technice mamy jednoznacznie wyznaczone kierunki składowych głównych, ale nie zwroty! Zmiana znaku kierunku głównego w żaden sposób nie wpływa na to, jak dobrze przybliża on dane, o czym najłatwiej się przekonać patrząc na ilustrację z początku tego skryptu. Konsekwencją tej zmiany znaku jest m.in. to, że na wykresie pierwszych dwóch składowych głównych otrzymanych funkcją princomp oś Y jest odbita symetrycznie.
Mimo że może się to wydawać oczywiste, to warto to dobrze zapamiętać, żeby w przyszłości nie dać się zmylić.

library(patchwork)  # przydatna biblioteka do łączenia wykresów z ggplot2
p1 <- ggplot(data.frame(pr.pca$x)) + geom_point(aes(x=PC1, y=PC2)) + theme_minimal() + ggtitle('', subtitle='Otrzymane funkcją prcomp')
p2 <- ggplot(data.frame(prin.pca$scores)) + geom_point(aes(x=Comp.1, y=Comp.2)) + theme_minimal() + ggtitle('', subtitle='Otrzymane funkcją princomp')
p1 + p2 + plot_annotation(title = 'Porównanie pierwszych dwóch składowych głównych')

knitr::opts_chunk$set(echo = F, include=F)

Zadanie 2. Jak przekonaliśmy się w poprzednim zadaniu, analiza składowych głównych może pomóc nam odkryć, że w danych znajdują się różne grupy obserwacji (tak jak na przykład sportowcy dwóch płci). Może również pomóc znaleźć nam kombinacje liniowe cech które dobrze róznicują te grupy, jeśli odpowiadają one za znaczną część wariancji w danych. W dodatku jest to tak zwana technika uczenia bez nadzoru, co oznacza, że nie musimy wiedzieć do jakich grup należą obserwacje - opiera się ona wyłącznie na wartościach zmiennych objaśniających.
W sytuacji gdy znamy te grupy, możemy wykorzystać PCA do znalezienia liniowej kombinacji cech które pozwolą te klasy rozróżnić (choć akurat do tego celu lepiej nadaje się liniowa analiza dyskryminacyjna, z którą spotkaliśmy się na jednym z poprzednich laboratoriów). Na przykład, wiemy, że pierwsza główna składowa z poprzedniego zadania może być wykorzystana do odróżnienia płci sportowca (choć akurat do tego celu również istnieją metody lepsze niż PCA).

W praktyce, nawet jeśli w danych mamy wyraźnie wyodrębnione grupy, istnieje kilka przypadków w których PCA zawodzi. Jednym z nich jest obecność w danych znacznego szumu - to znaczy wielu zmiennych, które nie mają żadnego wpłwu na interesujące nas grupy. Tego typu szum może ukryć obecność jakichkolwiek grup w danych, o czym przekonamy się w tym zadaniu.

Przeprowadź analizę składowych głównych na kolumnach ilościowych danych iris. Przedstaw dwie pierwsze składowe na wykresie punktowym, pokolorowanym w zależności od gatunku. Następnie utwórz dane złożone z tych kolumn ilościowych oraz stu dodatkowych kolumn z obserwacjami wylosowanymi ze standardowego rozkładu normalnego. Przeprowadź analizę składowych głównych na tak stworzonych danych, zwizualizuj piersze dwie składowe główne i sprawdź, w jaki sposób obecność nieinformatywnych kolumn wpłynęła na separację gatunków. Sprawdź co się stanie po zmniejszeniu odchylenia standardowego rozkładu normalnego z którego losujemy szum do wartości 0.1 oraz po zwiększeniu do 4. Skąd wynika obserwowane zjawisko?

W wielu praktycznych zastosowaniach istnieją różnorodne metody usuwania szumu z danych, oparte na przykład na różnicach w rozkładzie cech. Czasami na podstawie wiedzy dziedzinowej wiemy, że interesujące cechy mają rozkład inny niż normalny, i na tej podstawie możemy odrzucić kolumny które mogą zaburzać analizę.
Jeśli jest to możliwe, to warto wykorzystać tego typu metody przed przystąpieniem do eksploracyjnej analizy danych.

Zadanie 3. W tym zadaniu przekonamy się, że PCA jest również wrażliwe na obserwacje odstające (i to jak!).
Stwórz macierz zawierającą 100 obserwacji z 3-wymiarowego standardowego rozkładu normalnego (co oznacza po prostu to, że każda spośród 3 kolumn zawiera 100 obserwacji wylosowanych ze standardowego rozkładu normalnego).
Przemnóż pierwszą kolumnę przez 4, a drugą przez 2. Jaki kształt mają dane otrzymane w ten sposób? Jeśli chcesz, to możesz zweryfikować swoje przypuszczenia korzystając z interaktywnego wykresu 3d stworzonego za pomocą pakietu rgl.

Na podstawie kształtu danych, jakie kierunki główne spodziewasz się uzyskać? Przeprowadź analizę składowych głównych i zweryfikuj swoje podejrzenia patrząc na wektory macierzy kierunków głównych (na przykład macierzy rotation z wyniku funkcji prcomp). Przedstaw pierwsze dwie składowe główne na wykresie punktowym.

Zamień ostatnią obserwację na wektor c(40, 40, 40). Przeprowadź ponownie analizę składowych głównych. Na co wskazuje pierwszy kierunek główny? Przedstaw pierwsze dwie składowe główne na wykresie punktowym i porównaj go z poprzednim.

Regresja składowych głównych

Zadanie 4. Załaduj dane biopsy z pakietu MASS. Wybierz zmienne numeryczne i usuń brakujące obserwacje. Zrób regresję składowych głównych aby wyjaśnić zmienną V2. Wykorzystaj funkcję pcr z biblioteki pls. Odpowiedz na następujące pytania: