Eksploracja danych geoprzestrzennych

Redukcja wymiarowości

Author

Krzysztof Dyba

Wprowadzenie

W rzeczywistych zbiorach danych często spotykamy się z istotnie dużą liczbą zmiennych (cech), które opisują pewne obserwacje (obiekty). Jest to tak zwane zjawisko wielowymiarowości danych (high dimensionality). W praktyce, powoduje to problemy związane z ilością wymaganej pamięci do przetwarzania danych, dłuższym czasem obliczeń, trudnościami w interpretacji zależności i ich wizualizacją w przestrzeni dwuwymiarowej, a także zmniejszenie skuteczności modelu wraz ze wzrostem liczby zmiennych (w takiej sytuacji model zazwyczaj słabo generalizuje zależności).

Rozwiązaniem może być zastosowanie metody opartej na analizie głównych składowych (principal component analysis, PCA), która jest najczęściej stosowanym podejściem do redukcji wymiarowości. Przekształca ona zbiór danych ze skorelowanymi zmiennymi w nowy zbiór nieskorelowanych (ortogonalnych) zmiennych nazywanymi głównymi składowymi (principal components), zachowując przy tym jak najwięcej zmienności (informacji). Główne składowe są liniowymi kombinacjami oryginalnych zmiennych, uporządkowanymi według ilości wariancji, którą wyjaśniają. Oznacza to, że pierwsza główna składowa (PC1) wyjaśnia największą zmienność (wariancję), a każda kolejna składowa (PC2, PC3, itd.) mniej.

Analiza głównych składowych umożliwia przede wszystkim zmniejszenie liczby zmiennych w zbiorze danych i usunięcie współliniowości między nimi. Oprócz tego, przydatna jest do wizualizacji danych wielowymiarowych w dwóch wymiarach wykorzystując najczęściej dwie najistotniejsze główne składowe (PC1 i PC2).

Technika ta wymaga spełnienia kilku założeń:

  1. Relacje między zmiennymi są liniowe.
  2. Zmienne są typu numerycznego.
  3. Wszystkie obserwacje mają wartości (brak NA w zbiorze danych).
  4. Zmienne powinny mieć rozkład zbliżony do normalnego i brak wartości odstających.
  5. Zmienne muszą posiadać jednakowe skale (wymagana standaryzacja).

Niemniej, nie jest to jedyne podejście do redukcji wymiarowości. Wyróżnić można także skalowanie wielowymiarowe (multidimensional scaling), autoenkodery oparte o sieci neuronowe czy techniki służące wyłącznie do wizualizacji t-SNE i UMAP.

Wariancja

Przed przystąpieniem do analizy, należy wytłumaczyć czym jest właściwie wariancja. Jest to miara statystyczna, która opisuje ilościowo zróżnicowanie (rozrzut) wartości w zbiorze danych, określając jak bardzo odbiegają od wartości średniej tego zbioru. Wysoka wariancja wskazuje, że punkty są rozproszone, podczas gdy niska wariancja oznacza, że są one skupione blisko średniej. Zerowa wariancja oznacza, że wszystkie wartości są identyczne (brak zmienności).

# wariancja dla próby
sum((x - mean(x))^2 / (length(x) - 1))

Pamiętaj, że wariancja jest wyrażona w jednostkach kwadratowych (\(j^2\)) oryginalnych danych, co może sprawić, że interpretacja będzie mniej intuicyjna. Z tego powodu często preferowane jest odchylenie standardowe, tj. pierwiastek kwadratowy z wariancji.

Analiza głównych składowych zakłada, że główne składowe o największej wariancji zawierają najważniejsze informacje. Jeśli główna składowa nie wykazuje prawie żadnej wariancji, to ma ona niewielki wkład w reprezentację zbioru danych i można ją pominąć.

Analiza danych tabelarycznych

library("corrplot")

Przygotowanie danych

Analizę przeprowadzimy na przykładzie danych społeczno-ekonomicznych z World Bank Data. Dla uproszczenia przyjmijmy, że interesują nas wyłącznie kraje europejskie.

# wczytanie danych w formacie .csv
data = read.csv2("../dane/world_bank_data.csv")
# selekcja krajów europejskich
data = data[data$region == "EU" & !is.na(data$region), ]
# wyświetlenie struktury danych
str(data)
'data.frame':   47 obs. of  12 variables:
 $ country   : chr  "Albania" "Andorra" "Austria" "Belarus" ...
 $ region    : chr  "EU" "EU" "EU" "EU" ...
 $ code      : chr  "ALB" "AND" "AUT" "BLR" ...
 $ gdp       : num  3.27e+10 NA 4.12e+11 1.68e+11 4.96e+11 ...
 $ pop       : int  2889167 70473 8611088 9513000 11285721 3810416 7177991 163692 4224404 1165300 ...
 $ pop_growth: num  1.07 3.06 0.36 0.27 0.38 0.31 -0.17 0.74 0.04 1.3 ...
 $ internet  : num  63.3 96.9 83.9 62.2 85.1 65.1 56.7 NA 69.8 71.7 ...
 $ urban_pop : num  27.3 NA 30.9 26.3 18.5 21 23.1 NA 27.6 NA ...
 $ life_exp  : num  77.8 NA 81.3 73 80.6 76.4 75.4 80.6 77.3 80.1 ...
 $ literacy  : num  96.8 NA NA 99.7 NA ...
 $ export_gdp: num  0.271 NA 0.534 0.601 0.844 ...
 $ gdp_pc    : num  11305 NA 47824 17661 43992 ...

Jak możemy zauważyć powyżej, w naszej ramce danych mamy liczne wartości brakujące NA (w szczególności w kolumnie literacy), co uniemożliwi nam wykonanie analizy głównych składowych. Zatem musimy je usunąć.

# usuń całą kolumnę
data$literacy = NULL
# usuń wiersze z brakującymi wartościami
idx = complete.cases(data)
data = data[idx, ]

Poglądowo sprawdźmy również ogólną korelację liniową pomiędzy zmiennymi.

# macierz korelacji
cor_mat = cor(data[, 4:11], method = "pearson")
corrplot(cor_mat, method = "number", type = "lower", diag = FALSE)

Najsilniejsza korelacja występuję pomiędzy całkowitą populacją kraju a produktem krajowym brutto (PKB). Wynika to z prostego związku, iż PKB mierzy wartość wszystkich dóbr i usług wytworzonych w gospodarce w danym okresie, a większa liczba ludności zwiększa produkcję oraz konsumpcję. Kolejna silna zależność występuje pomiędzy produktem krajowym brutto na osobę a dostępem do internetu i spodziewaną długością życia. Pozostałe zmienne nie wykazuje aż tak silnych zależności.

Analiza głównych składowych jest wrażliwa na skalę zmiennych, dlatego ważne jest uprzednie ustandaryzowanie danych, tj. wyśrodkowanie (poprzez odjęcie średniej) i przeskalowanie (poprzez podzielenie przez odchylenie standardowe). Zapewnia to, że zmienne o większej skali nie zdominują pozostałych zmiennych. Do standaryzacji danych służy funkcja scale().

# skalowanie danych
data_scale = scale(data[4:11])

Po tej operacji średnia wartość każdej zmiennej wynosi 0, a odchylenie standardowe jest równe 1.

Analiza głównych składowych

Do analizy głównych składowych można wykorzystać między innymi funkcję prcomp(). Jeśli dokonaliśmy wcześniej standaryzacji danych, to koniecznie musimy nadać argumentom scale. i center wartość FALSE, aby uniknąć podwójnej standaryzacji.

pca = prcomp(data_scale, scale. = FALSE, center = FALSE)
summary(pca)
Importance of components:
                          PC1    PC2    PC3     PC4     PC5     PC6     PC7
Standard deviation     1.8263 1.4490 1.0349 0.84846 0.65724 0.46883 0.30846
Proportion of Variance 0.4169 0.2624 0.1339 0.08998 0.05399 0.02748 0.01189
Cumulative Proportion  0.4169 0.6794 0.8133 0.90323 0.95723 0.98470 0.99659
                           PC8
Standard deviation     0.16507
Proportion of Variance 0.00341
Cumulative Proportion  1.00000

W podsumowaniu wyników otrzymujemy:

  • Odchylenie standardowe (Standard deviation) – wskazuje istotność każdej głównej składowej (im większa wartość, tym większa istotność).
  • Proporcja wariancji (Proportion of Variance) – proporcja całkowitej wariancji wyjaśniona przez każdą główną składową.
  • Proporcja skumulowana (Cumulative Proportion) – skumulowana część wariancji wyjaśniona przez pierwsze \(k\) głównych składowych.

Dodatkowo, otrzymaliśmy dwa obiekty rotation i x. Ten pierwszy reprezentuje macierz ładunków (loadings), natomiast drugi zawiera nowe wartości obserwacji po transformacji rzutowane na główne składowe (można je później wykorzystać do modelowania).

# macierz ładunków PC1 i PC2
pca$rotation[, 1:2]
                  PC1        PC2
gdp         0.4263113  0.3668460
pop         0.3563355  0.4905988
pop_growth  0.2100404 -0.2323769
internet    0.3709802 -0.3990428
urban_pop  -0.3091641 -0.2715945
life_exp    0.4668035 -0.2207007
export_gdp -0.1106244 -0.3903555
gdp_pc      0.4287361 -0.3696659

Ładunki informują, w jakim stopniu każda zmienna wejściowa ma wpływ na każdą główną składową, tzn. wysoka wartość ładunku w pierwszej składowej głównej oznacza, że jest silnie skorelowana z tą składową. Natomiast znak wskazuje kierunek zależności.

Dla przykładu, zmienne life_exp, gdp oraz gdp_pc mają wysoki dodatni ładunek na pierwszej składowej (PC1), oznacza to, że ta składowa jest silnie zależna od tych zmiennych. Zobaczymy jak wygląda zależność pomiędzy pierwszą główną składową (PC1) a zmienną gdp na wykresie rozrzutu.

plot(data$gdp, pca$x[, 1], pch = 19, xlab = "PKB", ylab = "PC1")

Interpretacja

Najważniejszym krokiem jest zrozumienie, w jaki sposób interpretować otrzymane wyniki. W tym celu będą pomocne dwa wykresy.

Wykres piargowy

Wykres piargowy (scree plot) to wykres wyjaśnionej wariancji dla każdej głównej składowej. Pomaga określić, ile głównych składowych należy zachować do dalszej analizy. Punkt przegięcia, tzw. “łokieć” wykresu zazwyczaj wskazuje punkt, w którym dodanie większej liczby składowych w niewielkim stopniu przyczynia się do ogólnej wyjaśnionej wariancji. Aby go wyświetlić wystarczy użyć funkcji plot() na obiekcie klasy prcomp.

plot(pca, main = "Analiza głównych składowych", xlab = "Główna składowa")
abline(v = 3.7, lty = 2, col = "red")

Powyższy wykres sugeruje, że powinniśmy zachować pierwsze trzy główne składowe, które wyjaśniają ponad 80% zmienności w danych. Dodanie czwartej głównej składowej pozwoliłoby na wyjaśnienie ponad 90% zmienności. Pozostałe składowe reprezentują mniejszą ilość informacji (zasadniczo jest to szum w danych).

Biplot

Biplot to wykres punktowy przedstawiający rzutowane obserwacje na wybrane główne składowe oraz strzałki reprezentujące ładunki oryginalnych zmiennych, które tworzą pewne kąty. Jeśli wartość kąta pomiędzy strzałkami jest bliska:

  • 0° to zmienne są silnie dodatnio skorelowane.
  • 180° to zmienne są również silnie skorelowane, ale przeciwnie.
  • 90° to zmienne są niezależne (ortogonalne).

Biplot można stworzyć używając funkcji biplot() na obiekcie klasy prcomp. Argument choices pozwoli nam wybrać główne składowe.

biplot(pca, choices = 1:2, cex = 0.6, xlabs = data$country)

W przypadku większych zbiorów danych, składających się z minimum kilkuset obserwacji, odczytanie ich nazw nie będzie możliwe. W takiej sytuacji warto zamienić etykiety nazw na kropki, ustawiając argument xlabs = rep("·", nrow(data)). Alternatywnie, można całkowicie je pominąć i interpretować jedynie położenie głównych składowych.

Niemniej, interpretacja głównych składowych jest skomplikowanym procesem. Z tego powodu warto rozważyć najprostsze techniki redukcji wymiarowości oparte na usuwaniu cechy o bliskiej zeru wariancji (mała zmienność to mało informacji) oraz cech, które są ze sobą silnie skorelowane (np. korelacja powyżej 90%), aby zmniejszyć redundancję danych.

Analiza danych przestrzennych

library("terra")

W tej sekcji zastosujemy redukcję wymiarowości do danych przestrzennych, a wynik zaprezentujemy na mapie. Zacznijmy od wylistowania rastrów, które mamy dostępne w katalogu dane używając funkcji list.files(). Doprecyzujmy, że interesują nas wyłącznie pliki rastrowe z rozszerzeniem .tif.

# listowanie plików w folderze
files = list.files("../dane/", pattern = "\\.tif$", full.names = TRUE)
files = files[1:3]

Następnie używając funkcji rast() możemy wczytać wybrane trzy warstwy jako jeden raster wielokanałowy, tj. składający się z trzech warstw. Pamiętaj, żeby stworzyć raster wielokanałowy, to poszczególne warstwy muszą posiadać tę samą geometrię (liczba kolumn i wierszy), zakres przestrzenny oraz układ przestrzenny!

# wczytanie rastrów
r = rast(files)
r
class       : SpatRaster 
dimensions  : 270, 300, 3  (nrow, ncol, nlyr)
resolution  : 100, 100  (x, y)
extent      : 345000, 375000, 493000, 520000  (xmin, xmax, ymin, ymax)
coord. ref. : Transverse Mercator 
sources     : 01_roughness.tif  
              02_green.tif  
              03_distwaters.tif  
names       : xmod@urban, MEAN_green@urban, DIST_waters@urban

Wybrane rastry zostały wczytane jako obiekt r. Kolejne kroki analizy wykonamy na mniejszej próbie losowej. Takie podejście stosuje się w przypadku przetwarzania ogromnych zbiorów danych, ponieważ nie jest możliwe, aby wczytać je do pamięci komputera. Do losowania próby służy funkcja spatSample(), w której możemy zdefiniować jej wielkość oraz usuwanie brakujących wartości NA. Szczególnie ważne jest zadbanie o powtarzalność analizy, więc koniecznie jest wcześniejsze zdefiniowanie ziarna losowości set.seed().

# ziarno losowości
set.seed(1)
# losowanie próby
smp = spatSample(r, size = 5000, na.rm = TRUE, as.points = TRUE)

Teraz możemy wyświetlić wylosowane punkty na mapie.

plot(r[[1]])
plot(smp, add = TRUE, cex = 0.4, col = "red")

Większość narzędzi czy modeli operuje na typowych strukturach danych takich jak macierze czy ramki danych (niekoniecznie pliki rastrowe czy wektorowe). Z tego powodu dokonamy konwersji punktów (klasa SpatVector) na ramkę danych wykorzystując funkcję as.data.frame().

# konwersja punktów na ramkę danych
df_smp = as.data.frame(smp)

Następnie dane z próby należy zestandaryzować do wspólnej skali przy użyciu funkcji scale(). Istotną czynnością jest zapisane współczynników centrowania i skalowania z próby, ponieważ będzie trzeba je później zastosować dla całych rastrów (wszystkich komórek).

# standaryzacja próby
df_smp_scale = scale(df_smp, center = TRUE, scale = TRUE)
# współczynnik centrowania
v_center = attr(df_smp_scale, "scaled:center")
# współczynnik skalowania
v_scale = attr(df_smp_scale, "scaled:scale")

Po standaryzacji danych możemy obliczyć główne składowe za pomocą poznanej już funkcji prcomp().

pca = prcomp(df_smp_scale, center = FALSE, scale. = FALSE)

Po wyznaczeniu głównych składowych, przechodzimy do najtrudniejszego etapu, polegającego na ekstrapolacji wyliczonych głównych składowych z próby na cały raster. Zasadniczo, tę procedurę przeprowadza się sekwencyjnie z podziałem na mniejsze bloki rastra, jednakże nasze rastry są małe i możemy je wczytać całe do pamięci.

Ponownie wykorzystajmy funkcję as.data.frame(), aby dokonać konwersji rastrów (SpatRaster) do ramki danych. Ustawmy także argument na.rm = FALSE do zachowania pustych wartości w rastrze. Następnie, kluczowym krokiem jest wyznaczenie indeksów komórek, które posiadają wartości używając funkcji complete.cases(), ponieważ będzie to niezbędne do wstawienia wartości w odpowiednie miejsca w wektorze.

df = as.data.frame(r, na.rm = FALSE)
# indeksy komórek z wartościami
idx = complete.cases(df)
# usuń puste wartości
df = df[idx, ]

Tak jak wcześniej, standaryzujemy całe dane, z tą różnicą, iż definiujemy uprzednio wyliczone współczynniki centrowania (v_center) i skalowania (v_scale).

# standaryzacja całego zbioru
df = scale(df, center = v_center, scale = v_scale)

Mamy już przygotowane dane, teraz możemy wykonać ekstrapolację używając funkcji predict() i podając dwa wymagane argumenty. Pierwszy to obliczone składowe główne w obiekcie pca, a drugi obiekt to ramka danych df.

# PCA dla całego zbioru danych
pr = predict(pca, df)

Wynikiem powyższej predykcji jest macierz, w której w wierszach znajdują się komórki rastra (blisko 80 tysięcy), a w kolumnach składowe główne. Zauważ, że raster wejściowy składa się z 81 tysięcy komórek, natomiast z predykcji otrzymaliśmy ich mniej, co spowodowane jest wykluczeniem brakujących wartości. Rozwiązanie tego problemu jest proste. Należy wstawić wartości predykcji tylko w tych komórkach, które posiadały wartości.

Stwórzmy zatem nowy wektor zawierający wyłącznie brakujące wartości NA o takiej długości, ile raster posiada wszystkich komórek (270 * 300). Następnie używając indeksów z obiektu idx przypiszemy wartości pierwszej głównej składowej do wektora vec w odpowiednich miejscach.

# stwórz pusty wektor
vec = rep(NA, ncell(r))
# przypisz wartości z PC1
vec[idx] = pr[, 1]

Następnie tworzymy nowy raster o nazwie output kopiując metadane z rastra r oraz definiując jeden kanał i wartości na podstawie wektora vec.

# stwórz nowy raster
output = rast(r, nlyrs = 1, vals = vec)

Finalnie, dokonajmy wizualizacji pierwszej głównej składowej (PC1) na mapie.

plot(output, main = "Pierwsza składowa główna (PC1)")

Zadanie

Dokonaj redukcji wymiarowości rastrów znajdujących się w katalogu dane używając analizy głównych składowych i losowej próby. Wcześniej przygotuj zbiór danych do analizy, tj. usuń brakujące i odstające wartości oraz spróbuj transformować zmienne do rozkładu normalnego. Zaprezentuj pierwszą składową na mapie oraz zinterpretuj, co ona przedstawia. Otrzymane wyniki i wnioski przedstaw w formie raportu (Quarto lub Markdown) z wykresami oraz wyjaśnij, dlaczego podjąłeś określone decyzje.