Kriging
DIEGO CAMILO CELADA LOZADA
5/5/2021
Predicciones espaciales “Kriging”
1. Carga de libreria y datos
library(sp)
library(gstat)## Warning: package 'gstat' was built under R version 4.0.5library(nortest)Se carga la libreria SP para trabajar con la base de datos meuse y la libreria gstat para realizar el analisis
2. Analisis exploratorio
Se plantea el uso de la base de datos meuse, hacemos una exploración de las principales caracteristicas de esta.
data(meuse)
head(meuse)class(meuse)## [1] "data.frame"str(meuse)## 'data.frame': 155 obs. of 14 variables:
## $ x : num 181072 181025 181165 181298 181307 ...
## $ y : num 333611 333558 333537 333484 333330 ...
## $ cadmium: num 11.7 8.6 6.5 2.6 2.8 3 3.2 2.8 2.4 1.6 ...
## $ copper : num 85 81 68 81 48 61 31 29 37 24 ...
## $ lead : num 299 277 199 116 117 137 132 150 133 80 ...
## $ zinc : num 1022 1141 640 257 269 ...
## $ elev : num 7.91 6.98 7.8 7.66 7.48 ...
## $ dist : num 0.00136 0.01222 0.10303 0.19009 0.27709 ...
## $ om : num 13.6 14 13 8 8.7 7.8 9.2 9.5 10.6 6.3 ...
## $ ffreq : Factor w/ 3 levels "1","2","3": 1 1 1 1 1 1 1 1 1 1 ...
## $ soil : Factor w/ 3 levels "1","2","3": 1 1 1 2 2 2 2 1 1 2 ...
## $ lime : Factor w/ 2 levels "0","1": 2 2 2 1 1 1 1 1 1 1 ...
## $ landuse: Factor w/ 15 levels "Aa","Ab","Ag",..: 4 4 4 11 4 11 4 2 2 15 ...
## $ dist.m : num 50 30 150 270 380 470 240 120 240 420 ...Se realiza el analisis de la variable Zinc de la base de datos meuse, por lo cual analizamos los datos que contiene
head(meuse$zinc) # Obtenemos los seis primeros datos de la variable zinc## [1] 1022 1141 640 257 269 281hist(meuse$zinc) # Realizamos un histograma de la variable
boxplot(meuse$zinc) # Diagramamos una caja de bigotes de la variable
summary(meuse$zinc)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 113.0 198.0 326.0 469.7 674.5 1839.0Se determina que la variable tiene un sesgo positivo en la distribución de los datos, presenta asimetria hacia el lado derecho.
Realizamos la prueba de Kolmogorov - Smirnov
ks.test(meuse$zinc,mean(meuse$zinc),sd(meuse$zinc))## Warning in ks.test(meuse$zinc, mean(meuse$zinc), sd(meuse$zinc)): cannot compute
## exact p-value with ties##
## Two-sample Kolmogorov-Smirnov test
##
## data: meuse$zinc and mean(meuse$zinc)
## D = 0.6129, p-value = 0.8495
## alternative hypothesis: two-sidedlillie.test(meuse$zinc)# Comprobamos que se rechaza la hipotesis nula y los datos no proceden de una distribución normal##
## Lilliefors (Kolmogorov-Smirnov) normality test
##
## data: meuse$zinc
## D = 0.16643, p-value = 2.492e-11Realizamos una transformación de la variable zinc por medio de la funcion logaritmo, para buscar una distribución simetrica de los datos
meuse$logZn = log10(meuse$zinc)Realizamos un analisis de los datos con la nueva variable
head(meuse$logZn)## [1] 3.009451 3.057286 2.806180 2.409933 2.429752 2.448706hist(meuse$logZn)
boxplot(meuse$logZn)
3. Analisis de la estructura espacial
coordinates(meuse) = c('x', 'y')
plot(meuse, pch=1)
data(meuse.riv)
lines(meuse.riv)
plot(meuse, asp = 1, cex = 4*meuse$zinc/max(meuse$zinc), pch = 1)
lines(meuse.riv)
ve = variogram(logZn ~ 1, meuse, cutoff = 1300, width = 90)
veplot(ve, plot.numbers = T, asp = 1)
show.vgms()
Semivariograma Esferico
vt = vgm(psill = 0.12, model = 'Sph', range = 850, nugget = 0.01)
vtplot(ve, pl = T, model = vt)
va = fit.variogram(ve, vt)
vaplot(ve, pl = T, model = va)
Semivariograma Exponencial
vt1 = vgm(psill = 0.12, model = 'Exp', range = 850, nugget = 0.01)
vt1plot(ve, pl = T, model = vt1)
va1 = fit.variogram(ve, vt1)
va1plot(ve, pl = T, model = va1)
4. Predicciónes espaciales (Kriging Ordinario)
data(meuse.grid)
head(meuse.grid)coordinates(meuse.grid) = c('x', 'y')
class(meuse.grid)## [1] "SpatialPointsDataFrame"
## attr(,"package")
## [1] "sp"gridded(meuse.grid) = T
class(meuse.grid)## [1] "SpatialPixelsDataFrame"
## attr(,"package")
## [1] "sp"plot(meuse.grid)
Prediccion con el semivariograma esferico
OK = krige(logZn ~ 1, locations = meuse, newdata = meuse.grid, model = va)## [using ordinary kriging]names(OK)## [1] "var1.pred" "var1.var"head(OK)## Object of class SpatialPixelsDataFrame
## Object of class SpatialPixels
## Grid topology:
## cellcentre.offset cellsize cells.dim
## x 178460 40 78
## y 329620 40 104
## SpatialPoints:
## x y
## 1 181180 333740
## 2 181140 333700
## 3 181180 333700
## 4 181220 333700
## 5 181100 333660
## 6 181140 333660
## Coordinate Reference System (CRS) arguments: NA
##
## Data summary:
## var1.pred var1.var
## Min. :2.782 Min. :0.03355
## 1st Qu.:2.832 1st Qu.:0.04019
## Median :2.857 Median :0.04897
## Mean :2.859 Mean :0.04733
## 3rd Qu.:2.884 3rd Qu.:0.05396
## Max. :2.940 Max. :0.05958OK$pred = 10^(OK$var1.pred) # ELIMINAMOS EL LOGARITMO
pts.s = list('sp.points', meuse, col = 'white', pch=1, cex = 4*meuse$zinc/max(meuse$zinc))
pts.s## [[1]]
## [1] "sp.points"
##
## [[2]]
## coordinates cadmium copper lead zinc elev dist om ffreq soil
## 1 (181072, 333611) 11.7 85 299 1022 7.909 0.00135803 13.6 1 1
## 2 (181025, 333558) 8.6 81 277 1141 6.983 0.01222430 14.0 1 1
## 3 (181165, 333537) 6.5 68 199 640 7.800 0.10302900 13.0 1 1
## 4 (181298, 333484) 2.6 81 116 257 7.655 0.19009400 8.0 1 2
## 5 (181307, 333330) 2.8 48 117 269 7.480 0.27709000 8.7 1 2
## 6 (181390, 333260) 3.0 61 137 281 7.791 0.36406700 7.8 1 2
## 7 (181165, 333370) 3.2 31 132 346 8.217 0.19009400 9.2 1 2
## 8 (181027, 333363) 2.8 29 150 406 8.490 0.09215160 9.5 1 1
## 9 (181060, 333231) 2.4 37 133 347 8.668 0.18461400 10.6 1 1
## 10 (181232, 333168) 1.6 24 80 183 9.049 0.30970200 6.3 1 2
## 11 (181191, 333115) 1.4 25 86 189 9.015 0.31511600 6.4 1 2
## 12 (181032, 333031) 1.8 25 97 251 9.073 0.22812300 9.0 1 1
## 13 (180874, 333339) 11.2 93 285 1096 7.320 0.00000000 15.4 1 1
## 14 (180969, 333252) 2.5 31 183 504 8.815 0.11393200 8.4 1 1
## 15 (181011, 333161) 2.0 27 130 326 8.937 0.16833600 9.1 1 1
## 16 (180830, 333246) 9.5 86 240 1032 7.702 0.00000000 16.2 1 1
## 17 (180763, 333104) 7.0 74 133 606 7.160 0.01222430 16.0 1 1
## 18 (180694, 332972) 7.1 69 148 711 7.100 0.01222430 16.0 1 1
## 19 (180625, 332847) 8.7 69 207 735 7.020 0.00000000 13.7 1 1
## 20 (180555, 332707) 12.9 95 284 1052 6.860 0.00000000 14.8 1 1
## 21 (180642, 332708) 5.5 53 194 673 8.908 0.07034680 10.2 1 1
## 22 (180704, 332717) 2.8 35 123 402 8.990 0.09751360 7.2 1 1
## 23 (180704, 332664) 2.9 35 110 343 8.830 0.11393200 7.2 1 1
## 24 (181153, 332925) 1.7 24 85 218 9.020 0.34232100 7.0 1 2
## 25 (181147, 332823) 1.4 26 75 200 8.976 0.38580400 6.9 1 2
## 26 (181167, 332778) 1.5 22 76 194 8.973 0.42928900 6.3 1 2
## 27 (181008, 332777) 1.3 27 73 207 8.507 0.31511600 5.6 1 2
## 28 (180973, 332687) 1.3 24 67 180 8.743 0.32057400 4.4 1 2
## 29 (180916, 332753) 1.8 22 87 240 8.973 0.24986300 5.3 1 2
## 30 (181352, 332946) 1.5 21 65 180 9.043 0.48906400 4.8 1 2
## 31 (181133, 332570) 1.3 29 78 208 8.688 0.47277800 2.6 1 2
## 32 (180878, 332489) 1.3 21 64 198 8.727 0.28795700 1.0 1 2
## 33 (180829, 332450) 2.1 27 77 250 8.328 0.27162200 2.4 1 2
## 34 (180954, 332399) 1.2 26 80 192 7.971 0.38580700 1.9 1 2
## 35 (180956, 332318) 1.6 27 82 213 7.809 0.41841700 3.1 1 2
## 37 (180710, 332330) 3.0 32 97 321 6.986 0.24447400 1.6 1 2
## 38 (180632, 332445) 5.8 50 166 569 7.756 0.13570900 3.5 1 2
## 39 (180530, 332538) 7.9 67 217 833 7.784 0.04849650 8.1 1 1
## 40 (180478, 332578) 8.1 77 219 906 7.000 0.00000000 7.9 1 1
## 41 (180383, 332476) 14.1 108 405 1454 6.920 0.00135803 9.5 1 1
## 42 (180494, 332330) 2.4 32 102 298 7.516 0.13570900 1.4 1 2
## 43 (180561, 332193) 1.2 21 48 167 8.180 0.26622000 NA 1 2
## 44 (180451, 332175) 1.7 22 65 176 8.694 0.21184300 NA 1 2
## 45 (180410, 332031) 1.3 21 62 258 9.280 0.32057200 2.0 1 2
## 46 (180355, 332299) 4.2 51 281 746 7.940 0.08122200 5.1 1 2
## 47 (180292, 332157) 4.3 50 294 746 6.360 0.19008600 5.3 1 2
## 48 (180283, 332014) 3.1 38 211 464 7.780 0.28794100 4.5 1 2
## 49 (180282, 331861) 1.7 26 135 365 8.180 0.42382600 4.9 1 2
## 50 (180270, 331707) 1.7 24 112 282 9.420 0.55428900 4.5 1 2
## 51 (180199, 331591) 2.1 32 162 375 8.867 0.60322500 5.5 1 2
## 52 (180135, 331552) 1.7 24 94 222 8.292 0.61407100 3.4 1 2
## 53 (180237, 332351) 8.2 47 191 812 8.060 0.00135803 11.1 1 1
## 54 (180103, 332297) 17.0 128 405 1548 7.980 0.00000000 12.3 1 1
## 55 (179973, 332255) 12.0 117 654 1839 7.900 0.00543210 16.5 1 1
## 56 (179826, 332217) 9.4 104 482 1528 7.740 0.00543210 13.9 1 1
## 57 (179687, 332161) 8.2 76 276 933 7.552 0.00543210 8.1 1 1
## 58 (179792, 332035) 2.6 36 180 432 7.760 0.14657800 3.1 1 1
## 59 (179902, 332113) 3.5 34 207 550 6.740 0.13568400 5.8 1 1
## 60 (180100, 332213) 10.9 90 541 1571 6.680 0.07033330 10.2 1 1
## 61 (179604, 332059) 7.3 80 310 1190 7.400 0.04848310 12.0 1 1
## 62 (179526, 331936) 9.4 78 210 907 7.440 0.00543210 14.1 1 1
## 63 (179495, 331770) 8.3 77 158 761 7.360 0.00543210 14.5 1 1
## 64 (179489, 331633) 7.0 65 141 659 7.200 0.03166630 14.8 1 1
## 65 (179414, 331494) 6.8 66 144 643 7.220 0.01222430 13.3 1 1
## 66 (179334, 331366) 7.4 72 181 801 7.360 0.01222430 15.2 1 1
## 67 (179255, 331264) 6.6 75 173 784 5.180 0.03733950 11.4 1 1
## 69 (179470, 331125) 7.8 75 399 1060 5.800 0.21184600 9.0 1 1
## 75 (179692, 330933) 0.7 22 45 119 7.640 0.45103700 3.6 1 1
## 76 (179852, 330801) 3.4 55 325 778 6.320 0.57587700 6.9 1 1
## 79 (179140, 330955) 3.9 47 268 703 5.760 0.07568690 7.0 1 1
## 80 (179128, 330867) 3.5 46 252 676 6.480 0.12481000 6.2 1 1
## 81 (179065, 330864) 4.7 55 315 793 6.480 0.10302400 6.5 1 1
## 82 (179007, 330727) 3.9 49 260 685 6.320 0.15746900 5.7 1 1
## 83 (179110, 330758) 3.1 39 237 593 6.320 0.20097600 7.0 1 1
## 84 (179032, 330645) 2.9 45 228 549 6.160 0.20097600 7.3 1 1
## 85 (179095, 330636) 3.9 48 241 680 6.560 0.26622000 8.2 1 1
## 86 (179058, 330510) 2.7 36 201 539 6.900 0.29883500 4.3 1 1
## 87 (178810, 330666) 2.5 36 204 560 7.540 0.08122470 4.4 1 1
## 88 (178912, 330779) 5.6 68 429 1136 6.420 0.07035500 8.2 1 1
## 89 (178981, 330924) 9.4 88 462 1383 6.280 0.01222430 8.5 1 1
## 90 (179076, 331005) 10.8 85 333 1161 6.340 0.00000000 9.6 1 1
## 123 (180151, 330353) 18.1 76 464 1672 7.307 0.05377230 17.0 1 1
## 160 (179211, 331175) 6.3 63 159 765 5.700 0.05936620 12.8 1 1
## 163 (181118, 333214) 2.1 32 116 279 7.720 0.21184300 5.9 1 2
## 70 (179474, 331304) 1.8 25 81 241 7.932 0.12481000 2.9 2 2
## 71 (179559, 331423) 2.2 27 131 317 7.820 0.12481000 4.5 2 1
## 91 (179022, 330873) 2.8 36 216 545 8.575 0.09215160 10.7 2 1
## 92 (178953, 330742) 2.4 41 145 505 8.536 0.11394100 9.4 2 1
## 93 (178875, 330516) 2.6 33 163 420 8.504 0.17921600 9.0 2 1
## 94 (178803, 330349) 1.8 27 129 332 8.659 0.23359600 7.0 2 1
## 95 (179029, 330394) 2.0 38 148 400 7.633 0.33686100 6.5 2 1
## 96 (178605, 330406) 2.7 37 214 553 8.538 0.07035500 9.4 2 1
## 97 (178701, 330557) 2.7 34 226 577 7.680 0.05936620 10.2 2 1
## 98 (179547, 330245) 0.9 19 54 155 7.564 0.25534100 6.4 2 1
## 99 (179301, 330179) 0.9 22 70 224 7.760 0.36406700 7.6 2 1
## 100 (179405, 330567) 0.4 26 73 180 7.653 0.42929500 7.0 2 1
## 101 (179462, 330766) 0.8 25 87 226 7.951 0.38032800 5.6 2 1
## 102 (179293, 330797) 0.4 22 76 186 8.176 0.24987400 6.5 2 1
## 103 (179180, 330710) 0.4 24 81 198 8.468 0.26621200 6.6 2 1
## 104 (179206, 330398) 0.4 18 68 187 8.410 0.45103700 5.9 2 1
## 105 (179618, 330458) 0.8 23 66 199 7.610 0.30971000 6.5 2 1
## 106 (179782, 330540) 0.4 22 49 157 7.792 0.29335900 6.4 2 1
## 108 (179980, 330773) 0.4 23 63 203 8.760 0.53235100 7.2 2 2
## 109 (180067, 331185) 0.4 23 48 143 9.879 0.61951300 6.6 2 3
## 110 (180162, 331387) 0.2 23 51 136 9.097 0.68472500 4.3 2 2
## 111 (180451, 331473) 0.2 18 50 117 9.095 0.80974200 5.3 2 3
## 112 (180328, 331158) 0.4 20 39 113 9.717 0.88038900 4.1 2 3
## 113 (180276, 330963) 0.2 22 48 130 9.924 0.74959100 6.1 2 3
## 114 (180114, 330803) 0.2 27 64 192 9.404 0.57575200 7.5 2 3
## 115 (179881, 330912) 0.4 25 84 240 10.520 0.58148400 8.8 2 3
## 116 (179774, 330921) 0.2 30 67 221 8.840 0.49452000 5.7 2 3
## 117 (179657, 331150) 0.2 23 49 140 8.472 0.32058000 6.1 2 3
## 118 (179731, 331245) 0.2 24 48 128 9.634 0.33685100 7.1 2 3
## 119 (179717, 331441) 0.2 21 56 166 9.206 0.24985200 4.1 2 2
## 120 (179446, 331422) 0.2 24 65 191 8.470 0.07568690 6.0 2 1
## 121 (179524, 331565) 0.2 21 84 232 8.463 0.07568690 6.6 2 1
## 122 (179644, 331730) 0.2 23 75 203 9.691 0.16285300 6.8 2 1
## 124 (180321, 330366) 3.7 53 250 722 8.704 0.09749160 9.1 2 2
## 125 (180162, 331837) 0.2 33 81 210 9.420 0.44014200 5.9 2 2
## 126 (180029, 331720) 0.2 22 72 198 9.573 0.46190000 4.9 2 2
## 127 (179797, 331919) 0.2 23 86 139 9.555 0.22270100 7.1 2 1
## 128 (179642, 331955) 0.2 25 94 253 8.779 0.10302400 8.1 2 1
## 129 (179849, 332142) 1.2 30 244 703 8.540 0.09213530 8.3 2 1
## 130 (180265, 332297) 2.4 47 297 832 8.809 0.04848840 10.0 2 1
## 131 (180107, 332101) 0.2 31 96 262 9.523 0.16833100 5.9 2 1
## 132 (180462, 331947) 0.2 20 56 142 9.811 0.38581000 5.0 2 2
## 133 (180478, 331822) 0.2 16 49 119 9.604 0.48906400 4.5 2 2
## 134 (180347, 331700) 0.2 17 50 152 9.732 0.57602000 5.4 2 2
## 135 (180862, 333116) 0.4 26 148 415 9.518 0.08121940 2.3 2 1
## 136 (180700, 332882) 1.6 34 162 474 9.720 0.03733690 7.5 2 1
## 161 (180201, 331160) 0.8 18 37 126 9.036 0.77169800 4.6 2 3
## 162 (180173, 331923) 1.2 23 80 210 9.528 0.33682900 5.8 2 2
## 137 (180923, 332874) 0.2 20 80 220 9.155 0.22812300 4.4 3 1
## 138 (180467, 331694) 0.2 14 49 133 10.080 0.59776100 4.4 3 2
## 140 (179917, 331325) 0.8 46 42 141 9.970 0.44558000 4.5 3 2
## 141 (179822, 331242) 1.0 29 48 158 10.136 0.39667500 5.2 3 2
## 142 (179991, 331069) 0.8 19 41 129 10.320 0.58147800 4.6 3 3
## 143 (179120, 330578) 1.2 31 73 206 9.041 0.28796600 6.9 3 1
## 144 (179034, 330561) 2.0 27 146 451 7.860 0.23359600 7.0 3 1
## 145 (179085, 330433) 1.5 29 95 296 8.741 0.36406700 5.4 3 1
## 146 (179236, 330046) 1.1 22 72 189 7.822 0.33145400 6.2 3 1
## 147 (179456, 330072) 0.8 20 51 154 7.780 0.21184600 5.0 3 1
## 148 (179550, 329940) 0.8 20 54 169 8.121 0.10302900 5.1 3 1
## 149 (179445, 329807) 2.1 29 136 403 8.231 0.07035500 8.1 3 1
## 150 (179337, 329870) 2.5 38 170 471 8.351 0.14657600 8.0 3 1
## 151 (179245, 329714) 3.8 39 179 612 7.300 0.05377230 8.8 3 1
## 152 (179024, 329733) 3.2 35 200 601 7.536 0.11928600 9.3 3 1
## 153 (178786, 329822) 3.1 42 258 783 7.706 0.09214350 8.4 3 1
## 154 (179135, 329890) 1.5 24 93 258 8.070 0.24986300 7.7 3 1
## 155 (179030, 330082) 1.2 20 68 214 8.226 0.37494000 5.7 3 1
## 156 (179184, 330182) 0.8 20 49 166 8.128 0.42383700 4.7 3 1
## 157 (179085, 330292) 3.1 39 173 496 8.577 0.42383700 9.1 3 1
## 158 (178875, 330311) 2.1 31 119 342 8.429 0.27709000 6.5 3 1
## 159 (179466, 330381) 0.8 21 51 162 9.406 0.35860600 5.7 3 1
## 164 (180627, 330190) 2.7 27 124 375 8.261 0.01222430 5.5 3 3
## lime landuse dist.m logZn
## 1 1 Ah 50 3.009451
## 2 1 Ah 30 3.057286
## 3 1 Ah 150 2.806180
## 4 0 Ga 270 2.409933
## 5 0 Ah 380 2.429752
## 6 0 Ga 470 2.448706
## 7 0 Ah 240 2.539076
## 8 0 Ab 120 2.608526
## 9 0 Ab 240 2.540329
## 10 0 W 420 2.262451
## 11 0 Fh 400 2.276462
## 12 0 Ag 300 2.399674
## 13 1 W 20 3.039811
## 14 0 Ah 130 2.702431
## 15 0 Ah 220 2.513218
## 16 1 W 10 3.013680
## 17 1 W 10 2.782473
## 18 1 W 10 2.851870
## 19 1 W 10 2.866287
## 20 1 <NA> 10 3.022016
## 21 1 Am 80 2.828015
## 22 1 Am 140 2.604226
## 23 1 Ag 160 2.535294
## 24 0 Ah 440 2.338456
## 25 0 W 490 2.301030
## 26 0 W 530 2.287802
## 27 0 Ab 400 2.315970
## 28 0 Ag 400 2.255273
## 29 0 Ah 330 2.380211
## 30 0 Ag 630 2.255273
## 31 0 B 570 2.318063
## 32 0 Ag 390 2.296665
## 33 0 Ah 360 2.397940
## 34 0 B 500 2.283301
## 35 0 B 550 2.328380
## 37 0 Ab 340 2.506505
## 38 0 Ab 210 2.755112
## 39 1 Am 60 2.920645
## 40 1 W 10 2.957128
## 41 1 W 20 3.162564
## 42 0 Am 170 2.474216
## 43 0 Ga 320 2.222716
## 44 0 W 260 2.245513
## 45 0 Ah 360 2.411620
## 46 0 Ah 100 2.872739
## 47 0 Am 200 2.872739
## 48 0 Ah 320 2.666518
## 49 0 Ah 480 2.562293
## 50 0 Bw 660 2.450249
## 51 0 Bw 690 2.574031
## 52 0 Ab 710 2.346353
## 53 1 Ah 10 2.909556
## 54 1 W 10 3.189771
## 55 1 W 10 3.264582
## 56 1 W 10 3.184123
## 57 1 W 20 2.969882
## 58 0 Fw 200 2.635484
## 59 0 Fw 140 2.740363
## 60 1 Fw 70 3.196176
## 61 1 W 20 3.075547
## 62 1 W 10 2.957607
## 63 1 Fw 10 2.881385
## 64 1 W 20 2.818885
## 65 1 Ah 10 2.808211
## 66 1 W 20 2.903633
## 67 1 W 20 2.894316
## 69 0 Ah 270 3.025306
## 75 1 Fw 560 2.075547
## 76 0 Bw 750 2.890980
## 79 1 Ab 80 2.846955
## 80 1 Ab 130 2.829947
## 81 0 W 110 2.899273
## 82 0 W 200 2.835691
## 83 1 Ah 260 2.773055
## 84 0 W 270 2.739572
## 85 0 W 320 2.832509
## 86 0 Ah 360 2.731589
## 87 1 Am 80 2.748188
## 88 1 W 100 3.055378
## 89 1 W 70 3.140822
## 90 1 W 20 3.064832
## 123 1 W 50 3.223236
## 160 1 W 80 2.883661
## 163 0 W 290 2.445604
## 70 1 Ah 160 2.382017
## 71 0 W 160 2.501059
## 91 0 W 140 2.736397
## 92 0 W 150 2.703291
## 93 0 Ah 220 2.623249
## 94 0 Am 280 2.521138
## 95 1 Am 450 2.602060
## 96 1 Am 70 2.742725
## 97 0 Am 70 2.761176
## 98 0 W 340 2.190332
## 99 0 W 470 2.350248
## 100 0 Am 630 2.255273
## 101 0 Am 460 2.354108
## 102 0 Am 320 2.269513
## 103 0 Ah 320 2.296665
## 104 0 W 540 2.271842
## 105 0 W 420 2.298853
## 106 0 SPO 380 2.195900
## 108 0 W 500 2.307496
## 109 0 Am 760 2.155336
## 110 0 Ah 750 2.133539
## 111 0 Fw 1000 2.068186
## 112 0 Ah 860 2.053078
## 113 0 Ah 680 2.113943
## 114 0 Fw 500 2.283301
## 115 0 STA 650 2.380211
## 116 0 DEN 630 2.344392
## 117 0 Fw 410 2.146128
## 118 0 Ah 390 2.107210
## 119 0 Ah 310 2.220108
## 120 0 Ah 70 2.281033
## 121 0 W 70 2.365488
## 122 1 STA 150 2.307496
## 124 0 Bw 80 2.858537
## 125 0 Ah 450 2.322219
## 126 0 Aa 530 2.296665
## 127 0 W 240 2.143015
## 128 1 Tv 70 2.403121
## 129 0 Fw 70 2.846955
## 130 0 Ah 60 2.920123
## 131 0 Ah 190 2.418301
## 132 0 Ah 450 2.152288
## 133 0 Am 550 2.075547
## 134 0 Am 650 2.181844
## 135 0 Am 100 2.618048
## 136 0 W 170 2.675778
## 161 1 Ah 860 2.100371
## 162 0 W 410 2.322219
## 137 0 Aa 290 2.342423
## 138 0 Am 680 2.123852
## 140 0 Am 540 2.149219
## 141 0 Am 480 2.198657
## 142 0 W 720 2.110590
## 143 0 W 380 2.313867
## 144 0 W 310 2.654177
## 145 0 Ah 430 2.471292
## 146 0 Ah 370 2.276462
## 147 0 Fw 290 2.187521
## 148 0 W 150 2.227887
## 149 0 Bw 70 2.605305
## 150 0 Bw 220 2.673021
## 151 0 W 80 2.786751
## 152 0 W 120 2.778874
## 153 0 Ah 120 2.893762
## 154 0 Am 260 2.411620
## 155 0 Ah 440 2.330414
## 156 0 Am 540 2.220108
## 157 0 Ah 520 2.695482
## 158 0 Ah 350 2.534026
## 159 0 W 460 2.209515
## 164 0 W 40 2.574031
##
## $col
## [1] "white"
##
## $pch
## [1] 1
##
## $cex
## [1] 2.2229473 2.4817836 1.3920609 0.5589995 0.5851006 0.6112017 0.7525829
## [8] 0.8830886 0.7547580 0.3980424 0.4110930 0.5459489 2.3839043 1.0962480
## [15] 0.7090810 2.2446982 1.3181077 1.5464927 1.5986949 2.2882001 1.4638390
## [22] 0.8743883 0.7460576 0.4741707 0.4350190 0.4219685 0.4502447 0.3915171
## [29] 0.5220228 0.3915171 0.4524198 0.4306688 0.5437738 0.4176183 0.4632953
## [36] 0.6982055 1.2376291 1.8118543 1.9706362 3.1625884 0.6481784 0.3632409
## [43] 0.3828167 0.5611746 1.6226210 1.6226210 1.0092442 0.7939097 0.6133768
## [50] 0.8156607 0.4828711 1.7661773 3.3670473 4.0000000 3.3235454 2.0293638
## [57] 0.9396411 1.1963023 3.4170745 2.5883632 1.9728113 1.6552474 1.4333877
## [64] 1.3985862 1.7422512 1.7052746 2.3056009 0.2588363 1.6922240 1.5290919
## [71] 1.4703643 1.7248505 1.4899402 1.2898314 1.1941272 1.4790647 1.1723763
## [78] 1.2180533 2.4709081 3.0081566 2.5252855 3.6367591 1.6639478 0.6068515
## [85] 0.5241979 0.6895052 1.1854269 1.0984231 0.9135400 0.7221316 0.8700381
## [92] 1.2028276 1.2550299 0.3371397 0.4872213 0.3915171 0.4915715 0.4045677
## [99] 0.4306688 0.4067428 0.4328439 0.3414899 0.4415443 0.3110386 0.2958129
## [106] 0.2544861 0.2457858 0.2827624 0.4176183 0.5220228 0.4806960 0.3045133
## [113] 0.2784122 0.3610658 0.4154432 0.5046221 0.4415443 1.5704187 0.4567700
## [120] 0.4306688 0.3023382 0.5502991 1.5290919 1.8096792 0.5698749 0.3088635
## [127] 0.2588363 0.3306145 0.9026645 1.0309951 0.2740620 0.4567700 0.4785209
## [134] 0.2892877 0.3066884 0.3436650 0.2805873 0.4480696 0.9809679 0.6438282
## [141] 0.4110930 0.3349647 0.3675911 0.8765633 1.0244698 1.3311582 1.3072322
## [148] 1.7030995 0.5611746 0.4654704 0.3610658 1.0788472 0.7438825 0.3523654
## [155] 0.8156607print(spplot(OK, 'pred',asp = 1, col.regions = rev(bpy.colors(64)), main = 'Predicciones Esferico', sp.layout = list(pts.s)))
Prediccion con el semivariograma exponencial
OK1 = krige(logZn ~ 1, locations = meuse, newdata = meuse.grid, model = va1)## [using ordinary kriging]OK1$pred = 10^(OK1$var1.pred) # ELIMINAMOS EL LOGARITMO
pts.s = list('sp.points', meuse, col = 'white', pch=1, cex = 4*meuse$zinc/max(meuse$zinc))
print(spplot(OK1, 'pred',asp = 1, col.regions = rev(bpy.colors(64)), main = 'Predicciones Exponencial', sp.layout = list(pts.s)))
5. Validación cruzada definir el mejor modelo teorico
modelo esferico
valida = krige.cv(log(zinc) ~1, meuse, va, nmax = 40, nfold = 5)
head(valida)## coordinates var1.pred var1.var observed residual zscore fold
## 1 (181072, 333611) 6.781272 0.03439635 6.929517 0.148244366 0.79932241 1
## 2 (181025, 333558) 6.775021 0.03307464 7.039660 0.264639268 1.45514652 3
## 3 (181165, 333537) 6.313826 0.03421126 6.461468 0.147642015 0.79822516 5
## 4 (181298, 333484) 6.077401 0.04289458 5.549076 -0.528325009 -2.55093831 3
## 5 (181307, 333330) 5.590255 0.03284617 5.594711 0.004456147 0.02458766 5
## 6 (181390, 333260) 5.420330 0.04360783 5.638355 0.218025025 1.04405672 2bubble(valida, 'residual', main = 'Log(zinc): 5-fold cv residuales')
# Error medio
ME = round(mean(valida$residual),3)
ME## [1] 0.013# Mean Absolute Error
MAE = round(mean(abs(valida$residual)),3)
MAE## [1] 0.314# Root Mean Squre Error (RMSE)
RMSE = round(sqrt(mean(valida$residual^2)),3)
RMSE## [1] 0.425# Mean Squared Deviation Ratio (MSDR)
MSDR = mean(valida$residual^2/valida$var1.var)
MSDR## [1] 4.653911modelo exponencial
valida1 = krige.cv(log(zinc) ~1, meuse, va1, nmax = 40, nfold = 5)
bubble(valida1, 'residual', main = 'Log(zinc): 5-fold cv residuales')
# Error medio
ME1 = round(mean(valida1$residual),3)
ME1## [1] 0.026# Mean Absolute Error
MAE1 = round(mean(abs(valida1$residual)),3)
MAE1## [1] 0.313# Root Mean Squre Error (RMSE)
RMSE1 = round(sqrt(mean(valida1$residual^2)),3)
RMSE1## [1] 0.43# Mean Squared Deviation Ratio (MSDR)
MSDR1 = mean(valida1$residual^2/valida1$var1.var)
MSDR1## [1] 4.6249816. Comprobar los valores predichos y observados
head(valida)## coordinates var1.pred var1.var observed residual zscore fold
## 1 (181072, 333611) 6.781272 0.03439635 6.929517 0.148244366 0.79932241 1
## 2 (181025, 333558) 6.775021 0.03307464 7.039660 0.264639268 1.45514652 3
## 3 (181165, 333537) 6.313826 0.03421126 6.461468 0.147642015 0.79822516 5
## 4 (181298, 333484) 6.077401 0.04289458 5.549076 -0.528325009 -2.55093831 3
## 5 (181307, 333330) 5.590255 0.03284617 5.594711 0.004456147 0.02458766 5
## 6 (181390, 333260) 5.420330 0.04360783 5.638355 0.218025025 1.04405672 2cor.test(valida$var1.pred,valida$observed)##
## Pearson's product-moment correlation
##
## data: valida$var1.pred and valida$observed
## t = 16.896, df = 153, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.7438575 0.8556859
## sample estimates:
## cor
## 0.8068842plot(valida$observed,valida$var1.pred)