Análisis espacial con R
\[Z(s) = m(s) + \epsilon'(s)+ \epsilon''\]
{fig-align:“center”}
\[\hat z(s_0) = \hat m(s_0) + \hat e(s_0)\] \[\hat z(s_0) = \sum_{k=0 }^p\hat \beta_k\cdot q_k(s_0)+\sum_{i=1}^n \lambda_i\cdot [z(s_0)-q_k(s_0)\cdot \hat\beta_k]\]
\[\hat \sigma^2_E(s_0)=E\left[(\hat z(s_0)-z(s_0))\cdot (\hat z(s_0)-z(s_0))^T\right]\]
\[E[\hat z(s_0)-z(s_0)]=0\]
\[z(s)=q^T\cdot \beta+\epsilon(s)\] \[E[\epsilon(s)]=0\] \[E[\epsilon\cdot \epsilon^T(s)]=C\]
BLUP\[\hat z(s_0) = q_0^T\cdot\hat \beta + \hat \lambda_0^T\cdot(z-q\cdot \beta)\]
\[\hat \beta =(q^T\cdot C^{-1}\cdot q)^{-1}\cdot q^T\cdot C^{-1}\cdot z\] \[\hat \lambda_0=C^{-1}\cdot c_0\] O en forma simplificada
\[\hat z(s_0) = \mu + \hat \lambda^T\cdot(z-\mu)\]
Bajo los supuestos:
Son prácticamente lo mismo:
Se diferencian en el método de resolver el problema pero no en los resultados.
Predicción espacial se basa en la estimación de valores desconocido \(z(s_0)\) basado en datos conocidos \(Z(s_i)\)
\(Z_{s_i}\) estará compuesto por una parte determinística (tendencia, ej, OLS, GLS) y la variación de la correlación espacial (variograma)
Los métodos de predicción espacial geoestadísticos son derivados del BLUP.
Son casos particulares en donde se considera o no la parte determinística \((\mu)\).
Los modelos de predicción espacial se basan en la autocorrelación espacial y utilizan el modelo de variograma para el cálculo de la covarianza espacial.
Por decadas Kriging ha sido sinónimo de interpolación geoestadística.
Tiene su origen en la industria minera en los años cincuenta.
La idea original viene del ingeniero minero D.G. Krige y del estadístico H.S. Sichel.
Primera vez publicada en 1951.
Pero tomo casí una decada hasta que el matemático Frances G. Matheron derivara las formulas y estableciera todo el campo de geoestadística lineal.
Una versión estandar de Kriging es Kriging Ordinario
Un caso particular del modelo universal de variación donde solo considera la variación de la correlación espacial.
Por lo tanto no tenemos coeficientes de regresión y \(\beta = 0\)
\[\hat z(s_0) = q_0^T\cdot\hat \beta + \hat \lambda_0^T\cdot(z-q\cdot \beta)\]
\[\hat z_{OK}(s_0) = \hat \lambda_0^T\cdot z\]
De alguna manera podriamos ver a Kriging como una sofisticación del método de inverso de la distancia (IDW).
\(\hat \lambda_0\) son los pesos de kriging.
Los pesos reflejan la verdadera estructura de autocorrelación espacial (variograma).
\[\lambda_0=C^{-1}\cdot c_0; \: C(|h|=0)=C_0+C_1\]
Las salidas de cualquier modelo de predicción estadística son dos mapas:
El promedio de la varianza de todos los puntos es la varianza global
La varianza global nos sirve para evaluar la precisión
\[\hat \sigma^2_{OK}(s_0) = C_0+C_1-\sum_{i=1}^n w_i(s_0)\cdot C(s_0,s_i) + \phi\]
{gstat}library(gstat)
library(sp)
library(sf)
library(stars)
data(meuse)
meuse_sf <- st_as_sf(meuse,coords = c('x','y'),crs = 28992)
data(meuse.grid)
meuse_grid_st <- meuse.grid |> st_as_stars() |> st_set_crs(28992)
zinc.ok <- krige(log(zinc)~1, meuse_sf, meuse_grid_st, model=vgm(psill=0.714, "Exp", range=449, nugget=0))[using ordinary kriging]{gstat}
{gstat} utiliza simulaciones sequenciales por defectoEste es un algoritmo simple que recorre aleatoriamente las ubicaciones de predicción y en cada ubicación:
\[\hat z_{RK}(s_0) = \sum_{k=0 }^p\hat \beta_k\cdot q_k(s_0)+\sum_{i=1}^n \lambda_i\cdot [z(s_0)-q_k(s_0)\cdot \hat\beta_k]\]
.center[
]
Lo que modelamos geoestadísticamente son los residuos \(e(s_i) = z(s_0)-\hat z_{OLS/GLS}(s_0)\)
Obtenemos el variograma de los residuos.
Aplicamos Kriging ordinario para espacializar los residuos.
Luego sumamos la parte deterministica con los resiudos interpolados
{fig-align=‘center’)
Determinar un modelo de regresión lineal: \(z=6.64-0.195\cdot q\)
Derivar los residuos en todos los puntos observados: \(e^*(s_i)=z(s_i)-[b_0+b_1\cdot q(s_i)]\)
Model la estructura de la covarianz de los residuos (variograma):
(5. Sumar el residuo interpolado con el estimado de la regresión lineal:
\[\hat z(5,5)=\beta_0+\beta_1\cdot q_i+\sum_{i=1}^n \lambda_i(s_0)\cdot e(s_i)\]
\[\hat z(5,5)=6.68-0.199\cdot 12-.081=4.21\]
En este caso específico un poco diferente del resultado obtenido con \(\hat z_{OK}=4.30\)
En R
zinc.rk <- krige(log1p(zinc)~dist+soil, meuse_sf, meuse_grid_st,
model=vgm(psill=0.151, "Exp", range=374, nugget=0.055))[using universal kriging]Entonces, para predicción espacial hemos visto diferentes modelos de espacialización.
¿Cómo sabemos que modelo de predicción aplicar?
{fig-align=‘center’)
La varianza de estimación de OK o los modelos de regresión es la estimación estadística de la incerteza del modelo.
Pero el verdadero poder de predicción solo puede ser evaluado utilizando un set de datos de control independiente.
Así, el error de predicción es denominado como la precisión de la predicción.
La verdadera calidad de una predicción espacial puede ser evaluada de mejor forma al comparar valores estimado \((\hat z)\) con valores observados en punto de validación \((z^*)\).
Debido a que la recolección de datos independientes es muchas veces impracticable y costoso, la validación se hace mediante validación cruzada.
Validación cruzada: se separa el set de datos original en dos set de datos
k-fold: el set original de datos es separado en k grupos iguales y cada uno es usado para validación cruzada.
leave-one-out (LOO): cada punto es usado para validación cruzada.
Jackknifing: similar a LOO pero es utilizada para evaluar el sesgo del análisis estadístico y no de las predicciones.
LOOCV
K-FOLD
{fig-align=‘center’)
Generalmente calculamos las siguientes indicadores de calidad en CV
\[ME=\frac{1}{l}\cdot \sum_{j=1}^l\left[\hat z(s_j)-z^*(s_j)\right];\: E[ME]=0\] - La raiz cuadrada del error cuadrático medio (RMSE)
\[RMSE = \sqrt{\frac{1}{l}\cdot \sum_{j=1}^l \left[ \hat z(s_j)-z^*(s_j)\right]^2};\: E[RMSE]=0\]
O en forma estandarizada lo podemos escribir como:
\[RMSE_r = \frac{RMSE}{\sigma}\]
\[R^2=1-\frac{\sum_{j=1}^l [\hat z(s_j) - z^*(s_j)]^2}{\sum_{j=1}^l [z^*(s_j)-\bar z^*(s_j)]^2}\]
Ejemplo En R. Dejamos 100 datos para modelar y 55 para validar
Calculamos el \(R^2\) utilizando los valores cross-validados
En R {gstat} tiene la función krige.cv quwe permite realizar CV utilizando LOO o k-fold
Si vemos los puntos en el mapa
Si vemos los puntos en el mapa
Si vemos los puntos en el mapa
Si vemos los puntos en el mapa
