-
Notifications
You must be signed in to change notification settings - Fork 7
/
12-ocena.Rmd
410 lines (316 loc) · 15.4 KB
/
12-ocena.Rmd
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
# Ocena jakości estymacji {#ocena-jakosci-estymacji}
Odtworzenie obliczeń z tego rozdziału wymaga załączenia poniższych pakietów oraz wczytania poniższych danych:
```{r 12-ocena-1, message=FALSE, warning=FALSE}
library(sf)
library(stars)
library(gstat)
library(tmap)
library(rsample)
library(ggplot2)
library(geostatbook)
data(punkty)
data(siatka)
```
```{r 12-ocena-2, echo=FALSE}
par(mar = c(rep(0, 4)))
```
## Wizualizacja jakości estymacji
### Mapa
Do oceny przestrzennej jakości estymacji można zastosować mapę przestawiającą błędy estymacji (rycina \@ref(fig:mapa)).
Wyliczenie błędów estymacji odbywa się poprzez odjęcie od obserwacji wartości estymowanej.
```{r 12-ocena-3, eval=FALSE}
blad_estymacji = obserwowane - estymacja
```
```{r mapa, echo=FALSE, message=FALSE, fig.cap = "Mapa przedstawiająca rozkład przestrzenny błędów estymacji."}
set.seed(224)
punkty_podzial = initial_split(punkty, prop = 0.75, strata = temp)
train = training(punkty_podzial)
test = testing(punkty_podzial)
vario = variogram(temp ~ 1, data = train)
model = vgm(10, model = "Sph", range = 4000, nugget = 0.5)
fitted = fit.variogram(vario, model)
test_sk = krige(temp ~ 1, train, newdata = test,
model = fitted, beta = 15, debug.level = 0)
blad_estymacji_sk = test$temp - test_sk$var1.pred
test_sk$blad_estymacji_sk = blad_estymacji_sk
test_sk$true = test$temp
cuts = c(-5, -3, -1, 1, 3, 5)
tm_shape(test_sk) +
tm_symbols(col = "blad_estymacji_sk", breaks = cuts, title.col = "") +
tm_layout(main.title = "Błąd estymacji", legend.outside = TRUE)
```
### Histogram
Błąd estymacji można również przedstawić na wykresach, między innymi, na histogramie (rycina \@ref(fig:hist)).
```{r hist, echo=FALSE, fig.cap="Rozkład wartości błędów estymacji."}
ggplot(test_sk, aes(blad_estymacji_sk)) +
geom_histogram() +
xlab("Błąd estymacji") +
ylab("Liczebność")
```
### Wykres rozrzutu
Do porównania pomiędzy wartością estymowaną a obserwowaną może również posłużyć wykres rozrzutu (rycina \@ref(fig:point)).
```{r point, echo=FALSE, fig.cap="Wykres rozrzutu porównujący rzeczywiste (obserwowane) wartości i wartości estymowane."}
ggplot(test_sk, aes(var1.pred, true)) +
geom_point() +
xlab("Estymacja") +
ylab("Obserwacja")
```
## Statystyki jakości estymacji
### Podstawowe statystyki
W momencie, gdy trzeba określić jakość estymacji lub porównać wyniki pomiędzy estymacjami należy zastosować tzw. statystyki jakości estymacji.
Do podstawowych statystyk ocen jakości estymacji należą:
- Średni błąd estymacji (MPE, ang. *mean prediction error*).
- Pierwiastek błędu średniokwadratowego (RMSE, ang. *root mean square error*).
- Współczynnik determinacji (R^2^, ang. *coefficient of determination*).
Idealna estymacja dawałaby brak błędu oraz współczynnik determinacji pomiędzy pomiarami (całą populacją) i szacunkiem równy 1.
Należy jednak zdawać sobie sprawę, że dane wejściowe są obarczone błędem/niepewnością, dlatego też w praktyce idealna estymacja nie jest osiągana.
Wysokie, pojedyncze wartości błędu mogą świadczyć, np. o wystąpieniu wartości odstających.
### Średni błąd estymacji
Średni błąd estymacji (MPE) można wyliczyć korzystając z poniższego wzoru:
$$ MPE=\frac{\sum_{i=1}^{n}(v_i - \hat{v}_i)}{n} $$
, gdzie $v_i$ to wartość obserwowana a $\hat{v}_i$ to wartość estymowana.
Średni błąd estymacji można też wyliczyć używając funkcji `mean()` w R.
```{r 12-ocena-4, eval=FALSE}
MPE = mean(obserwowane - estymacja)
```
Optymalnie wartość średniego błędu estymacji powinna być jak najbliżej 0.
### Pierwiastek błędu średniokwadratowego
Pierwiastek błędu średniokwadratowego (RMSE) jest możliwy do wyliczenia poprzez wzór:
$$ RMSE=\sqrt{\frac{\sum_{i=1}^{n}(v_i-\hat{v}_i)^2}{n}} $$
, gdzie $v_i$ to wartość obserwowana a $\hat{v}_i$ to wartość estymowana.
RMSE można też wyliczyć w R.
```{r 12-ocena-5, eval=FALSE}
RMSE = sqrt(mean((obserwowane - estymacja) ^ 2))
```
Optymalnie wartość pierwiastka błędu średniokwadratowego powinna być jak najmniejsza.
### Współczynnik determinacji
Współczynnik determinacji (R^2^) jest możliwy do wyliczenia poprzez wzór:
$$ R^2 = 1 - \frac{\sum_{i=1}^{n} (\hat v_i - v_i)^2}{\sum_{i=1}^{n} (v_i - \overline{v_i})^2} $$
, gdzie $v_i$ to wartość obserwowana, $\hat{v}_i$ to wartość estymowana, a $\overline{v}$ średnia arytmetyczna wartości obserwowanych.
R^2^ można też wyliczyć w R.
```{r 12-ocena-6, eval=FALSE}
R2 = 1 - sum((estymacja - obserwowane) ^ 2) /
sum((obserwowane - mean(obserwowane)) ^ 2)
```
lub
```{r 12-ocena-7, eval=FALSE}
R2 = cor(obserwowane, estymacja) ^ 2
```
Współczynnik determinacji przyjmuje wartości od 0 do 1, gdzie model jest lepszy im wartość tego współczynnika jest bliższa jedności.
## Jakość wyników estymacji
### Walidacja wyników estymacji
Dokładne dopasowanie modelu do danych może w efekcie nie dawać najlepszych wyników.
Szczególnie będzie to widoczne w przypadku modelowania, w którym dane obarczone są znacznym szumem (zawierają wyraźny błąd) lub też posiadają kilka wartości odstających.
W efekcie ważne jest stosowanie metod pozwalających na wybranie optymalnego modelu.
Do takich metod należy, między innymi, walidacja podzbiorem (ang. *jackknifing*) oraz kroswalidacja (ang. *crossvalidation*).
### Walidacja podzbiorem
Walidacja podzbiorem polega na podziale zbioru danych na dwa podzbiory - treningowy i testowy.
Zbiór treningowy służy do stworzenia semiwariogramu empirycznego, zbudowania modelu oraz estymacji wartości.
Następnie wynik estymacji porównywany jest z rzeczywistymi wartościami ze zbioru testowego.
Zaletą tego podejścia jest stosowanie danych niezależnych od estymacji do oceny jakości modelu.
Wadą natomiast jest konieczność posiadania (relatywnie) dużego zbioru danych.
Na poniższym przykładzie zbiór danych dzielony jest używając funkcji `initial_split()` z pakietu **rsample**.
Użycie tej funkcji powoduje obiektu zawierającego wydzielony podzbiór treningowy i testowy.
Ważną zaletą funkcji `initial_split()` jest to, iż w zbiorze treningowym i testowym zachowane są podobne rozkłady wartości.
W przykładzie użyto argumentu `prop = 0.75`, który oznacza, że 75% danych będzie należało do zbioru treningowego, a 25% do zbioru testowego.
Następnie, korzystając ze stworzonego obiektu, budowane są dwa zbiory danych - treningowy (`train`) oraz testowy (`test`).
```{r 12-ocena-8}
set.seed(224)
punkty_podzial = initial_split(punkty, prop = 0.75, strata = temp)
train = training(punkty_podzial)
test = testing(punkty_podzial)
```
Dalszym krokiem jest stworzenie semiwariogramu empirycznego oraz jego modelowanie w oparciu o zbiór treningowy (rycina \@ref(fig:12-ocena-9)).
```{r 12-ocena-9, fig.cap = "Model semiwariogramu w oparciu o zbiór testowy."}
vario = variogram(temp ~ 1, locations = train)
model = vgm(model = "Sph", nugget = 0.5)
fitted = fit.variogram(vario, model)
plot(vario, model = fitted)
```
Do porównania wyników estymacji w stosunku do zbioru testowego posłuży funkcja `krige()`.
Wcześniej wymagała ona podania wzoru, zbioru punktowego, siatki oraz modelu.
W tym przypadku jednak chcemy porównać wynik estymacji i testowy zbiór punktowy.
Dlatego też, zamiast obiektu siatki definiujemy obiekt zawierający zbiór testowy (`test`).
```{r 12-ocena-10}
test_sk = krige(temp ~ 1,
locations = train,
newdata = test,
model = fitted,
beta = 15)
summary(test_sk)
```
Uzyskane w ten sposób wyniki możemy określić używając statystyk jakości estymacji lub też wykresów (ryciny \@ref(fig:12-ocena-12), \@ref(fig:12-ocena-13), \@ref(fig:12-ocena-14)).
```{r 12-ocena-11}
blad_estymacji_sk = test$temp - test_sk$var1.pred
summary(blad_estymacji_sk)
MPE = mean(test$temp - test_sk$var1.pred)
MPE
RMSE = sqrt(mean((test$temp - test_sk$var1.pred) ^ 2))
RMSE
R2 = cor(test$temp, test_sk$var1.pred) ^ 2
R2
```
```{r 12-ocena-12, fig.cap = "Mapa przedstawiająca rozkład przestrzenny błędów estymacji uzyskanych używając krigingu prostego."}
test_sk$blad_estymacji_sk = blad_estymacji_sk
tm_shape(test_sk) +
tm_symbols(col = "blad_estymacji_sk", title.col = "") +
tm_layout(main.title = "Błąd estymacji", legend.outside = TRUE)
```
```{r 12-ocena-13, fig.cap="Rozkład wartości błędów estymacji uzyskanych używając krigingu prostego."}
ggplot(test_sk, aes(blad_estymacji_sk)) +
geom_histogram() +
xlab("Błąd estymacji") +
ylab("Liczebność")
```
```{r 12-ocena-14, fig.cap="Wykres rozrzutu porównujący rzeczywiste (obserwowane) wartości i wartości estymowane uzyskane używając krigingu prostego."}
test_sk$true = test$temp
ggplot(test_sk, aes(var1.pred, true)) +
geom_point() +
xlab("Estymacja") +
ylab("Obserwacja")
```
W sytuacji, gdy uzyskany model jest wystarczająco dobry, możemy również uzyskać estymację dla całego obszaru z użyciem funkcji `krige()`, tym razem jednak podając obiekt siatki (rycina \@ref(fig:12-ocena-15)).
```{r 12-ocena-15, fig.cap="Estymacja uzyskana dla całego obszaru."}
test_sk = krige(temp ~ 1,
locations = train,
newdata = siatka,
model = fitted,
beta = 15)
tm_shape(test_sk) +
tm_raster("var1.pred", style = "cont", palette = "-Spectral")
```
### Kroswalidacja
W przypadku kroswalidacji te same dane wykorzystywane są do budowy modelu, estymacji, a następnie do oceny prognozy.
Procedura kroswalidacji LOO (ang. *leave-one-out cross-validation*) składa się z poniższych kroków:
1. Zbudowanie matematycznego modelu z dostępnych obserwacji.
2. Dla każdej znanej obserwacji następuje:
- Usunięcie jej ze zbioru danych.
- Użycie modelu do wykonania estymacji w miejscu tej obserwacji.
- Wyliczenie reszty (ang. *residual*), czyli różnicy pomiędzy znaną wartością a estymacją.
3. Podsumowanie otrzymanych wyników.
W pakiecie **gstat**, kroswalidacja LOO jest dostępna w funkcjach `krige.cv()` oraz `gstat.cv()`.
Działają one bardzo podobnie jak funkcje `krige()` oraz `gstat()`, jednak w przeciwieństwie do nich nie wymagają podania obiektu siatki.
```{r loovv, eval=FALSE}
vario = variogram(temp ~ 1, data = punkty)
model = vgm(model = "Sph", nugget = 0.5)
fitted = fit.variogram(vario, model)
cv_sk = krige.cv(temp ~ 1,
locations = punkty,
model = fitted,
beta = 15)
summary(cv_sk)
```
```{r trueloovv, echo=FALSE, cache=TRUE}
vario = variogram(temp ~ 1, data = punkty)
model = vgm(model = "Sph", nugget = 0.5)
fitted = fit.variogram(vario, model)
cv_sk = krige.cv(temp ~ 1,
locations = punkty,
model = fitted,
beta = 15,
verbose = FALSE)
summary(cv_sk)
```
Uzyskane w ten sposób wyniki możemy określić używając statystyk jakości estymacji lub też wykresów (ryciny \@ref(fig:12-ocena-17), \@ref(fig:12-ocena-18), \@ref(fig:12-ocena-19)).
```{r 12-ocena-16}
summary(cv_sk$residual)
MPE = mean(cv_sk$residual)
MPE
RMSE = sqrt(mean((cv_sk$residual) ^ 2))
RMSE
R2 = cor(cv_sk$observed, cv_sk$var1.pred) ^ 2
R2
```
```{r 12-ocena-17, fig.cap = "Mapa przedstawiająca rozkład przestrzenny błędów estymacji uzyskanych używając krigingu prostego i kroswalidacji."}
tm_shape(cv_sk) +
tm_symbols(col = "residual", title.col = "",
breaks = c(-10, -5, -3, -1, 1, 3, 5, 10)) +
tm_layout(main.title = "Błąd estymacji", legend.outside = TRUE)
```
```{r 12-ocena-18, fig.cap="Rozkład wartości błędów estymacji uzyskanych używając krigingu prostego i kroswalidacji."}
ggplot(cv_sk, aes(residual)) +
geom_histogram() +
xlab("Błąd estymacji") +
ylab("Liczebność")
```
```{r 12-ocena-19, fig.cap="Wykres rozrzutu porównujący rzeczywiste (obserwowane) wartości i wartości estymowane uzyskane używając krigingu prostego i kroswalidacji."}
ggplot(cv_sk, aes(var1.pred, observed)) +
geom_point() +
xlab("Estymacja") +
ylab("Obserwacja")
```
Podobnie jak w walidacji podzbiorem, gdy uzyskany model jest wystarczająco dobry, estymację dla całego obszaru uzyskuje się z użyciem funkcji `krige()` (rycina \@ref(fig:12-ocena-20)).
```{r 12-ocena-20, fig.cap="Estymacja uzyskana dla całego obszaru po oceneniu modelu używając kroswalidacji."}
cv_skk = krige(temp ~ 1,
locations = punkty,
newdata = siatka,
model = fitted,
beta = 15)
tm_shape(cv_skk) +
tm_raster("var1.pred", style = "cont", palette = "-Spectral")
```
<!--
```{r 12-ocena-21, eval=FALSE}
# ok_loocv = krige.cv(temp ~ 1, punkty, model=model_zl2)
# summary(ok_loocv)
```
- Tutaj inne przykłady
- Wykresy z loocv
- wykresy porównujące
```{r 12-ocena-22, eval=FALSE}
# ok_fit = gstat(formula=temp ~ 1, data=punkty, model=model_zl2)
# ok_loocv = gstat.cv(OK_fit, debug.level=0, random=FALSE)
# spplot(pe[6])
```
##
- prezentacja 5 Ani
- spatinter folder
- AIC
## Walidacja wyników estymacji
### Walidacja wyników estymacji | Kriging zwykły - LOO crossvalidation
krige.cv
```{r 12-ocena-23, eval=FALSE }
# OK_fit = gstat(id="OK_fit", formula=temp ~ 1, data=punkty, model=fitted)
# pe = gstat.cv(OK_fit, debug.level=0, random=FALSE)
# spplot(pe[6])
#
# z = predict(OK_fit, newdata = grid, debug.level = 0)
# grid2 = grid
# grid2$OK_pred = z$OK_fit.pred
# grid2$OK_se = sqrt(z$OK_fit.var)
# library('rasterVis')
# spplot(grid2, 'OK_pred')
# spplot(grid2, 'OK_se')
```
### Walidacja wyników estymacji | K Kriging uniwersalny - LOO crossvalidation
```{r 12-ocena-24, eval=FALSE }
# KU_fit = gstat(id="KU_fit", formula=temp~odl_od_morza, data=punkty, model=fitted_ku)
# pe = gstat.cv(KU_fit, debug.level=0, random=FALSE)
# spplot(pe[6])
# dodanie odległości od morza do siatki !!
# z_KU = predict(KU_fit, newdata = grid, debug.level = 0)
# grid$KU_pred = z$KU_fit.pred
# grid$KU_se = sqrt(z$KU_fit.var)
# library('rasterVis')
# spplot(grid, 'KU_pred')
# spplot(grid, 'KU_se')
```
-->
## Zadania {#z12}
1. Wydziel obiekt `punkty` w taki sposób aby 70% danych należało zbioru treningowego, a 30 % danych do zbioru testowego.
Zwizualizuj oba nowe zbiory danych.
2. Stwórz optymalne modele zmiennej `temp` ze zbioru treningowego używając krigingu zwykłego, kokrigingu oraz krigingu uniwersalnego.
3. Wykonaj estymacje korzystając z powyższych modeli.
Porównaj uzyskane estymacje korzystając ze statystyk jakości estymacji oraz wizualizacji jakości estymacji i używając zbioru testowego.
Który ze stworzonych modeli można uznać za najlepszy?
Dlaczego?
4. Porównaj uzyskane modele używając kroswalidacji.
Jak wygląda rozkład reszt z uzyskanych estymacji?
Który model można uznać za najlepszy oglądając rozkłady reszt?
5. Stwórz optymalne modele zmiennej `temp` ze zbioru treningowego używając krigingu zwykłego, kokrigingu oraz krigingu uniwersalnego, ale tym razem wcześniej transformując wartości tej zmiennej używając metod opisanych w sekcji \@ref(trans-danych).
Wykonaj estymacje korzystając z powyższych modeli.
Porównaj uzyskane estymacje korzystając ze statystyk jakości estymacji oraz wizualizacji jakości estymacji i używając zbioru testowego.
Który ze stworzonych modeli można uznać za najlepszy?
Dlaczego?
Czy któryś z uzyskanych modeli dla zmiennej po transformacji jest lepszy niż dla zmiennej przed transformacją?
<!-- jedno zadanie używające wcześniej używanych technik estymacji -->