Predicción Espacial (Kriging)

Librería gstat:

Modelado, Predicción y Simulación Geoestadística Espacial y Espacio-Temporal
Modelado de variogramas; (co)kriging simple, ordinario y universal de punto o bloque; kriging espacio-temporal; (co)simulación secuencial gaussiana o de indicadores; funciones de utilidad de trazado de mapas de variogramas y variogramas; admite sf y estrellas.

library(sp) # Para el manejo de datos espaciales
library(gstat) 

Análisis y descripción de los datos

Datos: Meuse

conjunto de datos de tipo “DataFrame” proporcionado por el paquete sp que comprende cuatro metales pesados medidos en la capa superior del suelo en una llanura aluvial a lo largo del río Meuse.

data(meuse)
data(meuse.riv)  # Rivera del río 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 ...

coordinates(meuse) = c("x", "y") # Convierte a datos espaciales
# Grafico
plot(meuse, asp = 1, pch = 1, main="Río Meuse")
lines (meuse.riv) 

Exploración de la variable Zinc

En el siguiente diagrama de barras y el diagrama de cajas y bigotes, observamos que la distribución de la variable zinc no es simétrica, pues está sesgada a la derecha. En el resúmen, se observa que la media es mayor que la mediana y en el diagrama de cajas se calculan algunos valores atípicos.

hist(meuse$zinc) 

boxplot(meuse$zinc, horizontal=T)

summary(meuse$zinc)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  113.0   198.0   326.0   469.7   674.5  1839.0 

# Grafico
plot(meuse, asp = 1, cex = 4*meuse$zinc/max(meuse$zinc), 
     pch = 1,  main="Río Meuse - Zinc")
lines (meuse.riv) 

Calculando logaritmo

De acuerdo a lo anterior, al ser una distribución no simétrica, se aconseja aplicar logaritmo para transformar los valores y obtener una distribución simétrica (normal). Además, esto reduce los posibles outliers o valores atípicos.

meuse$logZn = log10(meuse$zinc)
summary(meuse$logZn)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  2.053   2.297   2.513   2.556   2.829   3.265 

boxplot(meuse$logZn, horizontal=T)

Análisis estructural

La primera etapa en el desarrollo de un análisis geoestadístico es la determinación de la dependencia espacial entre los datos medidos de una variable. Esta fase es también conocida como análisis estructural. Para llevarla a cabo, con base en la información muestral, se usan tres funciones: El semivariograma, el covariograma y el correlograma, aunque en la práctica se emplea el semivariograma y no las otras dos funciones.

Para interpretar el semivariograma experimental se parte del criterio de que a menor distancia entre los sitios mayor similitud o correlación espacial entre las observaciones. Por ello en presencia de autocorrelación se espera que para valores de h pequeños el semivariograma experimental tenga magnitudes menores a las que este toma cuando las distancias h se incrementan.

En nuestro caso, se concluye que hay autocorrelación espacial.

• variogram: Calcula el variograma muestral o residual a partir de los datos, o en caso de que se proporcione un modelo lineal, para los residuos, con opciones de variograma direccional, robusto y agrupado, y para intervalos de distancia irregulares.
  Sus parámetros logZn~1 indican ausencia de regresores, meuse los datos, cutoff=1300 la distancia de separación espacial hasta la cual los pares de puntos se incluyen en las estimaciones de semivarianza y width=90, el ancho de los intervalos de distancia.
  Este comando retorna un objeto con los siguientes atributos: np: el número de pares de puntos para esta estimación; dist: la distancia promedio de todos los pares de puntos considerados para esta estimación y gamma: la estimación real del variograma de la muestra. Además indica la dirección horizontal, vertical y el id del par combinado.

ve = variogram(logZn ~ 1, meuse, cutoff = 1300, width = 90)
ve

    np       dist      gamma dir.hor dir.ver   id
1   41   72.24836 0.02649954       0       0 var1
2  212  142.88031 0.03242411       0       0 var1
3  320  227.32202 0.04818895       0       0 var1
4  371  315.85549 0.06543093       0       0 var1
5  423  406.44801 0.08025949       0       0 var1
6  458  496.09401 0.09509850       0       0 var1
7  455  586.78634 0.10656591       0       0 var1
8  466  677.39566 0.10333481       0       0 var1
9  503  764.55712 0.11461332       0       0 var1
10 480  856.69422 0.12924402       0       0 var1
11 468  944.02864 0.12290106       0       0 var1
12 460 1033.62277 0.12820318       0       0 var1
13 422 1125.63214 0.13206510       0       0 var1
14 408 1212.62350 0.11591294       0       0 var1
15 173 1280.65364 0.11719960       0       0 var1

plot(ve, plot.numbers = T)

Todos estos modelos tienen tres parámetros:

• Pepita: Se denota por \(C_0\) y representa una discontinuidad puntual del semivariograma en el origen.
  

• Meseta: Es la cota superior del semivariograma. Se denota por \(C_0+C_1\).

• Rango: La separación o distancia entre pares de puntos en la cual ya no hay dependencia espacial.

Método para seleccionar el mejor ajuste del modelo teórico

1. Conocer el proceso espacial
2. Ajuste visual
3. Ajuste automático y comparar la bondad del ajuste

A continuación se muestra un gráfico (Trellis) de los modelos diferentes teóricos de variogramas mediante la función show.vgms()

show.vgms()

Ajuste visual (Modelo esférico)

En nuestro ejemplo, visualmente se observa:

• Rango: Aproximadamente 850 m.
  

• Pepita = Aproximadamente 0.01
  

• Meseta = 0.13

• vgm: Genera un modelo de variograma o agrega a un modelo existente ingresando los parámetros: rango, silla sin \(C_0\), pepita \(C_0\), y el modelo que para este caso, es esférico (Sph).

vt = vgm(psill=0.12, range=850, nugget=0.01,
         model="Sph")
vt

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

Ajuste automático

• fit.variogram: Ajusta rangos y/o umbrales de un modelo de variograma simple o anidado a un variograma de muestra. Cuyos parámetros son el variograma y el variograma ajustado.

va = fit.variogram(ve, vt)
va

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

Interpolación: Kriging ordinario

Lo primero es crear una cuadrícula dónde se predecirán la concentración de Zinc.

data(meuse.grid) # Cuadrícula
coordinates(meuse.grid) = c("x", "y")
plot(meuse.grid)

# Promueve una estructura no cuadriculada a una cuadrícula
gridded(meuse.grid) = T

• krige: Función para kriging simple, ordinario o universal (a veces llamado kriging de deriva externa), función para la interpolación ponderada de distancia inversa. (Para nuesto caso kriging ordinario)

OK = krige(logZn~1, locations=meuse, newdata=meuse.grid, model=va)

[using ordinary kriging]

names(OK)

[1] "var1.pred" "var1.var" 

# Retomando los valores originales del Zinc.
OK$pred = 10^(OK$var1.pred) 

# Graficando
pts.s = list("sp.points", meuse, col="white", 
             pch=1, cex = 4*meuse$zinc/max(meuse$zinc))
print(spplot(OK, "pred", asp=1, 
             col.regions=rev(bpy.colors(64)),
             main = "Predicciones modelo Esférico",
             sp.layout=list(pts.s)))

Ejercicio:

Al variograma experimental calculado ve, realizar el ajuste e interpolación con el modelo teórico Exponencial

# Ajuste visual:
vt1 = vgm(psill=0.12, range=850, nugget=0.01, model="Exp")
plot(ve, pl = T, model = vt1)

vt1

  model psill range
1   Nug  0.01     0
2   Exp  0.12   850

# Ajuste automatico:
va1 = fit.variogram(ve, vt1)
plot(ve, pl=T, model=va1)

va1

  model       psill    range
1   Nug 0.004989641   0.0000
2   Exp 0.154851224 635.6953

# Interpolacion:
OK1 = krige(logZn~1, locations=meuse, newdata=meuse.grid, model=va1)

[using ordinary kriging]

OK1$pred = 10^(OK1$var1.pred) 
print(spplot(OK1, "pred", asp=1, 
             col.regions=rev(bpy.colors(64)),
             main = "Predicciones modelo Exponencial",
             sp.layout=list(pts.s)))

Validación cruzada

Consiste en excluir la observación de uno de los n puntos muestrales y con los n-1 valores restantes y el modelo de semivariograma escogido, predecir vía kriging el valor de la variable en estudio en la ubicación del punto que se excluyó.

krige.cv Funciones de validación cruzada para (co) kriging, kriging de puntos simples, ordinarios o universales en un vecindario local. Retorna un DataFrame que contiene las coordenadas de los datos o las de la primera variable en el objeto, y columnas de predicción y varianza de predicción de puntos de datos validados cruzados, valores observados, residuales, zscore (residuo dividido por el error estándar de kriging).

# Modelo esferico
valida = krige.cv(log(zinc)~1, meuse, va, nmax=40, nfold=5) # Esferico

# Modelo exponencial
valida1 = krige.cv(log(zinc)~1, meuse, va1, nmax=40, nfold=5)

bubble(valida, "residual", 
       main = "Log(zinc): 5-fold cv residuales \nModelo esférico    ")

Errores de predicción

Se piensa que si el modelo de semivarianza elegido describe bien la estructura de autocorrelación espacial, entonces la diferencia entre el valor observado y el valor predicho debe ser pequeña. Este procedimiento se realiza en forma secuencial con cada uno de los puntos muestrales y así se obtiene un conjunto de n “errores de predicción”.

Inicialmente podemos usar las siguientes medidas: error medio (ME), error absoluto medio (MAE), error cuadrático medio (RMSE) y relación de desviación cuadrática media (MSDR) de los residuos con varianza kriging.

Errores de predicción para el modelo esférico

# Mean Error (ME)
ME = round(mean(valida$residual),3)
ME

[1] -0.006

# Mean Absolute Error
MAE = round(mean(abs(valida$residual)),3)
MAE

[1] 0.309

# Root Mean Squre Error (RMSE)
RMSE = round(sqrt(mean(valida$residual^2)),3)
RMSE

[1] 0.408

# Mean Squared Deviation Ratio (MSDR)
MSDR = mean(valida$residual^2/valida$var1.var)
MSDR

[1] 4.513082

Errores de predicción para el modelo exponencial

# Mean Error (ME)
ME1 = round(mean(valida1$residual),3)
ME1

[1] 0.007

# 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.444

# Mean Squared Deviation Ratio (MSDR)
MSDR1 = mean(valida1$residual^2/valida1$var1.var)
MSDR1

[1] 4.642942

Se puede decir que aunque los errores son muy similares, el modelo esférico tiene un mejor comportamiento ya que sus errores son más pequeños por milésimas.

lm.cv Esta función proporciona medidas de validación interna y cruzada de precisión predictiva para la regresión lineal ordinaria

lm.cv = lm(valida$var1.pred ~ valida$observed)
summary(lm.cv)

Call:
lm(formula = valida$var1.pred ~ valida$observed)

Residuals:
     Min       1Q   Median       3Q      Max 
-1.09336 -0.22331 -0.01363  0.23209  0.72965 

Coefficients:
                Estimate Std. Error t value Pr(>|t|)    
(Intercept)      2.11730    0.22426   9.441   <2e-16 ***
valida$observed  0.63943    0.03782  16.907   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.3388 on 153 degrees of freedom
Multiple R-squared:  0.6514,    Adjusted R-squared:  0.6491 
F-statistic: 285.8 on 1 and 153 DF,  p-value: < 2.2e-16

plot(valida$var1.pred ~ valida$observed)
par(new=TRUE)

plot(1.96640+0.66687*valida$observed~valida$observed,
     type="l", ylab="valida$var1.pred", axes=FALSE)

Tarea: Distancia Inversa Ponderada (IDW)

La interpolación mediante distancia inversa ponderada determina los valores de celda a través de una combinación ponderada linealmente de un conjunto de puntos de muestra. La ponderación es una función de la distancia inversa. La superficie que se interpola debe ser la de una variable dependiente de la ubicación.

Vecindad de IDW del punto seleccionado

Este método presupone que la variable que se representa cartográficamente disminuye su influencia a mayor distancia desde su ubicación de muestra. Por ejemplo, al interpolar una superficie de poder adquisitivo de los consumidores para analizar las ventas minoristas de un sitio, el poder adquisitivo de una ubicación más distante tendrá menos influencia porque es más probable que las personas compren cerca de sus casas.

La Distancia Inversa Ponderada (IDW) es matemática (determinista) asumiendo que los valores más cercanos están más relacionados que otros con su función.

Referencias
[1] https://desktop.arcgis.com/es/arcmap/10.3/tools/spatial-analyst-toolbox/how-idw-works.htm
[2] https://acolita.com/interpolacion-con-la-distancia-inversa-ponderada-idw/