4. Autocorrelación Espacial

Importando librerías:

GWmodel(Modelos ponderados geográficamente): Paquete que introduce técnicas de una rama particular de las estadísticas espaciales, denominadas modelos ponderados geográficamente.
spdep: (Dependencia espacial: esquemas de ponderación, estadísticas y modelos):
- Es una colección de funciones para crear objetos de matriz de pesos espaciales para resumir estos objetos y analizar datos espaciales.
- Contiene una colección de pruebas para la autocorrelación espacial, incluida la: I global de Moran y la C de Geary, la I local de Moran, entre otras.
maptools (Herramientas para manejar objetos espaciales): Conjunto de herramientas para manipular datos geográficos.

library(GWmodel)
library(spdep)
library(maptools)

Leyendo y analizando datos

Es una base de datos (COUNTY_ATLANTIC.shp) que equivalen a 666 coordenadas de los condados correspondientes a las 13 colonias que tenía Reino Unido en la costa este de América del Norte, además se unen mediante la función merge con otra base de datos (data_atlantic_1998_2012.csv) que contiene información sobre la pobreza, el Pm25, entre otras en la misma región.

Nota: Los datos fueron tomados del anterior ejercicio: 3-TareaGWSS.R, allí podrá entender, y observar una descripción minuciosa de estos.

Condados = readShapeSpatial(file.choose()) # Leyendo datos espaciales
Datos = read.csv(file.choose(), header=T)
SPDF = merge(Condados, Datos, by = "FIPS") # Uniendo datos
# names (SPDF) # nombres del objeto SPDF
# dim(SPDF)   #--> [1] 666  20
# class(SPDF) #--> "SpatialPolygonsDataFrame" sp
#--------------------------------------------------------------------
# Distribución espacial de la tasa de pobreza.
spplot(SPDF, "POV", main="Distribución espacial de la tasa de pobreza \nen las 11 colonias")

Extraemos las coordenadas

Además, extraemos las coordenadas \((x,y)\) de la base de datos espacial y la almacenamos en la variable coords, para su posterior uso en el código.

# Extrayendo coordenadas desde la base de datos espaciales.
coords = coordinates(SPDF)

Creando la lista de vecinos usando el criterio de continuidad tipo reina (QUEEN), torre (ROOK) y K-Vecinos (KN-4)

poly2nb: (Construye una lista de vecinos a partir de una lista de polígonos): La función crea una lista de vecinos basada en regiones con límites contiguos, que comparten uno o más puntos de límites.

La condición de contiüidad se cumple cuando al menos un punto en el límite de un polígono está dentro de la distancia de ajuste de al menos un punto de su vecino. Esta relación viene dada por el argumento QUEEN = VERDADERO.
Si QUEEN = FALSO, al menos dos puntos límite deben estar dentro de la distancia de ataque entre sí, con el nombre convencional de una relación de “torre”.

knn2nb: (Lista de vecinos del objeto Knn): La función convierte un objeto knn devuelto por knearneigh en una lista de vecinos de clase nb con una lista de vectores enteros que contienen ID de números de regiones vecinas.

knearneig: (K Vecinos más cercanos para pesos espaciales): La función devuelve una matriz con los índices de puntos pertenecientes al conjunto de los k vecinos más cercanos entre sí.

# Construyendo lista de vecinos con criterio tipo reina
X_nb_queen = poly2nb(SPDF, queen=TRUE)  

# Construyendo lista de vecinos con criterio tipo torre
X_nb_rook = poly2nb(SPDF, queen=FALSE) 

# Construyendo lista de vecinos con k-vecinos más cercanos
IDs = row.names(as(SPDF, "data.frame")) # Id de las filas
X_kn4 = knn2nb(knearneigh(coords, k=4), row.names=IDs)

Mostrando resultados

X_nb_queen

## Neighbour list object:
## Number of regions: 666 
## Number of nonzero links: 3612 
## Percentage nonzero weights: 0.8143278 
## Average number of links: 5.423423

X_nb_rook

## Neighbour list object:
## Number of regions: 666 
## Number of nonzero links: 3526 
## Percentage nonzero weights: 0.7949391 
## Average number of links: 5.294294

X_kn4

## Neighbour list object:
## Number of regions: 666 
## Number of nonzero links: 2664 
## Percentage nonzero weights: 0.6006006 
## Average number of links: 4 
## Non-symmetric neighbours list

Graficando

par(mfrow=c(1,3))

plot(SPDF, main = "Criterio de contigüidad\n tipo reina")
plot(X_nb_queen, coords, add=T)

plot(SPDF, main = "Criterio de contigüidad\n tipo torre")
plot(X_nb_rook, coords, add=T)

plot(SPDF, main = "Criterio de contigüidad\n 4 vecinos más cercanos")
plot(X_kn4, coords, add=T)

Construcción de la matriz de ponderaciones espacial (asignar pesos a las áreas que están vinculadas)

Para considerar la estructura de vecindad y los efectos espaciales en la estimación de los parámetros se especifica una matriz cuadrada (W) denominada matriz de pesos o matriz de ponderaciones, cuyos elementos no negativos dan evidencia de la intensidad de la interdependencia existente entre cada par de unidades geográficas i y j.

nb2listw (Pesos espaciales para listas de vecinos): Complementa una lista de vecinos con pesos espaciales para el esquema de codificación elegido.

#   Matriz de pesos normalizados en fila 
# usando QUEEN(W) (Tipo reina):
X_nbq_w = nb2listw(X_nb_queen)
X_nbq_w

## Characteristics of weights list object:
## Neighbour list object:
## Number of regions: 666 
## Number of nonzero links: 3612 
## Percentage nonzero weights: 0.8143278 
## Average number of links: 5.423423 
## 
## Weights style: W 
## Weights constants summary:
##     n     nn  S0       S1       S2
## W 666 443556 666 270.2818 2754.166

#   Matriz de pesos normalizados en fila 
# usando ROOK(W) (Tipo torre):
X_nbr_w = nb2listw(X_nb_rook)
X_nbr_w

## Characteristics of weights list object:
## Neighbour list object:
## Number of regions: 666 
## Number of nonzero links: 3526 
## Percentage nonzero weights: 0.7949391 
## Average number of links: 5.294294 
## 
## Weights style: W 
## Weights constants summary:
##     n     nn  S0       S1       S2
## W 666 443556 666 275.3852 2754.204

#   Matriz de pesos normalizados en fila usando k = 4(W):
X_kn4_w = nb2listw(X_kn4)
X_kn4_w

## Characteristics of weights list object:
## Neighbour list object:
## Number of regions: 666 
## Number of nonzero links: 2664 
## Percentage nonzero weights: 0.6006006 
## Average number of links: 4 
## Non-symmetric neighbours list
## 
## Weights style: W 
## Weights constants summary:
##     n     nn  S0      S1       S2
## W 666 443556 666 305.125 2718.125

Resumen estadístico de la matriz de pesos para cada criterio

summary(unlist(X_nbq_w$weights))  # Criterio reina

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
## 0.09091 0.14286 0.16667 0.18439 0.20000 1.00000

summary(unlist(X_nbr_w$weights))  # Criterio torre

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.1000  0.1429  0.1667  0.1889  0.2000  1.0000

summary(unlist(X_kn4_w$weights))  # Criterio k-vecinos

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    0.25    0.25    0.25    0.25    0.25    0.25

Evaluando existencia de autocorrelación espacial

(a) Índice de Moran:

𝐼 toma valores entre -1 y 1. Valores cercanos a 1 implican dependencia espacial positiva (valores vecinos similares) y valores cercanos a -1 implican dependencia espacial negativa (valores vecinos diferentes). Valores cercanos a cero implican no autocorrelación espacial.

moran.test: (Prueba I de Moran para autocorrelación espacial) Utiliza una matriz de ponderaciones espaciales en forma de lista de ponderaciones. La opción alternativa especifica la hipótesis alternativa, debe ser mayor (predeterminado), menor o como en este caso, bilateral (two.side).

# Primero usamos tipo reina:
Ind_Moran_reina = moran.test(SPDF$POV, listw = X_nbq_w, 
                             alternative = "two.side")

# Luego tipo torre
Ind_Moran_torre = moran.test(SPDF$POV, listw = X_nbr_w, 
                             alternative = "two.side")

# Ahora, usamos k = 4:
Ind_Moran_k = moran.test(SPDF$POV, listw = X_kn4_w, 
                         alternative = "two.side")

Visualizamos resultados:

Ind_Moran_reina

## 
##  Moran I test under randomisation
## 
## data:  SPDF$POV  
## weights: X_nbq_w    
## 
## Moran I statistic standard deviate = 25.742, p-value < 2.2e-16
## alternative hypothesis: two.sided
## sample estimates:
## Moran I statistic       Expectation          Variance 
##      0.6314942529     -0.0015037594      0.0006046621

Ind_Moran_torre

## 
##  Moran I test under randomisation
## 
## data:  SPDF$POV  
## weights: X_nbr_w    
## 
## Moran I statistic standard deviate = 25.543, p-value < 2.2e-16
## alternative hypothesis: two.sided
## sample estimates:
## Moran I statistic       Expectation          Variance 
##      0.6325554088     -0.0015037594      0.0006161701

Ind_Moran_k

## 
##  Moran I test under randomisation
## 
## data:  SPDF$POV  
## weights: X_kn4_w    
## 
## Moran I statistic standard deviate = 22.688, p-value < 2.2e-16
## alternative hypothesis: two.sided
## sample estimates:
## Moran I statistic       Expectation          Variance 
##      0.5915958412     -0.0015037594      0.0006833555

Observamos el índice de Moran para:
* Criterio reina: I = 0.6315
* Criterio torre: I = 0.6326
* k-vecinos: I = 0.5916

Conclusión Por lo que los valores se acercan más a 1, significa la existencia de autocorrelación espacial positiva.

Recordemos La autocorrelación espacial positiva se presenta cuando valores altos [bajos] de la variable de estudio en una localización están asociados con valores altos [bajos] en las regiones vecinas.

De manera que la tasa de pobreza en los diferentes puntos se va a presentar de forma similar y esto afirma la posible existencia de clústeres.

Prueba del estadístico I de Moran

Se calcula un pseudo valor p determinando, la proporción de los valores I de Moran local generados a partir de permutaciones que muestran más agrupaciones que los datos originales. Si esta proporción (el pseudo valor p) es pequeña (menos de 0.05), puede concluir que sus datos muestran agrupaciones estadísticamente significativas.

moran.mc Esta es una prueba de permutación para el estadístico I de Moran calculado usando permutaciones aleatorias, para establecer el rango del estadístico observado en relación con los nsim valores simulados.

moran.mc(SPDF$POV, listw = X_nbq_w, nsim=999) # Reina

## 
##  Monte-Carlo simulation of Moran I
## 
## data:  SPDF$POV 
## weights: X_nbq_w  
## number of simulations + 1: 1000 
## 
## statistic = 0.63149, observed rank = 1000, p-value = 0.001
## alternative hypothesis: greater

moran.mc(SPDF$POV, listw = X_nbr_w, nsim=999) # Torre

## 
##  Monte-Carlo simulation of Moran I
## 
## data:  SPDF$POV 
## weights: X_nbr_w  
## number of simulations + 1: 1000 
## 
## statistic = 0.63256, observed rank = 1000, p-value = 0.001
## alternative hypothesis: greater

moran.mc(SPDF$POV, listw = X_kn4_w, nsim=999) # K-vecinos

## 
##  Monte-Carlo simulation of Moran I
## 
## data:  SPDF$POV 
## weights: X_kn4_w  
## number of simulations + 1: 1000 
## 
## statistic = 0.5916, observed rank = 1000, p-value = 0.001
## alternative hypothesis: greater

Dado que el p-value para cada método o criterio es igual a 0.001, y este valor es menor a 0.05, significa que se muestran agrupaciones estadísticamente significativas.

(b) Prueba C de Geary:

El índice de Geary (C) es un índice de comparaciones por pares entre las diferentes zonas; por lo general varía entre 0 y 2. Teóricamente, un valor de 1 indica ausencia de autocorrelación, es decir, que los valores de una zona no están relacionados con los valores de zonas cercanas. Los valores inferiores a 1 indican autocorrelación espacial positiva, mientras que valores superiores a 1 indican autocorrelación espacial negativa (Griffith, 1987).

geary.test: (Prueba C de Geary para autocorrelación espacial) Utiliza una matriz de ponderaciones espaciales en forma de lista de ponderaciones. La opción alternativa especifica la hipótesis alternativa, debe ser mayor (predeterminado), menor o como en este caso, bilateral (two.side).

# Primero usamos tipo reina:
C_Geary_reina = geary.test(SPDF$POV, listw = X_nbq_w, 
                             alternative = "two.side")

# Luego tipo torre
C_Geary_torre = geary.test(SPDF$POV, listw = X_nbr_w, 
                             alternative = "two.side")

# Ahora, usamos k = 4:
C_Geary_k = geary.test(SPDF$POV, listw = X_kn4_w, 
                         alternative = "two.side")

Visualizamos resultados:

C_Geary_reina

## 
##  Geary C test under randomisation
## 
## data:  SPDF$POV 
## weights: X_nbq_w 
## 
## Geary C statistic standard deviate = 23.5, p-value < 2.2e-16
## alternative hypothesis: two.sided
## sample estimates:
## Geary C statistic       Expectation          Variance 
##      0.3794742478      1.0000000000      0.0006972438

C_Geary_torre

## 
##  Geary C test under randomisation
## 
## data:  SPDF$POV 
## weights: X_nbr_w 
## 
## Geary C statistic standard deviate = 23.393, p-value < 2.2e-16
## alternative hypothesis: two.sided
## sample estimates:
## Geary C statistic       Expectation          Variance 
##      0.3772234532      1.0000000000      0.0007087578

C_Geary_k

## 
##  Geary C test under randomisation
## 
## data:  SPDF$POV 
## weights: X_kn4_w 
## 
## Geary C statistic standard deviate = 21.9, p-value < 2.2e-16
## alternative hypothesis: two.sided
## sample estimates:
## Geary C statistic       Expectation          Variance 
##      0.4050635405      1.0000000000      0.0007379694

Observamos el C de Geary para:
* Criterio reina: C =0.3795
* Criterio torre: C = 0.3772
* k-vecinos: C = 0.4074

Conclusión Por lo que los valores son menores que 1, confirma la existencia de autocorrelación espacial positiva.

Prueba del estadístico C de geary

geary.mc Esta es una prueba de permutación para el estadístico C de Geary calculado usando permutaciones aleatorias.

geary.mc(SPDF$POV, listw = X_nbq_w, nsim=999) # Reina

## 
##  Monte-Carlo simulation of Geary C
## 
## data:  SPDF$POV 
## weights: X_nbq_w 
## number of simulations + 1: 1000 
## 
## statistic = 0.37947, observed rank = 1, p-value = 0.001
## alternative hypothesis: greater

geary.mc(SPDF$POV, listw = X_nbr_w, nsim=999) # Torre

## 
##  Monte-Carlo simulation of Geary C
## 
## data:  SPDF$POV 
## weights: X_nbr_w 
## number of simulations + 1: 1000 
## 
## statistic = 0.37722, observed rank = 1, p-value = 0.001
## alternative hypothesis: greater

geary.mc(SPDF$POV, listw = X_kn4_w, nsim=999) # K-vecinos

## 
##  Monte-Carlo simulation of Geary C
## 
## data:  SPDF$POV 
## weights: X_kn4_w 
## number of simulations + 1: 1000 
## 
## statistic = 0.40506, observed rank = 1, p-value = 0.001
## alternative hypothesis: greater

Dado que el p-value para cada método o criterio es igual a 0.001, y este valor es menor a 0.05, significa que se muestran agrupaciones estadísticamente significativas.

(c) Prueba G(d) de Getis y Ord

globalG.test: La estadística G global para la autocorrelación espacial, que complementa las medidas locales de Gi LISA: localG..

# Ahora, usando k = 4:
G_Getis_k = globalG.test(SPDF$POV, listw = X_kn4_w, 
                         alternative = "two.side") 
G_Getis_k

## 
##  Getis-Ord global G statistic
## 
## data:  SPDF$POV 
## weights: X_kn4_w 
## 
## standard deviate = 15.199, p-value < 2.2e-16
## alternative hypothesis: two.sided
## sample estimates:
## Global G statistic        Expectation           Variance 
##       1.629446e-03       1.503759e-03       6.837882e-11

Representación gráfica de la dependencia espacial

Un gráfico de dispersión que permite analizar la autocorrelación espacial global, es el denominado Scatterplot de Moran univariante (Anselin, 1993). Este gráfico se construye en un plano cartesiano, en el eje de abscisas se ubican los valores de la variable de estudio estandarizada (POV), y en el eje de ordenadas, se ubican los valores del retardo espacial de la variable estandarizada. Entiéndase por retardo espacial de una unidad espacial i, al valor medio de todos los valores de la variable estandarizada correspondientes a las unidades espaciales vecinas de i.

Si la nube de puntos tiende a formarse sobre la diagonal principal del plano cartesiano (cuadrante I y III, en contra de las manecillas de reloj), esto nos indicará la presencia de autocorrelación espacial positiva.
Si la nube de puntos tiende a ubicarse sobre la diagonal secundaria del plano cartesiano (cuadrante II y IV, en contra de las manecillas de reloj) , será un indicio de existencia de autocorrelación espacial negativa.
De lo contrario, si la nube de puntos de distribuye sobre los cuatro cuadrantes, esto nos indicará ausencia de autocorrelación espacial.

moran.plot: (Diagrama de dispersión de Moran) Una gráfica de datos espaciales frente a sus valores espaciales rezagados, aumentada al informar el resumen de las medidas de influencia para la relación lineal entre los datos y el rezago.

scatterplot_moran = moran.plot(SPDF$POV, listw = X_kn4_w, main = "Scatterplot de Moran para la pobreza")

Dado que la nube de puntos tiende a formarse sobre la diagonal principal del plano cartesiano (cuadrante I y III), esto indica la presencia de autocorrelación espacial positiva.

Prueba de presencia de autocorrelación espacial (local):

Anselin (1995) se refirió a esto como una estadística “LISA”, para Indicador local de autocorrelación espacial.

localmoran: La estadística local espacial I de Moran se calcula para cada zona en función del objeto de ponderaciones espaciales.

lag.listw: Usando una representación dispersa de listw de una matriz de ponderaciones espaciales, calcula el vector de retardo.

LISA <- function(criterio) {
  
  P = localmoran(SPDF$POV, listw = criterio)
  #head(P) # Se muestran los primeros 6 datos.
  #dim(P) # La dimension de P
  cPOVE = SPDF$POV - mean(SPDF$POV) # Pobreza - Media(Pobreza)
  lagDV = lag.listw(X_nbq_w, SPDF$POV) # Calcula el retardo (promedios)
  clagDV = lagDV - mean(lagDV) # Retardo - Media(Retardo)
  p = P[,5] # Toma la columna: Pr(z > 0) de P

  # Se inicializa vector numerico de longitud filas de P (666)
  quadrant = vector(mode="numeric",length=nrow(P))+5
  quadrant[cPOVE>0 & clagDV>0 & p<= 0.05] = 1 # Alto-Alto
  quadrant[cPOVE<0 & clagDV<0 & p<= 0.05] = 2 # Bajo-Bajo
  quadrant[cPOVE<0 & clagDV>0 & p<= 0.05] = 3 # Bajo-Alto
  quadrant[cPOVE>0 & clagDV<0 & p<= 0.05] = 4 # Alto-Bajo

  # Grafico  
  brks = c(1,2,3,4,5)
  colors = c("red", "blue", "light blue", "pink", "white")
  plot(Condados, border ="lightgray", col=colors[findInterval(quadrant,brks,all.inside=FALSE)])
  legend("bottomright", legend = c("High-High", "Low-Low", "Low-High", "High-Low", "Insignificant"), fill = colors, bty="n", cex=0.7, y.intersp=1, x.intersp=1)
  box()
  title("LISA Cluster Map")
}

# Usando criterio reina
LISA(X_nbq_w)

# Usando K-vecinos, k=4
LISA(X_kn4_w)