1. Introducción

Dentro del análisis de los datos espaciales y la estadística espacial, uno de los conceptos más importantes a conocer es la autocorrelación espacial. Esta suele complicar las pruebas estadísticas y facilitar la interpolación espacial.

La autocorrelación (sea esta espacial o temporal) es una medida de similitud (correlación) entre observaciones cercanas. Es así que, a fines de entender mejor la autocorrelación espacial, revisaremos primero la autocorrelación temporal.

1.1 Autocorrelación temporal

Cuando medimos un fenómeno a lo largo del tiempo, es probable que dos observaciones que se encuentran cercanas en la dimensión tiempo sean similares. Por ejemplo, es probable que si una persona subió de 50 a 80 kg de peso, lo haya hecho de manera gradual, pasando de 50 a 60, de 60 a 70 y de 70 a 80, en lugar de, por ejemplo, pasar de 50 a 100, y luego de 100 a 80. Cuando un fenómeno se comporta de esta manera, decimos que existe autocorrelación temporal. Veamos un ejemplo con datos simulados, cargando primero las librerías necesarias.

# Librerías necesarias
library(tidyverse)
library(forecast)
library(ggpubr)
library(raster)
library(sf)
library(spdep)
library(ggrepel)
# Declaramos la función select
select <- dplyr::select
data_autocor_t <- tibble(peso = rnorm(100, 75, 6))
data_autocor_t %>% head()
## # A tibble: 6 x 1
##    peso
##   <dbl>
## 1  87.6
## 2  75.0
## 3  79.6
## 4  71.3
## 5  73.9
## 6  71.1

Grafiquemos el peso contra su rezago en el tiempo (ordenándolo de menor a mayor y sin ordenarlo):

autocor_peso_desorden <- data_autocor_t %>% 
  ggplot(aes(x=peso,y=lag(peso)))+
  geom_point()+
  labs(title="Sin autocorrelación \ntemporal",
       x="peso en t", y="peso en t-1")
autocor_peso_orden <- data_autocor_t %>% 
  arrange(peso) %>% 
  ggplot(aes(x=peso,y=lag(peso)))+
  geom_point()+
  labs(title="Con autocorrelación \ntemporal", 
       x="peso en t", y="peso en t-1")
ggarrange(autocor_peso_desorden, autocor_peso_orden, ncol=2)
## Warning: Removed 1 rows containing missing values (geom_point).

## Warning: Removed 1 rows containing missing values (geom_point).

Como se puede notar, si existiese autocorrelación temporal en los datos, se vería una clara tendencia entre los datos en el tiempo t versus los datos en el tiempo t-1. Para hacerlo más evidente, utilicemos un correlograma. En el análisis de series temporales, el correlograma, también conocido como un gráfico de autocorrelación, es una representación gráfica de las autocorrelaciones de la muestra versus el tiempo.

acf_desorden <- ggAcf(data_autocor_t$peso)+
  labs(title="Sin autocorrelación \ntemporal")
acf_orden <- ggAcf(sort(data_autocor_t$peso))+
  labs(title="Con autocorrelación \ntemporal")
ggarrange(acf_desorden, acf_orden, ncol=2)

Cuando existe autocorrelación temporal, los datos mostrarán un comportamiento como el del panel derecho, donde la altura de cada línea muestra la intensidad de la relación entre el periodo t y el periodo t-1.

1.2. Autocorrelación espacial

La autocorrelación espacial es muy similar a la autocorrelación temporal con un toque de dificultad adicional, dado que el tiempo se mueve en una sola dimensión, mientras que el espacio tiene al menos dos dimensiones. Así, las medidas de autocorrelación espacial muestran el grado al que las observaciones en ubicaciones espaciales (sean estas puntos, áreas o rásters) son similares entre sí. Entonces, para calcular la autocorrelación espacial necesitaremos dos ingredientes: observaciones y ubicaciones.

Cabe notar que la autocorrelación espacial en una variable puede ser exógena o endógena. Si es que la relación es de naturaleza exógena, la variable está causada por otra variable espacialmente autocorrelacionada. Si la relación de la naturaleza es endógena, la variable está causada por el mismo fenómeno.

Una de las formas más comunes de medir la autocorrelación espacial es a través de la I de Morán; sin embargo, existen otros índices como el C de Geary y el índice de recuento de combinaciones.

Analicemos ahora un ejemplo de autocorrelación espacial para Diekirch.

# Cargamos datos de ejemplo
diekirch <- shapefile(system.file("external/lux.shp", package = "raster"))
diekirch <- diekirch[diekirch$NAME_1=="Diekirch",]
# Creamos una variable adicional
diekirch$valor <- c(10,6,4,11,6)
diekirch@data
##   ID_1   NAME_1 ID_2   NAME_2 AREA valor
## 0    1 Diekirch    1 Clervaux  312    10
## 1    1 Diekirch    2 Diekirch  218     6
## 2    1 Diekirch    3  Redange  259     4
## 3    1 Diekirch    4  Vianden   76    11
## 4    1 Diekirch    5    Wiltz  263     6

Sobre estos datos, supongamos que nos interesa analizar la autocorrelación espacial sobre la variable AREA. Si es que dicha autocorrelación existiese, las regiones de tamaño similar se agruparían espacialmente. Veamos entonces los polígonos con un gráfico.

diekirch %>% 
  st_as_sf() %>% 
  ggplot()+
  geom_sf(aes(fill=AREA))+
  coord_sf()+
  scale_fill_gradientn(colors = rev(terrain.colors(3)))+
  geom_sf_label(aes(label=ID_2))
## Warning in st_point_on_surface.sfc(sf::st_zm(x)): st_point_on_surface may not
## give correct results for longitude/latitude data

En efecto podemos ver que los polígonos parecen estar ordenados por tamaño (4-2-3-5-1).

2. Cálculo de la I de Morán

El índice de autocorrelación espacial de Moran viene dado por:

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

Donde \(y\) representa los valores de la variable de análisis y \(w\) representa las entradas de la matriz de ponderaciones espaciales.

Esta expresión matemática podría causarnos temor, pero en realidad es una versión ampliada del cálculo del coeficiente de correlación, a la que se le agrega una matriz de ponderaciones espaciales.

Veamos entonces cómo realizar este cálculo. Para ello, determinaremos primero la matriz de ponderaciones espaciales.

vecinos <- poly2nb(diekirch, row.names = diekirch$Id)
vecinos
## Neighbour list object:
## Number of regions: 5 
## Number of nonzero links: 14 
## Percentage nonzero weights: 56 
## Average number of links: 2.8
matriz_pe <- nb2mat(vecinos, style = "B")
matriz_pe
##   [,1] [,2] [,3] [,4] [,5]
## 0    0    1    0    1    1
## 1    1    0    1    1    1
## 2    0    1    0    0    1
## 3    1    1    0    0    0
## 4    1    1    1    0    0
## attr(,"call")
## nb2mat(neighbours = vecinos, style = "B")

Similar a como lo hicimos en la clase anterior, observemos esta matriz de manera gráfica:

vecinos_sf <- as(object = nb2lines(nb = vecinos,
                                              coords = coordinates(diekirch)),
                            'sf') %>% 
  st_set_crs(st_crs(diekirch))
## Warning in CRS(proj4string): CRS: projargs should not be NULL; set to NA
diekirch %>% 
  st_as_sf() %>% 
  ggplot() +
  geom_sf() +
  geom_sf(data = vecinos_sf, color="red", size=1) +
  coord_sf()

Una vez que tenemos la matriz de ponderaciones calculadas, pasemos a la construcción de los términos incluidos en la I de Morán.

\(n\):

# Número de observaciones
n <- nrow(diekirch@data)
n
## [1] 5

\(y\) e \(\overline{y}\):

y <- diekirch$AREA
y_media <- mean(y)
y_media
## [1] 225.6

Creemos \((y_i-\overline{y})(y_j-\overline{y})\):

yiyj <- expand.grid(y-y_media, y-y_media)[,1]*expand.grid(y-y_media, y-y_media)[,2]
yiyj
##  [1]   7464.96   -656.64   2885.76 -12925.44   3231.36   -656.64     57.76
##  [8]   -253.84   1136.96   -284.24   2885.76   -253.84   1115.56  -4996.64
## [15]   1249.16 -12925.44   1136.96  -4996.64  22380.16  -5595.04   3231.36
## [22]   -284.24   1249.16  -5595.04   1398.76
yiyj <- rep(y-y_media, each = n)*rep(y-y_media)
yiyj
##  [1]   7464.96   -656.64   2885.76 -12925.44   3231.36   -656.64     57.76
##  [8]   -253.84   1136.96   -284.24   2885.76   -253.84   1115.56  -4996.64
## [15]   1249.16 -12925.44   1136.96  -4996.64  22380.16  -5595.04   3231.36
## [22]   -284.24   1249.16  -5595.04   1398.76

Con este dato calculado construimos una matriz, la cual debe ser multiplicada por la matriz de ponderación espacial, es decir: \(w_{ij}(y_i-\overline{y})(y_j-\overline{y})\).

wyiyj <- matrix(yiyj, ncol = n, nrow = n)*matriz_pe
wyiyj
##        [,1]    [,2]    [,3]      [,4]    [,5]
## 0      0.00 -656.64    0.00 -12925.44 3231.36
## 1   -656.64    0.00 -253.84   1136.96 -284.24
## 2      0.00 -253.84    0.00      0.00 1249.16
## 3 -12925.44 1136.96    0.00      0.00    0.00
## 4   3231.36 -284.24 1249.16      0.00    0.00
## attr(,"call")
## nb2mat(neighbours = vecinos, style = "B")

Luego sumamos los valores de esta matriz para obtener \(\sum_{i=1}^{n}\sum_{j=1}^{n}w_{ij}(y_i-\overline{y})(y_j-\overline{y})\):

sum_wyiyj <- sum(wyiyj)
sum_wyiyj
## [1] -17005.36

Creamos el término \(\sum_{i=1}^{n}\sum_{j=1}^{n}w_{ij}\):

sum_w <- sum(matriz_pe)
sum_w
## [1] 14

Y finalmente el término \(\sum_{i=1}^{n}(y_i-\overline{y})^2\):

sum_yi2 <- sum((y-y_media)^2)
sum_yi2
## [1] 32417.2

Con todos estos terminos podemos construir el indicador:

IM <- (n/sum_yi2)*(sum_wyiyj/sum_w)
IM
## [1] -0.1873494

Este indicador en sí mismo no tiene una interpretación hasta que sea comparado contra un valor esperado de la I de Morán, el cual corresponde al valor de \(\frac{-1}{(n-1)}\). Es decir:

IME <- -1/(n-1)
IME
## [1] -0.25

Realicemos la prueba de la I de Morán con el paquete spdep.

# Creamos una lista de vecinos, necesario para la función moran
lista_pe <- nb2listw(vecinos, style = "B")
lista_pe
## Characteristics of weights list object:
## Neighbour list object:
## Number of regions: 5 
## Number of nonzero links: 14 
## Percentage nonzero weights: 56 
## Average number of links: 2.8 
## 
## Weights style: B 
## Weights constants summary:
##   n nn S0 S1  S2
## B 5 25 14 28 168
moran(diekirch$AREA, lista_pe, n=length(lista_pe$neighbours), S0=Szero(lista_pe))
## $I
## [1] -0.1873494
## 
## $K
## [1] 2.663502

Como podemos ver, el resultado de I es igual al anterior, y lo que nos queda es realizar el test de Morán, que tiene como hipótesis nula que no existe autocorrelación espacial.

moran.test(diekirch$AREA, lista_pe, randomisation = F)
## 
##  Moran I test under normality
## 
## data:  diekirch$AREA  
## weights: lista_pe    
## 
## Moran I statistic standard deviate = 0.34626, p-value = 0.3646
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##        -0.1873494        -0.2500000         0.0327381

Bajo este test, vemos que el p-valor es 0.3646, lo cual nos indica que no rechacemos la hipótesis nula. Sin embargo, tenemos apenas 5 observaciones, lo cual podría llevarnos a un resultado sesgado. Por ello, es mejor utilizar el test con simulación de Monte Carlo, como se muestra a continuación.

moran.mc(diekirch$AREA, lista_pe, nsim=100)
## 
##  Monte-Carlo simulation of Moran I
## 
## data:  diekirch$AREA 
## weights: lista_pe  
## number of simulations + 1: 101 
## 
## statistic = -0.18735, observed rank = 68.5, p-value = 0.3218
## alternative hypothesis: greater

Bajo esta metodología, se asignan valores aleatorios a los polígonos en estudio un número de veces determinado para obtener una distribución del valor esperado de la I de Morán y realizar así el test.

Podemos observar que en este caso el test devuelve un p-valor de 0.3267, llevándonos a la misma conclusión anterior. ¿Qué resultado obtenemos si realizamos la prueba sobre la variable valor?

moran.mc(diekirch$valor, lista_pe, nsim=100)
## 
##  Monte-Carlo simulation of Moran I
## 
## data:  diekirch$valor 
## weights: lista_pe  
## number of simulations + 1: 101 
## 
## statistic = 0.17289, observed rank = 100.5, p-value = 0.009901
## alternative hypothesis: greater

Como podemos observar, el p-valor es de prácticamente 0, sugiriendo el rechazo de la hipótesis nula. Es decir, no rechazamos la hipótesis alternativa que indica que existe autocorrelación espacial positiva (estadístico > 0). Es decir, valores altos se agrupan espacialmente, y similar sucede con valores bajos. Si existiese autocorrelación espacial negativa, valores altos rechazarían valores altos, y existiese un proceso competitivo.

Finalmente, podemos construir un gráfico de dispersion de Morán para visualizar la autocorrelación espacial, como se muestra a continuación. Es decir, exploraremos la relación entre los datos y sus vecinos como un gráfico de dispersión.

moran.plot(diekirch$valor, nb2listw(vecinos, style = "W"), return_df = T, plot = F) %>% 
  head()
##    x       wx is_inf labels      dfb.1_      dfb.x      dffit        cov.r
## 1 10 7.666667   TRUE      0  0.09446936 -0.1582700 -0.2261335 3.4240315117
## 2  6 7.750000   TRUE      1  6.59965887 -4.1744043  8.9451521 0.0002190733
## 3  4 6.000000   TRUE      2 -1.14283699  0.9486403 -1.2033144 1.7607029403
## 4 11 8.000000   TRUE      3  0.06940487 -0.1002145 -0.1244925 5.1496907004
## 5  6 6.666667  FALSE      4 -0.16620251  0.1051261 -0.2252702 2.6213489041
##       cook.d       hat
## 1 0.03688963 0.3920455
## 2 0.51088070 0.2556818
## 3 0.65971144 0.5284091
## 4 0.01155573 0.5681818
## 5 0.03544209 0.2556818
moran.plot(diekirch$valor, nb2listw(vecinos, style = "W"), return_df = T, plot = F) %>% 
  ggplot(aes(x, wx))+
  geom_point()+
  geom_smooth(method='lm', formula= y~x)+
  geom_label_repel(aes(label=labels))+
  labs(x="valor", y="rezago espacial de valor")

Como podemos notar, si trazamos una línea de regresión entre los valores y sus rezagos espaciales existe un patrón lineal.

moran.plot(diekirch$AREA, nb2listw(vecinos, style = "W"), return_df = T, plot = F) %>% 
  ggplot(aes(x, wx))+
  geom_point()+
  geom_smooth(method='lm', formula= y~x)+
  geom_label_repel(aes(label=labels))+
  labs(x="AREA", y="rezago espacial de AREA")

Si revisamos la misma gráfica para la variable AREA veamos que el intervalo de confianza es más grande y la línea es casi plana, coincidiendo con el resultado del test de Morán.

¿Has entendido todos los conceptos? No dudes en hacer preguntas. ¡Nos vemos en el próximo capítulo!