Geostatystyczna analiza danych

Transformacja danych

Author

Anna Dmowska

1 Pakiety i fukcje pomocnicze

# Wczytanie niezbędnych bibliotek
library(sf)         # Dane wektorowe
library(tmap)       # Tworzenie map
library(gstat)      # Geostatystyka (kriging, variogramy)
library(terra)      # Dane rastrowe
library(tidyverse)  # Przetwarzanie danych i wykresy (ggplot2)
library(patchwork)  # Łączenie wykresów

# Funkcja pomocnicza do interpolacji gstat na obiektach SpatRaster
# Pozwala na użycie gstat bezpośrednio na siatce rastrowej z pakietu terra
interpolate_gstat = function(model, x, ...) {
  # Konwersja współrzędnych rastra na punkty (sf) z zachowaniem układu współrzędnych (CRS)
  v = st_as_sf(x, coords = c("x", "y"), crs = st_crs(model$data[[1]]$data))
  # Wykonanie predykcji (krigingu) dla tych punktów
  p = predict(model, v, ...)
  # Usunięcie geometrii, aby zwrócić czyste wartości tabelaryczne do rastra
  st_drop_geometry(p)
}

# Funkcja do transformacji odwrotnej dla logarytmu (Back-transformation) 
rev_trans = function(pred, var, obs) {
  # 1. Sprawdzenie, czy wejście to SpatRaster (z pakietu terra)
  is_raster = inherits(pred, "SpatRaster")
  
  # 2. Wyciągnięcie wartości do obliczeń (jeśli to raster)
  # Jeśli to punkty, używamy danych bezpośrednio
  p_vals = if(is_raster) terra::values(pred) else pred
  v_vals = if(is_raster) terra::values(var) else var
  
  # 3. Core: Logarytmiczna transformacja odwrotna (anti-log + bias correction)
  # Wzór: exp(pred + var/2)
  bt = exp(p_vals + (v_vals / 2))
  
  # 4. Obliczenie współczynnika k0 (Global Bias Correction)
  mu_bt = mean(bt, na.rm = TRUE)
  mu_original = mean(obs, na.rm = TRUE)
  k0 = mu_original / mu_bt
  
  # 5. Zastosowanie korekty
  btt_values = bt * k0
  
  # 6. Zwrócenie odpowiedniego formatu
  if (is_raster) {
    # Jeśli wejściem był raster, zwracamy nowy raster na bazie szablonu 'pred'
    res_rst = pred
    terra::values(res_rst) = btt_values
    return(res_rst)
  } else {
    # Jeśli wejściem były punkty, zwracamy wektor wartości
    return(btt_values)
  }
}

# Definicja siatki interpolacyjnej (obszar badań, rozdzielczość 100m, układ UTM 33N)
siatka = rast(xmin=616550, xmax=643850, ymin=5785850, ymax = 5822150, res = 100, crs = "EPSG:32633")

2 Przygotowanie danych

# Wczytanie danych z pliku CSV
df = read.csv("data/cw2_dane.csv")

# Konwersja ramki danych na obiekt przestrzenny (sf)
# coords: wskazujemy kolumny z długością (X) i szerokością (Y) geograficzną
# crs: definiujemy układ współrzędnych (EPSG:32633 to UTM strefa 33N)
punkty = st_as_sf(df, coords = c("Xutm", "Yutm"), crs = "EPSG:32633")

# Uśrednienie wartości dla punktów o identycznych współrzędnych
punkty = punkty  |> 
  group_by(geometry) |>  # Grupujemy dane po unikalnych lokalizacjach
  summarise(
    # Obliczamy średnią wartość parametru Ba (bar) dla każdej lokalizacji
    # na.rm = TRUE pozwala zignorować ewentualne braki danych
    Ba = mean(Ba, na.rm = TRUE), 
    
    # Zachowujemy pierwszy napotkany numer ID (Lp) dla porządku
    Lp = first(Lp) 
  ) |> 
  st_cast("POINT") # Upewniamy się, że wynik pozostaje w formacie punktowym

3 Transformacje danych

Transformacje danych są stosowane w geostatystyce głównie w celu poradzenia sobie z silną asymetrią rozkładu oraz obecnością wartości ekstremalnych. W obu przypadkach zastosowanie transformacji ułatwia proces dopasowania modelu struktury przestrzennej (wariogramu).

Większość metod interpolacyjnych, takich jak kriging, najlepiej radzi sobie z danymi o rozkładzie zbliżonym do normalnego.

Istnieje kilka metod transformacji danych, np. metoda Boxa-Coxa, logarytmizacja, normalizacja.

Procedura pracy z transformacją Transformacja danych składa się z trzech kluczowych etapów:

1, Transformacja bezpośrednia: Zamiana wartości oryginalnych na nową skalę (np. logarytmiczną).

Analiza geostatystyczna: Obliczenie semiwariogramu empirycznego i wykonanie estymacji (krigingu) na „nowych” danych.
Transformacja odwrotna (Back-transformation): Przywrócenie wyników estymacji do oryginalnych jednostek (np. z logarytmu z powrotem na mg/kg), aby mapa była interpretowalna.

4 Metody transformacji

4.1 Dane oryginalne

# Wizualizacja rozkładu wyjściowego
ggplot(punkty, aes(x = Ba)) + 
  geom_histogram() + 
  labs(title = "Dane oryginalne") + 
  theme_bw()

4.2 Transformacja logarytmiczna

Jest najstarszą i najczęściej stosowaną metodą, gdy dane mają rozkład log-normalny (rozkład asymetryczny prawostronny: wiele małych wartości, kilka bardzo dużych).

Wada: Wymaga transformacji wstecznej (tzw. back-transformation), która jest wrażliwa na wariancję (może zawyżać średnią). Jeśli dane zawierają zera, trzeba dodać stałą, np. $log(x + 1)$.

# Tworzymy nową kolumnę z logarytmem naturalnym
punkty$Ba_log = log(punkty$Ba)

4.3 Transformacja Boxa-Coxa

Transformacja Boxa-Coxa stanowi uogólnienie transformacji potęgowej. Można ją stosować wyłącznie dla dodatnich danych. Jesto to bardziej zaawansowana metoda, która automatycznie szuka najlepszego wariantu transformacji dobierając parametr $\lambda$, który najlepiej przybliża dane do rozkładu normalnego.

\[y = \frac{x^\lambda - 1}{\lambda}\]

(dla $\lambda \neq 0$) lub $log(x)$ (dla $\lambda = 0$).

Parametr $\lambda$	Rozdzaj transformacji
-1	1/x
0	Logarytm naturalny ($ln(x)$)
0.5	Pierwiastek kwadratowy ($\sqrt{x}$)
1	Brak zmian, dane mają rozkład normalny
2	Kwadrat ($x^2$)

library(bestNormalize)

# Budujemy model transformacji Box-Cox
# Funkcja automatycznie estymuje najlepsze lambda
bc_fit <- boxcox(punkty$Ba)

# Sprawdzamy wybrane parametry (np. wartość lambda)
print(bc_fit)

Standardized Box Cox Transformation with 978 nonmissing obs.:
 Estimated statistics:
 - lambda = -0.2876564 
 - mean (before standardization) = 2.12854 
 - sd (before standardization) = 0.2566151

# Przypisujemy przetransformowane wartości (x.t oznacza x-transformed)
punkty$Ba_bc <- bc_fit$x.t

4.4 Normalizacja (NST - Normal Score Transformation NST)

Polega na przypisaniu rang do danych empirycznych a następnie mapowaniu kwantyli rozkładu empirycznego na kwantyle rozkładu normalnego $N(0,1)$.

Zaleta: Zawsze daje idealny rozkład normalny. Bardzo dobrze radzi sobie z ekstremalnymi wartościami odstającymi.

# NST zawsze wymusza rozkład normalny, niezależnie od kształtu oryginału
# Używamy funkcji orderNorm (NST) z pakietu bestNormalize
nst_fit <- orderNorm(punkty$Ba)

# Wyciągnięcie wyników (Z-scores)
punkty$Ba_nst <- nst_fit$x.t

# Informacja o modelu (przydatna do późniejszej transformacji odwrotnej)
print(nst_fit)

orderNorm Transformation with 978 nonmissing obs and ties
 - 134 unique values 
 - Original quantiles:
  0%  25%  50%  75% 100% 
   4   19   27   39  940

4.5 Porównanie metod transformacji

Poniższe wykres porównują rozkład zmiennej oryginalnej z rozkładem po transformacji

# Przygotowanie wykresów do zestawienia
d1 = ggplot(punkty, aes(x = Ba)) + geom_histogram() + labs(title = "Dane oryginalne") + theme_bw()
d2 = ggplot(punkty, aes(x = Ba_nst)) + geom_histogram() + labs(title = "NST (Normal Score)") + theme_bw()
d3 = ggplot(punkty, aes(x = Ba_bc)) + geom_histogram() + labs(title = "Box-Cox") + theme_bw()
d4 = ggplot(punkty, aes(x = Ba_log)) + geom_histogram() + labs(title = "Logarytm") + theme_bw()

# Wyświetlenie wszystkich wykresów na jednym panelu
library(patchwork)
(d1 + d2) / (d3 + d4)

4.6 Wybór najlepszej metody transformacji

Pakiet bestNormalize pozwala na automatyczne przetestowanie wbudowanych metod transformacji i wybór tej, która najlepiej sprawdzi się w przypadku analizowanych danych.

# Automatyczny wybór najlepszej metody
best_fit <- bestNormalize(punkty$Ba)

# Najlepsza metoda (często jest to Box-Cox lub orderNorm/NST)
print(best_fit)

Best Normalizing transformation with 978 Observations
 Estimated Normality Statistics (Pearson P / df, lower => more normal):
 - arcsinh(x): 1.8341
 - Box-Cox: 1.422
 - Center+scale: 16.8453
 - Double Reversed Log_b(x+a): 20.3853
 - Log_b(x+a): 1.8426
 - orderNorm (ORQ): 1.2061
 - sqrt(x + a): 5.2524
 - Yeo-Johnson: 1.4523
Estimation method: Out-of-sample via CV with 10 folds and 5 repeats
 
Based off these, bestNormalize chose:
orderNorm Transformation with 978 nonmissing obs and ties
 - 134 unique values 
 - Original quantiles:
  0%  25%  50%  75% 100% 
   4   19   27   39  940

5 Przykład: Transformacja logarytmiczna

5.1 Dane oryginalne

Rozkład danych

ggplot(punkty, aes(x = Ba)) + 
  geom_histogram() + 
  labs(title = "Dane oryginalne", subtitle = "Rozkład wartości zmiennej Ba wraz z uwzględnieniem wartości anomalnych") + 
  theme_bw()

Semiwariogram empiryczny na danych oryginalnych

vario = variogram(Ba ~ 1, locations = punkty, cutoff = 15000, width = 400)
plot(vario)

5.2 Transformacja logarytmiczna

W pierwszym kroku tworzymy semiwariogram, jednakże zamiast używać bezpośrednio wartości zmiennej poddajemy ją logaritmizacji poprzez log(Ba). Możemy wcześniej do zbioru danych dodać nową zmienną zawierającą wartość zlogarytmizowaną lub zastosować logarytm bezpośrednio w funkcji tworzącej semiwarogram (formuła log(Ba)~1).

# Stworzenie semiwariogramu empirycznego z wykorzystaniem funkcji variogram(). 
vario_log = variogram(log(Ba) ~ 1, locations = punkty, cutoff = 15000, width = 400)
plot(vario_log)

W kolejnym kroku dopasowujemy model do utworzonego semiwariogramu empirycznego.

# Dopasownie modelu do semiwariogramu empirycznnego
model_log = vgm(model = "Sph", range = 1700, psill = 0.15, nugget = 0.2,
               add.to = vgm(model = "Gau", psill = 0.35, range = 10000))
  
plot(vario_log, model_log)

Następnie wykonujemy estymację wybierając jedną z metod krigingu. W poniższym przykładzie estymacja jest wykonywana metodą krigingu prostego. Kriging prosty wymaga podania wartości średniej arytmetycznej - w przypadku wykonywania semiwariogramu na wartościach zlogarytmizowanych, średnia także musi zostać obliczona dla zlogarytmizowanych wartości.

log_srednia = mean(log(punkty$Ba), na.rm = TRUE)

Aby móc porównać wartości zmiennej Ba przed, po transformacji logarytmicznej, oraz po transformacji odwrotnej wykonamy kriging dla zbioru punktów, a nie dla siatki. W ten sposób dla każdego punktu będziemy dysponować 4 wartościami:

oryginalną wartością zmiennej Ba (punkty$Ba),
zlogarytmizowaną wartością zmiennej Ba (punkty$log_Ba),
estymowaną wartością w skali logarytmicznej (sk$var1.pred),
estymowaną wartością po transformacji odwrotnej (sk$var1.pred2). Do wykonania estymacji na zbiorze punktowym służy funkcja predict(), która wymaga podania parametrów kriging (obiekt sk_param) oraz zbioru punktowego, dla którego wykonamy estymację (punkty).

## Zdefiniowanie parametrów krigingu - kriging prosty (Simple Krigig SK)
sklog_param = gstat(formula = log(Ba) ~ 1, 
            locations = punkty,
            model = model_log,
            beta = mean(punkty$Ba, na.rm = TRUE))   

## Wykonanie estymacji dla punktów      
sklog = predict(sklog_param, punkty)

[using simple kriging]

5.2.1 Transformacja odwrotna

Przywrócenie wartości do oryginalnej jednostki, np. po estymacji, wymaga zastosowania odpowiedniej metody transformacji odwrotnej. Dla transformacji logarytmicznej wzór na transformację odwrotną jest następujący:

\[y = k_0 \cdot exp[ln(\hat{y}_{OK})+\frac{\sigma^2_{OK}}{2}]\]

, gdzie

$k_0$ - współczynnik korekcyjny obliczany jako iloraz średniej oryginalnej i średniej po transformacji wstecznej
$exp$ - funkcja wykładnicza
$ln$ - logarytm naturalny
$\hat{y}_{OK}$ - estymacja krigingu zwykłego
$\sigma^2_{OK}$ - wariancja krigingu zwykłego

Funkcja rev_trans() stanowi implementację powyższego wzoru w R. Funckja przyjmuje 3 argumenty:

pred - estymowane wartości, których jednostką jest logarytm oryginalnej jednostki
var - wariancja kolejnych estymacji
obs - wartości pomiarów w punktach w oryginalnej jednostce

rev_trans = function(pred, var, obs) {
  # 1. Sprawdzenie, czy wejście to SpatRaster (z pakietu terra)
  is_raster = inherits(pred, "SpatRaster")
  
  # 2. Wyciągnięcie wartości do obliczeń (jeśli to raster)
  # Jeśli to punkty, używamy danych bezpośrednio
  p_vals = if(is_raster) terra::values(pred) else pred
  v_vals = if(is_raster) terra::values(var) else var
  
  # 3. Core: Logarytmiczna transformacja odwrotna (anti-log + bias correction)
  # Wzór: exp(pred + var/2)
  bt = exp(p_vals + (v_vals / 2))
  
  # 4. Obliczenie współczynnika k0 (Global Bias Correction)
  mu_bt = mean(bt, na.rm = TRUE)
  mu_original = mean(obs, na.rm = TRUE)
  k0 = mu_original / mu_bt
  
  # 5. Zastosowanie korekty
  btt_values = bt * k0
}

Do wykonania transformacji odwrotnej wykorzystujemy funkcję rev_trans(), która została powyżej zdefiniowana. Funkcja rev_trans() przyjmuje 3 argumenty: zlogaritmizowaną estymację, wariancję estymacji i oryginalne wartości pomiarów w punktach.

sklog$var1.pred2 = rev_trans(sklog$var1.pred, sklog$var1.var, punkty$Ba)

Poniższa rycina przedstawia porównanie wartości zmiennej oryginalnej, zlogarytmizowanej, oraz wartości estymowanych: zlogarytmizowanych i po transformacji odwrotnej.

punkty$log_Ba = log(punkty$Ba)
d1 = ggplot(punkty, aes(x = Ba)) + geom_density() + labs(title = "Ba - zmienna oryginalna") + theme_bw()
d2 = ggplot(sklog, aes(x = var1.pred2)) + geom_density() + labs(title = "Ba - wartości estymowane po transformacji odwrotnej") + theme_bw()
d3 = ggplot(punkty, aes(x = log_Ba)) + geom_density() + labs(title = "Ba - zmienna zlogarytmizowana (ln(Ba))") + theme_bw()
d4 = ggplot(sklog, aes(x = var1.pred)) + geom_density() + labs(title = "Ba - wartości estymowane zlogarytmizowane") + theme_bw()

d1/d2/d3/d4 + plot_layout(heights = unit(c(3, 3,3,3), c('cm', 'cm', 'cm', 'cm')))

Poniżej wykonamy także estymację z wykorzystaniem krigingu prostego w siatce.

siatka = rast(xmin=616550, xmax=643850, ymin=5785850, ymax = 5822150, res = 100, crs = "EPSG:32633")

#parametry krigingu SK              
sks_param = gstat(formula = log(Ba) ~ 1, 
            locations = punkty,
            model = model_log,
            beta = mean(log(punkty$Ba), na.rm = TRUE))

#interpolacja
sks = interpolate(siatka, sks_param, fun = interpolate_gstat)

[using simple kriging]
[using simple kriging]

#transformacja odwrotna 
sks$var1.pred2 = rev_trans(sks$var1.pred, sks$var1.var, punkty$Ba)

Poniżej wykonamy także estymację z wykorzystaniem krigingu zwykłego w siatce.

oklog_param = gstat(formula = log(Ba) ~ 1, 
            locations = punkty,
            model = model_log,
            nmax = 30)

#interpolacja         
oklog = interpolate(siatka, oklog_param, fun = interpolate_gstat)

[using ordinary kriging]
[using ordinary kriging]

#transformacja odwrotna 
oklog$var1.pred2 = rev_trans(oklog$var1.pred, oklog$var1.var, punkty$Ba)

6 Przykład: Normalizacja (NST - Normal Score Transformation)

W pierwszym kroku transformujemy dane wykorzystując funkcje pakietu bestNormalize.

library(bestNormalize)
# 1. Budowa modelu NST (orderNorm)
nst_fit <- orderNorm(punkty$Ba)
# 2. Wyciągnięcie przetransformowanych danych (Z-scores)
punkty$Ba_nst <- nst_fit$x.t
# Sprawdzenie wyniku
print(nst_fit)

orderNorm Transformation with 978 nonmissing obs and ties
 - 134 unique values 
 - Original quantiles:
  0%  25%  50%  75% 100% 
   4   19   27   39  940

Następnie tworzymy semiwariogram empiryczny dla danych po transformacji.

# Stworzenie semiwariogramu empirycznego z wykorzystaniem funkcji variogram(). 
vario_nst = variogram(Ba_nst ~ 1, locations = punkty, cutoff = 15000, width = 400)
plot(vario_nst)

W kolejnym kroku dopasowujemy model do utworzonego semiwariogramu empirycznego.

# Dopasownie modelu do semiwariogramu empirycznnego
model_nst = vgm(model = "Sph", range = 1700, psill = 0.55, nugget = 0.2,
               add.to = vgm(model = "Gau", psill = 0.55, range = 11000))
  
plot(vario_nst, model_nst)

Następnie wykonujemy estymację metodą krigingu zwykłego.

oknst_param = gstat(formula = Ba_nst ~ 1, 
            locations = punkty,
            model = model_nst,
            nmax = 30)

#interpolacja         
oknst = interpolate(siatka, oknst_param, fun = interpolate_gstat)

[using ordinary kriging]
[using ordinary kriging]

Kolejnym krokiem jest wykonanie transformacji odwrotnej wykorzystując w tym celu funkcję predict() z argumentem inverse = TRUE.

#transformacja odwrotna 
# 1. Pobieramy wartości z rastra
pred_nst <- values(oknst$var1.pred)

# 2. Transformacja wsteczna do oryginalnej skali
pred_original <- predict(nst_fit, newdata = pred_nst, inverse = TRUE)

# 3. Włożenie wartości z powrotem do rastra
oknst_2 <- oknst$var1.pred
values(oknst_2) <- pred_original

6.1 Porównanie semiwariogramów po transformacji

vario_log$wariant = "log"
vario_nst$wariant = "nst"
v_comb = rbind(vario_log, vario_nst)
ggplot(v_comb, aes(x = dist, y = gamma, color = wariant)) + 
  geom_point() + geom_line() +   
  expand_limits(y = 0) + 
  theme_bw()

Dla porównania poniżej przedstawiono semiwariogram empiryczny wykonany w oparciu o dane po wykluczeniu anomalii.

punkty_bez_anomalii <- punkty[!punkty$Lp%in%c(434, 557, 690, 731, 973,974,975,976,977),]
vario_bez_anomalii = variogram(Ba ~ 1, locations = punkty_bez_anomalii, width = 400, cutoff = 15000)

ggplot(vario_bez_anomalii, aes(x = dist, y = gamma)) + 
  geom_point() + geom_line() + 
  expand_limits(y = 0) + 
  theme_bw()

Możemy zwrócić uwagę, że po zastosowaniu transformacji logarytmicznej nadal mamy doczynienia z pierwszym odstającym punktem, który należy pominąć przy dopasowywaniu modelu. Semiwariogram empiryczny wykonany po wykonaniu normalizacji (NST) jest mniej chaotyczny i bardziej wygładzony, ułatwiając tym samym dopasowanie modelu.

7 Porównanie wyników krigingu

Kriging OK wykonany na podstawie danych po wykluczeniu anomalii.

punkty_bez_anomalii <- punkty[!punkty$Lp%in%c(434, 557, 690, 731, 973,974,975,976,977),]
vario_bez_anomalii = variogram(Ba ~ 1, locations = punkty_bez_anomalii, width = 400, cutoff = 15000)
model_zl = vgm(model = "Sph", range = 1500, psill = 1000, nugget = 30,
               add.to = vgm(model = "Gau", psill = 600, range = 10000))
plot(vario_bez_anomalii, model_zl)

ok_param = gstat(formula = Ba ~ 1, 
            locations = punkty_bez_anomalii,
            model = model_zl,
            nmax = 30)

#interpolacja         
ok = interpolate(siatka, ok_param, fun = interpolate_gstat)

[using ordinary kriging]
[using ordinary kriging]

# Mapa dla Krigingu Zwyczajnego (bez anomalii)
mapa_ok = tm_shape(ok) +
  tm_raster(col = "var1.pred", 
            col.scale = tm_scale_continuous(values = "-brewer.spectral", midpoint = NA)) +
  tm_title("OK: Bez anomalii")

# Mapa dla Krigingu Zwyczajnego po logarytmizacji
mapa_oklog = tm_shape(oklog) +
  tm_raster(col = "var1.pred2", 
            col.scale = tm_scale_continuous(values = "-Spectral")) +
  tm_title("OK: Logarytmizacja")

# Mapa dla Krigingu Zwyczajnego po transformacji NST
mapa_oknst = tm_shape(oknst_2) +
  tm_raster(col = "var1.pred", 
            col.scale = tm_scale_continuous(values = "-Spectral")) +
  tm_title("OK: Normal Score (NST)")

tmap_arrange(mapa_ok, mapa_oklog, mapa_oknst, ncol = 2)

Parametr $\lambda$	Rozdzaj transformacji
-1	1/x
0	Logarytm naturalny (\(ln(x)\))
0.5	Pierwiastek kwadratowy (\(\sqrt{x}\))
1	Brak zmian, dane mają rozkład normalny
2	Kwadrat (\(x^2\))