Introducción a la Econometría - Matrices de Peso Espacial

Objetivo

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

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: spData

## The legacy packages maptools, rgdal, and rgeos, underpinning the sp package,
## which was just loaded, will retire in October 2023.
## Please refer to R-spatial evolution reports for details, especially
## https://r-spatial.org/r/2023/05/15/evolution4.html.
## It may be desirable to make the sf package available;
## package maintainers should consider adding sf to Suggests:.
## The sp package is now running under evolution status 2
##      (status 2 uses the sf package in place of rgdal)

## 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

## Linking to GEOS 3.11.0, GDAL 3.5.3, PROJ 9.1.0; sf_use_s2() is TRUE

library(maptools)

## Loading required package: sp

## Please note that 'maptools' will be retired during October 2023,
## plan transition at your earliest convenience (see
## https://r-spatial.org/r/2023/05/15/evolution4.html and earlier blogs
## for guidance);some functionality will be moved to 'sp'.
##  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)
library(tmaptools)
library(rgdal)

## Please note that rgdal will be retired during October 2023,
## plan transition to sf/stars/terra functions using GDAL and PROJ
## at your earliest convenience.
## See https://r-spatial.org/r/2023/05/15/evolution4.html and https://github.com/r-spatial/evolution
## rgdal: version: 1.6-7, (SVN revision 1203)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 3.5.3, released 2022/10/21
## Path to GDAL shared files: /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/rgdal/gdal
##  GDAL does not use iconv for recoding strings.
## GDAL binary built with GEOS: TRUE 
## Loaded PROJ runtime: Rel. 9.1.0, September 1st, 2022, [PJ_VERSION: 910]
## Path to PROJ shared files: /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/rgdal/proj
## PROJ CDN enabled: FALSE
## Linking to sp version:1.6-1
## To mute warnings of possible GDAL/OSR exportToProj4() degradation,
## use options("rgdal_show_exportToProj4_warnings"="none") before loading sp or rgdal.

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:

nepal.shp<-readShapePoly('nepal.shp')

## Warning: shapelib support is provided by GDAL through the sf and terra packages
## among others

class(nepal.shp)

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

Podemos visualizar el mapa que usaremos en el laboratorio:

Mapa de Nepal

Análisis Exploratorio Preliminar de los Datos.

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

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

Mapa de quantiles

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

tm_shape(nepal.shp) +
tm_polygons(col = "pcinc", title = "Ingresos Per Capita", style = "quantile") +
tm_layout(legend.outside = TRUE)

Usando el paquete leaflet, podemos usar un mapa de quantiles

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

tmap_mode("view") #Para ver el mismo mapa de manera interactiva#

## tmap mode set to interactive viewing

tm_shape(nepal.shp) +
  tm_polygons(col = "pcinc", title = "Nivel de Ingresos per capita", style = "quantile") +
  tm_layout(legend.outside = TRUE)

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 ejemplo, 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(nepal.shp, row.names = nepal.shp$id, queen = FALSE)
summary(rook.fr)

## Neighbour list object:
## Number of regions: 75 
## Number of nonzero links: 366 
## Percentage nonzero weights: 6.506667 
## Average number of links: 4.88 
## Link number distribution:
## 
##  2  3  4  5  6  7  9 
##  3 15 15 15 15  9  3 
## 3 least connected regions:
## 31 42 43 with 2 links
## 3 most connected regions:
## 18 20 39 with 9 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(nepal.shp, border="grey")
plot(rook.fr,  coordinates(nepal.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  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 
##  6  3  7  5  3  6  7  7  4  4  3  6  4  3  3  6  3  9  3  9  3  4  6  5  7  5 
## 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 
##  5  4  4  3  2  5  3  4  5  3  3  6  9  4  4  2  2  5  5  6  7  5  3  4  4  7 
## 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 
##  7  6  3  4  6  5  6  5  5  6  7  5  6  4  5  5  7  3  6  6  4  6  4

rowSums(rook.fr.st)

##  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 
##  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  1 
## 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 
##  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  1 
## 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 
##  1  1  1  1  1  1  1  1  1  1  1  1  1  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(nepal.shp, row.names = nepal.shp$id, queen = TRUE)
summary(queen.fr)

## Neighbour list object:
## Number of regions: 75 
## Number of nonzero links: 366 
## Percentage nonzero weights: 6.506667 
## Average number of links: 4.88 
## Link number distribution:
## 
##  2  3  4  5  6  7  9 
##  3 15 15 15 15  9  3 
## 3 least connected regions:
## 31 42 43 with 2 links
## 3 most connected regions:
## 18 20 39 with 9 links

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

plot(nepal.shp, border="grey")
plot(queen.fr,  coordinates(nepal.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:

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

coords<-coordinates(nepal.shp)

2. 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(nepal.shp, border="grey")
plot(kn1,  coordinates(nepal.shp),  add=TRUE,  col="red")

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

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

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

plot(nepal.shp, border="grey")
plot(kn5,  coordinates(nepal.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. Sin embargo, su cálculo requiere mayor tiempo y recursos computacionales. Por favor, verificar si su computador cumple con los requerimientos para correr esta clase de matrices

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 = nepal.shp$id, style="W")
summary(dist.W)

## Characteristics of weights list object:
## Neighbour list object:
## Number of regions: 75 
## Number of nonzero links: 5550 
## Percentage nonzero weights: 98.66667 
## Average number of links: 74 
## Link number distribution:
## 
## 74 
## 75 
## 75 least connected regions:
## 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 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 73 74 75 with 74 links
## 75 most connected regions:
## 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 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 73 74 75 with 74 links
## 
## Weights style: W 
## Weights constants summary:
##    n   nn S0       S1       S2
## W 75 5625 75 3.917943 301.1786

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

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

Con estas herramientas, ya podemos realizar un análisis para detectar la existencia de autocorrelación espacial.