9 Estymacja lokalnego rozkładu prawdopodobieństwa

library("ggplot2")
library("gstat")
library("sp")
library("geostatbook")
data(punkty)
data(siatka)

9.1 Kriging danych kodowanych

9.1.1 Kriging danych kodowanych (ang. Indicator kriging)

Kriging danych kodowanych to metoda krigingu oparta o dane kategoryzowane lub też dane przetworzone z postaci ciągłej do binarnej. Jest ona zazwyczaj używana jest to oszacowania prawdopodobieństwa przekroczenia zdefiniowanej wartości progowej, może być również używana do szacowania wartości z całego rozkładu. Wartości danych wykorzystywane do krigingu danych kodowanych są określone jako 0 lub 1, co reprezentuje czy wartość danej zmiennej jest powyżej czy poniżej określonego progu.

9.1.2 Wady i zalety krigingu danych kodowanych

Zalety:

  • Możliwość zastosowania, gdy nie interesuje nas konkretna wartość, ale znalezienie obszarów o wartości przekraczającej dany poziom
  • Nie jest istotny kształt rozkładu

Wady:

  • Potencjalnie trudne do modelowania semiwariogramy (szczególnie skrajnych przedziałów)
  • Czasochłonność/pracochłonność

9.2 Kriging danych kodowanych | Przykłady

9.2.1 Binaryzacja danych

Pierwszym krokiem w krigingu danych kodowanych jest stworzenie zmiennej binarnej. Na poniższym przykładzie tworzona jest nowa zmienna temp_ind. Przyjmuje ona wartość TRUE (czyli 1) dla pomiarów temperatury wyższych niż 20 stopni Celsjusza, a dla pomiarów równych i niższych niż 20 stopni Celsjusza jej wartość wynosi FALSE (czyli 0).

summary(punkty$temp)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   8.706  13.284  15.309  15.950  18.273  26.139
punkty$temp_ind <- punkty$temp > 20
summary(punkty$temp_ind)
##    Mode   FALSE    TRUE 
## logical     206      44

W przykładzie, próg został wyznaczony arbitralnie. Istnieje oczywiście szereg innych możliwości wyznaczania progu. Można wykorzystać wiedzę zewnętrzną (np. toksyczne stężenie analizowanej substancji) lub też posłużyć się wykresem dystrybuanty do określenia istotnej zmiany wartości.

ggplot(punkty@data, aes(temp)) + stat_ecdf()

9.2.2 Kriging danych kodowanych (ang. Indicator kriging) | Modelowanie

Tworzenie i modelowanie semiwariogramu empirycznego w krigingu danych kodowanych wygląda tak samo jak, np. w przypadku krigingu zwykłego. Korzystając z funkcji variogram() tworzony jest semiwariogram empiryczny, używając vgm() tworzony jest model “ręczny”, który następnie jest dopasowywany z użyciem funkcji fit.variogram().

vario_ind <- variogram(temp_ind~1, punkty)
plot(vario_ind)

model_ind <- vgm(0.09, model = "Sph", range = 2000, nugget = 0.01)
plot(vario_ind, model = model_ind)

fitted_ind <- fit.variogram(vario_ind, model_ind)
fitted_ind
##   model      psill   range
## 1   Nug 0.01761559    0.00
## 2   Sph 0.08607730 1919.46
plot(vario_ind, model = fitted_ind)

9.2.3 Kriging danych kodowanych (ang. Indicator kriging) | Interpolacja

Ostatnim etapem jest stworzenie interpolacji geostatystycznej z pomocą funkcji krige. Wymaga ona czterech argumentów - wzoru (temp_ind~1), zbioru punktowego (punkty), siatki do interpolacji (siatka) oraz modelu (fitted_ind).

ik <- krige(temp_ind~1, punkty, siatka, model = fitted_ind)
## [using ordinary kriging]

W wyniku estymacji otrzymuje się mapę przestawiającą prawdopodobieństwo przekroczenia zadanej wartości (w tym wypadku jest to 20 stopni Celsjusza) oraz uzyskaną wariancję predykcji.

spplot(ik, "var1.pred")
spplot(ik, "var1.var")

9.2.4 Kriging danych kodowanych (ang. Indicator kriging)

Alternatywnie, zamiast tworzenia nowej zmiennej (takiej jak temp_ind), można wykorzystać funkcję I. Z jej użyciem można definiować przyjęte progi bezpośrednio do funkcji variogram i krige. Na poniższych przykładach w ten sposób ustalono trzy progi - poniżej 20, poniżej 16, oraz poniżej 12 stopni Celsjusza.

vario_ind20 <- variogram(I(temp < 20)~1, punkty)
fitted_ind20 <- fit.variogram(vario_ind20, vgm(0.08, "Sph", 3500, nugget = 0.01))
vario_ind16 <- variogram(I(temp < 16)~1, punkty)
fitted_ind16 <- fit.variogram(vario_ind16, vgm(0.18, "Sph", 3500, nugget = 0.03))
vario_ind12 <- variogram(I(temp < 12)~1, punkty)
fitted_ind12 <- fit.variogram(vario_ind12, vgm(0.13, "Sph", 2000, nugget = 0.03))

ik20 <- krige(I(temp < 20)~1, punkty, siatka, model = fitted_ind20, nmax = 30)
## [using ordinary kriging]
ik16 <- krige(I(temp < 16)~1, punkty, siatka, model = fitted_ind16, nmax = 30)
## [using ordinary kriging]
ik12 <- krige(I(temp < 12)~1, punkty, siatka, model = fitted_ind12, nmax = 30)
## [using ordinary kriging]