Aplicaciones de geoestadística

Maestría en Estadística Aplicada

Roberto Trespalacios

• Semivariogramas.
• Covariograma y correlograma.
• Estacionalidad fuerte y débil.
• Isotropía y anisotropía.

La existencia de la covarianza implica que la varianza existe, es finita y no depende de \( h \), es decir \( Var(Z(s_i)) =C(0)=\sigma^2 \) . Así mismo la estacionariedad de segundo orden implica la siguiente relación entre la función de semivarianza y la de autocovarianza:

\[ \begin{align*} \gamma(Z(s),Z(s+h))&=\frac{1}{2}E[Z(\boldsymbol{s}+\boldsymbol{h})-m+ Z(\boldsymbol{s})+m )]^2\\ &=\frac{1}{2} [E[Z(\boldsymbol{s}+\boldsymbol{h})+m]^2 + E[Z(\boldsymbol{s})-m]^2-2E[Z(\boldsymbol{s}+\boldsymbol{h})-m][Z(\boldsymbol{s})-m]]\\ &= \frac{1}{2} \sigma^2 +\frac{1}{2} \sigma^2-E[(Z(\boldsymbol{s}+h)-m)(Z(\boldsymbol{s})-m)]\\ &=\sigma^2-C(h) \end{align*} \]

• Cuando se definió la estacionariedad débil en el capítulo anterior se mencionó que se asumía que la varianza de los incrementos de la variable regionalizada era finita. A esta función denotada por \( 2\gamma(h) \) se le denomina variograma.
• Utilizando la definición teórica de la varianza en términos del valor esperado de una variable aleatoria, tenemos:

\[ \begin{align*} 2\gamma(h)=&Var[Z(\boldsymbol{s}+h]-Z(\boldsymbol{s})]\\ =&E[(Z(\boldsymbol{s}+h)-Z(\boldsymbol{s}))^2]-(E[Z(\boldsymbol{s}+h)-Z(\boldsymbol{s})])^2\\ =&E[(Z(\boldsymbol{s}+h)-Z(\boldsymbol{s}))^2] \end{align*} \]

• La mitad del variograma \( \gamma(h) \), se conoce como la función de semivarianza y caracteriza las propiedades de dependencia espacial del proceso.
• Dada una realización del fenómeno, la función de semivarianza es estimada, por el método de momentos, a través del semivariograma experimental. Esta, la definos como \( \hat{\gamma}(h) \)

\[ \hat{\gamma}(h)=\frac{\sum(Z(s+h)-Z(s))^2}{2n} \]

• La función de semivarianza se calcula para varias distancia \( h \).
• En la práctica, debido a irregularidad en el muestreo y por ende en las distancias entre los sitios, se toman intervalos de distancia \( \{[0,h], (h,2h],\dots,(2h,3h], \dots \} \)

• 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.
• A continuación se hace una revisión de los conceptos asociados a cada una de ellas y se describen sus bondades y limitaciones.

La función de covarianza muestral o el Covariograma entre parejas de observaciones que se encuentran a una distancia \( h \) se calcula, empleando la formula clásica de la covarianza muestral, por:

\[ \begin{align*} C(h)=&Cov[Z(\boldsymbol{s}+h]-Z(\boldsymbol{s})]\\ =&\frac{\sum_{i=1}^n(Z(\boldsymbol{s}+h)-m)(Z(\boldsymbol{s})-m)}{n}\\ =&\frac{\sum_{i=1}^n(Z(\boldsymbol{s}+h)Z(s))}{n} -m^2=C(h) \end{align*} \]

donde \( m \) representa el valor promedio en todo punto de la región de estudio y \( n \) es el número de parejas de puntos que se encuentran a una distancia \( h \).

El correlograma (coeficiente de correlación lineal) entre dos variables aleatorias espaciales \( Z(\boldsymbol{s}) \) y \( Z(\boldsymbol{s}+h) \)

\[ \begin{align*} \rho(\boldsymbol{h})=&Cor[Z(\boldsymbol{s}+h),Z(\boldsymbol{s})]\\ =&\frac{Cov[Z(\boldsymbol{s}+h),Z(\boldsymbol{s})]}{\sqrt{Var[Z(\boldsymbol{s}+h)]Var[Z(\boldsymbol{s})]}}\\ =&\frac{C(h)}{C(0)}\\ \end{align*} \]

Bajo el supuesto de estacionariedad cualquiera de las tres funciones de dependencia espacial mencionadas, es decir semivariograma, covariograma o correlograma, puede ser usada en la determinación de la relación espacial entre los datos.

Observación

Cuando existe superposición de los variogramas direccionales: el variograma es isotrópico; de lo contrario, diremos que el campo es anisotrópico.

El variograma exhibe una anisotropía más compleja: diferentes formas según la dirección.

• Análisis de isotropía espacial (en R): Usar el raster del taller 1, para hacer un análisis isotrópico espacial, de la variable temperatura (\( Z(s) \)) en cada punto del mapa. Gráfique para cada una de las diferentes posibilidades de desplazamiento del rango de ángulo \( [0, \pi] \), con incrementos de \( \pi/4 \).

• Cálculo de semivariograma-variograma(manual): Dado los 5 puntos de la variable lluvia, sobre \( D \subset \mathbb{R}^2 \), estime el variograma direccional este-oeste usar una tolerancia de 0 grados.

Punto	Coordenadas	Lluvia
1	(2,4)	33
2	(-3,0)	55
3	(1,1)	41
4	(3,-2)	36
5	(-1,2)	35

A continuación, veremos un ejemplo de autocorrelación temporal, que nos auyuda a entender la autocorrelación espacial.

set.seed(0)
d <- sample(100, 10)

Cálculo de la auto-correlación

a <- d[-length(d)]
b <- d[-1]
plot(a, b, xlab='t', ylab='t-1')

La corrrelación es

cor(a, b)

La autocorrelación calculada arriba es muy pequeña. Aunque esta es una muestra aleatoria, usted (casi) nunca obtiene un valor de cero. Calculamos la autocorrelación de “un retraso”, que es, comparamos cada valor con su vecino inmediato, y no con otros valores cercanos.

Después de los números en la autocorrelación se vuelve muy fuerte (como era de esperar).

d <- sort(d)
a <- d[-length(d)]
b <- d[-1]
plot(a, b, xlab='t', ylab='t-1')

La correlación es

cor(a, b)

La función acf muestra la autocorrelación calculada de una manera diferente (es 1 para cada punto en sí mismo, cuando se compara con el vecino más cercano, y luego se reduce gradualmente).

acf(d)

• El concepto de autocorrelación espacial es una extensión de la autocorrelación temporal.
• Aunque es un poco más complicado; ya que el tiempo es unidimensional, y solo va en una dirección, siempre adelante.

• Los objetos espaciales tienen (al menos) dos dimensiones y formas complejas, y puede no ser obvio cómo determinar qué está “cerca”.

• Las ecuaciones espaciales de autocorrelación espacial describen el grado en que son puntos, áreas o celdas de trama, son similares entre sí. Entonces necesitamos dos cosas: observaciones y ubicaciones.

• La autocorrelación espacial en una variable puede ser:

  • exógena (es causada por otra variable espacialmente autocorrelacionada, por ejemplo, lluvia)
  • endógena (es causada por el proceso en juego, por ejemplo, la diseminación de una enfermedad)

• Una estadística comúnmente utilizada que describe la autocorrelación espacial es el estadístico de Moran

• Otros índices incluyen C de Geary y, para datos binarios, el índice de recuento de uniones.

• El semi-variograma(del que hablamos anteriormente) también expresa la cantidad de autocorrelación espacial en un conjunto de datos.

El siguiente ejemplo se basa en un problema de autocorrelación de regiones.

Digamos que estamos interesados en la autocorrelación espacial en la variable “área”. Si hubiera autocorrelación espacial, las regiones de un tamaño similar estarían agrupadas espacialmente.

Aquí hay un diagrama de los polígonos. Utilizo las coordenadas para obtener los centroides de los polígonos para colocar las etiquetas.

library(raster)
p <- shapefile(system.file("external/lux.shp", package="raster"))
p <- p[p$NAME_1=="Diekirch", ]
p$value <- c(10, 6, 4, 11, 6)
data.frame(p)

Veamos la gráfica

par(mai=c(0,0,0,0))

xy <- coordinates(p)
plot(p, col=2:7)
points(xy, cex=6, pch=20, col='white')
text(p, 'ID_2', cex=1.5)

Ahora necesitamos determinar qué polígonos están “cerca” y cómo cuantificar eso. Aquí usaremos la adyacencia como criterio. Para encontrar polígonos adyacentes, podemos usar la librería spdep.

library(spdep)
w <- poly2nb(p, row.names=p$Id)
class(w)
summary(w)

• El summary(w) nos dice algo sobre el vecindario.
• El número promedio de vecinos es 2.8, 3 polígonos tienen 2 vecinos y 1 tiene 4 vecinos (¿cuál es ese?).

Para más detalles, podemos ver la estructura de w.

str(w)

Pregunta de tarea: ¿Qué el significan las primeras 5 líneas devueltas por str(w)

Veamos los enlaces entre los polígonos para cada centro.

plot(p, col='gray', border='blue', lwd=2)
plot(w, xy, col='red', lwd=2, add=TRUE)

Podemos transformar w en una matriz de ponderaciones espaciales. Una matriz de ponderaciones espaciales refleja la intensidad de la relación geográfica entre las observaciones.

wm <- nb2mat(w, style='B') #W

Ahora veamos el índice de autocorrelación espacial de Moran.

• El estadístico de Moran es una medida de la autocorrelación espacial desarrollada por Patrick Moran.
• La autocorrelación espacial se caracteriza por una correlación en una señal de ubicaciones cercanas en el espacio.
• Recoredemos, que la autocorrelación espacial es más compleja que la autocorrelación unidimensional porque la correlación espacial es multidimensional (es decir, espacio de 2 o 3 dimensiones) y multidireccional.

\[ I = \frac{n}{\sum_{i=1}^n (y_i - \bar{y})^2} \frac{\sum_{i=1}^n \sum_{j=1}^n w_{ij}(y_i - \bar{y})(y_j - \bar{y})}{\sum_{i=1}^n \sum_{j=1}^n w_{ij}} \]

donde \( n \) es el número de unidades espaciales indexados por \( i \) y \( j \); \( y \) es la variable de interés; \( \bar{y} \) es la media de \( y \); y \( w_{ij} \) es un elemento de una matriz de pesos espaciales.

El valor esperado del estadístico de Moran bajo la hipótesis nula de no autocorrelación espacial es

\[ E(I)=\frac{-1}{n-1} \]

• Después de hacer esto “a mano”“, ahora usemos el paquete spdep para calcular Moran y hacer una prueba de significación.
• Para hacer esto, necesitamos crear un objeto de pesos espaciales tipo 'listw' (en lugar de la matriz que utilizamos anteriormente).
• Para obtener el mismo valor, usamos "style = 'B'” para usar pesos de distancia binarios (TRUE/FALSE).

ww <- nb2listw (w, style = 'B')

• Ahora podemos usar la función moran.
• La función se define como moran(y, ww, n, Szero(ww)).

moran(p$value, ww, n=length(ww$neighbours), S0=Szero(ww))

Ahora podemos probar la significancia. Primero analíticamente, usando lógica basada en regresión lineal y suposiciones.

moran.test(p$value, ww, randomisation=FALSE)

• En lugar del enfoque anterior, debe usar la simulación de Monte Carlo.
• Ese es el método preferido (de hecho, el único método bueno).
• La idea es que los valores se asignen aleatoriamente a los polígonos, y se calcula el estadístico de Moran.
• Esto se repite varias veces para establecer una distribución de valores esperados.
• El valor observado del estadístico de Moran se puede comparar con la distribución simulada.

#Cual es el numero maximo de simulaciones
moran.mc(p$value, ww, nsim=ns)

Usar los datos del raster del mapa de Colombia para evaluar el estadístico de Moran (todo lo que se hizo en clase).(recuerde que debe excluir a San Andres y providencia)