Objetivos

El objetivo de este laboratorio será introducir a los estudiantes a los siguientes elementos clave que son utilizados en el análisis exploratorio de datos espaciales (ESDA) y modelos econométricos. Estos elementos son:

Exploración preliminar de datos espaciales.
Matrices de peso espacial
Medidas de autocorrelación espacial: global y local.

Paquetes que se utilizarán

Antes de comenzar, debemos incluir los paquetes estadísticos necesarios para llevar a cabo el análisis exploratorio espacial, la construcción de matrices de peso espacial y medidas de autocorrelación. Preliminarmente, estos paquetes son:

(Estos paquetes se incluyen de acuerdo al supuesto de que ya están instalados en R. En el caso de que no estén instalados, usar el comando install.package(“nombre del paquete”). También, es posible instalar los paquetes en el menú Herramientas (tools) –> instalar paquetes)

library(spdep)

## Loading required package: sp

## Loading required package: spData

## To access larger datasets in this package, install the spDataLarge
## package with: `install.packages('spDataLarge',
## repos='https://nowosad.github.io/drat/', type='source')`

## Loading required package: sf

## Warning: package 'sf' was built under R version 3.6.2

## Linking to GEOS 3.7.2, GDAL 2.4.2, PROJ 5.2.0

library(maptools)

## Checking rgeos availability: TRUE

library(RColorBrewer)
library(leaflet)
library(dplyr)

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

library(ggplot2)
library(tmap)

## Warning: package 'tmap' was built under R version 3.6.2

library(tmaptools)

## Warning: package 'tmaptools' was built under R version 3.6.2

También, fijamos el directorio de trabajo.

Información Espacial

Para agregar información espacial de tipo Shapes (en este caso, es un polígono) al ambiente global de trabajo debemos usar el siguiente código:

france.shp<-readShapePoly('guerry.shp')

## Warning: readShapePoly is deprecated; use rgdal::readOGR or sf::st_read

class(france.shp)

## [1] "SpatialPolygonsDataFrame"
## attr(,"package")
## [1] "sp"

Podemos visualizar el mapa que usaremos en el laboratorio:

Mapa de Francia: Guerry

Análisis Exploratorio Preliminar de los Datos.

Podemos usar otra interfaz para ver el mapa, a través del paquete leaflet

leaflet(france.shp) %>%
  addPolygons(stroke = FALSE, fillOpacity = 0.5, smoothFactor = 0.5) %>%
  addTiles() #adds a map title, the default is OpenStreetMap

Mapa de quantiles

spplot(france.shp, "suicids", at=quantile(france.shp@data$suicids, p = c(0, .25, .5, .75, 1), na.rm = TRUE), col.regions = brewer.pal(9, "Blues"))

spplot(france.shp, "crm_prp", at=quantile(france.shp@data$crm_prp, p = c(0, .25, .5, .75, 1), na.rm = TRUE), col.regions = brewer.pal(9, "Blues"))

Usando el paquete leaflet, podemos usar un mapa de quantiles

qpal<-colorQuantile("OrRd", france.shp@data$suicids, n=9) 
leaflet(france.shp) %>%
  addPolygons(stroke = FALSE, fillOpacity = .8, smoothFactor = 0.2, color = ~qpal(suicids)
  ) %>%
  addTiles()

Otra versión del mapa de quantiles puede ser a través del paquete tmap

tm_shape(france.shp) + 
  tm_fill("suicids",
          palette = "Reds", 
          style = "quantile", 
          title = "Number of suicids") +
  tm_borders(alpha=.4)

## Warning: Currect projection of shape france.shp unknown. Long-lat (WGS84)
## is assumed.

Matrices de Peso Espacial (Matrices W)

La matriz de peso espacial, o la también conocida como la matriz W, representa el proceso de influencia que tienen las observaciones. Esta matriz es definida de forma exógena por parte del (la) investigador/a. Por tanto, la matriz W representa la teoría que asume el(la) investigador/a acerca del proceso de influencia de las observaciones, para representar el proceso de contagio de una manera sencilla.

La matriz W es el corazón de la Econometría Espacial. La definición más simple es la siguiente:

\[\begin{align*} W = \begin{bmatrix} w_{11} & w_{12} & \cdots & w_{1n} \\ w_{21} & w_{22} & \cdots &w_{2n} \\ \vdots & \vdots & \ddots & \vdots \\ w_{n1} & w_{n2} & \cdots & w_{nn} \end{bmatrix} \end{align*} \] Donde cada elemento de esta matriz de conectividad está definido como:

\[ w_{ij}=\begin{cases} 1, \mbox{si } j \in N(i) \\ 0, \mbox{en el otro caso} \end{cases} \] Donde \(N(i)\) es el conjunto de vecinos de la localización i. De esta forma, la matriz W es una matriz cuadrada no estocástica cuyos elementos \(w_{ij}\) reflejan la intensidad de la interdependencia existente en cada par de observaciones. Por convención, asumimos que \(w_{ii}=0\)

Para definir la matriz W, existen diversos criterios para definir el conjunto de vecinos de una localización, \(N(i)\). Estos criterios dependerán de la forma en la cual el/la investigador/a cree que se lleva a cabo el proceso de interacción de las observaciones a analizar. Habitualmente, se recurre a criterios de contigüidad física para establecer la interacción (interdependencia) de las observaciones.

Es necesario recalcar que no solo existen criterios de contigüidad espacial. Algunas de estas matrices se basan en el criterio que el/la investigador/a considere prudente para representar las relaciones de interacción de las observaciones, de acuerdo a la evidencia empírica o teórica determinada para reflejar de manera realista la relación espacial. Estos criterios pueden ser por similaridad en ingresos, por tipo de población (por ejemplo, afroamericana), por intensidad en los flujo de comercio, etc.

En este laboratorio, revisaremos los criterios más usados para establecer la dependencia espacial de las observaciones. Este tipo de matrices son:

Matriz Tipo Torre (rook)

La matriz tipo torre (o rook) define como vecinos de una observación \(i\) si las otras observaciones comparten algún lado de \(i\). Este tipo de matriz se construye en base a los movimientos de la torre en un tablero de ajedrez.

Para construir la matriz W tipo rook, debemos usar la siguiente función:

rook.fr<-poly2nb(france.shp, row.names = france.shp$dept, queen = FALSE)
summary(rook.fr)

## Neighbour list object:
## Number of regions: 85 
## Number of nonzero links: 420 
## Percentage nonzero weights: 5.813149 
## Average number of links: 4.941176 
## Link number distribution:
## 
##  2  3  4  5  6  7  8 
##  7 11 12 16 30  7  2 
## 7 least connected regions:
## 25 29 57 62 66 68 75 with 2 links
## 2 most connected regions:
## 49 77 with 8 links

A través del resumen podemos revisar algunos detalles del vecindario, como por ejemplo, la cantidad de vecinos que tiene cada observación, etc. También, podemos revisar la red de interconexión de las observación a través del siguiente gráfico.

plot(france.shp, border="grey")
plot(rook.fr,  coordinates(france.shp),  add=TRUE,  col="red")

La Estandarización de la Matriz W

La estandarización de la matriz W corresponde a la siguiente operación por filas:

\[ w_{ij}^{*}=\dfrac{w_{ij}}{\sum_{j=1}^n w_{ij}} \] Donde \(w_{ij}^{*} \in W^{*}\). El resultado de la estandarización permite que cada fila de la matriz \(W^{*}\) sume la unidad.

A pesar de que estandarización no es necesaria, en algunas instancias esta estandarización es útil. Por ejemplo, a través de la matriz de peso espacial estandarizada, podemos definir lo siguiente:

\[L(y)=Wy\] Donde cada elemento es igual a:

\[ L(y_{i})=\sum_{j=1}^n w_{ij}^{*}y=\dfrac{w_{ij}y_{j}}{\sum_{j=1}^n w_{ij}} \]

Esta operación representa el promedio de la variable y observada en todas las localizaciones que son vecinas a la localización \(i\). Esto dependerá del criterio W escogido.

Lo anterior, también se llama rezago (o retardo) espacial.

El Rezago (retardo) Espacial

Esta operación representa al promedio ponderado de los valores de las localizaciones vecinas, con ponderaciones fijas y dadas de forma exógena. El operador de rezago espacial se obtiene como el producto de la matriz de peso espacial por el vector de observaciones de una variable \(y\). Es decir,

\[L(y)=Wy\] Cada elemento de una variable rezagada (o retardada) espacialmente es igual a:

\[ \sum_{j}^{n}w_{ij}y_{j} \]

Cada elemento del rezago espacial es igual a un promedio ponderado de los valores de la variable en el subgrupo de observaciones vecinas.

La Matriz Estandarizada por Filas en R

En R podemos distinguir ambas formas de construir la matriz de peso espacial a través de la siguiente función:

rook.fr.bi<-listw2mat(nb2listw(rook.fr, style ="B")) #Binary Matrix#
rook.fr.st<-listw2mat(nb2listw(rook.fr, style ="W")) #Row Standarized Matrix#

La función nos permite usar distintos estilos de la matriz W.

Volviendo a nuestro ejemplo, usando la matriz W tipo torre, podemos corroborar que se cumplen los estilos calculados, revisando la suma por filas de cada una de ellas.

rowSums(rook.fr.bi)

##  1  2  3  4  5  7  8  9 10 11 12 13 14 15 16 17 18 19 21 22 23 24 25 26 27 
##  4  6  6  4  3  7  3  3  5  5  7  3  3  6  5  5  6  6  7  3  6  7  2  5  6 
## 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 
##  6  2  6  6  6  4  4  6  6  5  6  5  4  6  7  5  4  7  6  6  5  8  4  6  6 
## 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 75 76 77 78 79 
##  5  4  5  4  2  6  3  6  6  2  6  3  3  2  4  2  4  5  7  6  2  3  8  6  5 
## 80 81 82 83 84 85 86 87 88 89 
##  5  5  6  3  6  4  6  6  6  5

rowSums(rook.fr.st)

##  1  2  3  4  5  7  8  9 10 11 12 13 14 15 16 17 18 19 21 22 23 24 25 26 27 
##  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1 
## 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 
##  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1 
## 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 75 76 77 78 79 
##  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1 
## 80 81 82 83 84 85 86 87 88 89 
##  1  1  1  1  1  1  1  1  1  1

Matriz tipo reina (queen)

La matriz tipo reina (queen) define como vecinas de una observación \(i\) a las observacionese que comparten algún lado o vértice con \(i\). La construcción de esta matriz se basa en los movimientos de la reina en un tablero de ajedrez.

Para construir la matriz W tipo reina, debemos usar la siguiente función:

queen.fr<-poly2nb(france.shp, row.names = france.shp$dept, queen = TRUE)
summary(queen.fr)

## Neighbour list object:
## Number of regions: 85 
## Number of nonzero links: 420 
## Percentage nonzero weights: 5.813149 
## Average number of links: 4.941176 
## Link number distribution:
## 
##  2  3  4  5  6  7  8 
##  7 11 12 16 30  7  2 
## 7 least connected regions:
## 25 29 57 62 66 68 75 with 2 links
## 2 most connected regions:
## 49 77 with 8 links

Tal como lo hicimos en la matriz anterior, podemos graficar su vecindario:

plot(france.shp, border="grey")
plot(queen.fr,  coordinates(france.shp),  add=TRUE,  col="red")

Matriz basada en K-vecinos

Otro tipo de matriz de peso espacial está basada en definir una cantidad de vecinos fija. Este tipo de matriz de contiguidad se puede calcular a través del siguiente procedimiento:

Recuperamos los centroides de cada polígono (ej: departamento) de nuestro archivo vectorial (shp)

coords<-coordinates(france.shp)

Usamos la siguiente función para calcular la matriz de peso espacial. Por ejemplo, para 1, 2, 3, 4, 5 vecinos:

kn1<-knn2nb(knearneigh(coords,k=1))
kn2<-knn2nb(knearneigh(coords,k=2))
kn3<-knn2nb(knearneigh(coords,k=3))
kn4<-knn2nb(knearneigh(coords,k=4))
kn5<-knn2nb(knearneigh(coords,k=5))

A través del siguiente loop, podemos calcular matrices de peso espacial para una cantidad mayor de vecinos;

n<-20
for(i in 1:n) {
  w<-knn2nb(knearneigh(coords,k=i))
  sum<-rowSums(listw2mat(nb2listw(w)))
  print(sum)
}

Para visualizar el vecindario que se construye a través de este criterio, podemos graficar estas matrices de la siguiente forma:

plot(france.shp, border="grey")
plot(kn1,  coordinates(france.shp),  add=TRUE,  col="red")

plot(france.shp, border="grey")
plot(kn2,  coordinates(france.shp),  add=TRUE,  col="red")

plot(france.shp, border="grey")
plot(kn3,  coordinates(france.shp),  add=TRUE,  col="red")

plot(france.shp, border="grey")
plot(kn4,  coordinates(france.shp),  add=TRUE,  col="red")

plot(france.shp, border="grey")
plot(kn5,  coordinates(france.shp),  add=TRUE,  col="red")

W basada en distancias

En la literatura, una de las matrices W basada en distancias más utilizadas en la matriz basada en distancias euclideanas.

Para construir matrices basadas en distancias euclideanas, podemos usar la siguiente función:

dist.mat<-as.matrix(dist(coords, method="euclidean"))
dist.mat.inv<-1/dist.mat
diag(dist.mat.inv)<-0
dist.mat.inv

Para obtener la matriz W, basada en distancias usamos la siguiente función:

dist.W<-mat2listw(dist.mat.inv, row.names = france.shp$dept, style="W")
summary(dist.W)

## Characteristics of weights list object:
## Neighbour list object:
## Number of regions: 85 
## Number of nonzero links: 7140 
## Percentage nonzero weights: 98.82353 
## Average number of links: 84 
## Link number distribution:
## 
## 84 
## 85 
## 85 least connected regions:
## 1 2 3 4 5 7 8 9 10 11 12 13 14 15 16 17 18 19 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 with 84 links
## 85 most connected regions:
## 1 2 3 4 5 7 8 9 10 11 12 13 14 15 16 17 18 19 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 with 84 links
## 
## Weights style: W 
## Weights constants summary:
##    n   nn S0       S1      S2
## W 85 7225 85 3.206611 341.364

Comprobamos que la suma de las filas de la matriz sean iguales a 1:

dist.W<-listw2mat(dist.W)
rowSums(dist.W)

Medidas de Autocorrelación Espacial Global

El Diagrama de Moran: El Camino Largo

El diagrama de Moran (Moran scatterplot) es una forma rápida para medir la estructura espacial de las observaciones de una variable bajo estudio. Este diagrama se construye centrando la variable bajo estudio de acuerdo a su media (aunque no es necesario) y el promedio de los valores de las observaciones vecinas, \(Wy\). En este caso, la matriz \(W\) está estandarizada (normalizada) por filas. Al centrar la variable de interés a su media, podemos delinear los cuadrantes del diagrama de acuerdo a \(y=0\) y \(Wy=0\).

También, podemos dibujar el grado de ajuste de la relación entre la variable espacial bajo estudio y su rezago espacial a través de una línea de regresión. Con ello podemos ver con mayor claridad la fuerza de la relación entre la variable de interés y su rezago espacial.

Si las observaciones están distribuidas de manera aleatoria en el espacio, no se podría observar una relación particular entre \(y\) y \(Wy\). En este caso, la pendiente de la relación debería ser igual a cero.

En el caso contrario, el parámetro de regresión es distinto de cero. De esta forma, cada cuadrante se ajustaría a un tipo de asociación espacial específica.

Para construir manualmente este diagrama, podemos seguir el siguiente procedimiento:

Crear una matriz W estandarizada por filas (no estocástica por filas)

queen.v2<-nb2listw(poly2nb(france.shp, row.names = france.shp$dept, queen = TRUE),style ="W",zero.policy = TRUE)
summary(queen.v2)

## Characteristics of weights list object:
## Neighbour list object:
## Number of regions: 85 
## Number of nonzero links: 420 
## Percentage nonzero weights: 5.813149 
## Average number of links: 4.941176 
## Link number distribution:
## 
##  2  3  4  5  6  7  8 
##  7 11 12 16 30  7  2 
## 7 least connected regions:
## 25 29 57 62 66 68 75 with 2 links
## 2 most connected regions:
## 49 77 with 8 links
## 
## Weights style: W 
## Weights constants summary:
##    n   nn S0      S1       S2
## W 85 7225 85 37.2761 347.6683

queen.v2$weights[1]

## [[1]]
## [1] 0.25 0.25 0.25 0.25

Centramos la varible a estudiar a su media.

var<-scale(france.shp$suicids, center = TRUE, scale = FALSE)
var

Calculamos el rezago espacial de la variable bajo estudio.

var.lag<-lag.listw(queen.v2, var)
var.lag

Para estimar la línea de regresión para visualizar el grado de ajuste, calculamos la siguiente regresión.

mor<-lm(var.lag ~ var)
summary(mor)

## 
## Call:
## lm(formula = var.lag ~ var)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -27249 -11968  -4321   7655  61386 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -72.6279  1900.6986  -0.038     0.97    
## var            0.4017     0.0607   6.617 3.36e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 17520 on 83 degrees of freedom
## Multiple R-squared:  0.3454, Adjusted R-squared:  0.3375 
## F-statistic: 43.79 on 1 and 83 DF,  p-value: 3.361e-09

Finalmente, podemos graficar el diagrama de moran de la siguiente forma:

plot(var.lag ~ var,  pch=20, asp=1, las=1)
abline(mor, lwd=2)
abline(h=mean(var), lt=2)
abline(v=mean(var.lag), lt=2)

Donde la pendiente de la curva de regresión nos indica la relación líneal entre ambas variables.

El Índice de Moran

El índice de moran, o estadístico de moran mide la relación lineal entre la variable de interés y sus valores espacialmente rezagados. Es decir,

\[ corr(y, Wy)= \dfrac{Cov(y, Wy)}{\sqrt{Var(y)Var(Wy)}} \] De acuerdo a esta ecuación, usamos la expresión del rezago espacial. También, asumimos que bajo el supuesto de estacionaridad (el mismo de series de tiempo estacionarias):

\[ Var(y)=Var(Wy) \] Por tanto, este coeficiente de correlación es igual a:

\[ corr(y, Wy)= \dfrac{Cov(y, Wy)}{Var(y)} \] En términos empíricos, el índice de moran puede estimarse de la siguiente forma:

\[ I= \dfrac{\sum_{i=1}^{n} \sum_{j=1}^n w_{ij}z_{i}z_{j}}{\sum_{i=1}^{n}z_{i}^{2}} \] donde \(n\) es el número de observaciones, \(z_{i}\) es el valor de la variable para la observación \(i\), \(z_{j}\) es el valor de la variable para la observación \(j\). Tanto \(z_{i}\) y \(z_{j}\) están centradas en la media. Finalmente, \(w_{ij}\) es el elemento \(ij\) de la matriz de peso espacial estandarizada por filas.

En términos matriciales, esta formula puede expresarse de la siguiente forma: \[ I=\dfrac{z^{T}Wz}{z^{T}z} \] Donde \(W\) es la matriz de peso espacial estandarizada por filas.

Tal como puede notar, este índice es igual a los estimados que se obtienen a través de OLS. Por tanto, la pendiente de la regresión entre el rezago espacial de la variable bajo estudio (y) y la variable bajo estudio (x) es igual al índice de moran.

El índice de moran se encuentra en el intervalo entre \([-1,1]\), De esta forma, podemos notar tres situaciones:

Cuando la variable z está distribuida aleatoriamente, el índice de moran se encuentra cercano a 0.
Cuando el valor observado del índice de moran es positivo, indica la existencia de autocorrelación positiva.
Cundo el valor observado del índice moran es negativo, indica la existencia de autocorrelación negativa.

Contraste de Hipótesis para el Índice de Moran

El contraste de hipótesis para determinar la significancia del índice de moran se lleva a cabo teniendo en cuenta lo siguiente:

\[ H_{0}: \mbox{Las observaciones se distruyen aleatoriamente. No hay autocorrelación espacial} \\ H_{1}: \mbox{Existe autocorrelación espacial} \] Para llevar a cabo el contraste de hipótesis, podemos usar dos formas:

Bajo Aleatorización.
Bajo Normalidad.

Contraste de Hipótesis bajo Aleatorización

Este método no asume una distribución a priori de la variable bajo estudio. Por tanto, el estadístico de moran estimado usando los datos es comparado con la distribución de las variables derivadado de una reorganización de los datos de forma aleatoria. Si la hipótesis nula es cierta, todas las posibles combinaciones de datos son equiprobales (todos los valores son igualmente probales).

El procedimiento para realizar este constraste de hipótesis es el siguiente:

Creamos la distribución de índices de moran, el cual corresponde a cada uno de los valores permutados en el área de análisis. En otras palabras, por cada permutación de valores, obtenemos un índice de moran por cada permtuación. Esto lo podemos crear de la siguiente forma:

n <- 599   # Define the number of simulations
I.r <- vector(length=n)  # Create an empty vector

for (i in 1:n){
  # Randomly shuffle Suicids Values
  x <- sample(var, replace=FALSE)
  # Compute a new set of lagged values
  x.lag <- lag.listw(queen.v2, x)
  # Compute the regression slope and store its value
  M.r    <- lm(x.lag ~ x)
  I.r[i] <- coef(M.r)[2]
}

Graficamos la distribución obtenida y comparamos con el valor estimado con la relación espacial determinada a través de la matriz W.

hist(I.r, main=NULL, xlab="Moran's I", las=1, xlim =c(-0.2,0.5))
abline(v=coef(mor)[2], col="red")

Calculamos el p-valor.

Obtenemos el número de valores del índice de moran simulados que son más grandes que nuestro valor observado (un lado):

N.greater <- sum(coef(mor)[2] > I.r)
p <- min(N.greater + 1, n + 1 - N.greater) / (n + 1)
p

## [1] 0.001666667

En este caso, rechazamos la hipótesis nula en favor de la alternativa. Por tanto, existe autocorrelación espacial.

Contraste de Hipótesis bajo Normalidad

Para calcular esta hipótesis bajo normalidad, usamos

moran.test(france.shp$suicids, listw = queen.v2,randomisation = FALSE)

## 
##  Moran I test under normality
## 
## data:  france.shp$suicids  
## weights: queen.v2    
## 
## Moran I statistic standard deviate = 5.9281, p-value = 1.532e-09
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       0.401681223      -0.011904762       0.004867397

El Diagrama de Moran: El Camino Corto

plot_mor<-moran.plot(france.shp$suicids, listw = queen.v2, xlab="Suicids", ylab = "Lag.Suicids", pch=19, labels = FALSE)

El Índice de Moran bajo Aleatorización

moran.test(france.shp$suicids, listw = queen.v2,  randomisation = TRUE)

## 
##  Moran I test under randomisation
## 
## data:  france.shp$suicids  
## weights: queen.v2    
## 
## Moran I statistic standard deviate = 6.0884, p-value = 5.704e-10
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       0.401681223      -0.011904762       0.004614565

mc<-moran.mc(france.shp$suicids, queen.v2, nsim=599)
plot(mc)

El Índice de Moran bajo Normalidad

moran.test(france.shp$suicids, listw = queen.v2,randomisation = FALSE)

## 
##  Moran I test under normality
## 
## data:  france.shp$suicids  
## weights: queen.v2    
## 
## Moran I statistic standard deviate = 5.9281, p-value = 1.532e-09
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       0.401681223      -0.011904762       0.004867397

El Índice de Moran Local

El diagrama de moran que presentamos anteriormente, no entrega indicaciones acerca de la significancia de los clústeres espaciales detectados. Por tanto, necesitamos calcular el estadístico LISA para cada observación. De esta forma, obtenemos un indicador de cuales son los clústeres espaciales de valores similares significativos.

El estadístico LISA es usado para la detección de clústeres espaciales, como también para la inestabilidad local, outliers significativos y regímenes espaciales. El estadístico LISA puede ser calculado de la siguiente forma:

\[ I_{i}=\dfrac{z_{i}}{m_{0}}\sum_{j=1}^{n}w_{ij}z_{j} \] Donde

\[ m_{0}=\sum_{i=1}^{n}\dfrac{z_{i}^2}{n} \] Donde \(z_{i}\) es el valor de la observación \(i\), \(z_{i}\) es el valor de la observación \(j\). Ambas variables están centradas de acuerdo a su media. Finalmente, \(w_{ij}\) es la matriz estandarizada por filas.

La suma del estadístico de LISA para cada observación es proporcional al índice de moran global. En el caso de la matriz W sea estandarizada por filas, la media del LISA es igual al índice de moran global.

Para calcular el LISA, podemos usar la siguiente función:

local.lisa<-localmoran(france.shp$suicids, listw = queen.v2, zero.policy= TRUE)
local.lisa

##               Ii        E.Ii    Var.Ii         Z.Ii    Pr(z > 0)
## 1   0.0136137856 -0.01190476 0.2240847  0.053907578 4.785044e-01
## 2   0.5820875241 -0.01190476 0.1459577  1.554774093 5.999994e-02
## 3   1.0460471078 -0.01190476 0.1459577  2.769187745 2.809812e-03
## 4   0.4185274169 -0.01190476 0.2240847  0.909282006 1.816006e-01
## 5   0.2436267784 -0.01190476 0.3022118  0.464824307 3.210286e-01
## 7   0.2043366826 -0.01190476 0.1236357  0.614988505 2.692812e-01
## 8   0.2628357962 -0.01190476 0.3022118  0.499766445 3.086198e-01
## 9   1.5082782515 -0.01190476 0.3022118  2.765286878 2.843638e-03
## 10  0.6167048149 -0.01190476 0.1772085  1.493271238 6.768310e-02
## 11  0.8248384858 -0.01190476 0.1772085  1.987695816 2.342267e-02
## 12  0.6779861082 -0.01190476 0.1236357  1.962042734 2.487875e-02
## 13  0.5684480090 -0.01190476 0.3022118  1.055689932 1.455549e-01
## 14  0.0102867616 -0.01190476 0.3022118  0.040367461 4.839001e-01
## 15  2.1278842798 -0.01190476 0.1459577  5.600895240 1.066238e-08
## 16  0.1098563322 -0.01190476 0.1772085  0.289245259 3.861968e-01
## 17  0.0366995421 -0.01190476 0.1772085  0.115460235 4.540402e-01
## 18 -0.1546905050 -0.01190476 0.1459577 -0.373741511 6.457017e-01
## 19  0.2655846818 -0.01190476 0.1459577  0.726328285 2.338188e-01
## 21  0.2569612259 -0.01190476 0.1236357  0.764652180 2.222393e-01
## 22 -0.0659227290 -0.01190476 0.3022118 -0.098261311 5.391376e-01
## 23  0.6927128883 -0.01190476 0.1459577  1.844335851 3.256710e-02
## 24  0.0018394710 -0.01190476 0.1236357  0.039088461 4.844099e-01
## 25  0.0024610032 -0.01190476 0.4584658  0.021216575 4.915364e-01
## 26  0.1148454777 -0.01190476 0.1772085  0.301097047 3.816702e-01
## 27  0.4666447540 -0.01190476 0.1459577  1.252602782 1.051752e-01
## 28  0.4114399850 -0.01190476 0.1459577  1.108104365 1.339084e-01
## 29 -0.2100727612 -0.01190476 0.4584658 -0.292671236 6.151133e-01
## 30 -0.0595094163 -0.01190476 0.1459577 -0.124605126 5.495819e-01
## 31  0.9935909160 -0.01190476 0.1459577  2.631883726 4.245646e-03
## 32  0.7361212523 -0.01190476 0.1459577  1.957957191 2.511752e-02
## 33  0.0854244807 -0.01190476 0.2240847  0.205606675 4.185491e-01
## 34 -0.1791072281 -0.01190476 0.2240847 -0.353212890 6.380356e-01
## 35  0.0510900516 -0.01190476 0.1459577  0.164888849 4.345157e-01
## 36  0.0718457823 -0.01190476 0.1459577  0.219216949 4.132405e-01
## 37  0.2542255561 -0.01190476 0.1772085  0.632196461 2.636293e-01
## 38  0.0001299173 -0.01190476 0.1459577  0.031500758 4.874351e-01
## 39  0.0120012066 -0.01190476 0.1772085  0.056788978 4.773566e-01
## 40 -0.0114005109 -0.01190476 0.2240847  0.001065223 4.995750e-01
## 41  0.3877242110 -0.01190476 0.1459577  1.046028355 1.477740e-01
## 42  1.1571288785 -0.01190476 0.1236357  3.324719980 4.425369e-04
## 43  3.0478014262 -0.01190476 0.1772085  7.268376775 1.819164e-13
## 44 -0.0814803820 -0.01190476 0.2240847 -0.146977532 5.584251e-01
## 45  0.5533853532 -0.01190476 0.1236357  1.607679433 5.395270e-02
## 46  0.3236681009 -0.01190476 0.1459577  0.878361564 1.898738e-01
## 47  0.0101618381 -0.01190476 0.1459577  0.057759299 4.769702e-01
## 48 -1.3251097079 -0.01190476 0.1772085 -3.119537545 9.990943e-01
## 49  0.0130621989 -0.01190476 0.1068942  0.076363994 4.695648e-01
## 50 -0.0324433362 -0.01190476 0.2240847 -0.043387453 5.173037e-01
## 51  0.6166514571 -0.01190476 0.1459577  1.645245146 4.995963e-02
## 52  0.2804434701 -0.01190476 0.1459577  0.765221145 2.220699e-01
## 53 -0.0248244342 -0.01190476 0.1772085 -0.030690870 5.122420e-01
## 54  0.2946277479 -0.01190476 0.2240847  0.647545674 2.586394e-01
## 55  0.3751963956 -0.01190476 0.1772085  0.919564458 1.789002e-01
## 56 -0.0157411810 -0.01190476 0.2240847 -0.008104382 5.032331e-01
## 57  0.2163270396 -0.01190476 0.4584658  0.337071997 3.680313e-01
## 58  0.0273756713 -0.01190476 0.1459577  0.102816487 4.590543e-01
## 59  0.5277156850 -0.01190476 0.3022118  0.981595852 1.631495e-01
## 60  0.8282994030 -0.01190476 0.1459577  2.199233390 1.393067e-02
## 61  0.0189782980 -0.01190476 0.1459577  0.080836372 4.677860e-01
## 62  0.4990923144 -0.01190476 0.4584658  0.754683633 2.252194e-01
## 63  2.4221251746 -0.01190476 0.1459577  6.371070427 9.385671e-11
## 64  1.3566836029 -0.01190476 0.3022118  2.489528835 6.395627e-03
## 65  2.8092573300 -0.01190476 0.3022118  5.131831131 1.434684e-07
## 66  0.0791883876 -0.01190476 0.4584658  0.134534055 4.464901e-01
## 67  0.2307791702 -0.01190476 0.2240847  0.512666440 3.040923e-01
## 68  0.1666504315 -0.01190476 0.4584658  0.263705388 3.960035e-01
## 69 -0.0935156848 -0.01190476 0.2240847 -0.172401942 5.684392e-01
## 70 -0.0252228997 -0.01190476 0.1772085 -0.031637431 5.126194e-01
## 71 -0.1292555647 -0.01190476 0.1236357 -0.333744509 6.307138e-01
## 72  0.0967379966 -0.01190476 0.1459577  0.284372290 3.880625e-01
## 75  1.0440696013 -0.01190476 0.4584658  1.559552111 5.943288e-02
## 76  0.7087271602 -0.01190476 0.3022118  1.310864534 9.495177e-02
## 77  0.8275783894 -0.01190476 0.1068942  2.567644767 5.119601e-03
## 78  0.9094045041 -0.01190476 0.1459577  2.411525894 7.942962e-03
## 79  0.0412973257 -0.01190476 0.1772085  0.126382337 4.497146e-01
## 80  0.6037388894 -0.01190476 0.1772085  1.462470493 7.180616e-02
## 81  0.8999965039 -0.01190476 0.1772085  2.166234786 1.514663e-02
## 82  0.3439751914 -0.01190476 0.1459577  0.931515349 1.757935e-01
## 83  0.5363062321 -0.01190476 0.3022118  0.997222476 1.593283e-01
## 84  0.2638143537 -0.01190476 0.1459577  0.721694453 2.352412e-01
## 85 -0.3535129201 -0.01190476 0.2240847 -0.721642495 7.647428e-01
## 86  0.1538340325 -0.01190476 0.1459577  0.433821095 3.322092e-01
## 87 -0.0076029488 -0.01190476 0.1459577  0.011259991 4.955080e-01
## 88  0.0538509752 -0.01190476 0.1459577  0.172115563 4.316733e-01
## 89  0.5176850622 -0.01190476 0.1772085  1.258048368 1.041871e-01
## attr(,"call")
## localmoran(x = france.shp$suicids, listw = queen.v2, zero.policy = TRUE)
## attr(,"class")
## [1] "localmoran" "matrix"

Comprobamos que la media del índice LISA es igual al índice de moran global.

i_moran<-mean(local.lisa[,1])
i_moran

## [1] 0.4016812

Índice LISA (índice de moran local)

moran.map <- cbind(france.shp, local.lisa) #agregamos la variable a la base de datos
tm_shape(moran.map) +
  tm_fill(col = "Ii",
          style = "quantile",
          title = "Local Moran Statistic")

## Warning: Currect projection of shape moran.map unknown. Long-lat (WGS84) is
## assumed.

## Variable(s) "Ii" contains positive and negative values, so midpoint is set to 0. Set midpoint = NA to show the full spectrum of the color palette.

Identificación de Clusters Espaciales

quadrant <- vector(mode="numeric",length=nrow(local.lisa))

Centramos la variable alrededor de su media.

m.suicids <- france.shp$suicids - mean(france.shp$suicids)

Centramos el Moran Local alteredor de su media

m.local <- local.lisa[,1] - mean(local.lisa[,1])

Determinamos el grado de significancia (threshold)

signif <- 0.1

Se construye el cuadrante

quadrant[m.suicids >0 & m.local>0] <- 4  
quadrant[m.suicids <0 & m.local<0] <- 1      
quadrant[m.suicids <0 & m.local>0] <- 2
quadrant[m.suicids >0 & m.local<0] <- 3
quadrant[local.lisa[,5]>signif] <- 0

Graficamos

brks <- c(0,1,2,3,4)
colors <- c("white","blue",rgb(0,0,1,alpha=0.4),rgb(1,0,0,alpha=0.4),"red")
plot(france.shp,border="lightgray",col=colors[findInterval(quadrant,brks,all.inside=FALSE)])
box()
legend("bottomleft", legend = c("insignificant","low-low","low-high","high-low","high-high"),
       fill=colors,bty="n")

Laboratorio 1: Introducción a la Econometría Espacial