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.280 15.310 15.950 18.270 26.140
punkty$temp_ind <- punkty$temp > 20
summary(punkty$temp_ind)
## Mode FALSE TRUE NA's
## logical 206 44 0
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]