Geoestadistica

“Se juzga a una sociedad por sus ruidos, su arte, sus fiestas, no por sus estadisticas” - Atalli

0.1 Geoestadistica

Desarrollada originalmente para la industria minera, la geoestadística cubre el análisis de datos de medición basados en la ubicación. Permite la interpolación basada en modelos de mediciones con estimación de la incertidumbre.

Aunque el termino puede resultar confuso (Estadísticas en la tierra), pero todos los ejemplos que hemos visto hasta ahora son en la tierra. En este contexto, Geoestadística hace alusión a un tipo de datos especificos y en consecuencia de análisis. Desarrollada originalmente por empresas mineras para obtener información de yacimientos a partir me mediciones tomadas en unas pocas localizaciones. Las preguntas que resuelve la estadística espacial varian, por ejemplo es estimar el total de material vegetal en un campo, la polución en una ciudad. Como todo modelo estadístico, lo hace cuantificando la incertidembre.

A diferencia del estudio de patrones puntuales, la ubicación de los datos no es estadisticamente enteresante,las localizaciones han sido elegidas para recolectarkja los datos. Hay toda una rama de la estadistica para el muestreo de datos en el espacio.

0.1.1 Naturaleva de la variable respuesta

Continua (Particulas PM10 en el aire)
Fracción o porcentaje (Proporcion de enfermos)
Conteo (numero insectos de broca en un cultivo de cafe)
Binaria (presencia/ausenica insectos de broca en un cultivo de cafe)
Categórica (tipos de rocas encontradas en una muestra.)

Entender la naturaleza de la variable es crucial, pues la mayoría de operaciones en geoestadistidca van a usar esa variable como respuesta de un modelo. Del mismo modo que no debería usar datos de conteos en un modelo lineal, se debe usar un modelo lineal generalizado Poisson. Tampoco debería usar un modelo geoestadistico lineal si su variable respuesta es un conteo.

Como con cualquier tipo de datos, el primer paso es la exploración; las graficas permiten encontrar rápidamente patrones. Si los datos se ven influenciados por fronteras, direcciones, éstas deben estar incorporadas en el modelo.

0.1.1.1 Tendencias de gran escala

0.1.1.2 Anisotropia

varían según la dirección en que son examinadas

Generalmente cuando el variograma es calculado en distintas direcciones presenta distintos comportamientos con la variación de la distancia.

0.1.1.3 Discontinuidades

0.1.2 Canadian geochemical survey data

Su trabajo es estudiar la acidez (pH) de algunos datos de encuestas canadienses

library(sp)
library(raster)
library(gstat)
ca_geo<-readRDS("D:/DriveW7/2019/Sabana/Estadistica aplicada Salud Publica/Bibliografía/DataCamp Spatial Statistics in R/datasets/ca_geo.rds", refhook = NULL)



# ca_geo has been pre-defined
str(ca_geo, 1)

## Formal class 'SpatialPointsDataFrame' [package "sp"] with 5 slots

# See what measurements are at each location
names(ca_geo)

##  [1] "ID"   "Elev" "pH"   "Zn"   "Cu"   "Pb"   "Ni"   "Co"   "Ag"   "Mn"  
## [11] "Fe"   "Mo"   "U"    "W"    "Sn"   "Hg"   "As"   "Sb"   "Ba"   "Cd"  
## [21] "V"    "Bi"   "Cr"   "LoI"  "F"    "Au"

# Get a summary of the acidity (pH) values
summary(ca_geo$pH)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
##   3.900   6.100   6.600   6.531   7.000   8.700      33

# Look at the distribution
hist(ca_geo$pH)

# Make a vector that is TRUE for the missing data
miss <- is.na(ca_geo$pH)
table(miss)

## miss
## FALSE  TRUE 
##  1107    33

# Plot a map of acidity
spplot(ca_geo[!miss, ], "pH")

Es un poco difícil ver las tendencias espaciales a partir de estos datos. Exploraremos otras maneras de hacer que los patrones sean más claros.

0.1.3 Ajuste de una superficie de tendencia

Los datos de acidez muestran que el pH aumenta ampliamente de noreste a suroeste. Ajustar un modelo lineal con las coordenadas como covariables interpolará un plano plano a través de los valores.

# ca_geo has been pre-defined
str(ca_geo, 1)

## Formal class 'SpatialPointsDataFrame' [package "sp"] with 5 slots

# Are they called lat-long, up-down, or what?
coordnames(ca_geo)

## [1] "x" "y"

# Complete the formula
m_trend <- lm(pH ~ x + y, as.data.frame(ca_geo))

# Check the coefficients
summary(m_trend)

## 
## Call:
## lm(formula = pH ~ x + y, data = as.data.frame(ca_geo))
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.83561 -0.32091 -0.00761  0.33188  2.06249 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  8.358e+01  3.002e+00   27.84   <2e-16 ***
## x           -5.691e-06  3.483e-07  -16.34   <2e-16 ***
## y           -1.313e-05  5.319e-07  -24.68   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.5299 on 1104 degrees of freedom
##   (33 observations deleted due to missingness)
## Multiple R-squared:  0.4237, Adjusted R-squared:  0.4227 
## F-statistic: 405.9 on 2 and 1104 DF,  p-value: < 2.2e-16

0.1.4 Predicción a partir de una superficie de tendencia

Nuestra próxima tarea es calcular el pH en las ubicaciones que tienen datos faltantes en la fuente. Para ello puede utilizar la función predict() en el modelo ajustado del ejercicio anterior.

# ca_geo, miss, m_trend have been pre-defined
ls.str()

## ca_geo : Formal class 'SpatialPointsDataFrame' [package "sp"] with 5 slots
## m_trend : List of 13
##  $ coefficients : Named num [1:3] 8.36e+01 -5.69e-06 -1.31e-05
##  $ residuals    : Named num [1:1107] 1.206 0.5765 -0.2202 -0.1243 0.0415 ...
##  $ effects      : Named num [1:1107] -217.2874 7.5432 -13.0775 -0.1542 0.0126 ...
##  $ rank         : int 3
##  $ fitted.values: Named num [1:1107] 6.79 6.92 6.92 6.92 6.96 ...
##  $ assign       : int [1:3] 0 1 2
##  $ qr           :List of 5
##  $ df.residual  : int 1104
##  $ na.action    : 'omit' Named int [1:33] 15 42 43 244 340 360 460 479 480 494 ...
##  $ xlevels      : Named list()
##  $ call         : language lm(formula = pH ~ x + y, data = as.data.frame(ca_geo))
##  $ terms        :Classes 'terms', 'formula'  language pH ~ x + y
##  $ model        :'data.frame':   1107 obs. of  3 variables:
## miss :  logi [1:1140] FALSE FALSE FALSE FALSE FALSE FALSE ...

# Make a vector that is TRUE for the missing data
miss <- is.na(ca_geo$pH)


# Create a data frame of missing data
ca_geo_miss <- as.data.frame(ca_geo)[miss, ]

# Predict pH for the missing data
predictions <- predict(m_trend, newdata = ca_geo_miss, se.fit = TRUE)

# Compute the exceedence probability
pAlkaline <- 1 - pnorm(7, mean = predictions$fit, sd = predictions$se.fit)
hist(pAlkaline)

0.2 Variagrama

La primera ley de la geografía, o principio de autocorrelación espacial, fue formulado por el geógrafo Waldo Tobler de esta forma:

Todas las cosas están relacionadas entre sí, pero las cosas más próximas en el espacio tienen una relación mayor que las distantes.

Pero como entenderemos estar relacionado (varianza) y proximidad (distancia)?

Que tan similares son dos muestras?

\[\Delta =(A-B)^2 \] Se calcula la distancia para cada par de puntos a una distancia dada.

El variograma o semivariograma (semi por que en realidad no es una varianza) es una herramienta que permite analizar el comportamiento espacial de una variable sobre un área definida. También permite conocer el alcance, es decir, la distancia máxima a la que una muestra tiene influencia sobre otra muestra, una aplicación de esta información es conocer la vecindad en que se pueden buscar muestras para estimar el valor de un punto especifico.

0.2.0.0.1 Nube Variograma

choose(200,2)

## [1] 19900

0.2.0.0.2 Construcción de Bins

0.2.0.0.3 Promedio en Bins

Usualmente un variograma es muy denso como para tener información útil, por lo cual es necesario segmentarlo (bins) y resumirlo (promedio).

0.2.0.0.4 Variograma

A veces llamado confusamente semivariograma, pues el eje y representa una semivarianza.

Como se interpreta? De la primera ley de la geografia sabemos que la relación decrece en función de la distancia, valores más cercanos espacialmente tendrán menor varianza y aquellos más alejados mayor.

La parte izquierda da información de las muestras que estan juntas una de otra.

El valor a la distancia cero es la diferencia estimada de dos muestras en la misma ubicación. A esta diferencia se le conoce como Nugget (pepita), pues si la primera muestra extrae una pepita de oro, la segunda no debería encontrarla, es decir varianza!

Puede ocurrir que para distancias cercanas a cero el valor del variograma no se aproxima a cero

Range: corresponde a la distancia a partir de la cual dos observaciones son independientes. El rango se interpreta como la zona de inﬂuencia.

Un paso fundamental para geoestadistica es estimar el variograma, para lo cual necesitamos tres parámetros, nugget, sill y range.

A continuación se presentan algunas formas usuales de variogramas.

show.vgms()

0.3 Estimación del variograma

Puede utilizar el paquete gstat para trazar nubes de variogramas y los variogramas a partir de los datos. Recuerde que:

La nube de variograma muestra las diferencias de las mediciones contra la distancia para todos los pares de puntos de datos.
El variograma divide la nube en ubicaciones de distancia y calcula la diferencia promedio dentro de cada ubicación.

# install.packages("gstat")
library(gstat)
# ca_geo, miss have been pre-defined


# Make a cloud from the non-missing data up to 10km
plot(variogram(pH ~ 1, ca_geo[!miss, ], cloud = TRUE, cutoff = 10000))

# Make a variogram of the non-missing data
plot(variogram(pH ~ 1, ca_geo[!miss, ]))

0.3.1 Variograma con tendencia espacial

Resulta natural suporque que si si el suelo en un punto particular es alcalino, entonces el suelo a un metro de distancia probablemente será alcalino también. Pero, ¿puede decir lo mismo de la tierra a un kilómetro, o a diez kilómetros, o a cien kilómetros?

La forma del variograma anterior indica que hay una tendencia a gran escala en los datos. Puede ajustar un variograma considerando esta tendencia con gstat. Este variograma debería aplanarse, indicando que no hay más correlación espacial después de cierta distancia.

ph_vgm <- variogram(pH ~ x + y, ca_geo[!miss, ])
plot(ph_vgm)

0.3.2 Ajuste del modelo del variograma

A continuación, ajustará un modelo a su variograma usando la función gstat fit.variogram(). Necesita darle algunos valores iniciales como punto de partida para que el algoritmo de optimización se ajuste a un mejor modelo.

# parametros iniciales estimados a ojo
nugget <- 0.14
psill <- 0.14
range <- 10000

# Fit the variogram
v_model <- fit.variogram(
  ph_vgm, 
  model = vgm(
    model = "Ste",
    nugget = nugget,
    psill = psill,
    range = range,
    kappa = 0.5
  )
)

# Show the fitted variogram on top of the binned variogram
plot(ph_vgm, model = v_model)

print(v_model)

##   model     psill    range kappa
## 1   Nug 0.1545119     0.00   0.0
## 2   Ste 0.1442006 14379.39   0.5

0.4 Predicciones por Krigging

0.4.1 Kriging ordinario

Usualmente kriging se utiliza para predecir los píxeles (o nodos) de una malla regular que cubre la zona de estudio.

la variable es modelada a partir de si misma;
la media espacial no es conocida a priori, sino estimada de los datos.

0.4.2 Llenando los huecos

La parte final de la estimación geoestática es el kriging mismo. Esta es la aplicación del variograma junto con los puntos de datos de la muestra para producir estimaciones e incertidumbre en nuevas ubicaciones.

El cálculo de estimaciones e incertidumbre, junto con la suposición de una respuesta normal (Gaussiana) significa que se puede calcular cualquier función de las estimaciones - por ejemplo, la probabilidad de que una nueva ubicación tenga suelo alcalino.

0.4.3 Creación de una tabla de predicción

vamos a elaborar un mapa de la probabilidad de que el area de estudio sea alcalina. Para elllo, usaremos kriging a través de la función krige(). Esto requiere un objeto SpatialPixels que necesitará un poco de manipulación de datos para su creación. Comienza por definir una cuadrícula, crear puntos en esa cuadrícula, recortar a la región de estudio, y finalmente convertir a SpatialPixels. En el camino, se incorporarán algunas nuevas funciones.

GridTopology() define una cuadrícula rectangular. Toma como entradas tres vectores de longitud dos. La primera especifica la posición de la esquina inferior izquierda de la cuadrícula. El segundo especifica la anchura y la altura de cada rectángulo de la cuadrícula, y el tercero especifica el número de rectángulos en cada dirección.

Para asegurar que la cuadrícula y el área de estudio tengan las mismas coordenadas, es necesario realizar algunas tareas. SpatialPoints() convierte los puntos en un sistema de referencia de coordenadas (CRS), o proyección (diferentes paquetes usan terminología diferente para el mismo concepto). El CRS se crea envolviendo el área de estudio en projection(), luego en CRS(). Para el propósito de este ejercicio, por ahora no nos preocuparemos por lo que estas funciones hacen exactamente, sólo que esta manipulación de datos es necesaria para alinear la cuadrícula y el área de estudio.

Ahora que tiene esa alineación, crop(), como el nombre sugiere, corta la cuadrícula al área de estudio.

Finalmente, SpatialPixels() convierte los puntos de cuadrícula recortados de la trama en el objeto sp equivalente.

# Plot the polygon and points
plot(geo_bounds); points(ca_geo)

# Find the corners of the boundary
bbox(geo_bounds)

# Define a 2.5km square grid over the polygon extent. The first parameter is
# the bottom left corner.
grid <- GridTopology(c(537853, 5536290), c(2500, 2500), c(72, 48))

# Create points with the same coordinate system as the boundary
gridpoints <- SpatialPoints(grid, proj4string = CRS(projection(geo_bounds)))
plot(gridpoints)

# Crop out the points outside the boundary
cropped_gridpoints <- crop(gridpoints, geo_bounds)
plot(cropped_gridpoints)

# Convert to SpatialPixels
spgrid <- SpatialPixels(cropped_gridpoints)
coordnames(spgrid) <- c("x", "y")
plot(spgrid)

0.4.4 Predicciones en cuadrícula

Construir la red es la parte difícil. Ahora puede calcular estimaciones por kriggin sobre la cuadrícula usando el modelo de variogramas de antes (v_model) y la cuadrícula de SpatialPixels.

# Do kriging predictions over the grid
ph_grid <- krige(pH ~ x + y, ca_geo[!miss, ], newdata = spgrid, model = v_model)

# Calc the probability of pH exceeding 7
ph_grid$pAlkaline <- 1 - pnorm(7, mean = ph_grid$var1.pred, sd = sqrt(ph_grid$var1.var))

# Map the probability of alkaline samples
spplot(ph_grid, zcol = "pAlkaline")

0.5 Auto Kriging

# Kriging with linear trend, predicting over the missing points
# install.packages("automap")
library(automap)

ph_auto <- autoKrige(
  pH ~ x + y, 
  input_data = ca_geo[!miss, ], 
  new_data = ca_geo[miss, ], 
  model = "Mat"
)

## [using universal kriging]

# Plot the variogram, predictions, and standard error
plot(ph_auto)

0.5.1 Autokriging en ubicaciones puntuales

La función autoKrige() en el paquete automap calcula los variogramas, ajusta los modelos, realiza la selección de modelos y realiza el kriging haciendo múltiples llamadas a las funciones gstat que usó anteriormente. Puede ser un gran ahorro de tiempo, no obstante, siempre debe comprobar los resultados cuidadosamente.

En este ejemplo obtendrá predicciones en las ubicaciones de datos faltantes.

autoKrige() puede probar varios tipos de modelos de variogramas. En el ejemplo, usted usará un modelo de variograma de Matern, el cual es comúnmente usado en análisis de suelo y silvicultura. Puede ver una lista completa de los modelos disponibles llamando a vgm() sin argumentos.

# Kriging with linear trend, predicting over the missing points

ph_auto <- autoKrige(
  pH ~ x + y, 
  input_data = ca_geo[!miss, ], 
  new_data = ca_geo[miss, ], 
  model = "Mat"
)

## [using universal kriging]

# Plot the variogram, predictions, and standard error
plot(ph_auto)

0.5.2 Auto-kriging sobre una red

También puede usar autoKrige() sobre la cuadrícula de spiders del ejercicio anterior. Esto reúne todos los conceptos que ha aprendido en el capítulo. Es decir, el kriging es ideal para predecir datos faltantes, graficar cosas en una cuadrícula es mucho más claro que graficar puntos individuales, y el kriging automático que es menos dispendioso que el manual.

# Auto-run the kriging
ph_auto_grid <- autoKrige(pH ~ x + y, input_data = ca_geo[!miss, ], new_data = spgrid)

# Remember predictions from manual kriging
plot(ph_grid)

# Plot predictions and variogram fit
plot(ph_auto_grid)

# Compare the variogram model to the earlier one
v_model
ph_auto_grid$var_model

0.6 F-Otros métodos de interpolación no-geoestadísticos

Son empiricos
Son geometricos
No consideran la presencia de anisotropiasw es decir direcciones en las cuales la variación es privilegiada.
No tienen en cuenta el error asociado a la estimación, solo la media.

0.6.1 Polígonos de Thiessen:

Se le asigna a cada punto del espacio el valor del dato más próximo.

Se supone el hecho de contar con un conjunto de establecimientos que se desean colocar sobre una cierta región geográfica de tal manera que las ubicaciones sean lo más rentables posible. Por tanto, se debe hallar una configuración que permita que el número de clientes atraídos sea el más factible. La suposición lógica indica que los clientes irían al establecimiento más cercano a su domicilio y no a aquellos que sean muy lejanos. Con base en esto, los polgonos de Thiessen otorgan la configuración deseada por los establecimientos.

Cambios abruptos de bordes.
Solo utiliza un punto para cada predicción.
Dificil de estimar en tres dimensiones.

thiessen = krige(log10(zinc) ~ 1, meuse, meuse.grid, nmax = 1)
## [inverse distance weighted interpolation]
pts.s <- list("sp.points", meuse, col="white",pch=20)
spplot(thiessen, "var1.pred", asp=1, col.regions=rev(heat.colors(50)),
       sp.layout = list(pts.s),main="Thiessen")

0.6.2 Inverso de la distancia:

Se refiere a algunos fenómenos físicos cuya intensidad es inversamente proporcional al cuadrado de la distancia al centro donde se originan (sonido y luz)

Se le asigna mayor peso a las muestras cercanas y menor peso a las muestras alejadas del bloque.

No hay una manera objetiva de seleccionar el peso (inverso o inveso cuadrado…).
No hay una manera objetiva de seleccionar el radio de interpolación. En todos los casos, la distribución no regular puede sobre-enfatizar áreas (o lo contrario).

0.7 Taller

Con el conjunto de datos “Meuse” de polución del suelo con muestras de concentración de metales pesados

library(sp)
library(gstat)
data(meuse)
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 ...

0.7.1 Explore los datos

hist(meuse$zinc, breaks = 16)

Distribución no simétrica, sesgada hacia la derecha

Al ser una distribución no simétrica, se aplica logaritmo para transformar los valores y obtener una distribución simetrica (normal). Esto, además, reduce los posibles outliers.

meuse$logZn <- log10(meuse$zinc)
hist(meuse$logZn, breaks = 16)

Luego para transformar a los valores originales se utiliza antilog. antilog (x)=10^{x head(10}(meuse$logZn))

0.7.2 visualización

coordinates(meuse) <- c("x", "y")
class(meuse)

## [1] "SpatialPointsDataFrame"
## attr(,"package")
## [1] "sp"

plot(meuse, asp = 1, pch = 1) #asp=1 las dos escalas son iguales
#cargar otro conjunto de datos llamado meuse.riv que contiene lineas de las márgenes del río.
data(meuse.riv)
lines(meuse.riv)

0.7.2.1 Otra visualización con símbolos proporcionales

plot(meuse, asp = 1, cex = 4 * meuse$zinc/max(meuse$zinc),pch = 1) 
#cex tamaño de circulos proporcionales al valor
lines(meuse.riv)

0.7.3 Variograma empirico

ve <- variogram(logZn ~ 1, meuse, cutoff = 1300,width = 90)
plot(ve, plot.numbers = T, asp=1)

0.7.4 Ajuste automatico Variograma

vt <- vgm(psill = 0.12, model = "Sph", range = 850,nugget = 0.01) ## ajuste visual a partir de los valores a ojo y las posibles formas del variograma, mostradas anteriormente

va <- fit.variogram(ve, vt)

plot(ve, pl = T, model = va)

0.7.5 Interpolación - Kriging

data(meuse.grid) #malla de 40m x 40m, disponible con el dataset meuse.
coordinates(meuse.grid) <- c("x", "y")
gridded(meuse.grid) <- T #indica que el conjunto de datos es un raster

ok <- krige(logZn ~ 1, locations = meuse, newdata = meuse.grid, model = va)

## [using ordinary kriging]

ok$pred <- 10^(ok$var1.pred)#volver a valores originales

par(mfrow=c(2,1))
pts.s <- list("sp.points", meuse, col="white",pch=1, cex=4*meuse$zinc/max(meuse$zinc))
print(spplot(ok, "var1.pred", asp=1, col.regions=rev(heat.colors(50)),
             main="Predicción OK, log-ppm Zn",sp.layout = list(pts.s)), 
             split=c(1,1,2,1), more=TRUE)
pts.s <- list("sp.points", meuse, col="black", pch=20)
print(spplot(ok, zcol="var1.var",col.regions=rev(gray(seq(0,1,.01))), asp=1,
             main="Varianza OK, log-ppm Zn^2",sp.layout = list(pts.s)), 
             split=c(2,1,2,1), more=FALSE)

0.7.6 estimación por otros metodos no-geoestadisticos

0.7.6.1 Thiessen

thiessen = krige(log10(zinc) ~ 1, meuse, meuse.grid, nmax = 1)

## [inverse distance weighted interpolation]

pts.s <- list("sp.points", meuse, col="white",pch=20)
spplot(thiessen, "var1.pred", asp=1, col.regions=rev(heat.colors(50)),
       sp.layout = list(pts.s),main="Thiessen")

0.7.6.2 Inverso de la distancia (idw)

idw = idw(log10(zinc) ~ 1, meuse, meuse.grid)

## [inverse distance weighted interpolation]

spplot(idw, "var1.pred", asp=1, col.regions=rev(heat.colors(50)),
       sp.layout = list(pts.s),main="IDW")

#Cambiando el número de vecinos
idw$vecino<- krige(log10(zinc) ~ 1, meuse, meuse.grid,nmax=6)

## [inverse distance weighted interpolation]

## [inverse distance weighted interpolation]
#Cambiando pesos
idw$idp05 = idw(log10(zinc) ~ 1, meuse, meuse.grid, idp = 0.5)$var1.pred

## [inverse distance weighted interpolation]

## [inverse distance weighted interpolation]
idw$idp5 = idw(log10(zinc) ~ 1, meuse, meuse.grid, idp = 5)$var1.pred

## [inverse distance weighted interpolation]

## [inverse distance weighted interpolation]
idw$idp10 = idw(log10(zinc) ~ 1, meuse, meuse.grid, idp = 10)$var1.pred

## [inverse distance weighted interpolation]

## [inverse distance weighted interpolation]
spplot(idw, c("var1.pred","idp05", "idp5", "idp10"),  col.regions=rev(heat.colors(50))
       ,main="IDW")

0.7.7 Validación

thiessen.cv.a <- krige.cv(log10(zinc) ~ 1, meuse,nmax = 1)
idw.cv.a <- krige.cv(log10(zinc) ~ 1, meuse)
ok.cv.a <- krige.cv(log10(zinc) ~ 1, locations = meuse, model = va)

thiessen.cv.a[1:5,]

## class       : SpatialPointsDataFrame 
## features    : 5 
## extent      : 181025, 181307, 333330, 333611  (xmin, xmax, ymin, ymax)
## coord. ref. : NA 
## variables   : 6
## names       :        var1.pred, var1.var,         observed,           residual, zscore, fold 
## min values  : 2.44870631990508,      Inf, 2.40993312333129, -0.396246850652592,    Inf,    1 
## max values  : 3.05728564441821,     -Inf, 3.05728564441821, 0.0478347486195205,   -Inf,    5

ok.cv.a[1:5,]

## class       : SpatialPointsDataFrame 
## features    : 5 
## extent      : 181025, 181307, 333330, 333611  (xmin, xmax, ymin, ymax)
## coord. ref. : NA 
## variables   : 6
## names       :        var1.pred,           var1.var,         observed,          residual,            zscore, fold 
## min values  : 2.42583454100418, 0.0328122692709205, 2.40993312333129, -0.21711690771744, -1.05142115439224,    1 
## max values  : 2.94104620305062, 0.0426416369705371, 3.05728564441821, 0.118833193371807, 0.654094974314019,    5

que diferencia nota entre ambas predicciones?

0.7.7.1 Correlación

par(mfrow=c(1,3))
print(plot(var1.pred~observed,thiessen.cv.a, main="Thiessen"), split=c(1,1,3,1), more=TRUE)

## NULL

print(plot(var1.pred~observed,idw.cv.a, main="IDW"), split=c(2,1,3,1), more=TRUE)

## NULL

print(plot(var1.pred~observed,ok.cv.a, main="OK"), split=c(3,1,2,1), more=FALSE)

## NULL

cor(thiessen.cv.a$var1.pred,thiessen.cv.a$observed)

## [1] 0.685839

calcular las demas correlaciones y comparar.

0.7.7.2 pronósticos

boxplot(thiessen.cv.a$var1.pred, idw.cv.a$var1.pred, ok.cv.a$var1.pred, main="Thiessen, IDW, OK")

0.7.8 Residuales

boxplot(thiessen.cv.a$residual, idw.cv.a$residual, ok.cv.a$residual, main="Thiessen, IDW, OK")

“An expert is one who knows more and more about less and less” - Nicholas Murray Butler