7 Pruebas de autocorrelación

Tomado de Kopczewska (2020), sección 4.5.

Librerías usadas:

  • rgdal
  • spdep
  • gstat

7.1 Introducción

Los métodos de estadística espacial también se denominan análisis exploratorio de datos espaciales (ESDA). Las estadística espacial comprende métodos para probar la existencia de procesos de autocorrelación espacial. Las medidas de autocorrelación espacial se utilizan para evaluar la correlación de variables con respecto a la ubicación espacial.

La autocorrelación espacial significa que las observaciones geográficamente cercanas son más similares entre sí que las distantes. Esto hace posible predecir valores en un punto dado, conociendo los valores en otras ubicaciones.

La medida de autocorrelación espacial se puede utilizar de varias formas: para probar la existencia de autocorrelación espacial y características de la población; para determinar el grado de autocorrelación, por ejemplo, la distancia por encima de la cual las observaciones son independientes; determinar el modelo teórico apropiado para la estructura espacial observada o concluir sobre el proceso espacial.

Con base en las conclusiones de la ESDA, se puede rechazar la hipótesis de la aleatoriedad espacial, lo que abre el camino a la búsqueda de conglomerados espaciales (Anselin, 1999). Las conclusiones de la ESDA también son la base para la especificación espacial de los modelos.

En estadística espacial se utilizan dos tipos de medidas: medidas globales y locales.

Por un lado, las medidas globales se plasman en un indicador de autocorrelación espacial o similitud general de regiones. Su ventaja es que es un indicador sintética y su desventaja es promediar.

Por otro lado, gracias a las estadísticas locales, que se determinan para cada región, se puede responder a la pregunta de si la región está rodeada de regiones con valores altos o bajos o si es similar/diferente a las regiones vecinas. Los valores de las estadísticas locales no se generalizan para toda el área analizada y se trabaja con información sobre la posición de cada región en relación con los vecinos.

El uso de técnicas ESDA es el primer paso en el análisis de datos espaciales y permite la detección de patrones de autocorrelación espacial globales y locales en los datos.

La autocorrelación global se puede medir mediante las estadísticas I de Moran global (Moran, 1950; Cliff & Ord, 1981) y C de Geary.

Las medidas locales más utilizadas incluyen las estadísticas locales Moran’s II y las estadísticas locales Geary’s Ci, que forman parte del llamado indicador local de asociación espacial (LISA) (Anselin, 1995) y las estadísticas locales G Getis-Ord (Ord & Getis, 1995).

Las estadísticas globales pueden identificar conglomerados y relaciones espaciales solo para todo el sistema. Sin embargo, estas estadísticas se pueden desagregar en estadísticas locales, gracias a las cuales se pueden detectar patrones de relaciones espaciales locales entre las regiones y sus vecinos.

El cálculo de LISA permite detectar cuál de las regiones influye fuertemente en la formación de conglomerados locales y examinar los conglomerados locales no aleatorios, los denominados puntos calientes/fríos (hot spots o cold spots), inestabilidad local (heterocedasticidad), observaciones atípicas significativas (Ministeri, 2003). La siguiente tabla resume algunas de las características listadas:

Estadísticos espaciales globales y locales. Fuente: Kopczewska (2020). (\(\alpha\): nivel de significancia).
Estadístico Tipo Tipo de medición Interpretación Funciones en R
Moran’s I Global Una medida de autocorrelación espacial \(I> 0\): autocorrelación espacial positiva, los valores de observación a la distancia \(d\) son similares
\(I <0\): autocorrelación espacial negativa, los valores de observación a la distancia \(d\) son diferentes
\(I = 0\): los valores de observación a la distancia d se distribuyen aleatoriamente
moran ()
moran.test ()
moran.mc ()
moran.plot ()
Geary’s C Global Una medida de autocorrelación espacial \(0 <C <1\): autocorrelación espacial positiva, los valores de observación a la distancia \(d\) son similares
\(1 <C <2\): autocorrelación espacial negativa, los valores de observación a la distancia \(d\) son diferentes
geary ()
geary.test ()
geary.mc ()
Join-count Global en grupos Una medida de autocorrelación espacial \(H_0\) asume que no hay autocorrelación espacial, el valor \(p>\alpha\) confirma esta hipótesis
\(H_1\) asume la existencia de autocorrelación espacial, valor \(p <\alpha\) confirma esta hipótesis
joincount.test ()
joincount.mc ()
joincount.multi ()
Local Moran’s \(I_i\) Local Junto con Local Geary, crea LISA, una medida de relaciones espaciales. valor \(p <\alpha\): rodeado por valores relativamente altos, un grupo de valores altos
valor \(p> (1-\alpha)\): un área rodeada por valores relativamente bajos, un grupo de valores bajos
local.moran ()
localmoran.sad ()
Local Geary’s \(C_i\) Local Junto con Local Moran, crea LISA, una medida de diversidad espacial valor \(p> (1-\alpha)\): positivo, región rodeada de observaciones similares valor \(p <\alpha\): negativo, región rodeada por diferentes observaciones Existe el comando usdm::lisa() que permite al usuario calcular este estadístico, pero para rásteres.
Local \(G_i\), \(G^*\) Local Medida de relaciones espaciales \(G_i> 0\): área rodeada por valores relativamente altos, clúster de alto valor
\(G_i <0\): área rodeada por valores relativamente bajos, clúster de valor bajo
localG ()
Local \(H_i\), LOSH Local Medida de relaciones espaciales Interpretado junto con las estadísticas \(G^*\). Valores altos indican un entorno heterogéneo; valores bajos indican homogeneidad local. LOSH()
Empirical Bayes index Local Medida de relaciones espaciales Interpretado como las estadísticas I de Moran, utilizadas por interés (por ejemplo, número de casos / población) EBImoran.mc()

Supuestos:

La estadística espacial se basa en varios supuestos. El más importante de ellos es la estacionariedad, lo que significa que los datos analizados deben tener una distribución normal con promedio y varianza constantes.

Sin embargo, las estructuras de autocorrelación espacial no siempre son constantes en toda el área de estudio, por lo que es suficiente satisfacer una estacionariedad débil, donde la media y la varianza son constantes, mientras que la autocorrelación depende solo de la distancia entre los objetos probados.

Los patrones espaciales deben ser isotrópicos; es decir, no deben depender de la dirección del estudio. Además, la forma del área estudiada afecta la intensidad, alcance y tipo de proceso espacial estimado, así como su importancia.

También existen los llamados efectos de borde (Fortin, Dale y Ver Hoef, 2002). Este problema se revela durante las observaciones en el borde del área de estudio. Estas regiones tienen menos vecinos que los objetos ubicados en el medio. Como resultado de este efecto, pueden aparecer diferencias en las estimaciones dependiendo de la matriz de vecindad adoptada (por ejemplo, de primer o segundo orden) (Cressie, 1993b).

Un elemento común en econometría espacial es el llamado retardo/rezago espacial, que se determina para la variable examinada. El operador de rezago espacial es un promedio ponderado del valor de la variable en las regiones vecinas, de acuerdo con la matriz de ponderaciones espaciales declarada \(W\).

Cuando se utiliza la matriz de vecindad de primer orden para el cálculo según el criterio de contigüidad, el rezago espacial será el promedio de los valores en las regiones limítrofes con el área examinada. Si se adoptara la matriz de segundo orden, el rezago espacial se determinaría a partir de los valores de los vecinos de estos vecinos. Una interpretación análoga ocurre cuando se usa la matriz de distancia o la matriz \(k\) de los vecinos más cercanos.

En las siguientes secciones usaremos datos de Valor agregado anual por cantones del Ecuador tomados del Banco Central del Ecuador y el indicador de Necesidades Básicas Insatisfechas por personas.

# leo los datos
uu <- "https://raw.githubusercontent.com/vmoprojs/DataLectures/master/SpatialData/VABNoPetroleroCantones2007-2019.csv"
datos <- read.csv(uu,sep = ",")
datos$COD_CANT[nchar(datos$COD_CANT)==3] = paste("0",datos$COD_CANT[nchar(datos$COD_CANT)==3],sep = "")

datos <- subset(datos,datos$YEAR == 2019)

uu <- "https://raw.githubusercontent.com/vmoprojs/DataLectures/master/SpatialData/proyeccion_cantonal_total_2010-2020.csv"

poblacion <- read.csv(uu,sep = ";")

poblacion$CODIGO[nchar(poblacion$CODIGO)==3] = paste("0",poblacion$CODIGO[nchar(poblacion$CODIGO)==3],sep = "")

uu <- "https://raw.githubusercontent.com/vmoprojs/DataLectures/master/SpatialData/NBI_PER_CANT.csv"
NBI <- read.csv(uu, sep = ";")

NBI$CODIGO[nchar(NBI$CODIGO)==3] = paste("0",NBI$CODIGO[nchar(NBI$CODIGO)==3],sep = "")

datos <- merge(datos,poblacion[,c("CODIGO","A_2019")],by.x = "COD_CANT",by.y = "CODIGO",all.x = TRUE)

datos <- merge(datos,NBI[,c("CODIGO","POBRES_P")],by.x = "COD_CANT",by.y = "CODIGO",all.x = TRUE)

datos$VAB_PC <- datos$VAB/datos$A_2019 #miles de USD / pob



library(rgdal)
uu <- "https://github.com/vmoprojs/DataLectures/raw/master/SpatialData/2012_nxcantones.zip"
poligonos <- read_git_shp(uu)
## OGR data source with driver: ESRI Shapefile 
## Source: "/private/var/folders/0p/n_r_hl095sv7nktfp_8n9n_80000gn/T/RtmpeHuEkZ/file26b25979a69/nxcantones.shp", layer: "nxcantones"
## with 224 features
## It has 6 fields
poligonos <- poligonos[poligonos$DPA_PROVIN!="90",]
poligonos <- poligonos[poligonos$DPA_PROVIN!="20",]
poligonos@data = merge(poligonos@data,datos,by.x  ="DPA_CANTON",by.y = "COD_CANT",all.x=TRUE)



# Me quedo con los datos que no son perdidos:
poligonos <- poligonos[!is.na(poligonos@data$POBRES_P),]

library(spdep)
set.ZeroPolicyOption(TRUE)
## [1] FALSE
# Se construye la lista de vecinos
cont.nb <- poly2nb(poligonos)
# Se construye la matriz de pesos espaciales
cont.listw <- nb2listw(cont.nb, style="W",zero.policy = TRUE)

7.2 Estadísticos Globales

Las medidas globales en estadística espacial permiten determinar la tendencia espacial general en la región. Su análisis permite una visión sintética de los datos. También proporcionan una buena base para una evaluación posterior más detallada del fenómeno.

7.2.1 Estadístico global: \(I\) de Moran

El estadístico global \(I\) de Moran es probablemente la medida espacial más antigua, introducida por Moran en 1950 (Moran (1950)) y se usa para realizar pruebas de hipótesis de autocorrelación global.

Es una prueba basada en la covarianza, similar a una prueba de Pearson y el estadístico tiene la siguiente fórmula:

\[ I = \frac{\sum_i\sum_jw_{ij}z_iz_j/S_0}{\sum_iz_i^2/n} \]

donde \(z_i = x_i-\bar{x}\), \(S_0=\sum_i\sum_jw_{ij}\). \(x_i\) es la observación en la unidad espacial \(i\), \(\bar{x}\) es el promedio de todas las unidades espaciales, \(n\) es el total de unidades espaciales y \(w_{ij}\) es un elemento de la matriz de pesos espaciales \(W\). La matriz de pesos espaciales \(W\) se suele estandarizar por fila (las filas suman \(1\)).

7.2.1.1 Inferencia

La significancia de estadístico \(I\) de Moran puede evaluarse desde un enfoque paramétrico o no paramétrico; en este texto se adopta el enfoque no paramétrico. Para este objetivo se obtiene una distribución de la permutación del estadístico basado en una simulación de Monte Carlo de 10000 permutaciones.

La hipótesis nula asumida es la ausencia de autocorrelación, es decir, una distribución aleatoria de los valores.

moran(poligonos$POBRES_P, cont.listw, length(cont.nb), Szero(cont.listw))
## $I
## [1] 0.08994426
## 
## $K
## [1] 3.410684
MIran <- moran.test(poligonos$POBRES_P, cont.listw,zero.policy = TRUE)
MIran
## 
##  Moran I test under randomisation
## 
## data:  poligonos$POBRES_P  
## weights: cont.listw  
## n reduced by no-neighbour observations  
## 
## Moran I statistic standard deviate = 1.7443, p-value = 0.04056
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       0.084924120      -0.004950495       0.002654920
(MIran$estimate[1]-MIran$estimate[2])/sqrt(MIran$estimate[3])
## Moran I statistic 
##           1.74426
MInorm <- moran.test(poligonos$POBRES_P, cont.listw, randomisation=FALSE)
MInorm
## 
##  Moran I test under normality
## 
## data:  poligonos$POBRES_P  
## weights: cont.listw  
## n reduced by no-neighbour observations  
## 
## Moran I statistic standard deviate = 1.7423, p-value = 0.04073
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       0.084924120      -0.004950495       0.002660820
(MInorm$estimate[1]-MInorm$estimate[2])/sqrt(MInorm$estimate[3])
## Moran I statistic 
##          1.742326
#pvalue para Moran's I (distribucion normal)
pval.norm <- 1-pnorm(MInorm$statistic, mean=0, sd=1)
pval.norm
## Moran I statistic standard deviate 
##                         0.04072574
MImc <- moran.mc(poligonos$POBRES_P, cont.listw, 999)
MImc
## 
##  Monte-Carlo simulation of Moran I
## 
## data:  poligonos$POBRES_P 
## weights: cont.listw  
## number of simulations + 1: 1000 
## 
## statistic = 0.084924, observed rank = 956, p-value = 0.044
## alternative hypothesis: greater
(Zmi <- (MImc$statistic-mean(MImc$res))/sd(MImc$res))
## statistic 
##    1.8169
Sol <- t(rbind(c(MIran$estimate,MIran$statistic,MIran$p.value),
c(MInorm$estimate,MInorm$statistic,MInorm$p.value),
c(MImc$statistic,mean(MImc$res),var(MImc$res),Zmi,MImc$p.value)))
rownames(Sol) <- c("MI","E[MI]","V[MI]","z-value","p-value")
colnames(Sol) <- c("Randomization","Normal","Monte Carlo")
Sol
##         Randomization       Normal  Monte Carlo
## MI        0.084924120  0.084924120  0.084924120
## E[MI]    -0.004950495 -0.004950495 -0.003218683
## V[MI]     0.002654920  0.002660820  0.002353486
## z-value   1.744260255  1.742325632  1.816900233
## p-value   0.040556859  0.040725740  0.044000000

7.2.2 Interpretación de MI

  • Su media teórica es \(\bar{I} = \frac{-1}{(n-1)}\) y tiende a \(0\).
  • Signo positivo y significativo: conglomerados de valores similares. Es decir, puede ser valores altos, bajos, o ambos.
  • Signo negativo y significativo: valores que alternan. Puede deberse a ouliers espaciales, puede deberse a heterogeneidad espacial (un patrón parecido a un tablero de ajedrez).

7.2.3 Comparando MI’s

  • El estadístico de moran depende de la matriz de pesos espaciales \(W\).
  • Se puede comparar variables con la misma matriz de pesos espaciales, pero no se puede comprar variables (incluso la misma variable) con distintas matrices de pesos espaciales.
  • Para comparabilidad, incluso cuando usas diferentes pesos espaciales, puedes usar valores estandarizados como \(z-value\) del ejemplo anterior.

Recuerda que el estadístico de Moran es global, es decir que da cuenta de un patrón general en los datos, pero no da cuenta de la ubicación de los patrones. Por ejemplo, si es positivo indica que existe autocorrelación y que las unidades espaciales se agrupan en en valores similares, pero no indica dónde están esos grupos.

La existencia de correlación positiva indica la existencia de conglomerados, pero no el por qué se dan estos conglomerados. El por qué puede deberse a dos motivos:

  • contagio verdadero: evidencia de conglomerados debido a interacción espacial. Por ejemplo, si una persona contagia a otras de gripe en un lugar cerrado.
  • contagio aparente: evidencia de conglomerados debido a heterogeneidad espacial. Por ejemplo, si un grupo de personas se reúne en un lugar cerrado que ya estaba contaminado.

7.2.4 Gráfico de Moran

El gráfico de Moran nos permite empezar a conectar la idea de estadísiticos globales y locales. También es una guía para reconocer el tipo de autocorrelación con la que estamos tratando.

variable <- poligonos$POBRES_P
variable.std <- ((variable-mean(variable))/sd(variable))
moran.plot(variable.std, cont.listw, labels=as.character(poligonos$DPA_CANTON),pch=19,quiet=F) # 

## Potentially influential observations of
##   lm(formula = wx ~ x) :
## 
##      dfb.1_ dfb.x dffit   cov.r   cook.d hat    
## 0101 -0.06   0.16 -0.17    1.05_*  0.01   0.04_*
## 0601  0.04  -0.10  0.10    1.03_*  0.01   0.03  
## 0802 -0.21  -0.27 -0.34_*  0.94_*  0.06   0.01  
## 1001  0.08  -0.21  0.22    1.04_*  0.02   0.04_*
## 1101  0.08  -0.20  0.22    1.03_*  0.02   0.03_*
## 1314  0.16   0.08  0.17    0.97_*  0.01   0.01  
## 1410  0.16   0.13  0.21    0.97_*  0.02   0.01  
## 1602 -0.16   0.16 -0.23    0.97_*  0.02   0.01  
## 1603  0.16   0.09  0.19    0.96_*  0.02   0.01  
## 1701  0.07  -0.23  0.24    1.06_*  0.03   0.06_*
## 1705  0.00   0.01 -0.01    1.07_*  0.00   0.06_*
## 1708 -0.19   0.05 -0.19    0.95_*  0.02   0.00  
## 1802  0.05  -0.12  0.14    1.04_*  0.01   0.03_*
## 1805 -0.18   0.03 -0.18    0.95_*  0.02   0.00  
## 2401  0.15   0.05  0.16    0.97_*  0.01   0.01

Se puede tipificar los cuadrantes:

  • Cuadrante I & III:

  • Autocorrelación positiva.

  • Conglomerados de valores similares.

  • Unidades espaciales que son similares a sus vecinos.

  • Cuadrante II & IV:

  • Autocorrelación negativa.

  • outliers espaciales.

  • Unidades espaciales que son diferentes a sus vecinos.

Otra forma de entender los cuadrantes es mediante el tipo de autocorrelación:

  • Cuadrante I: High-High, unidades espaciales que están sobre su media, rodeadas de vecinas que están sobre su media.
  • Cuadrante II: Low-High, unidades espaciales que están debajo de su media, rodeadas de vecinas que están debajo de su media. Este cuadrante marca unidades espaciales candidatas a outliers espaciales.
  • Cuadrante III: Low-Low, unidades espaciales que están debajo su media, rodeadas de vecinas que están debajo de su media.
  • Cuadrante IV: High-Low, unidades espaciales que están sobre su media, rodeadas de vecinas que están debajo de su media. Este cuadrante marca unidades espaciales candidatas a outliers espaciales.

¿Cuáles de estas señales de los cuadrantes son significativas? Esto se aborda con estadísticos locales.

Veamos el gráfico de Moran como una regresión lineal:

x <- poligonos$POBRES_P# extraemos la variable
zx <- scale(x) # se estandariza la variable
mean(zx) # controlamos que el promedio sea correcto
## [1] -4.67039e-17
sd(zx) # controlamos que la desviación sea correcta
## [1] 1
wzx <- lag.listw(cont.listw, zx) # rezago espacial de x
morlm <- lm(wzx~zx) # regression
slope <- morlm$coefficients[2] # pendiente
intercept <- morlm$coefficients[1] # término constante
par(pty="s") # ventana cuadrada
plot(zx, wzx, xlab="zx",ylab="Rezago espacial de zx", pch="*")
abline(intercept, slope) # linea de regresion
abline(h=0, lty=2) # linea horizontal en y = 0
abline(v=0, lty=2) # linea vertical en x = 0

Graficamos las unidades espaciales según el cuadrante al que pertenecen en el gráfico de Moran:

# map of belonging to the quarters of the Moran scatterplot
# creating a variable for analysis
x <- poligonos$POBRES_P# creating a variable for analysis
zx <- scale(x) # variable standardisation
wzx <- lag.listw(cont.listw, zx) # spatial lag of x
cond1 <- ifelse(zx>=0 & wzx>=0, 1,0) # I quarter
cond2 <- ifelse(zx>=0 & wzx<0, 2,0) # II quarter
cond3 <- ifelse(zx<0 & wzx<0, 3,0) # III quarter
cond4 <- ifelse(zx<0 & wzx>=0, 4,0) # IV quarter
cond.all <- as.data.frame(cond1+cond2+cond3+cond4)
# chart - colour map, Figure 10b
brks <- c(1,2,3,4)
cols <- c("royalblue", "springgreen", "yellow", "red")

par(mar=c(5.1,1,4.1,1))
sp::plot(poligonos, col=cols[findInterval(cond.all$V1, brks)])
legend("bottomleft", legend=c("I - HH - high rodeado por high",
                              "II - LH - low rodeado por high", "III - LL - low
rodeado por low", "IV exercise - HL - high rodeado por low"),
       fill=cols, bty="n", cex=0.8)
title(main="Regiones según cuadrantes del gráfico de Moran")

7.2.5 Estadístico Global de Geary

Propuesto por Geary en 1954, se basa en la noción de disimilaridad:

\[ C = \frac{(n-1)}{2S_0}\frac{\sum_i\sum_jw_{ij}(x_i-x_j)^2}{\sum_iz_i^2} \]

donde \(z_i = x_i-\bar{x}\), \(S_0=\sum_i\sum_jw_{ij}\). \(x_i\) es la observación en la unidad espacial \(i\), \(x_j\) es la observación en la unidad espacial \(j\), \(\bar{x}\) es el promedio de todas las unidades espaciales, \(n\) es el total de unidades espaciales y \(w_{ij}\) es un elemento de la matriz de pesos espaciales \(W\).

  • El valor de \(C\) está comprendido entre \(0\) y \(2\), asintóticamente sigue una distribución normal.
  • Asumiendo ausencia de autocorrelación espacial, el valor esperado de \(C\) es \(1\).
  • Si \(0< C\leq 1\) indica regiones similares (autocorrelación positiva). Esto es equivalente a que su valor \(z\) sea negativo, \(z<0\).
  • Si \(1< C\leq 2\) indica regiones con valores alternantes (autocorrelación negativa). Esto es equivalente a que su valor \(z\) sea positivo, \(z>0\).
  • El estadístico \(C\) asume una matriz de pesos espaciales simétrica, en el caso de asimetría (como en el caso de los vecinos más cercanos), el estadísitico puede no ser de utilidad.

7.2.5.1 Inferencia

La significancia de estadístico \(C\) de Geary puede evaluarse desde un enfoque paramétrico o no paramétrico; en este texto se adopta principalemente el enfoque no paramétrico. Para este objetivo se obtiene una distribución de la permutación del estadístico basado en una simulación de Monte Carlo de 10000 permutaciones.

La hipótesis nula asumida es la ausencia de autocorrelación, es decir, una distribución aleatoria de los valores.

geary(poligonos$POBRES_P, cont.listw, length(cont.nb), length(cont.nb)-1,Szero(cont.listw))
## $C
## [1] 0.9040488
## 
## $K
## [1] 3.410684
GCran <- geary.test(poligonos$POBRES_P, cont.listw)
GCran
## 
##  Geary C test under randomisation
## 
## data:  poligonos$POBRES_P 
## weights: cont.listw  
## n reduced by no-neighbour observations 
## 
## Geary C statistic standard deviate = 2.6134, p-value = 0.004482
## alternative hypothesis: Expectation greater than statistic
## sample estimates:
## Geary C statistic       Expectation          Variance 
##       0.853354493       1.000000000       0.003148575
GCnorm <- geary.test(poligonos$POBRES_P, cont.listw,randomisation = FALSE)
GCnorm
## 
##  Geary C test under normality
## 
## data:  poligonos$POBRES_P 
## weights: cont.listw  
## n reduced by no-neighbour observations 
## 
## Geary C statistic standard deviate = 2.6389, p-value = 0.004158
## alternative hypothesis: Expectation greater than statistic
## sample estimates:
## Geary C statistic       Expectation          Variance 
##       0.853354493       1.000000000       0.003088003
GCmc <- geary.mc(poligonos$POBRES_P, cont.listw, 999)
GCmc
## 
##  Monte-Carlo simulation of Geary C
## 
## data:  poligonos$POBRES_P 
## weights: cont.listw  
## number of simulations + 1: 1000 
## 
## statistic = 0.85335, observed rank = 69, p-value = 0.069
## alternative hypothesis: greater
(Zgc <- (GCmc$statistic-mean(GCmc$res))/sd(GCmc$res))
## statistic 
## -1.500496

Un resumen de los resultados

Sol <- t(rbind(c(GCran$estimate,GCran$statistic*-1,GCran$p.value),
c(GCnorm$estimate,GCnorm$statistic*-1,GCnorm$p.value),
c(GCmc$statistic,mean(GCmc$res),var(GCmc$res),Zgc,GCmc$p.value)))
rownames(Sol) <- c("C","E[C]","V[C]","z-value","p-value")
colnames(Sol) <- c("Randomization","Normal","Monte Carlo")
Sol
##         Randomization       Normal  Monte Carlo
## C         0.853354493  0.853354493  0.853354493
## E[C]      1.000000000  1.000000000  0.942909025
## V[C]      0.003148575  0.003088003  0.003562097
## z-value  -2.613435234 -2.638942541 -1.500495560
## p-value   0.004481854  0.004158254  0.069000000

7.2.6 Variograma

El semi-variogama es un diagrama bi-dimensional donde en el eje \(x\) se muestra la distancia y en el eje \(y\) se tiene:

\[\begin{equation} \hat\gamma(\mathbf{h}) = \frac{1}{2\#N(\mathbf{h})}\sum(Z(s_i+\mathbf{h})-Z(s_i))^2 \end{equation}\]

# Variograma
library(gstat)
pr.v <- variogram(POBRES_P~ 1,  poligonos)
plot(pr.v)

7.3 Estadísticos Locales

Anselin (1995) propuso indicadores locales de relaciones espaciales (local indicators of spatial association - LISA). LISA incluye estadísticos locales: \(I_i\) de Moran y \(C_i\) de Geary.

El estadístico de Moran local permite la identificación de efectos de aglomeración espacial (similar al estadístico \(G\) Getis-Ord), mientras que el Geary local muestra similitudes y diferencias espaciales.

Las estadísticas locales de Moran brindan información sobre grupos de valor bajo o alto, mientras que los estadísticos locales de Geary muestran diferencias promedio entre el objeto y los vecinos, lo que ayuda a encontrar valores atípicos y patrones de similitud/diferencia.

Estos estadísticos están relacionadas con sus contrapartes globales y se pueden utilizar para estimar el impacto de las estadísticas individuales en sus contrapartes globales.

Gracias a LISA, es posible identificar los denominados puntos calientes, es decir, centros con valores altos rodeados de valores bajos, así como clusters locales en ausencia de autocorrelación global. Las islas de alto valor pueden interpretarse no solo como puntos calientes, sino también como valores atípicos. Entonces, los estadísticos locales son un indicador de inestabilidad local y desviaciones locales del patrón de autocorrelación global.

En resumen, los LISA nos presentan dos componentes:

  • un estadístico para cada unidad espacial que puede ser sujeto a pruebas de significancia.

  • cuantifica la relación entre el estadístico global y el local: la suma de los estadísticos locales es proporcional al estadístico global.

7.3.1 Estadístico local: \(I\) de Moran

\[ I_i = \frac{(x_i-\bar{x})\sum_{j=1}^{n}w_{ij}(x_j-\bar{x})}{\sum_{i=1}^{n}(x_i-\bar{x})^2/n} \] Para evaluar la significancia se evalúa los valores

  • \(p-\text{value}<\alpha\): las unidades espaciales que cumplen con esta condición indican tienen un índice local de Moran positivo y significativo. Esto significa que éstas unidades están rodeadas de valores altos de la variable analizada.

  • \(p-\text{value}>1-\alpha\): las unidades espaciales que cumplen con esta condición indican tienen un índice local de Moran negativo y significativo. Esto significa que éstas unidades están rodeadas de valores bajos de la variable analizada.

locM <- localmoran(poligonos$POBRES_P, cont.listw,alternative = "greater")
# locM <- localmoran_perm(poligonos$POBRES_P, cont.listw)
oid1 <- order(poligonos$DPA_CANTON)
locMorMat <- printCoefmat(data.frame(locM[oid1,], row.names=poligonos$DPA_CANTON[oid1]), check.names=FALSE)
##               Ii        E.Ii      Var.Ii        Z.Ii Pr.z...E.Ii..
## 0101  1.8138e+00 -3.6921e-02  1.0614e+00  1.7964e+00        0.0362
## 0102  1.7414e-01 -4.3472e-03  1.8262e-01  4.1767e-01        0.3381
## 0103  4.1489e-01 -1.3237e-03  7.0056e-02  1.5725e+00        0.0579
## 0104 -1.4863e-01 -3.1669e-03  3.3777e-01 -2.5030e-01        0.5988
## 0105  4.5277e-01 -1.0482e-03  3.6641e-02  2.3708e+00        0.0089
## 0106  8.9486e-01 -5.5214e-03  1.1805e+00  8.2868e-01        0.2036
## 0107 -1.6218e+00 -9.7991e-03  1.0382e+00 -1.5821e+00        0.9432
## 0108  1.3674e-01 -2.5409e-03  8.8685e-02  4.6769e-01        0.3200
## 0109  1.0303e-01 -2.9613e-04  1.0359e-02  1.0152e+00        0.1550
## 0110 -2.6998e-01 -1.0390e-03  5.5005e-02 -1.1467e+00        0.8742
## 0111 -8.5020e-03 -2.5916e-03  1.0906e-01 -1.7897e-02        0.5071
## 0112 -2.3728e-01 -5.2482e-03  2.2027e-01 -4.9438e-01        0.6895
## 0113  3.4386e-02 -8.0482e-04  5.7092e-02  1.4728e-01        0.4415
## 0114 -1.7998e-01 -4.7526e-03  1.4119e-01 -4.6636e-01        0.6795
## 0115 -2.3526e-03 -6.7624e-07  4.8009e-05 -3.3944e-01        0.6329
## 0201 -7.1932e-02 -4.4763e-05  3.1777e-03 -1.2752e+00        0.8989
## 0202 -5.0106e-01 -1.7852e-03  4.6317e-02 -2.3199e+00        0.9898
## 0203  1.0876e-01 -1.7806e-04  7.5113e-03  1.2569e+00        0.1044
## 0204  1.2012e-01 -2.2469e-03  1.5916e-01  3.0672e-01        0.3795
## 0205  3.3598e-02 -7.8076e-05  2.7318e-03  6.4430e-01        0.2597
## 0206  9.8631e-01 -3.7731e-03  1.3153e-01  2.7300e+00        0.0032
## 0207  2.8961e-02 -1.8708e-03  4.2933e-02  1.4880e-01        0.4409
## 0301  1.7788e+00 -1.1883e-02  4.9540e-01  2.5441e+00        0.0055
## 0302  6.3483e-02 -7.4899e-04  1.6091e-01  1.6012e-01        0.4364
## 0303 -6.1636e-02 -1.8734e-05  2.0045e-03 -1.3763e+00        0.9156
## 0304  1.3739e-01 -8.6266e-04  4.5676e-02  6.4688e-01        0.2589
## 0305  4.6022e-01 -1.6314e-03  5.6991e-02  1.9346e+00        0.0265
## 0306  1.2467e-01 -7.8076e-05  1.6785e-02  9.6292e-01        0.1678
## 0307 -3.2744e-01 -1.5403e-03  6.4887e-02 -1.2794e+00        0.8996
## 0401 -1.6499e+00 -2.1483e-02  1.4924e+00 -1.3330e+00        0.9087
## 0402 -8.1725e-02 -9.4006e-05  4.9812e-03 -1.1566e+00        0.8763
## 0403  5.4137e-01 -4.1514e-03  2.1908e-01  1.1655e+00        0.1219
## 0404  1.6404e-02 -8.7178e-05  4.6194e-03  2.4264e-01        0.4041
## 0405 -6.2487e-03 -4.6154e-03  1.1941e-01 -4.7265e-03        0.5019
## 0406  5.6747e-01 -3.8972e-03  1.3584e-01  1.5502e+00        0.0605
## 0501 -2.2810e-03 -3.5310e-03  3.7647e-01  2.0373e-03        0.4992
## 0502  1.8251e-01 -4.3500e-04  3.0869e-02  1.0413e+00        0.1489
## 0503 -2.9675e-01 -4.5285e-03  1.1717e-01 -8.5369e-01        0.8034
## 0504 -1.5030e-01 -3.1108e-03  1.6434e-01 -3.6309e-01        0.6417
## 0505  4.7170e-02 -3.3999e-05  1.1897e-03  1.3686e+00        0.0856
## 0506 -1.5424e-01 -1.5012e-03  5.2451e-02 -6.6694e-01        0.7476
## 0507  6.0308e-01 -7.3661e-03  3.0850e-01  1.0990e+00        0.1359
## 0601 -4.4155e-01 -2.2674e-02  1.1743e+00 -3.8654e-01        0.6505
## 0602 -9.7725e-02 -2.7846e-03  5.1727e-02 -4.1744e-01        0.6618
## 0603 -2.9725e-02 -7.0262e-03  2.4413e-01 -4.5940e-02        0.5183
## 0604 -8.1137e-02 -1.0809e-03  3.7782e-02 -4.1186e-01        0.6598
## 0605 -3.4670e-02 -5.3780e-05  2.8498e-03 -6.4845e-01        0.7417
## 0606  1.0710e+00 -8.9949e-03  1.9165e+00  7.8014e-01        0.2177
## 0607 -2.7191e-03 -4.9360e-04  1.7264e-02 -1.6938e-02        0.5068
## 0608 -1.3600e-01 -1.7431e-03  9.2211e-02 -4.4213e-01        0.6708
## 0609 -6.0517e-02 -2.0579e-04  1.0903e-02 -5.7759e-01        0.7182
## 0610 -2.7352e-01 -1.4353e-03  6.0473e-02 -1.1064e+00        0.8657
## 0701 -3.0388e-01 -9.7001e-03  2.8673e-01 -5.4939e-01        0.7086
## 0702  6.8565e-02 -5.6938e-04  1.9913e-02  4.8993e-01        0.3121
## 0703  6.3926e-01 -7.4768e-03  5.2684e-01  8.9102e-01        0.1865
## 0704  2.5289e-01 -1.2167e-03  1.3002e-01  7.0471e-01        0.2405
## 0705 -2.0614e-01 -1.3870e-03  1.4820e-01 -5.3188e-01        0.7026
## 0706  3.2383e-02 -1.2860e-04  2.7646e-02  1.9553e-01        0.4225
## 0707  5.6854e-03 -1.7806e-04  1.2639e-02  5.2156e-02        0.4792
## 0708  1.0795e+00 -1.4048e-02  3.1845e-01  1.9379e+00        0.0263
## 0709 -8.9810e-02 -7.7390e-03  2.6871e-01 -1.5833e-01        0.5629
## 0710  0.0000e+00  0.0000e+00  0.0000e+00         NaN           NaN
## 0711  6.3109e-01 -7.3046e-03  3.8427e-01  1.0299e+00        0.1515
## 0712  8.4284e-01 -9.9984e-03  5.2455e-01  1.1775e+00        0.1195
## 0713  7.3873e-01 -4.7526e-03  3.3580e-01  1.2830e+00        0.0997
## 0714  3.5345e-01 -5.9354e-04  1.7706e-02  2.6607e+00        0.0039
## 0801  8.5326e-01 -9.6017e-03  4.0123e-01  1.3622e+00        0.0866
## 0802 -1.9608e+00 -8.0699e-03  8.5647e-01 -2.1100e+00        0.9826
## 0803 -9.2086e-02 -1.1852e-02  6.2061e-01 -1.0185e-01        0.5406
## 0804  3.8871e-01 -5.2277e-03  2.1942e-01  8.4099e-01        0.2002
## 0805 -1.4728e-01 -1.6199e-03  1.7304e-01 -3.5017e-01        0.6369
## 0806  1.0033e-01 -2.6266e-04  1.1079e-02  9.5569e-01        0.1696
## 0807  9.1053e-02 -1.1206e-02  3.8774e-01  1.6422e-01        0.4348
## 0901 -1.5926e+00 -2.0468e-02  1.0624e+00 -1.5252e+00        0.9364
## 0902  4.6350e-01 -4.2629e-03  1.1033e-01  1.4083e+00        0.0795
## 0903 -2.2747e-01 -1.8277e-03  9.6681e-02 -7.2570e-01        0.7660
## 0904  9.6141e-01 -5.0127e-03  3.5408e-01  1.6241e+00        0.0522
## 0905  5.6268e-01 -6.6126e-03  2.7716e-01  1.0814e+00        0.1398
## 0906 -4.0467e-02 -6.1377e-05  2.1476e-03 -8.7190e-01        0.8084
## 0907  1.1765e-01 -1.9269e-03  3.5824e-02  6.3176e-01        0.2638
## 0908 -5.7075e-01 -5.3735e-03  2.8323e-01 -1.0623e+00        0.8560
## 0909 -6.8788e-02 -6.6168e-04  1.0320e-02 -6.7063e-01        0.7488
## 0910 -1.9001e-01 -8.6266e-04  3.6367e-02 -9.9187e-01        0.8394
## 0911 -2.3265e-03 -1.1066e-05  3.8723e-04 -1.1766e-01        0.5468
## 0912 -1.6926e-01 -1.5122e-03  8.0016e-02 -5.9301e-01        0.7234
## 0913 -7.1286e-02 -9.4443e-04  3.9810e-02 -3.5255e-01        0.6378
## 0914 -8.0135e-01 -7.8022e-03  5.4959e-01 -1.0704e+00        0.8578
## 0916 -1.2494e+00 -1.5260e-02  3.9059e-01 -1.9747e+00        0.9758
## 0918  1.6310e-01 -5.7470e-03  1.7056e-01  4.0884e-01        0.3413
## 0919  1.5585e+00 -7.8909e-03  5.5579e-01  2.1011e+00        0.0178
## 0920  7.7080e-02 -5.8663e-04  4.1623e-02  3.8069e-01        0.3517
## 0921  9.0311e-02 -1.5514e-03  8.2088e-02  3.2063e-01        0.3742
## 0922 -3.6900e-01 -4.7330e-03  1.9875e-01 -8.1708e-01        0.7931
## 0923  0.0000e+00  0.0000e+00  0.0000e+00         NaN           NaN
## 0924 -2.3176e-01 -3.5142e-03  1.2254e-01 -6.5203e-01        0.7428
## 0925 -3.9241e-02 -9.6781e-05  2.8886e-03 -7.2834e-01        0.7668
## 0927 -1.8305e-01 -1.5514e-03  8.2088e-02 -6.3347e-01        0.7368
## 0928 -5.4057e-01 -8.0699e-03  2.8010e-01 -1.0061e+00        0.8428
## 1001 -1.1452e+00 -3.3906e-02  2.3255e+00 -7.2873e-01        0.7669
## 1002  0.0000e+00  0.0000e+00  0.0000e+00         NaN           NaN
## 1003  3.1608e-02 -3.5544e-04  2.5225e-02  2.0125e-01        0.4203
## 1004  0.0000e+00  0.0000e+00  0.0000e+00         NaN           NaN
## 1005  1.6949e-01 -6.9517e-04  4.9318e-02  7.6633e-01        0.2217
## 1006 -8.3144e-02 -4.1436e-04  2.1949e-02 -5.5841e-01        0.7117
## 1101 -1.2156e+00 -2.7259e-02  1.8825e+00 -8.6611e-01        0.8068
## 1102  0.0000e+00  0.0000e+00  0.0000e+00         NaN           NaN
## 1103 -4.4154e-01 -4.2815e-03  1.7987e-01 -1.0310e+00        0.8487
## 1104  1.8948e-02 -1.8734e-05  7.9044e-04  6.7460e-01        0.2500
## 1105 -1.7707e-01 -1.8277e-03  6.3840e-02 -6.9357e-01        0.7560
## 1106 -1.7062e-01 -4.6643e-03  1.0674e-01 -5.0797e-01        0.6943
## 1107  5.8188e-01 -2.0481e-03  8.6239e-02  1.9884e+00        0.0234
## 1108 -4.7402e-02 -4.5476e-03  1.9100e-01 -9.8057e-02        0.5391
## 1109  1.5432e-01 -1.0716e-03  2.7823e-02  9.3161e-01        0.1758
## 1110  1.3434e-02 -1.0141e-05  5.3740e-04  5.7994e-01        0.2810
## 1111  4.8552e-01 -2.6283e-03  1.1060e-01  1.4678e+00        0.0711
## 1112  0.0000e+00  0.0000e+00  0.0000e+00         NaN           NaN
## 1113  0.0000e+00  0.0000e+00  0.0000e+00         NaN           NaN
## 1114 -5.5274e-01 -2.2811e-03  1.6157e-01 -1.3694e+00        0.9146
## 1115 -3.0379e-01 -1.5798e-03  1.1198e-01 -9.0310e-01        0.8168
## 1116  9.3564e-01 -4.4613e-03  3.1531e-01  1.6742e+00        0.0470
## 1201  0.0000e+00  0.0000e+00  0.0000e+00         NaN           NaN
## 1202  7.0375e-02 -6.6943e-03  4.7207e-01  1.1217e-01        0.4553
## 1203 -3.2249e-02 -1.3582e-05  1.4532e-03 -8.4562e-01        0.8011
## 1204  4.6677e-02 -1.2418e-03  1.3270e-01  1.3154e-01        0.4477
## 1205  4.6042e-02 -8.0482e-04  4.2616e-02  2.2693e-01        0.4102
## 1206  8.7022e-01 -2.7846e-03  2.9711e-01  1.6016e+00        0.0546
## 1207 -2.7583e-02 -1.5012e-03  3.2227e-01 -4.5943e-02        0.5183
## 1208 -2.6870e-01 -1.1382e-03  8.0711e-02 -9.4180e-01        0.8269
## 1209  9.9237e-01 -5.9770e-03  2.5068e-01  1.9940e+00        0.0231
## 1210 -2.2472e-02 -1.1066e-05  7.8564e-04 -8.0134e-01        0.7885
## 1211  6.6531e-01 -4.5962e-03  1.9303e-01  1.5247e+00        0.0637
## 1212  4.1312e-01 -4.7330e-03  3.3442e-01  7.2256e-01        0.2350
## 1213  0.0000e+00  0.0000e+00  0.0000e+00         NaN           NaN
## 1301 -4.8593e-01 -3.1829e-03  3.3947e-01 -8.2854e-01        0.7963
## 1302 -2.7212e-03 -4.1977e-03  2.2152e-01  3.1371e-03        0.4987
## 1303  7.9436e-02 -4.9360e-04  3.5025e-02  4.2709e-01        0.3347
## 1304  1.1802e-01 -3.8859e-04  4.1562e-02  5.8082e-01        0.2807
## 1305 -1.8488e-01 -7.9802e-03  5.6202e-01 -2.3596e-01        0.5933
## 1306  9.2901e-02 -1.1382e-03  8.0711e-02  3.3101e-01        0.3703
## 1307  5.7140e-01 -2.2335e-03  5.1237e-02  2.5342e+00        0.0056
## 1308 -1.2230e+00 -1.1774e-02  4.9091e-01 -1.7287e+00        0.9581
## 1309  6.7361e-01 -4.8022e-03  1.6723e-01  1.6590e+00        0.0486
## 1310  1.1926e+00 -8.0699e-03  2.8010e-01  2.2686e+00        0.0116
## 1311  8.4294e-01 -6.8593e-03  2.8742e-01  1.5851e+00        0.0565
## 1312  5.5499e-02 -7.5584e-05  3.1889e-03  9.8414e-01        0.1625
## 1313  2.0994e-01 -1.7852e-03  7.5186e-02  7.7215e-01        0.2200
## 1314  6.4685e-01 -1.0716e-03  7.5994e-02  2.3503e+00        0.0094
## 1315  2.0272e-01 -1.5403e-03  2.6131e-02  1.2636e+00        0.1032
## 1316  9.2879e-01 -8.4339e-03  3.5284e-01  1.5778e+00        0.0573
## 1317  7.1162e-01 -7.3661e-03  2.1825e-01  1.5390e+00        0.0619
## 1318  2.9762e-01 -1.0890e-02  5.7083e-01  4.0834e-01        0.3415
## 1319  1.0440e+00 -6.8593e-03  3.6100e-01  1.7490e+00        0.0401
## 1320 -3.7915e-02 -4.7330e-03  2.4963e-01 -6.6414e-02        0.5265
## 1321  2.1440e-01 -2.4668e-04  1.7508e-02  1.6222e+00        0.0524
## 1322 -1.3699e-01 -2.0481e-03  4.2088e-02 -6.5778e-01        0.7447
## 1401 -6.3482e-01 -2.9066e-03  1.5358e-01 -1.6125e+00        0.9466
## 1402  5.6799e-02 -1.1478e-03  8.1391e-02  2.0312e-01        0.4195
## 1403 -9.0823e-04 -1.7428e-04  7.3522e-03 -8.5597e-03        0.5034
## 1404 -4.5577e-01 -2.5409e-03  1.7993e-01 -1.0685e+00        0.8573
## 1405  1.1122e-02 -4.3500e-04  3.0869e-02  6.5778e-02        0.4738
## 1406 -1.2620e-01 -4.3472e-03  2.2937e-01 -2.5443e-01        0.6004
## 1407  3.6022e-03 -8.8059e-03  3.6827e-01  2.0447e-02        0.4918
## 1408 -8.0543e-02 -3.5544e-04  3.8017e-02 -4.1126e-01        0.6596
## 1409  1.8647e-01 -1.2292e-02  8.6192e-01  2.1409e-01        0.4152
## 1410  1.1248e+00 -2.9454e-03  6.3140e-01  1.4192e+00        0.0779
## 1411 -2.4232e-02 -6.1377e-05  6.5666e-03 -2.9828e-01        0.6173
## 1412 -4.6339e-01 -9.3789e-03  9.9409e-01 -4.5536e-01        0.6756
## 1501 -3.1470e-02 -2.1969e-05  7.6872e-04 -1.1342e+00        0.8717
## 1503 -3.3250e-01 -3.5738e-03  1.8871e-01 -7.5718e-01        0.7755
## 1504  2.1250e-01 -3.1829e-03  8.2466e-02  7.5107e-01        0.2263
## 1601  3.6014e-01 -2.1997e-03  6.5515e-02  1.4156e+00        0.0784
## 1602  1.3005e+00 -4.4135e-03  9.4472e-01  1.3426e+00        0.0897
## 1603  7.5359e-01 -1.3500e-03  2.8985e-01  1.4022e+00        0.0804
## 1604 -3.7412e-02 -9.1859e-03  6.4615e-01 -3.5114e-02        0.5140
## 1701 -9.3484e-01 -5.5095e-02  2.1965e+00 -5.9359e-01        0.7236
## 1702 -1.3779e-01 -2.3429e-03  9.8622e-02 -4.3129e-01        0.6669
## 1703  1.3005e+00 -8.3681e-03  1.7841e+00  9.7994e-01        0.1636
## 1704  1.3431e-01 -5.4573e-04  2.8904e-02  7.9323e-01        0.2138
## 1705  8.3016e-01 -5.0718e-02  2.5514e+00  5.5148e-01        0.2907
## 1707  4.8159e-02 -1.4650e-05  7.7631e-04  1.7290e+00        0.0419
## 1708  3.6658e-01 -2.8391e-04  2.0150e-02  2.5844e+00        0.0049
## 1709 -8.0894e-01 -7.1951e-03  5.0713e-01 -1.1258e+00        0.8699
## 1801 -1.3855e+00 -1.8375e-02  3.8779e+00 -6.9422e-01        0.7562
## 1802 -6.2389e-01 -2.5630e-02  2.6720e+00 -3.6599e-01        0.6428
## 1803  0.0000e+00  0.0000e+00  0.0000e+00         NaN           NaN
## 1804  0.0000e+00  0.0000e+00  0.0000e+00         NaN           NaN
## 1805  2.7171e-01 -1.7806e-04  1.2639e-02  2.4184e+00        0.0078
## 1806  6.7863e-02 -1.5798e-03  4.7082e-02  3.2004e-01        0.3745
## 1807 -3.9940e-02 -1.4650e-05  4.3727e-04 -1.9093e+00        0.9719
## 1808 -5.9691e-02 -7.7666e-04  2.0171e-02 -4.1481e-01        0.6609
## 1809  1.5976e-02 -3.6910e-04  1.5567e-02  1.3100e-01        0.4479
## 1901 -2.4332e-01 -6.4736e-03  1.9198e-01 -5.4056e-01        0.7056
## 1902 -3.4691e-03 -3.8859e-04  2.0585e-02 -2.1471e-02        0.5086
## 1903  1.4519e-02 -1.7525e-05  9.2866e-04  4.7701e-01        0.3167
## 1904 -2.2653e-01 -4.4613e-03  4.7521e-01 -3.2214e-01        0.6263
## 1905 -9.3063e-02 -2.6729e-04  1.4160e-02 -7.7981e-01        0.7822
## 1906 -6.0048e-03 -3.3999e-05  1.8017e-03 -1.4067e-01        0.5559
## 1907  6.3483e-02 -1.1749e-04  2.5258e-02  4.0019e-01        0.3445
## 1908 -1.5634e-01 -1.6199e-03  6.8237e-02 -5.9231e-01        0.7232
## 1909  8.3866e-01 -4.3285e-03  1.2864e-01  2.3503e+00        0.0094
##  [ reached getOption("max.print") -- omitted 15 rows ]
# Veamos los valores significativos <0.05
la <- na.omit(locMorMat[which(locMorMat$Pr.z...E.Ii..<0.05),])
la
##              Ii          E.Ii       Var.Ii     Z.Ii Pr.z...E.Ii..
## 0101 1.81382857 -3.692141e-02 1.0613807692 1.796438   0.036212471
## 0105 0.45277446 -1.048224e-03 0.0366411848 2.370836   0.008873961
## 0206 0.98631443 -3.773148e-03 0.1315324756 2.729965   0.003167050
## 0301 1.77875705 -1.188259e-02 0.4953987395 2.544080   0.005478301
## 0305 0.46022347 -1.631351e-03 0.0569913970 1.934643   0.026517061
## 0708 1.07950301 -1.404784e-02 0.3184458131 1.937855   0.026320476
## 0714 0.35345024 -5.935355e-04 0.0177059821 2.660706   0.003898846
## 0919 1.55853767 -7.890937e-03 0.5557866004 2.101148   0.017814004
## 1107 0.58188495 -2.048146e-03 0.0862393824 1.988430   0.023382086
## 1116 0.93563596 -4.461341e-03 0.3153142454 1.674176   0.047048047
## 1209 0.99237245 -5.977030e-03 0.2506784373 1.993995   0.023076290
## 1307 0.57139777 -2.233496e-03 0.0512370617 2.534199   0.005635236
## 1309 0.67361429 -4.802151e-03 0.1672307534 1.658969   0.048561012
## 1310 1.19257518 -8.069907e-03 0.2801047910 2.268581   0.011646895
## 1314 0.64684782 -1.071583e-03 0.0759941456 2.350341   0.009378104
## 1319 1.04400621 -6.859264e-03 0.3609993832 1.749016   0.040144089
## 1707 0.04815897 -1.464956e-05 0.0007763121 1.728985   0.041905926
## 1708 0.36657521 -2.839137e-04 0.0201503733 2.584388   0.004877594
## 1805 0.27170746 -1.780555e-04 0.0126385767 2.418449   0.007793408
## 1909 0.83866198 -4.328529e-03 0.1286434126 2.350330   0.009378378
## 2103 1.01185808 -5.671269e-03 0.1973245702 2.290637   0.010992195
## 2401 0.40026800 -4.500908e-04 0.0319392449 2.242212   0.012473844
# Veamos los valores significativos > 0.95
lb <- na.omit(locMorMat[locMorMat$Pr.z...E.Ii..>0.95,])
lb
##               Ii          E.Ii       Var.Ii      Z.Ii Pr.z...E.Ii..
## 0202 -0.50105637 -1.785173e-03 0.0463170021 -2.319884     0.9898264
## 0802 -1.96077729 -8.069907e-03 0.8564742647 -2.109988     0.9825703
## 0916 -1.24935621 -1.526019e-02 0.3905867534 -1.974650     0.9758461
## 1308 -1.22299700 -1.177365e-02 0.4909108951 -1.728713     0.9580698
## 1807 -0.03994032 -1.464956e-05 0.0004372697 -1.909315     0.9718892
## 2202 -0.88838041 -6.942477e-03 0.2412458247 -1.794576     0.9636394

Veamos el mapa de significancia del estadístico de Moran local de las unidades espaciales.

names(locMorMat)[5]<-"Prob"
brks <- c(min(locMorMat[,5],na.rm = TRUE), 0.05000, 0.95000, max(locMorMat[,5],na.rm = TRUE))
cols <- c("grey30", "grey90", "grey60")
cols <- c("royalblue", "springgreen", "yellow")
sp::plot(poligonos, col=cols[findInterval((locMorMat[,5]), brks)])
legend("bottomleft", legend=c("rodeado de relativamente altos valores,
locM>0", "no significativo", "rodeado de relativamente bajos valores, locM<0"),
fill = cols, bty="n", cex=0.75)
title(main=" Estadístico local de Moran ", cex=0.7)

7.3.2 Estadístico Local de Geary

El estadístico local de Geary se define como

\[ C_i = \frac{1}{m_2} \sum_j w_{ij}(x_i-x_j) \] donde \(m_2 = n^{-1}\sum_i(x_i-\bar{x})^2\). El valor esperado del Geary local se obtiene de la siguiente manera bajo aleatorización:

\[ E(C_i) = \frac{2nw_i}{(n-1)} \]

7.3.3 Estadístico Local: Getis-Ord \(G_i\) y \(G_i^*\)

Los estadísticos \(G_i\) y \(G_i^*\) también son LISA y están dados por:

\[ G_i = \frac{\sum_{j\neq i}w_{ij}x_j}{\sum_{j\neq i}x_j} \]

\[ G_i^* = \frac{\sum_{j}w_{ij}x_j}{\sum_{j}x_j} \] donde \(x_i\) es la observación en la unidad espacial \(i\) y \(w_{ij}\) es un elemento de la matriz de pesos espaciales \(W\). Note que \(G_i^*\) incluye a la unidad espacial analizada, mientras \(G_i\) no. Tanto \(G_i\) como \(G_i^*\) requieren que \(x\) sea positivo.

Inferencia

Para la inferencia se asume una distribución \(t\) para \(G_i\) y \(G_i^*\). La hipótesis nula es que \(G_i=0\).

  • Valores positivos y significativos en la unidad espacial significan una agrupación de unidad espacial de alto valor (unidad analizada y vecinas), el llamado cluster de alto valor.

  • Valores negativos del estadístico \(G_i\) indican la existencia de un grupo de regiones con valores bajos de la variable examinada, es decir, la unidad espacial \(i\) está rodeada por unidades espaciales similares de valor bajo.

  • Valores positivos y significativas de Getis-Ord tienen un valor \(p\) por debajo de \(0.05\), mientras que las estadísticas negativas y significativas tienen un valor \(p\) por encima de \(0.95\).

Veamos su cálculo:

# local Gi y Gi^*
locG <- localG(poligonos$POBRES_P, cont.listw,zero.policy = TRUE)
locGstar<-localG(poligonos$POBRES_P, nb2listw(include.self(cont.nb)))
summary(locG)
##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max.     NA's 
## -2.72996 -0.68409  0.09806  0.08202  0.86900  2.53420       12
summary(locGstar)
##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
## -2.96275 -0.65507  0.07007  0.05473  0.77959  2.62669

Revisamos la significancia:

# significancia con la prueba t para n=100
sig <- ifelse(locG<=-3.289 | locG>=3.289, "*", " ")
which(sig=="*")
## integer(0)
sig <- ifelse(locGstar<=-3.289 | locG>=3.289, "*", " ")
which(sig=="*")
## integer(0)

En este caso no tenemos significancia.

val <- qt(0.975,df=length(poligonos$POBRES_P)-1)
# significancia con la prueba t para n=100
sig <- ifelse(locG<=-val | locG>=val, "*", " ")
which(sig=="*")
##  [1]   5  17  21  23  66  68  88  90 111 129 140 143 147 181 187 200 203 213
sig <- ifelse(locGstar<=-val | locG>=val, "*", " ")
which(sig=="*")
##  [1]   1   5  17  21  23  27  60  66  88  90 111 129 140 143 147 181 187 200 203 213

En este caso si tenemos significancia. Veamos el mapa de las unidades espaciales:

# graph of significant locG statistics
brks<-c(-10,-val, val, 10)
cols<-c("grey40", "grey85", "grey55")
cols<-c("red", "grey85", "green")
sp::plot(poligonos, col=cols[findInterval(locG, brks)])
legend("bottomleft", legend=c("G significativo - rodeado de valores relativamente bajos", "G no significativo", "G significativo - rodeado por valores relativamente altos"), fill=cols, bty="n", cex=0.75)
title(main="Valores G significativos")

7.3.4 Heterocedasticidad espacial local

Es posible extender el análisis de similitud de vecindades mediante la estabilidad de una variable en el espacio utilizando la medida de heterocedasticidad espacial local (LOSH) (\(H_i\)) (heterocedasticidad espacial local), complementaria a la \(G_i\).

El análisis combinado de ambos indicadores permite medir la similitud e inestabilidad de una variable en el espacio. El estadístico \(H_i\) mide la varianza de una variable en la vecindad, indicando áreas con variabilidad uniforme y no uniforme, de manera análoga a la medida de autocorrelación en el indicador \(G_i\).

La medida local de heterocedasticidad espacial LOSH (Hi) viene dada por la fórmula:

\[ H_i = \frac{\sum_jw_{ij}|e_j|^\alpha}{\sum_jw_{ij}} \] donde \(e_j = x_j - \bar{x_j}\), \(w_{ij}\) es un elemento de la matriz de pesos espaciales \(W\). Si \(\alpha=1\), \(H_i\) corresponde a una medida de desviación absoluta, si \(\alpha=2\), \(H_i\) corresponde a una medida de varianza.

La interpretación combinada se la puede apreciar en la siguiente tabla:

Efectos espaciales usando \(H_i\) y \(G_i\). Fuente: Yamagata and Seya (2019).
\(H_i\) alto \(H_i\) bajo
\(G_i\) alto Un punto caliente con condiciones
locales heterogéneas
Un punto caliente con unidades
espaciales similares
\(G_i\) bajo Condiciones locales heterogéneas
pero a un nivel medio bajo (poco frecuente)
Condiciones locales homogéneas
y un nivel medio bajo

Veamos el cálculo del estadístico:

# Gi local
locG <- localG(poligonos$POBRES_P, cont.listw)
a <- summary(locG)
a
##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max.     NA's 
## -2.72996 -0.68409  0.09806  0.08202  0.86900  2.53420       12

Notemos que la siguiente función presenta un problema por valores perdidos:

#Hi LOSH local 
locH <- LOSH(poligonos$POBRES_P, cont.listw, a=2, var_hi=TRUE, zero.policy=TRUE,na.action=na.exclude)

Creamos LOSH1 para arreglar el problema:

LOSH1 <- function (x, listw, a = 2, var_hi = TRUE, zero.policy = NULL, 
          na.action = na.fail, spChk = NULL) 
{
  if (is.null(zero.policy)) 
    zero.policy <- get("zeroPolicy", envir = .spdepOptions)
  stopifnot(is.logical(zero.policy))
  n <- length(listw$neighbours)
  if (n != length(x)) 
    stop("Different numbers of observations")
  NAOK <- deparse(substitute(na.action)) == "na.pass"
  NAOK <- TRUE
  if (var_hi) {
    res <- matrix(nrow = n, ncol = 6)
    colnames(res) <- c("Hi", "E.Hi", "Var.Hi", "Z.Hi", "x_bar_i", 
                       "ei")
  }
  else {
    res <- matrix(nrow = n, ncol = 3)
    colnames(res) <- c("Hi", "x_bar_i", "ei")
  }
  Wi <- vapply(listw$weights, sum, FUN.VALUE = 0)
  res[, "x_bar_i"] <- lag.listw(listw, x, zero.policy = zero.policy, 
                                NAOK = NAOK)/Wi
  res[, "ei"] <- abs(x - res[, "x_bar_i"])^a
  denom_hi <- mean(res[, "ei"],na.rm = TRUE) * Wi
  res[, "Hi"] <- lag.listw(listw, res[, "ei"], zero.policy = zero.policy, 
                           NAOK = NAOK)/denom_hi
  if (var_hi) {
    var_ei <- (sum(res[, "ei"]^2,na.rm = TRUE)/n) - mean(res[, "ei"],na.rm = TRUE)^2
    res[, "Var.Hi"] <- (n - 1)^(-1) * denom_hi^(-2) * var_ei * 
      (n * vapply(listw$weights, function(y) sum(y^2), 
                  FUN.VALUE = 0) - Wi^2)
    res[, "Z.Hi"] <- (2 * res[, "Hi"])/res[, "Var.Hi"]
    res[, "E.Hi"] <- 1
  }
  res
}
#Hi LOSH local 
locH <- LOSH1(poligonos$POBRES_P, cont.listw, a=2, var_hi=TRUE, zero.policy=TRUE,na.action=na.exclude)
summary(locH)
##        Hi               E.Hi       Var.Hi            Z.Hi             x_bar_i      
##  Min.   :0.02379   Min.   :1   Min.   :0.2074   Min.   : 0.01665   Min.   :0.5620  
##  1st Qu.:0.38869   1st Qu.:1   1st Qu.:0.4650   1st Qu.: 1.09681   1st Qu.:0.7215  
##  Median :0.63771   Median :1   Median :0.7041   Median : 2.08861   Median :0.7719  
##  Mean   :0.96706   Mean   :1   Mean   :0.8607   Mean   : 2.98045   Mean   :0.7720  
##  3rd Qu.:1.22984   3rd Qu.:1   3rd Qu.:0.9433   3rd Qu.: 4.16863   3rd Qu.:0.8244  
##  Max.   :5.38932   Max.   :1   Max.   :2.8565   Max.   :17.54567   Max.   :0.9590  
##  NA's   :12                    NA's   :12       NA's   :12         NA's   :12      
##        ei          
##  Min.   :0.000002  
##  1st Qu.:0.001702  
##  Median :0.007519  
##  Mean   :0.021598  
##  3rd Qu.:0.022525  
##  Max.   :0.255834  
##  NA's   :12

Guardamos los valores relevantes del estadístico \(H_i\):

b <- summary(locH[,"Hi"])
b
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
## 0.02379 0.38869 0.63771 0.96706 1.22984 5.38932      12

Veamos el gráfico conjunto de \(G_i\) y \(H_i\):

par(mfrow = c(1,2))
# mapa de Gi 
brks <- c(a[1], a[2], a[3], a[5], a[6])
colfunc <- colorRampPalette(c("royalblue", "springgreen", "yellow", "red"))
coli <- colfunc(5)
sp::plot(poligonos, col=coli[findInterval(locG, brks)], main=" Gi statistics")
legend("bottomright", legend=c("Very low", "Low", "Medium", "High", "Very
high"), fill=coli, bty="n")


# mapa de LOSH 
brks <- c(b[1], b[2], b[3], b[5], b[6])
sp::plot(poligonos, col=coli[findInterval(locH[,"Hi"], brks)])
legend("bottomright", legend=c("Very low", "Low", "Medium", "High", "Very
high"), fill=coli, bty="n")
title(main="Hi statistics")

par(mfrow = c(1,1))

7.4 Cross correlación espacial entre dos variables

Global

Introducido por Chen (2015), el índice de cross-correlación espacial global (GSCI por sus siglas en inglés) está dado por:

\[ R_c = x^TWy \] y es un estadístico escalado entre \([-1,1]\); \(x\) y \(y\) son las variables estudiadas y \(W\) es la matriz de pesos espaciales.

Este es una medida indirecta que calcula correlación relacionada con distancia.

También se tiene una relación entre correlación tradicional de dos variables y \(R_c\):

\[ R_0 = R_c + R_p \]

donde \(R_0\) es un índice de correlación de Pearson, donde la dependencia consiste de relación entre valores de la variable y relaciones en el espacio. \(R_c\) es el coeficiente GSCI y \(R_p\) es una cross-correlación parcial, que expresa la influencia directa de las variables sin efectos espaciales.

Local

También se tiene medidas locales de cross-correlación, el índice de cross-correlación espacial local (LSCI por sus siglas en inglés):

\[ R_i^{(xy)} = x_i\sum_jw_{ij}y_j \] \[ R_j^{(yx)} = y_i\sum_jw_{ij}x_j \] donde \(i\) y \(j\) son diferentes unidades espaciales y \(w_{ij}\) son elementos de la matriz de pesos espaciales \(W\).

El GSCI refleja la suma de correlaciones cruzadas entre dos (todos) elementos cualesquiera, mientras que el LSCI mide la correlación cruzada entre un elemento dado y todos los demás en el sistema.

library(spatialEco) # crossCorrelation()
library(GISTools) # graphics, choropleth()
# cont.mat <- spdep::nb2mat(cont.nb,zero.policy=TRUE) # W matrix convertida a tipo matriz
cont.mat <- as.matrix(dist(coordinates(poligonos)))
x1 <- poligonos$VAB_PC 
x2 <- poligonos$POBRES_P
# versión con W
# ii <- crossCorrelation(x1, x2, w=cont.mat)
ii <- crossCorrelation(x1, x2, w=cont.mat, type=c("LSCI", "GSCI"), k=999,
scale.xy=TRUE, scale.partial=TRUE, scale.matrix=TRUE, alpha=0.05)
ii
## Moran's-I under randomization assumptions... 
##   First-order Moran's-I:  0.00121441 
##   First-order p-value:  817 
## Chen's SCI under randomization assumptions... 
## 
##  Summary statistics of local partial cross-correlation [xy] 
##      Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
## -1.000000 -0.051547 -0.005889 -0.006475  0.044144  1.000000 
##  
##   p-value based on 2-tailed t-test:  0.02102102 
##   p-value based on 2-tailed t-test observations above/below CI:  0.04204204 
## 
##  Counts of cluster types
## High.High  High.Low  Low.High   Low.Low 
##        28        43        61        83
  • scale.partial: estandariza el estadístico entre -1 y 1.
  • scale.matrix: estandariza \(W\) tal que sus filas sumen 1.
  • alpha: nivel de significancia.
  • k número de simulaciones.
  • scale.xy estandarización de las variables analizadas.
  • LSCI: local spatial cross-correlation index.
  • GSCI: global spatial cross-correlation index.

Los resultados muestran que el GSCI es 0.0012144 y su p-valor es 817, lo que indica que las dos variables no están correlacionadas en el espacio.

Veamos la correlación de pearson:

cor(x1,x2)
## [1] -0.3547163

Según las expresiones anteriores, tenemos que \(R_0=-0.3547\), \(R_c=0.00006\) y \(R_p = R_0-R_c\).

Respecto a la correlación local LSCI, vemos que con un \(\alpha = .1\) si es significativo porque su p-valor es 0.021021. Veamos algunos resultados:

head(ii$SCI)
##       lsci.xy    lsci.yx
## 1 -0.33685997 0.78317737
## 2  0.05493543 0.24457315
## 3  0.06117472 0.26364556
## 4  0.20225330 0.07121259
## 5  0.05273068 0.23528730
## 6  0.22108706 0.02374638

Las estadísticas locales lsci.xy deben interpretarse como la influencia de \(x\) sobre \(y\).

A continuación se calcula la correlación basándonos en las distancias en los centroides:

# version with coordinates and distance
X <- coordinates(poligonos)
# colnames(X) <- c("x","y")
# X <- as.matrix(dist(X))
iii <- crossCorrelation(x1, x2, coords=X, type = c("LSCI","GSCI"), k=999, dist.function="inv.power", scale.xy=TRUE, scale.partial=TRUE, scale.matrix=FALSE, alpha=0.05)
iii
## Moran's-I under randomization assumptions... 
##   First-order Moran's-I:  0.00121441 
##   First-order p-value:  810 
## Chen's SCI under randomization assumptions... 
## 
##  Summary statistics of local partial cross-correlation [xy] 
##      Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
## -1.000000 -0.051547 -0.005889 -0.006475  0.044144  1.000000 
##  
##   p-value based on 2-tailed t-test:  0.02602603 
##   p-value based on 2-tailed t-test observations above/below CI:  0.05205205 
## 
##  Counts of cluster types
## High.High  High.Low  Low.High   Low.Low 
##        28        43        61        83
caja <- bbox(poligonos)
# mapping of result ii
choropleth(poligonos, ii$SCI[,"lsci.xy"])
shades <- auto.shading(ii$SCI[,"lsci.xy"])
choro.legend(1100000,9445216+9445216*0.02, shades, cex=0.65, bty="n")
title(main="Matriz de contigüidad")

# mapping of result iii
choropleth(poligonos, iii$SCI[,"lsci.xy"])
shades <- auto.shading(iii$SCI[,"lsci.xy"])
choro.legend(1100000,9445216+9445216*0.02, shades, cex=0.65, bty="n")
title(main="Matriz de distancia inversa al cuadrado")

En los gráficos son visibles agrupaciones de valores altos y bajos. Las áreas claras (correlación cruzada espacial negativa fuerte) significan una diferenciación espacial mucho más fuerte (relación de las variables y efectos espaciales) que los campos más oscuros (correlación cruzada espacial negativa media).