Variograma

Análisis espacial con R

Dr. Francisco Zambrano

El Variograma

El método de inverso de la distancia utiliza \(\frac{1/D^\beta}{\sum 1/D^\beta}\) como peso para la influencia de la variable a diferentes distancias.
Quizas podemos encontrar algo mejor que el inverso de la distancia para esto (ej. Moran).
Podemos obtener los pesos de forma objetiva de manera que los pesos reflejen la verdadera autocorrelación.
Matheron (1962) propone el cálculo de la denóminada semivarianza entre valores vecinos.

\[ \gamma(h) = \frac{1}{2}E[(z(s_i)-z(s_i+h))^2] \]

\(Z_(s_i)\): valor de la variable en el punto \(s_i\)
\(Z_(s_i+h)\): valor de la variable en el punto \(s_i\) más una distancia \(h\)

Si tenemos \(n\) puntos generamos \(n\cdot (n-1)/2\) pares entre ellos.

El Variograma

Si consideramos un procesos isotrópico

\[\gamma(\tilde h_j) = \frac{1}{2\cdot N_h} \sum_{i=1}^{N_h} [(z(s_i)-z(s_i+h))^2]\]

La nuve de puntos que se genera con el cálculo de las semivarianza no es simple de graficar.
Por lo que se promedia la semivarianza dentro de una distancia de retraso (lag)
Lo que esperamos encontrar es que las semivarianzas son valores bajos a cortas distancias y que aumentan y se estabiliza a una cierta distancia.
Esto se puede interpretar como: las variables respuesta son más similares a cortas distancia hasta aumentar en donde en cierto punto donde la diferencia entre los pares de puntos es más menos igual a la varianza global.

El Variograma

El variograma es ajustado por mínimos cuadrados.
Los pesos son determinados por \(\frac{N_j}{j_j^2}\). Donde \(N_j\) es la cantidad de puntos a cierta distancia de retraso y \(h_j\) es la distancia.

El Variograma

Ejemplo: con los datos `meuse`

library(sf)
library(sp)
library(gstat)
data(meuse)
v <- variogram(log(zinc)~1,st_as_sf(meuse,coords = c('x','y')),cloud=TRUE)

El Variograma

Ejemplo: con los datos `meuse`

Para \(n\cdot (n-1)/2\) pares de puntos = 11.935

plot(v)

El Variograma

Ejemplo: con los datos `meuse`

Sin embargo el variograma experimental se va ajustando a distintos intervalos \(h_{ij}=|| s_i-s_j ||\)

hscat(log(zinc)~1,st_as_sf(meuse,coords=c('x','y')),(0:9)*100)

El Variograma

Ejemplo: con los datos `meuse`

Variograma experimental

plot(variogram(log(zinc)~1,st_as_sf(meuse,coords = c('x','y'))),plot.numbers =TRUE)

Cutoff, Lag width, dependencia en la dirección

El comando:

plot(variogram(log(zinc)~1,st_as_sf(meuse,coords = c('x','y'))))

calcula y grafíca el variograma
hace varias decisiones por nosotros (default)

Cutoff, Lag width, dependencia en la dirección

Dirección

Decide que la dirección es ignorada
Los pares de puntos son elegidos de acuerdo a distancia de seperación, no dirección.
Podemos elegir buscar por dirección.

plot(variogram(log(zinc)~1,st_as_sf(meuse,coords = c('x','y'))),alpha=c(0,45,90,135))

Cutoff, Lag width, dependencia en la dirección

Dirección

Cutoff, Lag width, dependencia en la dirección

Cutoff

Distancia máxima dentro de la cual los puntos son considerados
Por defecto es el 1/3 de la diagonal de la extensión total.

Lag width

\[\frac{cutoff}{15}\]

Estos valores resultan en un ajuste inicial
Puede que no sean apropiados, por lo que pueden ser sobreescritos

Cutoff, Lag width, dependencia en la dirección

plot(variogram(log(zinc)~1,st_as_sf(meuse,coords = c('x','y')),
               cutoff =1000,width=50),
               plot.numbers = TRUE)

Cutoff, Lag width, dependencia en la dirección

El lag width no tiene que tener un valor fijo y pueden ser intervalos

plot(variogram(log(zinc)~1,st_as_sf(meuse,coords = c('x','y')),
               boundaries=c(0,50,100,seq(250,1500,250))),
     plot.numbers = TRUE)

Modelado del variograma

El variograma es usado para realizar predicciones espaciales.
Así como el método inverso de la distancia define el peso (influencia), otros métodos utilizan el variograma para determinar la influencia.
La forma tradicional de ajustar un variograma es ajustar un modelo paramétrico
Una vez calculado el variograma experimental debemos hacer un ajuste utilizando alguno de los diferentes modelos de variograma autorizados:
- Lineal
- esférico
- exponencial
- circular
- gausiano
- bessel
- power

Modelado del variograma

En R

show.vgms()

Modelado del variograma

En {gstat} los modelos de variogramas son construidos usando un modelo o una combinación de modelos.

(v1 <- vgm(1,'Sph',300))

  model psill range
1   Sph     1   300

(v2 <- vgm(1,'Sph',300,0.5))

  model psill range
1   Nug   0.5     0
2   Sph   1.0   300

Modelado del variograma

Lista de modelos de variogramas

vgm()

   short                                      long
1    Nug                              Nug (nugget)
2    Exp                         Exp (exponential)
3    Sph                           Sph (spherical)
4    Gau                            Gau (gaussian)
5    Exc        Exclass (Exponential class/stable)
6    Mat                              Mat (Matern)
7    Ste Mat (Matern, M. Stein's parameterization)
8    Cir                            Cir (circular)
9    Lin                              Lin (linear)
10   Bes                              Bes (bessel)
11   Pen                      Pen (pentaspherical)
12   Per                            Per (periodic)
13   Wav                                Wav (wave)
14   Hol                                Hol (hole)
15   Log                         Log (logarithmic)
16   Pow                               Pow (power)
17   Spl                              Spl (spline)
18   Leg                            Leg (Legendre)
19   Err                   Err (Measurement error)
20   Int                           Int (Intercept)

Modelado del variograma

Hay muchos modelos para ajustar al variograma experimental
La mayoría de estudios han usado: exponencial, esférico, gausiano, Matérn o power.
Para ajustar uno de estos modelos mediante mínimos cuadrados, debemos hacer varios pasos:
1. Elegir un modelo apropiado
2. Elegir valores iniciales para los parámetros del modelo (sill, range, nugget)
3. Ajustar el modelo utilizando uno de los criterios de ajuste

Modelado del variograma

Para el variograma de log(zinc)

v <- variogram(log(zinc)~1,st_as_sf(meuse,coords = c('x','y')))
plot(v)

Modelado del variograma

El modelo esférico y quizas el exponencial lucen como buenas alternativas
Necesitamos valores iniciales
Revisemos los parametros de la función vgm para definir un modelo de variograma y sus parametros iniciales
La función fit.variogram aplica mínimos cuadrados para ajustar los parámetros del modelo de variograma

fit.variogram(v,vgm(psill=1,model='Sph',range=800,nugget=1))

  model      psill    range
1   Nug 0.05065923   0.0000
2   Sph 0.59060463 896.9976

Modelado del variograma

Pero que pasa si elegimos valores muy lejos de valores rasonables.

fit.variogram(v,vgm(psill=1,model='Sph',range=10,nugget=1))

  model psill range
1   Nug     1     0
2   Sph     1    10

Modelado del variograma

Como quedó el ajuste de nuestro modelo de variograma

vfit <- fit.variogram(v,vgm(psill=1,model='Sph',range=800,nugget=1)) 
plot(v,vfit)

Modelado del variograma

Anisotropía

Análizar la autocorrelación espacial en diferentes direcciones

v.dir <- variogram(log(zinc)~1,st_as_sf(meuse,coords = c('x','y')),alpha=(0:3)*45)
plot(v.dir)

Modelado del variograma

Anisotropía

El variograma de la anisotropía en 2D puede ser modelado utilizando una elipse en vez de una circunferencia-

Modelado del variograma

Anisotropía

v.anis <- vgm(0.6,"Sph",1600,0.05,anis = c(45,0.3))
plot(v.dir,v.anis)

Modelado del variograma

Anisotropía

El modelo ajustado tiene un rango de 1600 en el eje principal (45° NE)
Y de \(0.3 \cdot 1600 = 480\) en la dirección menor (135°)
El paquete {gstat} no proporciona ajuste automático del variograma con anisotropía.

Modelado del variograma

Mapa del variograma

Cuando tenemos mucha información uno puede querer ver un variograma en forma de mapa.

plot(variogram(log(zinc)~1,st_as_sf(meuse,coords = c('x','y')), map = TRUE, cutoff=1000, width=100))

Modelado del variograma

Multiples variables

Variables colocadas

\[\gamma(\tilde h_ij) = \frac{1}{2\cdot N_h} \sum_{i=1}^{N_h} [z_i(s)-z_i(s+h)]\cdot [z_j(s)-z_j(s+h)]\]

Variables no colocadas

\[\gamma(\tilde h_ij) = \frac{1}{2\cdot N_h} \sum_{i=1}^{N_h} [z_i(s)-m_i]\cdot [z_j(s)-m_j]\]

Modelado del variograma

Multiples variables

Cuando tenemos más de una variable observada en cada posición \(s(x,y)\)
En {gstat} ordenamos y copiamos las variables
En el caso de los datos meuse para el caso de 5 parámetros de concentraciones zinc, Cadmium,Cupper y lead

meuse.sf <- st_as_sf(meuse,coords = c('x','y'))
g <- gstat(NULL,"logCd",log(cadmium)~1,meuse.sf)
g <- gstat(g, "logCu",log(copper)~1,meuse.sf)
g <- gstat(g, "logPb",log(lead)~1,meuse.sf)
g <- gstat(g, "logZn",log(zinc)~1,meuse.sf)

Modelado del variograma

Multiples variables

data:
logCd : formula = log(cadmium)`~`1 ; data dim = 155 x 12
logCu : formula = log(copper)`~`1 ; data dim = 155 x 12
logPb : formula = log(lead)`~`1 ; data dim = 155 x 12
logZn : formula = log(zinc)`~`1 ; data dim = 155 x 12

El variograma experimental

vm <- variogram(g)

Modelado del variograma

Multiples variables

Para ajustar con múltiples variable utilizamos fit.lmc

vm.fit <- fit.lmc(vm,g,vgm(1,"Sph",800,1))
plot(vm,vm.fit)

Modelado del variograma

Multiples variables

Como se puede observarlas variables tienen una fuerte correlación cruzada.

Podemos ver sus correlaciones.

cor(as.data.frame(meuse)[c('cadmium','copper','lead','zinc')])

          cadmium    copper      lead      zinc
cadmium 1.0000000 0.9254499 0.7989466 0.9162139
copper  0.9254499 1.0000000 0.8183069 0.9082695
lead    0.7989466 0.8183069 1.0000000 0.9546913
zinc    0.9162139 0.9082695 0.9546913 1.0000000

Modelado del variograma

Variograma de los residuos

El variograma de los residuos es utilizado cuando, por ejemplo, tengo el ajuste de un modelo de regresión.
En el caso de los datos meuse

\[log[\hat Z(s)] = \beta_0 + \beta_1\cdot \sqrt{D(s)} + \epsilon\] Los residuso

\[log[Z(s)]-log[\hat Z(s)]\]

Modelado del variograma

Variograma de los residuos

plot(variogram(log(zinc)~sqrt(dist),meuse.sf))

Modelado del variograma

Variograma de los residuos

Para el modelado de los residuos {gstat} usa OLS
OLS considera que las observaciones son independientes.
Para considerar la estructura dependiente presente (espacial), pueden considerarse Generalized Least Square (GLS)

Modelado del variograma

Variograma de los residuos

Comparemos entre OLS y GLS

OLS

vt <- variogram(log(zinc)~sqrt(dist),meuse.sf)
vt.fit <- fit.variogram(vt, vgm(1,"Exp",300,1))
plot(vt,vt.fit)

Modelado del variograma

Variograma de los residuos

Comparemos entre OLS y GLS

GLS

g.wls <- gstat(NULL,"log-zinc",log(zinc)~sqrt(dist),meuse.sf,model=vt.fit,set=list(gls=1))
(variogram(g.wls)$gamma - vt$gamma)/vt$gamma

 [1]  2.465425e-05 -9.643664e-05 -2.069750e-04 -3.034732e-04 -5.839567e-05
 [6] -1.249649e-04  2.066984e-04  1.135141e-04  3.844557e-07 -4.096026e-05
[11] -2.551952e-04  2.041244e-04 -1.660468e-04  1.996696e-04  3.291614e-05

Variograma

El Variograma

El Variograma

El Variograma

El Variograma

Ejemplo: con los datos meuse

El Variograma

Ejemplo: con los datos meuse

El Variograma

Ejemplo: con los datos meuse

El Variograma

Ejemplo: con los datos meuse

Cutoff, Lag width, dependencia en la dirección

Cutoff, Lag width, dependencia en la dirección

Dirección

Cutoff, Lag width, dependencia en la dirección

Dirección

Cutoff, Lag width, dependencia en la dirección

Cutoff

Lag width

Cutoff, Lag width, dependencia en la dirección

Cutoff, Lag width, dependencia en la dirección

Modelado del variograma

Modelado del variograma

Modelado del variograma

Modelado del variograma

Modelado del variograma

Modelado del variograma

Modelado del variograma

Modelado del variograma

Modelado del variograma

Modelado del variograma

Modelado del variograma

Anisotropía

Modelado del variograma

Anisotropía

Modelado del variograma

Anisotropía

Modelado del variograma

Anisotropía

Modelado del variograma

Mapa del variograma

Modelado del variograma

Multiples variables

Modelado del variograma

Multiples variables

Modelado del variograma

Multiples variables

Modelado del variograma

Multiples variables

Modelado del variograma

Multiples variables

Modelado del variograma

Variograma de los residuos

Modelado del variograma

Variograma de los residuos

Modelado del variograma

Variograma de los residuos

Modelado del variograma

Variograma de los residuos

Modelado del variograma

Variograma de los residuos

Ejemplo: con los datos `meuse`

Ejemplo: con los datos `meuse`

Ejemplo: con los datos `meuse`

Ejemplo: con los datos `meuse`