-
Notifications
You must be signed in to change notification settings - Fork 2
/
tutorial.Rmd
831 lines (673 loc) · 33.8 KB
/
tutorial.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
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
---
title: "Geostatistics in R"
author: "Ellie White"
date: "December 9th, 2016"
output:
html_document:
df_print: paged
---
This tutorial was written for UC Davis course HYD 273: Introduction to Geostatistics: Statistical treatment of spatial data with hydrologic emphasis. Topics: theory of random functions, variogram analysis, Kriging/co-Kriging, indicator geostatistics, and stochastic simulation of spatial variability. Geostatistical software use.
Disclaimer: It will most definitely contain mistakes. Please let me know them at: [email protected]
Credit: Hijmans, R. (2016) Interpolation [Source code] http://www.rspatial.org for cleaned up example data and starter code.
# Contents
1.0 Spatial Data
2.0 Blank Interpolation Zone
3.0 Evaluation Metrics
4.0 The Average Model
5.0 Inverse Distance Weighted Model
6.0 Variogram
.... 6.1 Anisotropy
7.0 Kriging Models
.... 7.1 Simple Kriging
.... 7.2 Ordinary Kriging
.... 7.3 Universal Kriging
.... 7.4 Indicator Kriging
.... 7.5 Local Kriging
.... 7.6 CoKriging
.... 7.7 Block Kriging
8.0 Sequential Gaussian Simulation
9.0 Comparing Models
.... 9.1 Comparing RMSE
.... 9.2 Comparing Plots
10.0 Ensemble Model
11.0 Transition Probability Markov Chain
# Libraries
library(sp) --for geostatistical analysis
library(rgdal) --for spatial data transformations
library(raster) --for spatial data manipulation and analysis
library(gstat) --for geostatistical analysis
library(dismo) --for k fold cross-validation
library(spMC) --for transition probability markov chains
**note:** install each of these libraries by typing the following in your Console: install.packages("[name of library]")
```{r setup, include=FALSE}
library(knitr)
library(formatR)
opts_chunk$set(
fig.width =7.5,
fig.height=7.5,
collapse =TRUE,
tidy =TRUE,
tidy.opts=list(width.cutoff=120)
)
```
# 1.0 Spatial Data
**What:** Ozone Concentrations
**Units:** ppb (parts per billion)
**Time:** averages for 1980-2009
**Time resolution:** 30 yr average
**Data type:** .csv
**Data Source:** https://www.arb.ca.gov/aqd/aqdcd/aqdcddld.htm, or download the clean version from: http://rspatial.org/analysis/data/airqual.csv
**Provided by:** Cal EPA Air Resources Board
```{r data}
# download the clean version from: http://rspatial.org/analysis/data/airqual.csv, read in the data from your working directory
aq_data <- read.csv("inputdata/airqual.csv")
head(aq_data)
# change the OZDLYAV column to ppm (parts per million) so it is easier to read. This is the column of data we will be working with
aq_data$OZDLYAV <- aq_data$OZDLYAV*1000
aqd <- aq_data
# change data frame to a spatial points data frame
library(sp)
coordinates(aqd) <- ~LONGITUDE + LATITUDE
proj4string(aqd) <- CRS('+proj=longlat +datum=NAD83')
# transform data into a suitable projection, for California: Albers and Teale Albers are good
TA <- CRS("+proj=aea +lat_1=34 +lat_2=40.5 +lat_0=0 +lon_0=-120 +x_0=0 +y_0=-4000000 +datum=NAD83 +units=km +ellps=GRS80")
library(rgdal)
aq_ta <- spTransform(aqd, TA)
# for visualization purposes download California counties from link in: http://rspatial.org/analysis/rst/4-interpolation.html
library(raster)
ca <- shapefile("inputdata/counties/counties_2000.shp")
ca_ta <- spTransform(ca, TA)
# visualization of the data
plot(ca_ta, border='gray', axes=TRUE, main="Ozone Concentration (ppm)")
text(aq_ta, labels=round(aq_ta$OZDLYAV, 2), cex=.5, col='navy')
ucd_colors <- colorRampPalette(c("goldenrod", "navy"))
spplot(aq_ta, "OZDLYAV", col.regions=ucd_colors(6), key.space="right", main="Ozone Concentration (ppm)", sp.layout=list("sp.polygons", ca_ta, col="grey"))
hist(aq_data$OZDLYAV, main="Histogram", col="navy", border="goldenrod", xlab="OZDLYAV (ppm)")
hist(log(aq_data$OZDLYAV), main="Histogram", col="navy", border="goldenrod", xlab="log(OZDLYAV)")
```
# 2.0 Blank Interpolation Zone
This section creates a spatial grid over california (projected to Teale Albers) to interpolate or make predictions to, i.e. a blank slate to be later given values with each model developed.
```{r interpolation_zone}
r <- raster(ca_ta)
res(r) <- 10 # 10 km if your Coordinate Reference System's units are in km, the Teal Albers system is in km
# coerce r into a spatial grid, not very different from a raster
g <- as(r, 'SpatialGrid')
```
# 3.0 Evaluation Metrics
This section defines RMSE, R Squared, and NSE as functions to be used when evaluating and comparing models.
```{r rmse}
RMSE <- function(observed, predicted) {
sqrt(mean((predicted - observed)^2, na.rm=TRUE))
}
R2 <- function(observed, predicted) {
sum((predicted - mean(observed, na.rm=TRUE))^2, na.rm=TRUE)/sum((observed - mean(observed, na.rm=TRUE))^2, na.rm=TRUE)
}
# use this if you have hydrologic data, I have left it out of this analysis
NSE <- function(observed, predicted) {
1-(sum((predicted - observed)^2, na.rm=TRUE)/sum((mean(observed) - observed)^2, na.rm=TRUE))
}
```
# 4.0 The Average Model
This model is the arithmatic average of all the observations.
```{r average_model}
am <- mean(aqd$OZDLYAV)
# make a raster that has this value everywhere, we may need this for comparison to other models later.
am_raster <- r
values(am_raster) <- am
# evaluation of this model
am_rmse <- RMSE(am, aq_ta$OZDLYAV)
am_rmse
```
# 5.0 Inverse Distance Weighted Model
This model is a weighted combination of all observations. The weights will be purely based on distance.
```{r idw}
library(gstat)
# gstat objects hold all the information necessary for univariate or multivariate geostatistical prediction
gs <- gstat(formula=OZDLYAV~1, locations=aq_ta)
idw <- interpolate(r, gs)
idwr <- mask(idw, ca_ta)
# plot
plot(idwr, main="IDW")
plot(ca_ta, border='gray', add=TRUE)
contour(idwr, add=TRUE, nlevels=10)
```
## 5.1 Optimal Inverse Distance Weighted Model
In the IDW model there are two parameters that can be optimized: 1) the number of nearest observations that should be used for a prediction, and 2) the distance decay rate.
```{r opt_idw}
# optimization objective: minimize the RMSE
f1 <- function(x, test, train) {
nmx <- x[1]
idp <- x[2]
if (nmx < 1) return(Inf)
if (idp < .001) return(Inf)
m <- gstat(formula=OZDLYAV~1, locations=train, nmax=nmx, set=list(idp=idp))
p <- predict(m, newdata=test, debug.level=0)$var1.pred
RMSE(test$OZDLYAV, p)
}
# 5-fold cross validation
set.seed(20161209)
i <- sample(nrow(aq_ta), 0.2*nrow(aq_ta))
tst <- aq_ta[i,]
trn <- aq_ta[-i,]
# optimization model
opt <- optim(c(8, .5), f1, test=tst, train=trn)
opt
# optimal IDW model
# opt$par[1]: the number of nearest observations that should be used for a prediction
# opt$par[2]: the distance decay rate. the power on the distance in the demoninator
idw_opt <- gstat(formula=OZDLYAV~1, locations=aq_ta, nmax=opt$par[1], set=list(idp=opt$par[2]))
idw_b <- interpolate(r, idw_opt)
idw_best <- mask(idw_b, ca_ta)
# plot
plot(idw_best, main="IDW_BEST")
plot(ca_ta, border='gray', add=TRUE)
contour(idw_best, add=TRUE, nlevels=10)
```
# 6.0 Variogram
A variogram is a description of the spatial continuity of the data. The experimental variogram is a function fitted to the variogram. This function measures the variability between pairs of points at various distances and is used in Kriging models.
```{r variogram}
gs <- gstat(formula=OZDLYAV~1, locations=aq_ta)
v <- variogram(gs, width=20) # width is the width of the bins for grouping data
head(v)
plot(v, main="Variogram (Bins=20 km)")
# set cloud=TRUE to see all the pairs of observations
vc <- variogram(gs, width=20, cloud=TRUE)
head(vc)
plot(vc, main="Variogram Cloud")
# now to fit a model to the variogram. Here are all the types of fits we can try
show.vgms()
# fit a variogram model, best way is to eye ball the sill, range, and nugget
# here is an exponential fit
fve <- fit.variogram(v, vgm(85, "Exp", 75, 20))
fve
plot(variogramLine(fve, 400), type='l', ylim=c(0,120), main="Exponential Fit")
points(v[,2:3], pch=20, col='navy', cex=1.5)
# here is a spherical fit
fvs <- fit.variogram(v, vgm(85, "Sph", 75, 20))
fvs
plot(variogramLine(fvs, 400), type='l', ylim=c(0,120), main="Spherical Fit")
points(v[,2:3], pch=20, col='navy', cex=1.5)
```
## 6.1 Anisotropy
Anisotropy is when the range (geometric anisotropy) or the sill (zonal anisotropy) of the semivariogram varies in different directions. Both types of anisotropy can be present.
```{r anisotropy}
# use alpha to set the direction
vgm_aniso <- variogram(gs, alpha=c(0, 45, 90, 135))
plot(vgm_aniso)
# eye ball a model fit, eyeball the sill, range, and nugget
aniso_fit <- vgm(200, "Sph", 200, 0, anis=c(90, 1/3))
plot(vgm_aniso, aniso_fit, as.table=TRUE)
# that wasn't very good. let R try to fit a model through regression
raniso_fit <- fit.variogram(vgm_aniso, model=aniso_fit)
plot(vgm_aniso, raniso_fit)
# can also subset and eyeball a model fit in the 45 Degree direction
vgm_45 <- subset(vgm_aniso, vgm_aniso$dir.hor == 45)
plot(vgm_45)
vgm_45_sph <- vgm(150, "Sph", 200, 0)
plot(vgm_45, vgm_45_sph, main="Variogram (Dir=45 Degrees)") # notice this is another way of plotting variograms
```
# 7.0 Kriging Models
Originating in geostatistics, kriging or Gaussian process regression is a method of interpolation for which the interpolated values are modeled by a Gaussian process governed by prior covariances. We will look at a few different types of kriging below.
## 7.1 Simple Kriging --known mean
For simple kriging (and simulation based on simple kriging) you need to define a vector with the trend coefficients (including intercept); if no independent variables are defined the model only contains an intercept and beta should be the simple kriging mean.
```{r sk}
skp <- krige(formula=OZDLYAV~1, locations=aq_ta, g, model=fve, beta=mean(aq_data$OZDLYAV))
spplot(skp)
# plot nicely
sk <- brick(skp)
sk <- mask(sk, ca_ta)
names(sk) <- c('prediction', 'variance')
plot(sk)
title(main="Simple Kriging", line=6)
plot(sk$prediction, main="Simple Kriging")
plot(ca_ta, border='gray', add=TRUE)
contour(sk$prediction, add=TRUE, nlevels=10)
```
## 7.2 Ordinary Kriging --no known mean
Leave beta out of simple kriging to get an ordinary kriging estimate.
```{r ok}
okp <- krige(formula=OZDLYAV~1, locations=aq_ta, g, model=fve)
spplot(okp)
# plot nicely
ok <- brick(okp)
ok <- mask(ok, ca_ta)
names(ok) <- c('prediction', 'variance')
plot(ok)
title(main="Ordinary Kriging", line=6)
plot(ok$prediction, main="Ordinary Kriging")
plot(ca_ta, border='gray', add=TRUE)
contour(ok$prediction, add=TRUE, nlevels=10)
```
## 7.3 Universal Kriging --kriging with external drift
Universal kriging is useful when the variable to map is related to other spatially known variables. For universal kriging, suppose z is a trending surface either linearly or quadratically dependent on x and y. This section shows you how to do both.
```{r uk}
# set up the x and y coordinates
aq_ta$x <- coordinates(aq_ta)[, 1]
aq_ta$y <- coordinates(aq_ta)[, 2]
# cannot use g, a spatial grid, anymore. The grid here needs to be a spatial pixels data frame!
gpix <- as(r, 'SpatialPixelsDataFrame')
# linear trend surface
ukpl <- krige(formula=OZDLYAV~x+y, locations=aq_ta, gpix, model=fve)
spplot(ukpl)
# plot nicely
ukl <- brick(ukpl)
ukl <- mask(ukl, ca_ta)
names(ukl) <- c('prediction', 'variance')
plot(ukl)
title(main="Universal Kriging-Linear Trend Surface", line=6)
plot(ukl$prediction, main="Universal Kriging-Linear Trend Surface")
plot(ca_ta, border='gray', add=TRUE)
contour(ukl$prediction, add=TRUE, nlevels=10)
# quadratic trend surface
ukpq <- krige(formula=OZDLYAV~ x + y + I(x*y) + I(x^2) + I(y^2), locations=aq_ta, gpix, model=fve)
spplot(ukpq)
# plot nicely
ukq <- brick(ukpq)
ukq <- mask(ukq, ca_ta)
names(ukq) <- c('prediction', 'variance')
plot(ukq)
title(main="Universal Kriging-Quadratic Trend Surface", line=6)
plot(ukq$prediction, main="Universal Kriging-Quadratic Trend Surface")
plot(ca_ta, border='gray', add=TRUE)
contour(ukq$prediction, add=TRUE, nlevels=10)
```
## 7.4 Indicator Kriging
Use indicator kriging to map the probability of ozone exceeding a critical threshold.
```{r ik}
# pick a threshold, lets look at our variable again
summary(aq_data$OZDLYAV)
# develop a variogram for indicator kriging
gsi <- gstat(formula=I(OZDLYAV>30)~1, locations=aq_ta)
vi <- variogram(gsi, width=20)
fvis <- fit.variogram(vi, vgm(85, "Exp", 75, 20))
plot(variogramLine(fvis, 400), type='l', ylim=c(0,0.3), main="Exponential Fit")
points(vi[,2:3], pch=20, col='navy', cex=1.5)
# indicator kriging
ikp <- krige(formula=I(OZDLYAV>30)~1, locations=aq_ta, g, model=fvis, nmax=30, indicators=TRUE)
spplot(ikp)
# plot nicely
ik <- brick(ikp)
ik <- mask(ik, ca_ta)
names(ik) <- c('prediction', 'variance')
plot(ik)
title(main="Indicator Kriging (i>30ppm)", line=6)
```
## 7.5 Local Kriging
In kriging, by default, all observations are used. But we can use nmin, nmax and or maxdist to restrict the number of observations the kriging estimator can use.
**nmax:** the max number of nearest observations that should be used for a kriging prediction or simulation, where nearest is defined in terms of the space of the spatial locations.
**maxdist:** only observations within a distance of maxdist from the prediction location are used for prediction or simulation; if combined with nmax, both criteria apply.
**nmin:** if the number of nearest observations within distance maxdist is less than nmin, a missing value will be generated; see maxdist.
```{r lk}
# remember the dimensions of the data
dim(aq_data)
# and the clustering of data points
plot(ca_ta, border='gray', main="Data Points")
points(aq_ta, col='navy')
# now pick an appropriate nmin, nmax, and or maxdist
lkp <- krige(formula=OZDLYAV~1, locations=aq_ta, g, model=fve, nmin=20, nmax=100, maxdist=100)
spplot(lkp)
# plot nicely
lk <- brick(lkp)
lk <- mask(lk, ca_ta)
names(lk) <- c('prediction', 'variance')
plot(lk)
title(main="Local Kriging", line=6)
plot(lk$prediction, main="Local Kriging")
plot(ca_ta, border='gray', add=TRUE)
contour(lk$prediction, add=TRUE, nlevels=10)
```
## 7.6 CoKriging
Use cokriging, or making use of other data sets, to make a better prediction. If the new data is grounded in pysicall processes that influence the variable of interest, the variance in the error prediction can be reduced.
```{r cok}
# find potentially useful data sets, use scatterplot matrix
# this function will just print the r and pvalue in our scatterplot matrix below
panel.cor <- function(x, y, digits=2, cex.cor, ...){
usr <- par("usr"); on.exit(par(usr))
par(usr=c(0, 1, 0, 1))
# correlation coefficient
r <- cor(x, y)
txt <- format(c(r, 0.123456789), digits=digits)[1]
txt <- paste("r= ", txt, sep="")
text(0.5, 0.6, txt)
# p-value calculation
p <- cor.test(x, y)$p.value
txt2 <- format(c(p, 0.123456789), digits=digits)[1]
txt2 <- paste("p= ", txt2, sep="")
if(p<0.01) txt2 <- paste("p= ", "<0.01", sep="")
text(0.5, 0.4, txt2)
}
pairs(aq_data[,22:28], upper.panel=panel.cor)
# based on the r values in the scatter plot matrix above, I will be using the following data sets
summary(aq_ta$OZMAX8O)
summary(aq_ta$OZMX8ST)
# now to create a gstat object and the variogram model
gsc1 <- gstat(id="OZDLYAV", formula=OZDLYAV~1, locations=aq_ta)
gsc2 <- gstat(gsc1, id="OZMAX8O", formula=OZMAX8O~1, locations=aq_ta)
gsc3 <- gstat(gsc2, id="OZMX8ST", formula=OZMX8ST~1, locations=aq_ta)
# plot the variograms and cross-variograms
plot(variogram(gsc3))
# fit a model variogram to the target variogram
fvsc <- fit.variogram(variogram(gsc1), vgm(85, "Sph", 75, 20))
# fit a model variogram to all the variograms
vmc <- variogram(gsc3)
fvmc <- fit.lmc(vmc, gsc3, model=fvsc)
# plot the fitted variograms to all the sample variograms:
plot(variogram(gsc3), fvmc)
# co-kriging, WARNING: this may take a few minutes!
ckp <- predict(fvmc, g)
spplot(ckp)
# plot nicely
ck <- brick(ckp)
ck <- mask(ck, ca_ta)
plot(ck)
plot(ck$OZDLYAV.pred, main="Co-Kriging")
plot(ca_ta, border='gray', add=TRUE)
contour(ck$OZDLYAV.pred, add=TRUE, nlevels=10)
```
## 7.7 Block Kriging --for rectangular or irregular blocks
Use block kriging to smooth out the surface, and have a less computationally expensive model.
**block size:** a vector with 1, 2 or 3 values containing the size of a rectangular in x-, y- and z-dimension respectively (0 if not set), or a data frame with 1, 2 or 3 columns, containing the points that discretize the block in the x-, y- and z-dimension to define irregular blocks relative to (0,0) or (0,0,0)—see also the details section of predict.gstat.
```{r bk}
# let's try smaller blocks first
bkp_small <- krige(formula=OZDLYAV~1, locations=aq_ta, g, model=fve, block=c(50, 50))
spplot(bkp_small)
# plot nicely
bk_small <- brick(bkp_small)
bk_small <- mask(bk_small, ca_ta)
names(bk_small) <- c('prediction', 'variance')
plot(bk_small)
title(main="Block Kriging (b=50km)", line=6)
plot(bk_small$prediction, main="Block Kriging (b=50km)")
plot(ca_ta, border='gray', add=TRUE)
contour(bk_small$prediction, add=TRUE, nlevels=10)
# let's try larger blocks
bkp_large <- krige(formula=OZDLYAV~1, locations=aq_ta, g, model=fve, block=c(400, 400))
spplot(bkp_large)
# plot nicely
bk_large <- brick(bkp_large)
bk_large <- mask(bk_large, ca_ta)
names(bk_large) <- c('prediction', 'variance')
plot(bk_large)
title(main="Block Kriging (b=400km)", line=6)
plot(bk_large$prediction, main="Block Kriging (b=400km)")
plot(ca_ta, border='gray', add=TRUE)
contour(bk_large$prediction, add=TRUE, nlevels=10)
```
# 8.0 Sequential Gaussian Simulation
Sequential Gaussian simulation (SGS) simulates continuous variables, such as our ozone concentration data. Sequential indicator simulation (SIS) simulates discrete variables, such as an ozone indicator variable that splits the ozone concentration to high and low values based on a threshold value.
**nsim:** integer; if set to a non-zero value, conditional simulation is used instead of kriging interpolation. For this, sequential Gaussian or indicator simulation is used (depending on the value of indicators), following a single random path through the data.
```{r conditional_sim}
# SGS
condsimp <- krige(formula=OZDLYAV~1, locations=aq_ta, g, model=fve, nmax=30, nsim=4)
spplot(condsimp, main="four conditional simulations")
# plot nicely
condsim <- brick(condsimp)
condsim <- mask(condsim, ca_ta)
plot(condsim)
# SIS, note that the model should be set to the fit to the indicator variogram!
condsimp2 <- krige(formula=I(OZDLYAV>30)~1, locations=aq_ta, g, model=fvis, nmax=30, nsim=4, indicators=TRUE)
spplot(condsimp2, main="four conditional simulations")
# plot nicely
condsim2 <- brick(condsimp2)
condsim2 <- mask(condsim2, ca_ta)
plot(condsim2)
```
# 9.0 Comparing Models
To recap, the models developed are: Average Model(am), Inverse Distance Weighted (idw), Simple Kriging (sk), Ordinary Kriging (ok), Universal Kriging with a Linear Trend (ukl), Universal Kriging with a Quadratic Trend (ukq), Indicator Kriging (ik), Co-Kriging (ck), Block Kriging (bk).
note: I left co-kriging out of the cross validation, because the code would take a looooong time to run.
## 9.1 Comparing RMSE and NSE
Use k-fold crossvalidation to evaluate the models.
```{r compare_rmse}
library(dismo)
nfolds <- 5
k <- kfold(aq_ta, nfolds)
# define vectors to be populated later
ensrmse <- idwrmse <- skrmse <- okrmse <- uklrmse <- ukqrmse <- ikrmse <- lkrmse <- bk_lrmse <- rep(NA, 5)
ensrsq <- idwrsq <- skrsq <- okrsq <- uklrsq <- ukqrsq <- ikrsq <- lkrsq <- bk_lrsq <- rep(NA, 5)
# cross validation loop
for (i in 1:nfolds) {
test <- aq_ta[k!=i,]
train <- aq_ta[k==i,]
# idw_best
m <- gstat(formula=OZDLYAV~1, locations=train, nmax=opt$par[1], set=list(idp=opt$par[2]))
p1 <- predict(m, newdata=test, debug.level=0)$var1.pred
idwrmse[i] <- RMSE(test$OZDLYAV, p1)
idwrsq[i] <- R2(test$OZDLYAV, p1)
# sk
m <- gstat(formula=OZDLYAV~1, locations=train, model=fve, beta=mean(aq_ta$OZDLYAV))
p2 <- predict(m, newdata=test, debug.level=0)$var1.pred
skrmse[i] <- RMSE(test$OZDLYAV, p2)
skrsq[i] <- R2(test$OZDLYAV, p2)
# ok
m <- gstat(formula=OZDLYAV~1, locations=train, model=fve)
p3 <- predict(m, newdata=test, debug.level=0)$var1.pred
okrmse[i] <- RMSE(test$OZDLYAV, p3)
okrsq[i] <- R2(test$OZDLYAV, p3)
# ukl
m <- gstat(formula=OZDLYAV~x+y, locations=aq_ta, model=fve)
p4 <- predict(m, newdata=test, debug.level=0)$var1.pred
uklrmse[i] <- RMSE(test$OZDLYAV, p4)
uklrsq[i] <- R2(test$OZDLYAV, p4)
# ukq
m <- gstat(formula=OZDLYAV~x+y+I(x*y)+I(x^2)+I(y^2), locations=aq_ta, model=fve)
p5 <- predict(m, newdata=test, debug.level=0)$var1.pred
ukqrmse[i] <- RMSE(test$OZDLYAV, p5)
ukqrsq[i] <- R2(test$OZDLYAV, p5)
# ik
m <- gstat(formula=I(OZDLYAV>30)~1, locations=aq_ta, model=fvis, nmax=30)
p6 <- predict(m, newdata=test, debug.level=0, indicators=TRUE)$var1.pred
test$indicator <- ifelse(test$OZDLYAV>30, 1, 0)
ikrmse[i] <- RMSE(test$indicator, p6)
ikrsq[i] <- R2(test$indicator, p6)
# lk
m <- gstat(formula=OZDLYAV~1, locations=aq_ta, model=fve, nmin=20, nmax=100, maxdist=100)
p7 <- predict(m, newdata=test, debug.level=0)$var1.pred
lkrmse[i] <- RMSE(test$OZDLYAV, p7)
lkrsq[i] <- R2(test$OZDLYAV, p7)
# bk
m <- gstat(formula=OZDLYAV~1, locations=aq_ta, model=fve)
p8 <- predict(m, newdata=test, debug.level=0, block=c(400, 400))$var1.pred
bk_lrmse[i] <- RMSE(test$OZDLYAV, p8)
bk_lrsq[i] <- R2(test$OZDLYAV, p8)
w <- c(idwrsq[i], skrsq[i], okrsq[i], uklrsq[i], ukqrsq[i], ikrsq[i], lkrsq[i], bk_lrsq[i])
weights <- w / sum(w)
# ensemble model
ensemble <- p1*weights[1]+ p2*weights[2]+ p3*weights[3]+ p4*weights[4]+ p5*weights[5]+ p6*weights[6]+ p7*weights[7]+ p8*weights[8]
ensrmse[i] <- RMSE(test$OZDLYAV, ensemble)
ensrsq[i] <- R2(test$OZDLYAV, ensemble)
}
idw_rmse <- mean(idwrmse)
sk_rmse <- mean(skrmse)
ok_rmse <- mean(okrmse)
ukl_rmse <- mean(uklrmse)
ukq_rmse <- mean(ukqrmse)
ik_rmse <- mean(ikrmse)
lk_rmse <- mean(lkrmse)
bk_lrmse <- mean(bk_lrmse)
idw_rsq <- mean(idwrsq)
sk_rsq <- mean(skrsq)
ok_rsq <- mean(okrsq)
ukl_rsq <- mean(uklrsq)
ukq_rsq <- mean(ukqrsq)
ik_rsq <- mean(ikrsq)
lk_rsq <- mean(lkrsq)
bk_lrsq <- mean(bk_lrsq)
rms <- c(idw_rmse, sk_rmse, ok_rmse, ukl_rmse, ukq_rmse, ik_rmse, lk_rmse, bk_lrmse)
rsq <- c(idw_rsq, sk_rsq, ok_rsq, ukl_rsq, ukq_rsq, ik_rsq, lk_rsq, bk_lrsq)
names(rms) <- names(rsq) <- c('IDW_Best', 'SK', 'OK', 'UK_Lin', 'UK_Quad', 'IK', 'LK', 'BK_Large')
rms
rsq
# barplot the RMSE, store and print the rmse on the bars
ylim <- c(0, 1.1*max(rms))
xx <- barplot(rms, main="RMSE Comparison", xlab="Geostatistical Methods", ylim=ylim, col="goldenrod")
text(x=xx, y=rms, label=round(rms, digits=2), pos=3, cex=0.8, col="navy")
# barplot the RSQ, store and print the rmse on the bars
ylim <- c(0, 1.1*max(rsq))
xx <- barplot(rsq, main="R^2 Comparison", xlab="Geostatistical Methods", ylim=ylim, col="goldenrod")
text(x=xx, y=rsq, label=round(rsq, digits=2), pos=3, cex=0.8, col="navy")
```
## 9.2 Comparing Plots
Use plots to identify spatial differences and maybe go back to tweek model parameters.
```{r comparing_plots}
# compare prediction plots
s <- stack(idw_best, sk[[1]], ok[[1]], ukl[[1]], ukq[[1]], ik[[1]], lk[[1]], ck[[1]], bk_large[[1]])
names(s) <- c('IDW_Best', 'SK', 'OK', 'UK_Lin', 'UK_Quad', 'IK', 'LK', 'Co_K', 'BK_Large')
plot(s)
# compare errors plots
s_stdev <- stack(sk[[2]], ok[[2]], ukl[[2]], ukq[[2]], ik[[2]], lk[[2]], ck[[2]], bk_large[[2]])
names(s_stdev) <- c('SK_error', 'OK_error', 'UK_Lin_error', 'UK_Quad_error', 'IK_error', 'LK_error', 'Co_K_error', 'BK_Large_error')
plot(s_stdev)
```
# 10.0 Ensemble Model
We can combine any or all models by taking a weighted average of all the model predictions with R^2 as the weights.
```{r ensemble}
# from cross validation we have:
rmse_e <- mean(ensrmse)
rmse_e
rsq_e <- mean(ensrsq)
rsq_e
# weights were calculated in the cv loop
names(weights) <- c('IDW_Best', 'SK', 'OK', 'UK_Lin', 'UK_Quad', 'IK', 'LK', 'BK_Large')
weights
models <- stack(idw_best, sk[[1]], ok[[1]], ukl[[1]], ukq[[1]], ik[[1]], lk[[1]], bk_large[[1]])
ensemble_model <- sum(models*weights)
plot(ensemble_model, main="Ensemble Model")
plot(ca_ta, border='gray', add=TRUE)
contour(ensemble_model, add=TRUE, nlevels=10)
```
# 11.0 Transition Probability Markov Chain
Much like the TPROGS software, the spMC lirary is used for modelling of transition probability matrices, and implementing efficient algorithms for improving predictions and simulations of categorical random fields. It deals with both one-dimensional and multidimensional continuous lag models.
Three different fitting methods were implemented in the package. The first is based on the estimates of the main features that characterise the process, the second focuses on the minimisation of the discrepancies between the empirical and theoretical transition probabilities, and the third follows the maximum entropy approach. Once the model parameters are properly estimated, transition probabilities are calculated through the matrixvalued exponential function.
For the estimations of one-dimensional continuous lag models use:
**transiogram:** Empirical transition probabilities estimation
**tpfit:** One-dimensional model parameters estimation
**tpfit_ils:** Iterated least squares method for one-dimensional model parameters estimation
**tpfit_me:** Maximum entropy method for one-dimensional model parameters estimation
**tpfit_ml:** Mean length method for one-dimensional model parameters estimation
For the estimations of multidimensional continuous lag models use:
**pemt:** Pseudo-empirical multidimensional transiograms estimation
**multi_tpfit:** Multidimensional model parameters estimation
**multi_tpfit_ils:** Iterated least squares method for multidimensional model parameters estimation
**multi_tpfit_me:** Maximum entropy method for multidimensional model parameters estimation
**multi_tpfit_ml:** Mean length method for multidimensional model parameters estimation
These transition probabilities are then combined to predict the category in an unsampled position. Three algorithms are used to simulate spatial random fields; those based on the kriging techniques (Carle and Fogg, 1996), those using fixed and random path methods (Li, 2007a; Li and Zhang, 2007), or those using multinomial categorical simulation proposed by Allard et al. (2011).
For categorical spatial random field simulation and prediction use:
**sim:** Random field simulation
**sim_ck:** Conditional simulation based on indicator cokriging
**sim_ik:** Conditional simulation based on indicator kriging
**sim_mcs:** Multinomial categorical simulation
**sim_path:** Conditional simulation based on path algorithms
For more information see the spMC package (Sartore, 2016).
```{r transition_probability}
# change the ozone concentration data into categorical data
# make the groups, you can see the interval categories now in the data set
aq_ta$group <- cut(aq_ta$OZDLYAV, breaks=6)
table(aq_ta$group)
# we can lump high and very high, because we don't have enough data, but I will leave it for now
# make the groupings into categories: verylow, low, etc.
aq_ta$categories <- cut(aq_ta$OZDLYAV, breaks=6, labels=c("verylow", "low", "medlow", "medhigh", "high", "veryhigh"))
table(aq_ta$categories)
# plot the categorical variable now
library(RColorBrewer)
brwCol <- brewer.pal(nlevels(aq_ta$categories), "YlGnBu")
spplot(aq_ta, "categories", col.regions=brwCol, key.space="right", main="Ozone Concentration Categories", sp.layout=list("sp.polygons", ca_ta, col="grey"), pch=15)
# now we can look at the spatial transition probabilities in these categories
library(spMC)
# make two transiograms to estimate two transition rate matrices, i.e. the model coefficients, along the axial directions
# first, all couples of points along axial directions are chosen such that their lag-length is less than three km. After, we calculate the empirical transition probabilities for twenty points within the maximum distance
Trg <- list()
Trg[[1]] <- transiogram(aq_ta$categories, aq_ta@coords, max.dist=30, mpoints=50, direction=c(1, 0))
Trg[[2]] <- transiogram(aq_ta$categories, aq_ta@coords, max.dist=30, mpoints=50, direction=c(0, 1))
# if we want to compare these probabilities with the theoretical one, we first need to estimate two transition rate matrices, i.e. the model coefficients, along the axial directions
# estimate the transition rates through the mean lengths; estimate the parameters of the MC model
RTm <- list()
RTm[[1]] <- tpfit_ml(aq_ta$categories, aq_ta@coords, direction=c(1, 0))
RTm[[2]] <- tpfit_ml(aq_ta$categories, aq_ta@coords, direction=c(0, 1))
RTm
# estimate the transition rate matrices through the maximum entropy approach
ETm <- list()
ETm[[1]] <- tpfit_me(aq_ta$categories, aq_ta@coords, direction=c(1, 0))
ETm[[2]] <- tpfit_me(aq_ta$categories, aq_ta@coords, direction=c(0, 1))
ETm
# given the model coefficients, the transition probabilities for some specific lags are calculated
RTr <- list()
ETr <- list()
for (i in 1:2) {
RTr[[i]] <- predict(RTm[[i]], lags=Trg[[i]]$lags)
ETr[[i]] <- predict(ETm[[i]], lags=Trg[[i]]$lags)
}
# Since these probabilities are calculated with respect to some fixed directions, i.e. by considering a one-dimensional perspective, they can be graphically compared
# a transition probability matrix plot, it shows the probability dynamic related to one-dimensional lags along the specified direction
for (i in 1:2){
mixplot(list(Trg[[i]], RTr[[i]], ETr[[i]]), type=c("p", "l", "l"), pch="+", col=c(3, 1, 2), legend=FALSE, main=paste("One-dimensional transiograms", c("(X-axis)", "(Y-axis)")[i]))
}
# The function pemt() can be considered as another tool to check the anisotropy of the process. It estimates the transition rate matrix for each multidimensional lag direction and computes the transition probabilities. At the same time the function calculates other probabilities through the transition rates. Then the probabilities are drawn.
psEmpTr <- pemt(aq_ta$categories, aq_ta@coords, mpoints=40, max.dist=c(50, 50))
# dev.off() # uncomment and run this line if plot is not running
image(psEmpTr, col=rev(heat.colors(500)), useRaster=TRUE, breaks=c(0:500)/500, contour=FALSE, mar=c(0.1, 0.1, 0.1, 0.1))
# plot again with contour lines that are displayed for both the pseudo-empirical and the theoretical probabilities
# dev.off()
contour.pemt(psEmpTr, mar=c(0.1, 0.1, 0.1, 0.1), col=c("goldenrod", "navy"))
# methods to estimate transition rate matrices along axial directions. These functions implement algorithms based on the mean lengths and maximum entropy
MTr <- list()
MTr$average <- multi_tpfit_ml(aq_ta$categories, aq_ta@coords)
MTr$entropy <- multi_tpfit_me(aq_ta$categories, aq_ta@coords)
# with the output of these functions, we can draw the theoretical transition probability maps
# dev.off()
image(MTr$average, 40, max.dist=50, col=rev(heat.colors(500)), nlevels=5, breaks=0:500/500)
# dev.off()
image(MTr$entropy, 40, max.dist=50, col=rev(heat.colors(500)), nlevels=5, breaks=0:500/500)
# now, the simulation
# plot the map of observations again
X <- aq_ta$x
Y <- aq_ta$y
par(mfrow=c(1, 1), mar=c(5, 4, 4, 2))
plot(X, Y, col=brwCol[aq_ta$categories], pch=15, main="Ozone Concentration Categories")
legend("topright", legend=levels(aq_ta$categories), col=brwCol, pch=15)
# sample observations
set.seed(20161209)
smp <- sample(length(aq_ta):1, 100)
# plot the sample
plot(X, Y, type="n", main="Sample of 100 observations", xlab="X", ylab="Y")
points(X[smp], Y[smp], pch=19, col=brwCol[aq_ta$categories[smp]])
legend("topright", legend=levels(aq_ta$categories), col=brwCol, pch=19)
# make a grid
grid <- expand.grid(X,Y)
# the kriging algorithm will approximate the conditional probabilities by considering the twelve nearest neighbors for all points in the simulation grid.
iks <- sim_ik(MTr$average, data=aq_ta$categories[smp], coords=aq_ta@coords[smp,], grid, knn=12)
# Both fixed and random path simulation methods are performed by considering those nearest points along the axial directions within a radius of length ten.
fpth <- sim_path(MTr$average, data=aq_ta$categories[smp], coords=aq_ta@coords[smp,], grid, radius=1, TRUE)
rpth <- sim_path(MTr$average, data=aq_ta$categories[smp], coords=aq_ta@coords[smp,], grid, radius=1)
# The multinomial categorical simulation method will approximate the prediction probabilities by considering all sample points.
mcs <- sim_mcs(MTr$average, data=aq_ta$categories[smp], coords=aq_ta@coords[smp, ], grid)
posCol <- as.integer(iks$Prediction)
plot(X, Y, pch=15, col=brwCol[posCol], main="Kriging prediction map")
legend("topright", legend=levels(aq_ta$categories), col=brwCol, pch=15)
posCol <- as.integer(fpth$Prediction)
plot(X, Y, pch=15, col=brwCol[posCol], main="Fixed path prediction map")
legend("topright", legend=levels(aq_ta$categories), col=brwCol, pch=15)
posCol <- as.integer(rpth$Prediction)
plot(X, Y, pch=15, col=brwCol[posCol], main="Random path prediction map")
legend("topright", legend=levels(aq_ta$categories), col=brwCol, pch=15)
posCol <- as.integer(mcs$Prediction)
plot(X, Y, pch=15, col=brwCol[posCol],
main="Multinomial categorical prediction map")
legend("topright", legend=levels(aq_ta$categories), col=brwCol, pch=15)
# to compute the number of matches, we calculate the contingency table
ikTb <- table(aq_ta$categories[smp], iks$Prediction[smp])
fpTb <- table(aq_ta$categories[smp], fpth$Prediction[smp])
rpTb <- table(aq_ta$categories[smp], rpth$Prediction[smp])
mcTb <- table(aq_ta$categories[smp], mcs$Prediction[smp])
# relative frequencies of matches
ikPr <- sum(diag(ikTb))/length(aq_ta$categories)
fpPr <- sum(diag(fpTb))/length(aq_ta$categories)
rpPr <- sum(diag(rpTb))/length(aq_ta$categories)
mcPr <- sum(diag(mcTb))/length(aq_ta$categories)
# print in a table
probcov <- c(ikPr, fpPr, rpPr, mcPr)
names(probcov) <- "Probability of coverage"
methods <- c("Kriging", "Fixed path", "Random path", "Multinomial")
comparison_methods <- data.frame(methods, probcov)
comparison_methods
# in order to establish which is the best predictor, one should perform these simulations more than once. At each time, another 100 observations must be randomly selected.
```