Construcción de un mapa de predicción por medio de Geoestadística

Carlos Eduardo Caceres Gonzalez ccaceresg@unal.edu.co

Junio 03 del 2021

—————————————————————————————————————————————————————————————————

Paquetes requeridos

#setwd("C:/Academic/GEOESTADISTICA/Clase 03062021/R")
source("C:/Academic/GEOESTADISTICA/Clase 03062021/R/programas estadistica en R.txt")
require(car)

## Loading required package: car

## Warning: package 'car' was built under R version 4.0.4

## Loading required package: carData

require(MASS)

## Loading required package: MASS

require(akima)

## Loading required package: akima

## Warning: package 'akima' was built under R version 4.0.5

require(gstat)

## Loading required package: gstat

## Warning: package 'gstat' was built under R version 4.0.4

require(geoR)

## Loading required package: geoR

## Warning: package 'geoR' was built under R version 4.0.5

## --------------------------------------------------------------
##  Analysis of Geostatistical Data
##  For an Introduction to geoR go to http://www.leg.ufpr.br/geoR
##  geoR version 1.8-1 (built on 2020-02-08) is now loaded
## --------------------------------------------------------------

require(lattice)

## Loading required package: lattice

require(maptools)

## Loading required package: maptools

## Warning: package 'maptools' was built under R version 4.0.4

## Loading required package: sp

## Warning: package 'sp' was built under R version 4.0.4

## Checking rgeos availability: TRUE

require(rgdal)

## Loading required package: rgdal

## Warning: package 'rgdal' was built under R version 4.0.4

## rgdal: version: 1.5-23, (SVN revision 1121)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 3.2.1, released 2020/12/29
## Path to GDAL shared files: C:/Users/GEOMATICS/Documents/R/win-library/4.0/rgdal/gdal
## GDAL binary built with GEOS: TRUE 
## Loaded PROJ runtime: Rel. 7.2.1, January 1st, 2021, [PJ_VERSION: 721]
## Path to PROJ shared files: C:/Users/GEOMATICS/Documents/R/win-library/4.0/rgdal/proj
## PROJ CDN enabled: FALSE
## Linking to sp version:1.4-5
## To mute warnings of possible GDAL/OSR exportToProj4() degradation,
## use options("rgdal_show_exportToProj4_warnings"="none") before loading rgdal.
## Overwritten PROJ_LIB was C:/Users/GEOMATICS/Documents/R/win-library/4.0/rgdal/proj

require(sf)

## Loading required package: sf

## Warning: package 'sf' was built under R version 4.0.5

## Linking to GEOS 3.9.0, GDAL 3.2.1, PROJ 7.2.1

—————————————————————————————————————————————————————————————————

Mapa de ubiacion

#install.packages("leaflet")
library(leaflet)

## Warning: package 'leaflet' was built under R version 4.0.5

library(rgdal)
library(sp)

Samples <- read_sf("C:/Academic/GEOESTADISTICA/Clase 03062021/R/ESTACIONES.shp")

Area_Study <- read_sf("C:/Academic/GEOESTADISTICA/Clase 03062021/R/CUNDINAMARCA_WGS84.shp")

leaflet() %>% addTiles() %>% addProviderTiles(providers$OpenStreetMap) %>%   
  addPolygons(data = Samples, fill = TRUE, stroke = TRUE, color = "#03F", 
              popup = paste0("<b>",Samples$prec,"</b>"), group = "Estaciones") %>% 
  addPolygons(data = Area_Study, fill = TRUE, stroke = TRUE, color = "#f93", 
              popup = "<b>Area de Estudio</b>", group = "Area de Estudio") %>% 
  # add a legend
  addLegend("bottomright", colors = c("#03F", "#f93"), labels = c("Estaciones","Area de Estudio")) %>%   
  # add layers control
  addLayersControl(
    overlayGroups = c("Area de Estudio", "Estaciones"),
    options = layersControlOptions(collapsed = FALSE)
  )

\(\text{El mapa muestra la distribucion espacial de las estaciones en el Departamento de Cundinamarca.}\ \) \(\text{Se observa una distribucion aleatoria.}\ \)

—————————————————————————————————————————————————————————————————

Datos

trellis.par.set(sp.theme())
cundinamarca = readShapePoly("C:/Academic/GEOESTADISTICA/Clase 03062021/R/CUNDINAMARCA.shp")

## Warning: readShapePoly is deprecated; use rgdal::readOGR or sf::st_read

poligonos = polygons(cundinamarca)
datos = read.table("C:/Academic/GEOESTADISTICA/Clase 03062021/R/precipitacion.txt", sep = "\t", dec = ",", header = T)
datos

##             x         y prec
## 1    981235.0  990896.0  376
## 2   1037799.9 1081782.0  446
## 3   1051832.5 1067956.6  454
## 4   1025594.0 1086730.6  461
## 5   1031043.7 1069912.7  523
## 6   1016916.6  984470.8  535
## 7   1016379.2 1007482.9  561
## 8   1001576.7 1042491.9  578
## 9   1029281.6 1094108.2  601
## 10  1024490.7 1072608.4  625
## 11  1031618.4 1073916.2  645
## 12  1020812.9 1063329.1  646
## 13   976188.5 1001270.1  647
## 14   964271.5  951600.3  648
## 15  1018979.8 1076753.6  694
## 16  1020064.8 1042494.7  697
## 17  1045913.1 1083058.7  706
## 18  1004415.5 1048917.0  711
## 19   970456.7 1015061.8  717
## 20  1051459.5 1077533.3  731
## 21  1017237.6 1083885.9  733
## 22   987588.7 1025739.5  735
## 23  1027698.7 1062004.4  736
## 24  1014522.0 1029599.3  749
## 25  1017129.1 1074596.7  753
## 26   967829.1 1019497.3  767
## 27   973843.1 1022226.5  769
## 28   980345.3 1027289.2  786
## 29   947223.3 1029992.5  787
## 30  1015638.6  993184.4  794
## 31  1017924.9 1083443.7  801
## 32   979380.2  985367.3  802
## 33  1010830.8 1011164.4  813
## 34   977462.0 1030109.9  825
## 35  1013127.5 1073058.7  838
## 36   991571.6 1030206.6  842
## 37  1008287.4 1034640.9  843
## 38  1040382.6 1062784.6  844
## 39  1008971.4 1064609.3  856
## 40  1038811.9 1076076.3  859
## 41  1005638.7 1004628.4  868
## 42  1007720.4 1047114.7  868
## 43  1037259.5 1096910.1  876
## 44  1017562.9 1068912.7  879
## 45  1022362.6 1071070.5  880
## 46  1017527.6 1076897.0  888
## 47   954102.0  932562.5  892
## 48   983095.2 1022223.8  895
## 49  1023048.5 1041268.1  898
## 50   977006.8  991339.4  911
## 51   994190.5 1044338.9  919
## 52  1002965.3  984347.4  920
## 53  1025589.0 1099636.0  923
## 54  1030336.5 1089652.0  928
## 55  1030071.1 1062790.5  936
## 56  1012929.3 1065516.8  962
## 57   983099.3 1040658.2  969
## 58   988680.6  982158.6  979
## 59   936058.3  967939.0 1005
## 60  1017495.9 1070814.8 1008
## 61  1019412.7 1074752.1 1009
## 62  1041779.9 1043432.8 1011
## 63   954999.9 1066613.7 1019
## 64  1018200.7 1088575.0 1031
## 65  1018213.1 1040658.5 1056
## 66   965645.1 1022992.8 1060
## 67  1006270.1 1017334.4 1065
## 68  1038547.1 1053560.7 1086
## 69   936837.5  985389.2 1116
## 70  1002851.7 1053229.7 1148
## 71   994188.9 1024068.9 1160
## 72  1042238.7 1055409.6 1163
## 73  1020384.0 1051087.2 1185
## 74  1053538.6 1054952.7 1189
## 75  1033010.1 1036970.1 1221
## 76   950887.5  923352.5 1269
## 77  1059333.5 1047205.0 1274
## 78   917871.3  977167.0 1275
## 79  1037649.4 1052609.2 1280
## 80  1051467.5 1019485.8 1286
## 81   928478.0 1065607.3 1288
## 82  1004780.0 1071365.7 1359
## 83  1042256.7 1024080.5 1377
## 84   916729.8  979015.0 1392
## 85   938724.3 1027776.1 1396
## 86  1021921.8 1025907.6 1400
## 87  1013837.7 1069619.6 1400
## 88   953581.7 1084441.4 1446
## 89   935019.1 1027779.2 1448
## 90  1045497.7 1021317.7 1490
## 91   959580.2  986237.3 1542
## 92   953823.0 1080758.7 1561
## 93   938933.8 1052934.9 1664
## 94  1003432.9  948496.8 1665
## 95  1015067.9 1071223.4 1700
## 96   970833.0 1041248.1 1709
## 97   973833.8  994580.5 1752
## 98   954851.7  943653.7 1752
## 99   940651.8 1127305.4 1773
## 100 1043209.6 1006265.7 1778
## 101 1041474.9 1053108.9 1785
## 102  957193.2 1005646.2 1814
## 103  964610.7 1040664.8 1900
## 104 1020244.9  990365.6 1909
## 105  971088.3  982284.6 1957
## 106 1038381.7 1030669.3 2022
## 107 1012001.9  969850.9 2049
## 108 1051495.2 1027768.8 2068
## 109 1060499.5 1018099.0 2070
## 110  968532.8  975385.1 2105
## 111  974815.0 1073526.5 2107
## 112 1034856.4 1025912.4 2147
## 113  956890.7  979504.1 2168
## 114  961471.3 1039527.3 2235
## 115 1060781.1 1012735.7 2264
## 116  965984.4 1083361.4 2278
## 117 1061790.2 1013411.1 2286
## 118  981250.0 1051717.1 2421
## 119  960917.4 1040666.6 2461
## 120 1006374.2  964111.1 2468
## 121  957244.5 1092268.7 2605
## 122 1053960.4 1039592.3 2664
## 123 1066300.0 1020404.0 2759
## 124 1045956.4 1014859.8 2924
## 125 1073019.5 1024380.4 3269
## 126 1065782.8 1015515.6 3408
## 127 1073725.1  981704.3 3466
## 128 1069102.8 1012731.6 3668
## 129 1065598.1 1010881.7 3710
## 130 1079543.0  991541.3 3750
## 131 1074979.0 1028551.6 3776
## 132  949034.2  983069.4 3974
## 133  991600.4 1085808.4 3986

#xy = SpatialPoints(datos[c("x", "y")])
#plot(poligonos, main = "Distribución de Puntos")
#points(xy, pch = 3, cex = 0.5, col = "red")
#datos

—————————————————————————————————————————————————————————————————

Analisis exploratorio

#tabla.estadisticas.cont(datos)

x <- datos$prec

nombres <- c("Mínimo","Máximo","Promedio",  "Mediana","Moda","Desviación estándar", "Desviación mediana", "Asimetría","Kurtosis","Coef. Var. Promedio(%)","Coef. Var. Mediana(%)")

tabla1<- rbind(min(x),max(x),mean(x,na.rm=TRUE),  median(x), moda(x), sd(x,na.rm=TRUE), meda(x), asimetria(x), kurtosis(x), 100*cvariacion(x), 100*cmvariacion(x))

rownames(tabla1) <- nombres

tabla1

##                               [,1]
## Mínimo                  376.000000
## Máximo                 3986.000000
## Promedio               1403.052632
## Mediana                1065.000000
## Moda                    868.000000
## Desviación estándar     837.770267
## Desviación mediana      335.000000
## Asimetría                 1.391488
## Kurtosis                  4.364208
## Coef. Var. Promedio(%)   59.710537
## Coef. Var. Mediana(%)    31.455399

op = par(mfrow = c(1, 3))
hist(datos$prec, freq = FALSE, main = "", xlab = "Precipitacion", ylab = "Frecuencia")
curve(dnorm(x, mean(datos$prec), sd(datos$prec)), add = T)
boxplot(datos$prec)
qqPlot(datos$prec, ylab = "Precipitacion")

## [1] 133 132

par(op)

\(\text{Se observa en el histograma contruido una asimetria en la distribucion con concentracion hacia los valores bajos.}\ \) \(\text{Se aprecia una heterogeneidad de los valores de precipitacion.}\ \)

breaks = c(min(datos$prec), quantile(datos$prec,probs = c(0.2, 0.4, 0.6, 0.8),type = 6), max(datos$prec))
datos1 = datos[, c(1, 2, 3)]
coordinates(datos1) = c("x", "y")
spplot(datos1, "prec", cuts = breaks)

spplot(datos1, "prec", cuts = breaks[c(1, 2, 6)], col.regions = c(1,0))

spplot(datos1, "prec", cuts = breaks[c(1, 3, 6)], col.regions = c(1,0))

spplot(datos1, "prec", cuts = breaks[c(1, 4, 6)], col.regions = c(1,0))

spplot(datos1, "prec", cuts = breaks[c(1, 5, 6)], col.regions = c(1,0))

int.prec = interp(x = datos$x, y = datos$y, z = datos$prec)
contour(int.prec)

print(levelplot(int.prec$z))

persp(int.prec$x, int.prec$y, int.prec$z, xlab = "Este", ylab = "Norte", zlab = "Precipitacion",phi = 30, theta = 20, col = "lightblue", expand = 0.5, ticktype = "detailed")

—————————————————————————————————————————————————————————————————

Transformacion Box cox

summary(powerTransform(datos$prec))

## bcPower Transformation to Normality 
##            Est Power Rounded Pwr Wald Lwr Bnd Wald Upr Bnd
## datos$prec   -0.3288        -0.5      -0.6297       -0.028
## 
## Likelihood ratio test that transformation parameter is equal to 0
##  (log transformation)
##                            LRT df     pval
## LR test, lambda = (0) 4.589895  1 0.032161
## 
## Likelihood ratio test that no transformation is needed
##                            LRT df       pval
## LR test, lambda = (1) 73.48677  1 < 2.22e-16

lambda = powerTransform(datos$prec)$lambda
datos2 = data.frame(datos[, 1:2], tprec = ((datos[, 3]^lambda) - 1)/lambda)

#tabla.estadisticas.cont(datos2)
op = par(mfrow = c(1, 3))

hist(datos2$tprec, freq=FALSE, main="", xlab="Transformacion de la precipitacion", ylab="Frecuencia")
curve(dnorm(x, mean(datos2$tprec), sd(datos2$tprec)), add = T)
boxplot(datos2$tprec)
qqPlot(datos2$tprec, ylab = "Precipitacion")

## [1] 1 2

par(op)

breaks = c(min(datos2$tprec), quantile(datos2$tprec,probs=c(0.2,0.4,0.6,0.8),type = 6), max(datos2$tprec))
datos3 = datos2[, c(1, 2, 3)]
coordinates(datos3) = c("x", "y")
print(spplot(datos3, "tprec", cuts = breaks))

for (i in 2:5) {
    fnfile <- paste("EjemploClase0", i, ".pdf", sep = "")
    pdf(file = fnfile, width = 6, height = 6)
    print(spplot(datos3, "tprec", cuts = breaks[c(1, i, 6)], col.regions = c(1, 0)))
    dev.off()
    cat("\\includegraphics{", fnfile, "}\n\n", sep = "")
    }

## \includegraphics{EjemploClase02.pdf}

## \includegraphics{EjemploClase03.pdf}

## \includegraphics{EjemploClase04.pdf}

## \includegraphics{EjemploClase05.pdf}

int.tprec = interp(x = datos2$x, y = datos2$y, z = datos2$tprec)
contour(int.tprec)

print(levelplot(int.tprec$z))

persp(int.tprec$x, int.tprec$y, int.tprec$z, xlab = "Este", ylab = "Norte", zlab = "Transformacion de la precipitacion", phi = 30, theta = 20, col = "lightblue", expand = 0.5, ticktype = "detailed")

—————————————————————————————————————————————————————————————————

Analisis de tendencia

scatterplot(tprec~x, reg.line=lm, smooth=TRUE, spread=TRUE, boxplots=FALSE, span=0.5, data=datos2)

## Warning in plot.window(...): "reg.line" is not a graphical parameter

## Warning in plot.window(...): "spread" is not a graphical parameter

## Warning in plot.window(...): "span" is not a graphical parameter

## Warning in plot.xy(xy, type, ...): "reg.line" is not a graphical parameter

## Warning in plot.xy(xy, type, ...): "spread" is not a graphical parameter

## Warning in plot.xy(xy, type, ...): "span" is not a graphical parameter

## Warning in axis(side = side, at = at, labels = labels, ...): "reg.line" is not a
## graphical parameter

## Warning in axis(side = side, at = at, labels = labels, ...): "spread" is not a
## graphical parameter

## Warning in axis(side = side, at = at, labels = labels, ...): "span" is not a
## graphical parameter

## Warning in axis(side = side, at = at, labels = labels, ...): "reg.line" is not a
## graphical parameter

## Warning in axis(side = side, at = at, labels = labels, ...): "spread" is not a
## graphical parameter

## Warning in axis(side = side, at = at, labels = labels, ...): "span" is not a
## graphical parameter

## Warning in box(...): "reg.line" is not a graphical parameter

## Warning in box(...): "spread" is not a graphical parameter

## Warning in box(...): "span" is not a graphical parameter

## Warning in title(...): "reg.line" is not a graphical parameter

## Warning in title(...): "spread" is not a graphical parameter

## Warning in title(...): "span" is not a graphical parameter

scatterplot(tprec~y, reg.line=lm, smooth=TRUE, spread=TRUE, boxplots=FALSE, span=0.5, data=datos2)

## Warning in plot.window(...): "reg.line" is not a graphical parameter

## Warning in plot.window(...): "spread" is not a graphical parameter

## Warning in plot.window(...): "span" is not a graphical parameter

## Warning in plot.xy(xy, type, ...): "reg.line" is not a graphical parameter

## Warning in plot.xy(xy, type, ...): "spread" is not a graphical parameter

## Warning in plot.xy(xy, type, ...): "span" is not a graphical parameter

## Warning in axis(side = side, at = at, labels = labels, ...): "reg.line" is not a
## graphical parameter

## Warning in axis(side = side, at = at, labels = labels, ...): "spread" is not a
## graphical parameter

## Warning in axis(side = side, at = at, labels = labels, ...): "span" is not a
## graphical parameter

## Warning in axis(side = side, at = at, labels = labels, ...): "reg.line" is not a
## graphical parameter

## Warning in axis(side = side, at = at, labels = labels, ...): "spread" is not a
## graphical parameter

## Warning in axis(side = side, at = at, labels = labels, ...): "span" is not a
## graphical parameter

## Warning in box(...): "reg.line" is not a graphical parameter

## Warning in box(...): "spread" is not a graphical parameter

## Warning in box(...): "span" is not a graphical parameter

## Warning in title(...): "reg.line" is not a graphical parameter

## Warning in title(...): "spread" is not a graphical parameter

## Warning in title(...): "span" is not a graphical parameter

mod1reg = lm(tprec ~ x + y + I(x * y) + I(x^2) + I(y^2), data = datos2)
summary(mod1reg)

## 
## Call:
## lm(formula = tprec ~ x + y + I(x * y) + I(x^2) + I(y^2), data = datos2)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.125469 -0.025059 -0.000915  0.024278  0.121931 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  3.493e+00  3.682e+00   0.949   0.3446    
## x           -1.176e-05  5.432e-06  -2.165   0.0322 *  
## y            1.018e-05  4.099e-06   2.484   0.0143 *  
## I(x * y)    -1.366e-11  2.511e-12  -5.440 2.64e-07 ***
## I(x^2)       1.295e-11  2.258e-12   5.734 6.76e-08 ***
## I(y^2)       1.537e-12  1.962e-12   0.783   0.4350    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.04042 on 127 degrees of freedom
## Multiple R-squared:  0.3958, Adjusted R-squared:  0.372 
## F-statistic: 16.64 on 5 and 127 DF,  p-value: 1.288e-12

mod2reg = stepAIC(mod1reg, scope = list(upper = mod1reg$formula, lower = ~1), direction = "both")

## Start:  AIC=-847.57
## tprec ~ x + y + I(x * y) + I(x^2) + I(y^2)
## 
##            Df Sum of Sq     RSS     AIC
## - I(y^2)    1  0.001002 0.20851 -848.93
## <none>                  0.20751 -847.57
## - x         1  0.007660 0.21517 -844.75
## - y         1  0.010083 0.21759 -843.26
## - I(x * y)  1  0.048348 0.25586 -821.71
## - I(x^2)    1  0.053731 0.26124 -818.95
## 
## Step:  AIC=-848.93
## tprec ~ x + y + I(x * y) + I(x^2)
## 
##            Df Sum of Sq     RSS     AIC
## <none>                  0.20851 -848.93
## - x         1  0.008919 0.21743 -845.36
## - y         1  0.046848 0.25536 -823.97
## - I(x * y)  1  0.048359 0.25687 -823.19
## - I(x^2)    1  0.054432 0.26294 -820.08

summary(mod2reg)

## 
## Call:
## lm(formula = tprec ~ x + y + I(x * y) + I(x^2), data = datos2)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.126152 -0.025600 -0.000092  0.027088  0.124705 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  2.531e+00  3.466e+00   0.730   0.4666    
## x           -1.250e-05  5.342e-06  -2.340   0.0208 *  
## y            1.279e-05  2.385e-06   5.363 3.70e-07 ***
## I(x * y)    -1.310e-11  2.405e-12  -5.449 2.51e-07 ***
## I(x^2)       1.302e-11  2.252e-12   5.781 5.38e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.04036 on 128 degrees of freedom
## Multiple R-squared:  0.3929, Adjusted R-squared:  0.3739 
## F-statistic: 20.71 on 4 and 128 DF,  p-value: 3.526e-13

—————————————————————————————————————————————————————————————————

Analisis de residuales

datos4 = data.frame(datos[, 1:2], res = mod2reg$residuals)
#tabla.estadisticas.cont(datos4)
op = par(mfrow = c(1, 3))
hist(datos4$res, freq = FALSE, main = "", xlab = "Residuales", ylab = "Frecuencia")
curve(dnorm(x, mean(datos4$res), sd(datos4$res)), add = T)
boxplot(datos4$res)
qqPlot(datos4$res, ylab = "Residuales")

## [1]   1 133

par(op)

breaks = c(min(datos4$res), quantile(datos4$res,probs=c(0.2,0.4,0.6,0.8), type = 6),max(datos4$res))
datos5 = datos4[, c(1, 2, 3)]
coordinates(datos5) = c("x", "y")
print(spplot(datos5, "res", cuts = breaks))

for (i in 2:5){
fnfile <- paste("EjemploClase", i, ".pdf", sep = "")
pdf(file = fnfile, width = 6, height = 6)
print(spplot(datos5, "res", cuts = breaks[c(1, i, 6)], col.regions = c(1,0)))
dev.off()
cat("\\includegraphics{", fnfile, "}\n\n", sep = "")
}

## \includegraphics{EjemploClase2.pdf}

## \includegraphics{EjemploClase3.pdf}

## \includegraphics{EjemploClase4.pdf}

## \includegraphics{EjemploClase5.pdf}

int.res = interp(x = datos4$x, y = datos4$y, z = datos4$res)
contour(int.res)

print(levelplot(int.res$z))

persp(int.res$x, int.res$y, int.res$z, xlab = "Este", ylab = "Norte", zlab = "Residuales",phi = 30, theta = 20, col = "lightblue", expand = 0.5, ticktype = "detailed")

scatterplot(res~x, reg.line=lm, smooth=TRUE, spread=TRUE, boxplots=FALSE, span=0.5, data=datos4)

## Warning in plot.window(...): "reg.line" is not a graphical parameter

## Warning in plot.window(...): "spread" is not a graphical parameter

## Warning in plot.window(...): "span" is not a graphical parameter

## Warning in plot.xy(xy, type, ...): "reg.line" is not a graphical parameter

## Warning in plot.xy(xy, type, ...): "spread" is not a graphical parameter

## Warning in plot.xy(xy, type, ...): "span" is not a graphical parameter

## Warning in axis(side = side, at = at, labels = labels, ...): "reg.line" is not a
## graphical parameter

## Warning in axis(side = side, at = at, labels = labels, ...): "spread" is not a
## graphical parameter

## Warning in axis(side = side, at = at, labels = labels, ...): "span" is not a
## graphical parameter

## Warning in axis(side = side, at = at, labels = labels, ...): "reg.line" is not a
## graphical parameter

## Warning in axis(side = side, at = at, labels = labels, ...): "spread" is not a
## graphical parameter

## Warning in axis(side = side, at = at, labels = labels, ...): "span" is not a
## graphical parameter

## Warning in box(...): "reg.line" is not a graphical parameter

## Warning in box(...): "spread" is not a graphical parameter

## Warning in box(...): "span" is not a graphical parameter

## Warning in title(...): "reg.line" is not a graphical parameter

## Warning in title(...): "spread" is not a graphical parameter

## Warning in title(...): "span" is not a graphical parameter

scatterplot(res~y, reg.line=lm, smooth=TRUE, spread=TRUE, boxplots=FALSE, span=0.5, data=datos4)

## Warning in plot.window(...): "reg.line" is not a graphical parameter

## Warning in plot.window(...): "spread" is not a graphical parameter

## Warning in plot.window(...): "span" is not a graphical parameter

## Warning in plot.xy(xy, type, ...): "reg.line" is not a graphical parameter

## Warning in plot.xy(xy, type, ...): "spread" is not a graphical parameter

## Warning in plot.xy(xy, type, ...): "span" is not a graphical parameter

## Warning in axis(side = side, at = at, labels = labels, ...): "reg.line" is not a
## graphical parameter

## Warning in axis(side = side, at = at, labels = labels, ...): "spread" is not a
## graphical parameter

## Warning in axis(side = side, at = at, labels = labels, ...): "span" is not a
## graphical parameter

## Warning in axis(side = side, at = at, labels = labels, ...): "reg.line" is not a
## graphical parameter

## Warning in axis(side = side, at = at, labels = labels, ...): "spread" is not a
## graphical parameter

## Warning in axis(side = side, at = at, labels = labels, ...): "span" is not a
## graphical parameter

## Warning in box(...): "reg.line" is not a graphical parameter

## Warning in box(...): "spread" is not a graphical parameter

## Warning in box(...): "span" is not a graphical parameter

## Warning in title(...): "reg.line" is not a graphical parameter

## Warning in title(...): "spread" is not a graphical parameter

## Warning in title(...): "span" is not a graphical parameter

mod1regr = lm(res ~ x + y + I(x * y) + I(x^2) + I(y^2), data = datos4)
summary(mod1regr)

## 
## Call:
## lm(formula = res ~ x + y + I(x * y) + I(x^2) + I(y^2), data = datos4)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.125469 -0.025059 -0.000915  0.024278  0.121931 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)
## (Intercept)  9.622e-01  3.682e+00   0.261    0.794
## x            7.379e-07  5.432e-06   0.136    0.892
## y           -2.609e-06  4.099e-06  -0.636    0.526
## I(x * y)    -5.563e-13  2.511e-12  -0.222    0.825
## I(x^2)      -7.313e-14  2.258e-12  -0.032    0.974
## I(y^2)       1.537e-12  1.962e-12   0.783    0.435
## 
## Residual standard error: 0.04042 on 127 degrees of freedom
## Multiple R-squared:  0.004807,   Adjusted R-squared:  -0.03437 
## F-statistic: 0.1227 on 5 and 127 DF,  p-value: 0.9871

mod2regr = stepAIC(mod1regr, scope = list(upper = mod1regr$formula,lower = ~1), direction = "both")

## Start:  AIC=-847.57
## res ~ x + y + I(x * y) + I(x^2) + I(y^2)
## 
##            Df  Sum of Sq     RSS     AIC
## - I(x^2)    1 0.00000171 0.20751 -849.57
## - x         1 0.00003015 0.20754 -849.55
## - I(x * y)  1 0.00008018 0.20759 -849.52
## - y         1 0.00066186 0.20817 -849.15
## - I(y^2)    1 0.00100229 0.20851 -848.93
## <none>                   0.20751 -847.57
## 
## Step:  AIC=-849.57
## res ~ x + y + I(x * y) + I(y^2)
## 
##            Df  Sum of Sq     RSS     AIC
## - I(x * y)  1 0.00007863 0.20759 -851.52
## - x         1 0.00008391 0.20759 -851.52
## - y         1 0.00066305 0.20817 -851.15
## - I(y^2)    1 0.00100058 0.20851 -850.93
## <none>                   0.20751 -849.57
## 
## Step:  AIC=-851.52
## res ~ x + y + I(y^2)
## 
##          Df  Sum of Sq     RSS     AIC
## - x       1 0.00006306 0.20765 -853.48
## - y       1 0.00092143 0.20851 -852.93
## - I(y^2)  1 0.00092195 0.20851 -852.93
## <none>                 0.20759 -851.52
## 
## Step:  AIC=-853.48
## res ~ y + I(y^2)
## 
##          Df  Sum of Sq     RSS     AIC
## - y       1 0.00085840 0.20851 -854.93
## - I(y^2)  1 0.00085889 0.20851 -854.93
## <none>                 0.20765 -853.48
## 
## Step:  AIC=-854.93
## res ~ I(y^2)
## 
##          Df  Sum of Sq     RSS     AIC
## - I(y^2)  1 4.9063e-07 0.20851 -856.93
## <none>                 0.20851 -854.93
## 
## Step:  AIC=-856.93
## res ~ 1

summary(mod2regr)

## 
## Call:
## lm(formula = res ~ 1, data = datos4)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.126152 -0.025600 -0.000092  0.027088  0.124705 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.129e-17  3.446e-03       0        1
## 
## Residual standard error: 0.03974 on 132 degrees of freedom

—————————————————————————————————————————————————————————————————

Analisis estructural

geo = as.geodata(datos4, coords.col = 1:2, data.col = 3)
plot(geo, scatter3d = TRUE)

## plot.geodata: the argument scatter3d=TRUE requires the package scatterplot3d
##  which is not available, will plot an histogram instead

var1 = variog(geo)

## variog: computing omnidirectional variogram

plot(var1)

var2 = variog(geo, max.dist = 125000)

## variog: computing omnidirectional variogram

plot(var2)

plot(variog4(geo, max.dist = 125000))

## variog: computing variogram for direction = 0 degrees (0 radians)
##         tolerance angle = 22.5 degrees (0.393 radians)
## variog: computing variogram for direction = 45 degrees (0.785 radians)
##         tolerance angle = 22.5 degrees (0.393 radians)
## variog: computing variogram for direction = 90 degrees (1.571 radians)
##         tolerance angle = 22.5 degrees (0.393 radians)
## variog: computing variogram for direction = 135 degrees (2.356 radians)
##         tolerance angle = 22.5 degrees (0.393 radians)
## variog: computing omnidirectional variogram

—————————————————————————————————————————————————————————————————

Ajuste modelo teórico de semivarianza

#ev=eyefit(var2)
#ev

mod1exp=variofit(var2,ini=c(9e-04,22737.3),nugget=8e-04,cov.model="exponential",weights="equal")

## variofit: covariance model used is exponential 
## variofit: weights used: equal 
## variofit: minimisation function used: optim

mod1exp

## variofit: model parameters estimated by OLS (ordinary least squares):
## covariance model is: exponential
## parameter estimates:
##      tausq    sigmasq        phi 
##     0.0008     0.0009 22737.3000 
## Practical Range with cor=0.05 for asymptotic range: 68114.86
## 
## variofit: minimised sum of squares = 0

mod2exp=variofit(var2,ini=c(9e-04, 22737.3),nugget=8e-04,cov.model="exponential",weights="npairs")

## variofit: covariance model used is exponential 
## variofit: weights used: npairs 
## variofit: minimisation function used: optim

mod2exp

## variofit: model parameters estimated by WLS (weighted least squares):
## covariance model is: exponential
## parameter estimates:
##      tausq    sigmasq        phi 
##     0.0008     0.0009 22737.3000 
## Practical Range with cor=0.05 for asymptotic range: 68114.86
## 
## variofit: minimised weighted sum of squares = 1e-04

mod3exp=variofit(var2,ini=c(9e-04,22737.3),nugget=8e-04,cov.model="exponential",weights="cressie")

## variofit: covariance model used is exponential 
## variofit: weights used: cressie 
## variofit: minimisation function used: optim

mod3exp

## variofit: model parameters estimated by WLS (weighted least squares):
## covariance model is: exponential
## parameter estimates:
##      tausq    sigmasq        phi 
##     0.0008     0.0009 22737.3000 
## Practical Range with cor=0.05 for asymptotic range: 68114.86
## 
## variofit: minimised weighted sum of squares = 35.757

mod4exp=likfit(geo,ini=c(9e-04,22737.3),nugget=8e-04,cov.model="exponential",lik.method="ML")

## kappa not used for the exponential correlation function
## ---------------------------------------------------------------
## likfit: likelihood maximisation using the function optim.
## likfit: Use control() to pass additional
##          arguments for the maximisation function.
##         For further details see documentation for optim.
## likfit: It is highly advisable to run this function several
##         times with different initial values for the parameters.
## likfit: WARNING: This step can be time demanding!
## ---------------------------------------------------------------
## likfit: end of numerical maximisation.

mod4exp

## likfit: estimated model parameters:
##        beta       tausq     sigmasq         phi 
## "1.000e-04" "5.000e-04" "1.600e-03" "2.274e+04" 
## Practical Range with cor=0.05 for asymptotic range: 68114.86
## 
## likfit: maximised log-likelihood = 252.1

mod5exp=likfit(geo,ini=c(9e-04,22737.3),nugget=8e-04,cov.model="exponential",lik.method="REML")

## kappa not used for the exponential correlation function
## ---------------------------------------------------------------
## likfit: likelihood maximisation using the function optim.
## likfit: Use control() to pass additional
##          arguments for the maximisation function.
##         For further details see documentation for optim.
## likfit: It is highly advisable to run this function several
##         times with different initial values for the parameters.
## likfit: WARNING: This step can be time demanding!
## ---------------------------------------------------------------
## likfit: end of numerical maximisation.

mod5exp

## likfit: estimated model parameters:
##        beta       tausq     sigmasq         phi 
## "1.000e-04" "5.000e-04" "1.700e-03" "2.274e+04" 
## Practical Range with cor=0.05 for asymptotic range: 68114.86
## 
## likfit: maximised log-likelihood = 251

plot(var2,main=expression(paste("Exponencial estimated ",tau^2)), ylim = c(0, 0.0021))
lines(mod1exp, max.dist = 125000, col = 1)
lines(mod2exp, max.dist = 125000, col = 2)
lines(mod3exp, max.dist = 125000, col = 3)
lines(mod4exp, max.dist = 125000, col = 4)
lines(mod5exp, max.dist = 125000, col = 5)
legend(60000, 5e-04, legend = c("OLS", "WLS - npairs", "WLS - cressie", "ML", "REML"),col = 1:5, lty = c(1, 1, 2, 2), lwd = c(1, 2, 1, 2), cex = 0.7)

mod1sph=variofit(var2,ini=c(9e-04,22737.3),nugget=8e-04,cov.model="spherical",weights="equal")

## variofit: covariance model used is spherical 
## variofit: weights used: equal 
## variofit: minimisation function used: optim

mod1sph

## variofit: model parameters estimated by OLS (ordinary least squares):
## covariance model is: spherical
## parameter estimates:
##      tausq    sigmasq        phi 
##     0.0008     0.0008 22737.3000 
## Practical Range with cor=0.05 for asymptotic range: 22737.3
## 
## variofit: minimised sum of squares = 0

mod2sph=variofit(var2,ini=c(9e-04,22737.3),nugget=8e-04,cov.model="spherical",weights="npairs")

## variofit: covariance model used is spherical 
## variofit: weights used: npairs 
## variofit: minimisation function used: optim

mod2sph

## variofit: model parameters estimated by WLS (weighted least squares):
## covariance model is: spherical
## parameter estimates:
##      tausq    sigmasq        phi 
##     0.0008     0.0009 22737.3000 
## Practical Range with cor=0.05 for asymptotic range: 22737.3
## 
## variofit: minimised weighted sum of squares = 2e-04

mod3sph=variofit(var2,ini=c(9e-04,22737.3),nugget=8e-04,cov.model="spherical",weights="cressie")

## variofit: covariance model used is spherical 
## variofit: weights used: cressie 
## variofit: minimisation function used: optim

mod3sph

## variofit: model parameters estimated by WLS (weighted least squares):
## covariance model is: spherical
## parameter estimates:
##      tausq    sigmasq        phi 
##     0.0008     0.0009 22737.3000 
## Practical Range with cor=0.05 for asymptotic range: 22737.3
## 
## variofit: minimised weighted sum of squares = 84.6133

mod4sph=likfit(geo,ini=c(9e-04,22737.3),nugget=8e-04,cov.model="spherical",lik.method="ML")

## kappa not used for the spherical correlation function
## ---------------------------------------------------------------
## likfit: likelihood maximisation using the function optim.
## likfit: Use control() to pass additional
##          arguments for the maximisation function.
##         For further details see documentation for optim.
## likfit: It is highly advisable to run this function several
##         times with different initial values for the parameters.
## likfit: WARNING: This step can be time demanding!
## ---------------------------------------------------------------
## likfit: end of numerical maximisation.

mod4sph

## likfit: estimated model parameters:
##        beta       tausq     sigmasq         phi 
## "1.800e-03" "4.000e-04" "1.300e-03" "2.274e+04" 
## Practical Range with cor=0.05 for asymptotic range: 22737.3
## 
## likfit: maximised log-likelihood = 251

mod5sph=likfit(geo,ini=c(9e-04,22737.3),nugget=8e-04,cov.model="spherical",lik.method="REML")

## kappa not used for the spherical correlation function
## ---------------------------------------------------------------
## likfit: likelihood maximisation using the function optim.
## likfit: Use control() to pass additional
##          arguments for the maximisation function.
##         For further details see documentation for optim.
## likfit: It is highly advisable to run this function several
##         times with different initial values for the parameters.
## likfit: WARNING: This step can be time demanding!
## ---------------------------------------------------------------
## likfit: end of numerical maximisation.

mod5sph

## likfit: estimated model parameters:
##        beta       tausq     sigmasq         phi 
## "1.900e-03" "4.000e-04" "1.300e-03" "2.274e+04" 
## Practical Range with cor=0.05 for asymptotic range: 22737.3
## 
## likfit: maximised log-likelihood = 249.2

plot(var2, main = expression(paste("Esferico estimated ",tau^2)), ylim = c(0, 0.0021))
lines(mod1sph, max.dist = 125000, col = 1)
lines(mod2sph, max.dist = 125000, col = 2)
lines(mod3sph, max.dist = 125000, col = 3)
lines(mod4sph, max.dist = 125000, col = 4)
lines(mod5sph, max.dist = 125000, col = 5)
legend(60000, 5e-04, legend = c("OLS", "WLS - npairs", "WLS - cressie", "ML", "REML"),col = 1:5, lty = c(1, 1, 2, 2), lwd = c(1, 2, 1, 2), cex = 0.7)

mod1gau=variofit(var2,ini=c(9e-04,22737.3),nugget=8e-04,cov.model="gaussian",weights="equal")

## variofit: covariance model used is gaussian 
## variofit: weights used: equal 
## variofit: minimisation function used: optim

mod1gau

## variofit: model parameters estimated by OLS (ordinary least squares):
## covariance model is: gaussian
## parameter estimates:
##      tausq    sigmasq        phi 
##     0.0008     0.0009 22737.3000 
## Practical Range with cor=0.05 for asymptotic range: 39354.14
## 
## variofit: minimised sum of squares = 0

mod2gau=variofit(var2,ini=c(9e-04,22737.3),nugget=8e-04,cov.model="gaussian",weights="npairs")

## variofit: covariance model used is gaussian 
## variofit: weights used: npairs 
## variofit: minimisation function used: optim

mod2gau

## variofit: model parameters estimated by WLS (weighted least squares):
## covariance model is: gaussian
## parameter estimates:
##      tausq    sigmasq        phi 
##     0.0008     0.0009 22737.3000 
## Practical Range with cor=0.05 for asymptotic range: 39354.14
## 
## variofit: minimised weighted sum of squares = 0

mod3gau=variofit(var2,ini=c(9e-04,22737.3),nugget=8e-04,cov.model="gaussian",weights="cressie")

## variofit: covariance model used is gaussian 
## variofit: weights used: cressie 
## variofit: minimisation function used: optim

mod3gau

## variofit: model parameters estimated by WLS (weighted least squares):
## covariance model is: gaussian
## parameter estimates:
##      tausq    sigmasq        phi 
##     0.0008     0.0009 22737.3000 
## Practical Range with cor=0.05 for asymptotic range: 39354.14
## 
## variofit: minimised weighted sum of squares = 16.922

mod4gau=likfit(geo,ini=c(9e-04,22737.3),nugget=8e-04,cov.model="gaussian",lik.method="ML")

## kappa not used for the gaussian correlation function
## ---------------------------------------------------------------
## likfit: likelihood maximisation using the function optim.
## likfit: Use control() to pass additional
##          arguments for the maximisation function.
##         For further details see documentation for optim.
## likfit: It is highly advisable to run this function several
##         times with different initial values for the parameters.
## likfit: WARNING: This step can be time demanding!
## ---------------------------------------------------------------
## likfit: end of numerical maximisation.

mod4gau

## likfit: estimated model parameters:
##        beta       tausq     sigmasq         phi 
## "1.800e-03" "9.000e-04" "1.000e-03" "2.274e+04" 
## Practical Range with cor=0.05 for asymptotic range: 39354.14
## 
## likfit: maximised log-likelihood = 252.2

mod5gau=likfit(geo,ini=c(9e-04,22737.3),nugget=8e-04,cov.model="gaussian",lik.method="REML")

## kappa not used for the gaussian correlation function
## ---------------------------------------------------------------
## likfit: likelihood maximisation using the function optim.
## likfit: Use control() to pass additional
##          arguments for the maximisation function.
##         For further details see documentation for optim.
## likfit: It is highly advisable to run this function several
##         times with different initial values for the parameters.
## likfit: WARNING: This step can be time demanding!
## ---------------------------------------------------------------
## likfit: end of numerical maximisation.

mod5gau

## likfit: estimated model parameters:
##        beta       tausq     sigmasq         phi 
## "1.800e-03" "9.000e-04" "1.000e-03" "2.274e+04" 
## Practical Range with cor=0.05 for asymptotic range: 39354.14
## 
## likfit: maximised log-likelihood = 250.8

plot(var2, main = expression(paste("Gaussiano estimated ", tau^2)), ylim = c(0, 0.0021))
lines(mod1gau, max.dist = 125000, col = 1)
lines(mod2gau, max.dist = 125000, col = 2)
lines(mod3gau, max.dist = 125000, col = 3)
lines(mod4gau, max.dist = 125000, col = 4)
lines(mod5gau, max.dist = 125000, col = 5)
legend(60000, 5e-04, legend = c("OLS", "WLS - npairs", "WLS - cressie", "ML", "REML"),col = 1:5, lty = c(1, 1, 2, 2), lwd = c(1, 2, 1, 2), cex = 0.7)

—————————————————————————————————————————————————————————————————

Validacion cruzada

ve.fit1exp <- as.vgm.variomodel(mod1exp)
ve.fit2exp <- as.vgm.variomodel(mod2exp)
ve.fit3exp <- as.vgm.variomodel(mod3exp)
ve.fit4exp <- as.vgm.variomodel(mod4exp)
ve.fit5exp <- as.vgm.variomodel(mod5exp)
ve.fit1sph <- as.vgm.variomodel(mod1sph)
ve.fit2sph <- as.vgm.variomodel(mod2sph)
ve.fit3sph <- as.vgm.variomodel(mod3sph)
ve.fit4sph <- as.vgm.variomodel(mod4sph)
ve.fit5sph <- as.vgm.variomodel(mod5sph)
ve.fit1gau <- as.vgm.variomodel(mod1gau)
ve.fit2gau <- as.vgm.variomodel(mod2gau)
ve.fit3gau <- as.vgm.variomodel(mod3gau)
ve.fit4gau <- as.vgm.variomodel(mod4gau)
ve.fit5gau <- as.vgm.variomodel(mod5gau)
cvolsexp <- krige.cv(res ~ 1, datos5, ve.fit1exp, maxdist = 125000)
cvwlsnpairsexp <- krige.cv(res ~ 1, datos5, ve.fit2exp, maxdist = 125000)
cvwlsncressieexp <- krige.cv(res ~ 1, datos5, ve.fit3exp, maxdist = 125000)
cvmlexp <- krige.cv(res ~ 1, datos5, ve.fit4exp, maxdist = 125000)

cvremlexp <- krige.cv(res ~ 1, datos5, ve.fit5exp, maxdist = 125000)
cvolsgau <- krige.cv(res ~ 1, datos5, ve.fit1gau, maxdist = 125000)
cvwlsnpairsgau <- krige.cv(res ~ 1, datos5, ve.fit2gau, maxdist = 125000)
cvwlsncressiegau <- krige.cv(res ~ 1, datos5, ve.fit3gau, maxdist = 125000)
cvmlgau <- krige.cv(res ~ 1, datos5, ve.fit4gau, maxdist = 125000)
cvremlgau <- krige.cv(res ~ 1, datos5, ve.fit5gau, maxdist = 125000)
cvolssph <- krige.cv(res ~ 1, datos5, ve.fit1sph, maxdist = 125000)
cvwlsnpairssph <- krige.cv(res ~ 1, datos5, ve.fit2sph, maxdist = 125000)
cvwlsncressiesph <- krige.cv(res ~ 1, datos5, ve.fit3sph, maxdist = 125000)
cvmlsph <- krige.cv(res ~ 1, datos5, ve.fit4sph, maxdist = 125000)
cvremlsph <- krige.cv(res ~ 1, datos5, ve.fit5sph, maxdist = 125000)
sqrt((round(mean(cvolsexp$residual), 5) - 0)^2 + (round(mean(cvolsexp$zscore),5) - 0)^2 + (round(var(cvolsexp$zscore), 5) - 1)^2)

## [1] 0.00616786

sqrt((round(mean(cvwlsnpairsexp$residual), 5) - 0)^2 + (round(mean(cvwlsnpairsexp$zscore),5) - 0)^2 + (round(var(cvwlsnpairsexp$zscore), 5) - 1)^2)

## [1] 0.00616786

sqrt((round(mean(cvwlsncressieexp$residual), 5) - 0)^2 + (round(mean(cvwlsncressieexp$zscore),5) - 0)^2 + (round(var(cvwlsncressieexp$zscore), 5) - 1)^2)

## [1] 0.03219662

sqrt((round(mean(cvmlexp$residual), 5) - 0)^2 + (round(mean(cvmlexp$zscore),5) - 0)^2 + (round(var(cvmlexp$zscore), 5) - 1)^2)

## [1] 0.05267475

sqrt((round(mean(cvremlexp$residual), 5) - 0)^2 + (round(mean(cvremlexp$zscore),5) - 0)^2 + (round(var(cvremlexp$zscore), 5) - 1)^2)

## [1] 0.05235032

sqrt((round(mean(cvolssph$residual), 5) - 0)^2 + (round(mean(cvolssph$zscore),5) - 0)^2 + (round(var(cvolssph$zscore), 5) - 1)^2)

## [1] 0.03543913

sqrt((round(mean(cvwlsnpairssph$residual), 5) - 0)^2 + (round(mean(cvwlsnpairssph$zscore),5) - 0)^2 + (round(var(cvwlsnpairssph$zscore), 5) - 1)^2)

## [1] 0.04127076

sqrt((round(mean(cvwlsncressiesph$residual), 5) - 0)^2 + (round(mean(cvwlsncressiesph$zscore),5) - 0)^2 + (round(var(cvwlsncressiesph$zscore), 5) - 1)^2)

## [1] 0.09148157

sqrt((round(mean(cvmlsph$residual), 5) - 0)^2 + (round(mean(cvmlsph$zscore),5) - 0)^2 + (round(var(cvmlsph$zscore), 5) - 1)^2)

## [1] 0.07420701

sqrt((round(mean(cvremlsph$residual), 5) - 0)^2 + (round(mean(cvremlsph$zscore),5) - 0)^2 + (round(var(cvremlsph$zscore), 5) - 1)^2)

## [1] 0.07180047

sqrt((round(mean(cvolsgau$residual), 5) - 0)^2 + (round(mean(cvolsgau$zscore), 5) - 0)^2 + (round(var(cvolsgau$zscore), 5) - 1)^2)

## [1] 0.1549355

sqrt((round(mean(cvwlsnpairsgau$residual), 5) - 0)^2 + (round(mean(cvwlsnpairsgau$zscore),5) - 0)^2 + (round(var(cvwlsnpairsgau$zscore), 5) - 1)^2)

## [1] 0.1554954

sqrt((round(mean(cvwlsncressiegau$residual), 5) - 0)^2 + (round(mean(cvwlsncressiegau$zscore),5) - 0)^2 + (round(var(cvwlsncressiegau$zscore), 5) - 1)^2)

## [1] 0.1436571

sqrt((round(mean(cvmlgau$residual), 5) - 0)^2 + (round(mean(cvmlgau$zscore), 5) - 0)^2 + (round(var(cvmlgau$zscore), 5) - 1)^2)

## [1] 0.03725241

sqrt((round(mean(cvremlgau$residual), 5) - 0)^2 + (round(mean(cvremlgau$zscore), 5) - 0)^2 + (round(var(cvremlgau$zscore), 5) - 1)^2)

## [1] 0.0357216

—————————————————————————————————————————————————————————————————

Prediccion espacial Kriging

muestra = spsample(poligonos, n = 10000, "regular")
muestra1 = as.data.frame(muestra)
names(muestra1) = c("x", "y")
gridded(muestra1) = c("x", "y")
plot(muestra1)

tprec.ok=krige(tprec ~ x+y+I(x*y)+I(x^2),datos3,muestra1, ve.fit1exp)

## [using universal kriging]

tprec.ok$var1.pred=ifelse(tprec.ok$var1.pred>=max(datos2$tprec),max(datos2$tprec),tprec.ok$var1.pred)
prec.ok = tprec.ok
prec.ok$var1.pred=(prec.ok$var1.pred * lambda + 1)^(1/lambda)
prec.ok$var1.var=prec.ok$var1.var/(((prec.ok$var1.pred)^(lambda-1))^2)
summary(prec.ok)

## Object of class SpatialPixelsDataFrame
## Coordinates:
##        min     max
## x 910266.2 1113772
## y 904553.1 1136022
## Is projected: NA 
## proj4string : [NA]
## Number of points: 9994
## Grid attributes:
##   cellcentre.offset cellsize cells.dim
## x          911042.9 1553.482       131
## y          905329.8 1553.482       149
## Data attributes:
##    var1.pred         var1.var       
##  Min.   : 621.7   Min.   :   28818  
##  1st Qu.:1069.8   1st Qu.:  137531  
##  Median :1394.3   Median :  305832  
##  Mean   :1582.9   Mean   :  980569  
##  3rd Qu.:1715.8   3rd Qu.:  534289  
##  Max.   :3986.0   Max.   :14511200

li = list("sp.polygons", poligonos)
pts = list("sp.points", datos3, pch = 3, col = "black", cex = 0.2)
spplot(prec.ok, c("var1.pred"), as.table = TRUE, main = " ", sp.layout = list(li, pts), contour = FALSE,labels = FALSE, pretty = TRUE, col = "black", col.regions = terrain.colors(100))

spplot(prec.ok, c("var1.var"), as.table = TRUE, main = " ", sp.layout = list(li, pts), contour = FALSE,labels = FALSE, pretty = TRUE, col = "black", col.regions = terrain.colors(100))

\(\text{Se presenta el mapa de prediccion para la precipitacion en el Departamento de Cundinamarca.}\ \) \(\text{De igual manera se presenta el mapa de la varianza de las predicciones.}\ \)

—————————————————————————————————————————————————————————————————