Eksploracja danych geoprzestrzennych

Analiza danych

Author

Krzysztof Dyba 

# ziarno losowości
set.seed(1)

Rozkład danych

Rozkład danych opisuje w jaki sposób wartości w zbiorze są zorganizowane. Zrozumienie rozkładów ma kluczowe znaczenie dla interpretacji danych, wyboru odpowiednich metod statystycznych i wyciągnięcia wniosków. Testy i modele statystyczne opierają się na założeniach związanych z rozkładem wartości analizowanych zmiennych. Użycie niewłaściwej metody dla danego rozkładu może prowadzić do błędnych wyników, a finalnie wniosków.

Rodzaje rozkładów

Istnieje kilka podstawowych rodzajów rozkładu danych, w tym:

  • Rozkład normalny (rozkład Gaussa) – kształt przypomina dzwon i jest symetryczny wokół średniej. Odchylenie standardowe określa rozrzut od średniej. Średnia, mediana i moda są równe.
  • Rozkład jednostajny – ma płaski kształt i wszystkie wartości mają identyczne prawdopodobieństwo wystąpienia. Średnia i mediana są równe.
  • Rozkład dwumodalny – posiada dwa odrębne szczyty (mody), co może wskazywać na występowanie dwóch różnych podgrup lub procesów. W tym przypadku, średnia i mediana nie są zbyt przydatne.
  • Rozkład skośny – rozkład asymetryczny z większością wartości skupioną po jednej stronie. Może być prawoskośny, wtedy: średnia > mediana > moda lub lewoskośny, wtedy: średnia < mediana < moda.

Rozkład normalny

data = rnorm(1000, mean = 0, sd = 1)

hist(data, main = "Rozkład normalny", xlab = NULL, ylab = "Częstość")
abline(v = 0, col = "red", lty = 1) # średnia
abline(v = c(-1, 1), col = "orange", lty = 2) # 1 odchylenie standardowe (68%)

Rozkład jednostajny

data = runif(1000, min = 0, max = 10)

hist(data, main = "Rozkład jednostajny", xlab = NULL, ylab = "Częstość")

Rozkład dwumodalny

data1 = rnorm(1000, mean = 5, sd = 1) # pierwsza grupa
data2 = rnorm(1000, mean = 10, sd = 1) # druga grupa
bimodal_data = c(data1, data2)

hist(bimodal_data, main = "Rozkład dwumodalny", xlab = NULL, ylab = "Częstość")
abline(v = 5, col = "red")
abline(v = 10, col = "red")

Rozkład skośny

# rozkład prawoskośny (rozkład wykładniczy)
right_skew_data = rexp(1000, rate = 0.5)
# rozkład lewoskośny (rozkład beta)
left_skew_data = rbeta(1000, 5, 2)

par(mfrow = c(1, 2))
hist(right_skew_data, main = "Rozkład prawoskośny", xlab = NULL, ylab = "Częstość")
hist(left_skew_data, main = "Rozkład lewoskośny", xlab = NULL, ylab = NULL)

Testowanie hipotez

Testowanie hipotez to metoda podejmowania decyzji lub wyciągania wniosków dotyczących populacji na podstawie mniejszej próby. Polega ono na weryfikacji pewnego założenia (hipotezy) w wyniku czego można je potwierdzić lub obalić.

Na początku zdefiniowane są dwie przeciwstawne hipotezy:

  1. Hipoteza zerowa (\(H_0\)) – domyślne założenie stwierdzające brak efektu, brak różnicy czy brak związku. Tę hipotezę staramy się obalić.
  2. Hipoteza alternatywna (\(H_1\)) – zaprzeczenie hipotezy zerowej, czyli stwierdzenie istnienia efektu, różnicy czy związku. Podejrzewamy, że ta hipoteza może być prawdziwa i próbujemy ją potwierdzić poprzez obalenie hipotezy zerowej.

Następnie musimy określić poziom istotności testu (\(\alpha\)), czyli wartość progową (zazwyczaj przyjmuje się wartość 0,05 lub 0,01), która definiuje prawdopodobieństwo odrzucenia hipotezy zerowej (\(H_0\)), gdy jest ona prawdziwa. Wiąże się to z pewnym ryzykiem, że odrzuciliśmy właściwą hipotezę.

W finalnym etapie wykonujemy test statystyczny, z którego zwracane jest prawdopodobieństwo testowe nazywane wartością \(p\) (p-value) zakładające, że hipoteza zerowa (\(H_0\)) jest prawdziwa. Najważniejszym zadaniem jest porównanie założonego poziomu istotności testu (\(\alpha\)) z otrzymaną wartością \(p\):

  • Jeśli wartość \(p \leq \alpha\), to odrzucamy hipotezę zerową (\(H_0\)). Oznacza to, że wynik jest istotny statystycznie, co sugeruje, że hipoteza alternatywna (\(H_1\)) jest bardziej prawdopodobna. Jednak nie dowodzi tego bezpośrednio.
  • Jeśli \(p > \alpha\), to nie możemy odrzucić hipotezy zerowej (\(H_0\)). Oznacza to, że wynik nie jest statystycznie istotny i nie ma wystarczających dowodów na jej obalenie.

Co ważne, wartość \(p\) nie mierzy wielkości ani praktycznego wpływu efektu. Efekt może być istotny statystycznie, ale mały.

Rozkład normalny

Rozkład normalny jest fundamentem statystyki, ponieważ jest symetryczny i ma dobrze zdefiniowane właściwości (tj. 68% danych mieści się w zakresie 1 odchylenia standardowego od średniej, 95% w zakresie 2 itd.), dzięki czemu przewidywanie prawdopodobieństwa i przedziałów ufności jest prostsze.

Co więcej, centralne twierdzenie graniczne głosi, iż rozkład średniej z próby wystarczająco dużej liczby niezależnych zmiennych losowych będzie zbliżony do rozkładu normalnego bez względu od rozkładu bazowego. Wiele testów parametrycznych (np. t-test, ANOVA) czy regresja liniowa, zakładają, że dane lub reszty mają rozkład normalny. Jeśli dane znacznie odbiegają od rozkładu normalnego, metody te mogą dawać błędne wyniki. Niemniej, istnieją także alternatywy nieparametryczne dla innych rozkładów.

Testowanie normalności rozkładu

Do określenia czy dane mają rozkład normalny, można użyć metod wizualnych (histogram, wykres gęstości czy wykres kwantyl-kwantyl) lub testów statystycznych. Możemy wykorzystać test Shapiro-Wilka sprawdzający hipotezę zerową czy dane mają rozkład normalny lub test Kołmogorowa-Smirnowa porównujący empiryczny rozkład danych z rozkładem normalnym. Metody wizualne powinny zawsze uzupełniać testy statystyczne, ponieważ są bardziej odporne na odchylenia od normy.

Test Shapiro-Wilka

Hipotezą zerową (\(H_0\)) testu Shapiro-Wilka jest założenie, że dane pochodzą z populacji o rozkładzie normalnym (tj. wartości mają rozkład normalny). To założeniem chcemy obalić, wtedy przyjmiemy hipotezę alternatywną (\(H_1\)), która sugeruje, że dane nie mają rozkładu normalnego. Ograniczeniem tego testu jest maksymalna liczba obserwacji, która wynosi 5000.

data = rnorm(1000, mean = 0, sd = 1)
shapiro.test(data)

    Shapiro-Wilk normality test

data:  data
W = 0.99825, p-value = 0.4028

Wartość \(p\) jest powyżej 0,05, więc nie możemy odrzucić hipotezy zerowej (\(H_0\)), że dane mają rozkład zbliżony do normalnego. W uproszczeniu mówiąc, wartość \(p\) większa od 0,05 sugeruje, że dane mogą mieć rozkład normalny.

data = runif(1000, min = 0, max = 10)
shapiro.test(data)

    Shapiro-Wilk normality test

data:  data
W = 0.95014, p-value < 2.2e-16

Wartość \(p\) jest mniejsza niż 0,05, więc możemy odrzucić hipotezę zerową (\(H_0\)), że dane mają rozkład zbliżony do normalnego, przyjmując hipotezę alternatywną (\(H_1\)). W powyższym przykładzie rozkład jest jednostajny.

Test Kołmogorowa-Smirnowa

Interpretacja testu Kołmogorowa-Smirnowa jest identyczna jak w przypadku testu Shapiro-Wilka. Hipoteza zerowa (\(H_0\)) zakłada, że próba pochodzi z określonego rozkładu, w tym przypadku rozkładu normalnego pnorm. Jeśli wartość \(p\) jest większa od 0,05, to prawdopodobnie mamy do czynienia z rozkładem zmiennej zbliżonym do normalnego.

data = rnorm(1000, mean = 0, sd = 1)
ks.test(data, "pnorm", mean = mean(data), sd = sd(data))

    Asymptotic one-sample Kolmogorov-Smirnov test

data:  data
D = 0.015393, p-value = 0.9718
alternative hypothesis: two-sided

Transformacja danych

Jeśli zmienna nie ma rozkładu zbliżonego do normalnego, to możemy podjąć próbę jej transformacji, co potencjalnie pomoże spełnić wymagania testów parametrycznych i modeli liniowych, zmniejszy wpływ skośności oraz wartości odstających, a także zwiększy liniowość relacji między zmiennymi. Jednak, należy mieć na uwadze, iż transformacja nie zawsze pomaga w przypadku trudnych rozkładów (w takiej sytuacji przydatne są testy nieparametryczne oraz modele nieliniowe czy uczenia maszynowego). Warto także wypróbować różne transformacje, aby sprawdzić, która z nich daje najlepsze wyniki.

Podczas interpretacji wyników należy dokonać transformacji wstecznej do oryginalnej skali wartości. Jeśli do transformacji został użyty logarytm, to operacją odwrotną jest funkcja wykładnicza.

Pierwiastek

Pierwiastkowania nie można stosować do wartości ujemnych.

data_sqrt = sqrt(right_skew_data)

par(mfrow = c(1, 2))
hist(right_skew_data, main = "Rozkład prawoskośny", xlab = NULL, ylab = "Częstość")
hist(data_sqrt, main = "Transformacja pierwiastkowa", xlab = NULL, ylab = NULL)

Logarytm

Logarytmowania nie można przeprowadzić dla danych z wartościami zerowymi lub ujemnymi. W przypadku zer, można dodać bardzo małą wartość przed obliczeniem logarytmu, np. 0,000001.

data_log = log(right_skew_data)

par(mfrow = c(1, 2))
hist(right_skew_data, main = "Rozkład prawoskośny", xlab = NULL, ylab = "Częstość")
hist(data_log, main = "Transformacja logarytmiczna", xlab = NULL, ylab = NULL)

Potęga

data_pow = left_skew_data^2

par(mfrow = c(1, 2))
hist(left_skew_data, main = "Rozkład lewoskośny", xlab = NULL, ylab = "Częstość")
hist(data_pow, main = "Transformacja potęgowa", xlab = NULL, ylab = NULL)

Box-Cox

Podobnie jak transformacja logarytmiczna, transformacja Boxa-Coxa nie działa w przypadku wartości poniżej i równych zero. Funkcja dostępna jest w wielu pakietach, np. w bestNormalize.

library("bestNormalize")
data_boxcox = boxcox(right_skew_data, standardize = FALSE)
par(mfrow = c(1, 2))
hist(right_skew_data, main = "Rozkład prawoskośny", xlab = NULL, ylab = "Częstość")
hist(data_boxcox$x.t, main = "Transformacja Boxa-Coxa", xlab = NULL, ylab = NULL)

Yeo-Johnson

Transformacja Yeo-Johnsona jest rozszerzeniem transformacji Boxa-Coxa dla liczb ujemnych i równych zero. Funkcja dostępna jest w wielu pakietach, np. w bestNormalize.

library("bestNormalize")
data_yeo = yeojohnson(right_skew_data, standardize = FALSE)
par(mfrow = c(1, 2))
hist(right_skew_data, main = "Rozkład prawoskośny", xlab = NULL, ylab = "Częstość")
hist(data_yeo$x.t, main = "Transformacja Yeo-Johnsona", xlab = NULL, ylab = NULL)

Zależności między zmiennymi

Określenie związku pomiędzy zmiennymi jest istotne, ponieważ pozwala podjąć decyzję, które zmienne uwzględnić w modelu ze względu na ich wzajemne powiązanie, jak i wynik modelowania (tj. które zmienne są najistotniejsze). Oprócz tego, możemy uniknąć redundancji informacji poprzez wykluczenie mocno skorelowanych ze sobą zmiennych, dzięki czemu interpretacja działania modelu będzie prostsza i zwiększy się wydajność obliczeń.

Związek może występować między zmiennymi numerycznymi oraz kategorycznymi. W zależności od ich typu stosuje się różne metody analityczne. W przypadku zmiennych numerycznych do wizualizacji zależności można wykorzystać wykres rozrzutu, który pokaże czy istnieje zależność liniowa lub nieliniowa oraz kierunek (dodatni lub ujemny). Natomiast, podejście ilościowe polega na wykorzystaniu testów korelacji, które umożliwiają wyliczenia współczynnika korelacji (siły) oraz istotności statystycznej.

Testy korelacji

Możemy wyróżnić dwa podstawowe rodzaje testów korelacji:

  1. Parametryczny – zakłada liniową zależność pomiędzy zmiennymi oraz rozkład normalny (test Pearsona). Jest wrażliwy na wartości odstające.
  2. Nieparametryczny – zakłada monotoniczną zależność pomiędzy zmiennymi niekoniecznie liniową (test Spearmana i Kendalla).

Użyjmy przykładowego zestawu danych iris, aby przećwiczyć testowanie korelacji w praktyce.

iris = iris

Analizę zależności pomiędzy zmiennymi możemy rozpocząć od eksploracji wizualnej. W tym celu wykorzystamy funkcję pairs(), która stworzy macierz wykresów rozrzutu wszystkich zmiennych.

pairs(iris[, 1:4], col = iris$Species, pch = 19)

Od razu możemy zauważyć, że wartości zmiennych Petal.Length oraz Petal.Width układają się w rosnącą linię, co może świadczyć o dodatniej korelacji. Wykonajmy odpowiedni test, który potwierdzi naszą hipotezę, używając funkcji cor.test(). W teście Pearsona hipoteza zerowa (\(H_0\)) zakłada, że nie ma liniowej korelacji między dwiema zmiennymi, podczas gdy hipoteza alternatywna (\(H_1\)) zakłada istnienie takiej relacji.

cor.test(iris$Petal.Length, iris$Petal.Width, method = "pearson")

    Pearson's product-moment correlation

data:  iris$Petal.Length and iris$Petal.Width
t = 43.387, df = 148, p-value < 2.2e-16
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 0.9490525 0.9729853
sample estimates:
      cor 
0.9628654

Współczynnik korelacji liniowej Pearsona (cor) wynosi około 0,96, co świadczy o wysokiej dodatniej korelacji. Natomiast, wartość \(p\) poniżej 0,05 sugeruje, że jest ona istotna statystycznie.

cor.test(iris$Sepal.Length, iris$Sepal.Width, method = "pearson")

    Pearson's product-moment correlation

data:  iris$Sepal.Length and iris$Sepal.Width
t = -1.4403, df = 148, p-value = 0.1519
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 -0.27269325  0.04351158
sample estimates:
       cor 
-0.1175698

W przypadku zmiennych Sepal.Length oraz Sepal.Width nie można zauważyć, żeby wartości tworzyły jakiś znaczący wzorzec. Współczynnik korelacji obliczony z testu ma niską wartość wynoszącą około -0,12 oraz nie jest istotny statystycznie, ponieważ wartość \(p\) jest większa od 0,05.

Pamiętaj, że korelacja nie oznacza związku przyczynowo-skutkowego, a istotność statystyczna nie zawsze przekłada się na praktyczne zastosowanie.

Macierz korelacji

Istnieje również możliwość przedstawienia współczynników korelacji w formie wizualnej jako macierz korelacji za pomocą pakietu corrplot.

library("corrplot")
cor_mat = cor(iris[, 1:4], method = "pearson")
corrplot(cor_mat, method = "number")

Rycinę możemy zmodyfikować tak, aby pominąć współczynniki korelacji nieistotne statystycznie i wyświetlić tylko dolną część macierzy.

testRes = cor.mtest(iris[, 1:4], conf.level = 0.95)
corrplot(cor_mat, p.mat = testRes$p, method = "circle", type = "lower",
         insig = "blank", addCoef.col = "black", diag = FALSE)

Zadanie

Dokonaj eksploracyjnej analizy danych dla plików rastrowych (*.tif) zamieszczonych w katalogu dane. Zaprezentuj wczytane zmienne na mapach, a następnie wykorzystaj poznane na zajęciach techniki wizualne oraz opisowe, aby sprawdzić:

  • statystyki opisowe,
  • rozkłady wartości,
  • wartości brakujące,
  • wartości odstające,
  • korelacje pomiędzy zmiennymi.

Obliczenia wykonaj na podstawie losowej próby i napisz jaka wielkość próby została wykorzystana. Pamiętaj, aby ustawić ziarno losowości! Do każdego etapu analizy obowiązkowo dołącz komentarz z wnioskami.

Jako rozwiązanie należy przesłać:

  1. Działający plik źródłowy (np. Quarto .qmd lub RMarkdown .Rmd).
  2. Wygenerowany raport z wynikami przeprowadzonej analizy, rycinami oraz wnioskami (np. .html lub .pdf). Uwaga: Jeśli generujesz raport w formacie .html, to pamiętaj osadzić wszystkie niezbędne pliki.