La Autocorrelación Espacial se relaciona con el grado en el que la similitud de valores entre observaciones de un conjunto de datos se conecta con la similitud en la localización de dichas observaciones. De forma muy semejante a como ocurre con la Correlación Tradicional (que determina cómo los valores de una variable cambian en función de los valores de otra) y de forma análoga a su homólogo en una serie temporal (que relaciona el valor de una variable en algún punto del tiempo con los valores en periodos anteriores), la Autocorrelación Espacial relaciona los valores de la variable de interés en un lugar dado con los valores de ella misma en sus alrededores.
Una idea clave dentro de este contexto es el de Aleatoriedad Espacial, esto es, una situación en la cual la ubicación de una observación no guarda relación alguna con el valor de una variable; en otras palabras, una variable es aleatoria espacialmente si la distribución que sigue en el espacio parece no tener algún patrón discernible. Como tal, la Autocorrelación Espacial puede ser definida formalmente como la "absencia de Aleatoriedad Espacial", y, al igual que con la correlación tradicional, puede presentarse en dos variaciones:
En esta práctica, se aprenderá a explorar la Autocorrelación Espacial de un conjunto de datos determinado, analizando la presencia, naturaleza y robustez de los datos utilizados. Para esto, se utilizará un conjunto de herramientas colectivamente conocidos como Análisis Exploratorio de Datos Espaciales (ESDA, por sus siglas en inglés), diseñado específicamente para este propósito. El alcance de los Métodos ESDA es bastante amplio y abarca desde aproximaciones sencilllas, como Coropletas (Choropleths) y Consultas a Bases de Datos, hasta metodologías más avanzadas y robustas, como la Inferencia Estadística y la Identificación de la Dimensión Geográfica de los Datos. El propósito de esta sesión es explorar más a detalle estas últimas técnicas.
Los Métodos ESDA comúnmente se dividen en dos grupos principales: herramientas para analizar a la autocorrelación espacial a nivel global, y otras a nivel local. Las primeras consideran el patrón general que sigue la localización de los valores, y hacen posible realizar conclusiones sobre el grado de Agrupamiento (Clustering) que siguen los datos; pueden responderse preguntas como ¿Los valores en general siguen un patrón en particular en su distribución geográfica? o ¿Los valores similares se encuentran más cerca de lo que se encontrarían en una distribución aleatoria?. Se analizará la Autocorrelación Espacial Global utilizando el estadístico de I de Moran.
Los métodos utilizados para la Autocorrelación Espacial Local se enfocan en la inestabilidad espacial, esto es, que tanto se alejan ciertas ubicaciones en particular del patrón general; la idea detrás de esto es que, aún cuando existe un patrón general en términos de la naturaleza y robustez de la asociación espacial, algunas áreas en particular pueden divergir sustancialmente de éste mismo. Sin importar el grado general de concentración de los valores, existen focos donde valores inusualmente altos (o bajos) se encuentran cercanos a otros valores altos (o bajos), zonas que son conocidas como Hotspots (o Coldspots, para valores bajos); además, también es posible observar valores muy altos rodeados de valores muy bajos (o viceversa), los cuales reciben el nombre de Outliers Espaciales (Valores Atípicos). La ténica principal que se utilizará en esta sesión para explorar la Autocorrelación Espacial Local serán los Indicadores Locales de Asociación Espacial (LISA, por sus siglas en inglés).
Nota: Para instalar las dependencias en Colab, copien las siguientes líneas en una celda y ejecútenla.
# Important library for many geopython libraries
!apt install gdal-bin python-gdal python3-gdal
# Install rtree - Geopandas requirment
!apt install python3-rtree
# Install Geopandas
!pip install git+git://github.com/geopandas/geopandas.git
# Install descartes - Geopandas requirment
!pip install descartes
# Install Folium for Geographic data visualization
!pip install folium
# Install plotlyExpress
!pip install mapclassify
!pip install palettable
!pip install pysal
!pip install esda
%matplotlib inline
import seaborn as sns
import pandas as pd
from libpysal.weights import W,Queen, Rook, KNN, DistanceBand, min_threshold_distance, block_weights, lag_spatial
import esda
import geopandas as gpd
import numpy as np
import matplotlib.pyplot as plt
import mapclassify
np.random.seed(123)
En esta práctica, se continuarán utilizando las AGEB's como unidad espacial, así como la variable de Personas sin Derechohabiencia a Servicios Públic de Salud colectada por el Consejo Nacional de Evaluación de Política del Desarrollo Social (CONEVAL) en 2010.
NOTA: Bajar datos para trabajar en Colab
url = "https://www.dropbox.com/s/idewfb8v9stbt7b/pob_sinderechohab.zip?dl=1"
r = requests.get(url, allow_redirects=True)
open('/content/pob_sinderechohab.zip', 'wb').write(r.content)
Como es de costumbre, se creará una variable que contenga la ruta donde se encuentran los datos a utilizar, si van a trabajar en Colab, ajusten el valor:
f = 'data/'
Como se ha hecho siempre, se utiliza GeoPandas
para importar e inspeccionar los datos espaciales de la práctica:
# Importar el ShapeFile
agebs = gpd.read_file(f + 'pob_sinderechohab.shp')
# Graficar el GeoDataFrame
agebs.plot(cmap = 'Set1', figsize = (10,10))
<AxesSubplot:>
Antes de continuar con la parte analítica, resulta pertinente crear un Mapa de Coropletas para, de entrada, tener una idea visual de cómo es que se distribuye la variable en el espacio. A continuación, se muestra el código necesario para construir este mapa, en el cual cabe resaltar el uso del argumento scheme
dentro de la función .plot()
para indicar cómo deben de ser clasificados los datos, para después asignarles el color correspondiente:
# Preparación de la figura y sus filas
fig , filas = plt.subplots(1, figsize = (10,10))
# Gráfica del GeoDataFrame, generando un Mapa de Coropletas a través de 'scheme'
agebs.plot(ax = filas, column = 'p_sderech', scheme = 'fisherjenks' , cmap = 'YlGn', legend = True)
# Adición de título a la gráfica
fig.suptitle('Población sin Derechohabiencia a Servicios de Salud Públicos en la CDMX', size = 20)
# Mostrar la gráfica
plt.show()
Crea un mapa similar al anterior, pero utilizando como Métodos de Clasificación los Cuantiles (quantiles
) e Intervalor Iguales (equalinterval
), y responde las preguntas: ¿En qué difieren los mapas generados? ¿Cómo varía la interpretación de lo observado en cada uno de los casos? ¿Cuál es la distribución de los datos? Confirma esta última respuesta generando un Diagrama de Densidad de Kernal o un Histograma.
Como se ha discutido anteriormente, una Matriz de Pesos Espaciales es la forma en la que el espacio geográfico es codificado de forma numérica, a modo de que sea fácil para una computadora (o un método estadístico) el interpretarlo. Se han estudiado algunas de las varias formas conceptuales en las que puede definirse una Matriz de Pesos Espaciales, como los criterios de Contigüidad, Basados en Distancia, o de Bloques.
En este caso, se utilizará uno de los criterios más comúnes en este tipo de análisis: la Matriz de Contigüidad de Reina:
m_reina = Queen.from_dataframe(agebs, idVariable='cvegeo')
m_reina.transform = 'R'
('WARNING: ', '090090015012A', ' is an island (no neighbors)')
/home/plablo/miniconda3/envs/covid/lib/python3.10/site-packages/libpysal/weights/weights.py:172: UserWarning: The weights matrix is not fully connected: There are 8 disconnected components. There is 1 island with id: 090090015012A. warnings.warn(message)
Cabe recordar que, al modificar el parámetro transform
de la matriz, los valores de los pesos asignados se alteran; por ejemplo, al modificar este parámetro por 'R'
, se han estandarizado los pesos de modo que la suma de éstos sea igual a uno:
m_reina['090150001056A']
{'0901500010517': 0.25, '0901500010610': 0.25, '0901500010574': 0.25, '0901500010502': 0.25}
También necesitamos quitar la isla que detectamos, más adelante esta sla nos va a impedir realizar algunos cálculos
agebs = agebs.set_index('cvegeo').drop('090090015012A').reset_index()
Volvemos a calcular nuestra matriz de vecindad
m_reina = Queen.from_dataframe(agebs, idVariable='cvegeo')
m_reina.transform = 'R'
/home/plablo/miniconda3/envs/covid/lib/python3.10/site-packages/libpysal/weights/weights.py:172: UserWarning: The weights matrix is not fully connected: There are 7 disconnected components. warnings.warn(message)
Una vez que se encuentran listos los datos y la Matriz de Pesos Espaciales, es posible calcular el Rezago Espacial de la variable p_sderech
. Importante recordar que el rezago espacial es el resultado del producto entre la Matriz de Pesos Espaciales y una variable dada y que, si $W$ se encuentra estandarizada en filas, los resultados equivalen al valor promedio de una variable en la vecidad de una observación dada.
En la siguiente línea, se calcula el Rezago Espacial para la variable p_sderech
y es almacenada directamente en la tabla original:
agebs['rez_espacial'] = lag_spatial(m_reina , agebs['p_sderech'])
Puede compararse entre sí la apariencia del rezago espacial calculado con respecto a su variable original:
agebs[['p_sderech' , 'rez_espacial']].head()
p_sderech | rez_espacial | |
---|---|---|
0 | 389.0 | 337.0 |
1 | 323.0 | 419.5 |
2 | 448.0 | 484.5 |
3 | 240.0 | 409.0 |
4 | 1017.0 | 500.0 |
La forma en la que se interpreta rez_espacial
para una observación dada es la siguiente: el AGEB con Clave '0901500010235' (la segunda en la lista anterior) tiene una Población sin Derechohabiencia de 323 personas, y se encuentra rodeada por otras AGEB's que, en promedio, tienen 419.5 Pesonas sin Derechohabiencia.
Para ejemplificar lo anterior correctamente, pueden obtenerse primero los vecinos de la AGEB mencionada:
m_reina.neighbors['0901500010235']
['090150001024A', '0901500010076', '0901500010377', '0901500010362', '0901500010381', '0901500010220']
Y después, verificar los valores de las variables para estos vecinos:
vecinos = agebs.set_index('cvegeo').loc[m_reina.neighbors['0901500010235'] , 'p_sderech']
vecinos
cvegeo 090150001024A 226.0 0901500010076 1428.0 0901500010377 226.0 0901500010362 324.0 0901500010381 267.0 0901500010220 146.0 Name: p_sderech, dtype: float64
De forma sencilla, puede calcularse el valor promedio de estos valores, en cuyo caso podemos comprobar que es muy similar al Rezago Espacial calculado para dicha observación:
vecinos.mean()
436.1666666666667
Para algunas técnicas que se estudiarán más adelante, tiene más sentido el operar con la versión estandarizada de la variable, en lugar de sus valores brutos; estandarizar significa sencillamente el restar el valor de la media de la variable y dividirlo por la desviación estándar de cada observacion de la columna. Esto se consigue rápidamente a través de los siguientes comandos:
agebs['p_sderech_std'] = (agebs['p_sderech'] - agebs['p_sderech'].mean()) / agebs['p_sderech'].std()
agebs.head()
cvegeo | alcaldia | p_sderech | p_nescu | p_hacin | p_analf | geometry | rez_espacial | p_sderech_std | |
---|---|---|---|---|---|---|---|---|---|
0 | 0900700013628 | 09007 | 389.0 | 475.0 | 0.0 | 9.0 | POLYGON ((2810132.372 824698.172, 2810169.540 ... | 337.0 | -0.959260 |
1 | 0901500010235 | 09015 | 323.0 | 481.0 | 6.0 | 7.0 | POLYGON ((2798881.634 831643.241, 2798843.076 ... | 419.5 | -1.034395 |
2 | 0900200010097 | 09002 | 448.0 | 780.0 | 5.0 | 37.0 | POLYGON ((2792415.239 836846.390, 2792356.808 ... | 484.5 | -0.892094 |
3 | 0900200011184 | 09002 | 240.0 | 389.0 | 0.0 | 11.0 | POLYGON ((2792260.139 836768.777, 2792333.695 ... | 409.0 | -1.128883 |
4 | 0900300011285 | 09003 | 1017.0 | 1291.0 | 0.0 | 23.0 | POLYGON ((2802121.600 817466.682, 2802124.157 ... | 500.0 | -0.244339 |
Y, de forma análoga, debe de calcularse el Rezago Espacial para estos valores estandarizados:
agebs['rez_espacial_std'] = ps.lib.weights.lag_spatial(m_reina , agebs['p_sderech_std'])
agebs.head()
cvegeo | alcaldia | p_sderech | p_nescu | p_hacin | p_analf | geometry | rez_espacial | p_sderech_std | rez_espacial_std | |
---|---|---|---|---|---|---|---|---|---|---|
0 | 0900700013628 | 09007 | 389.0 | 475.0 | 0.0 | 9.0 | POLYGON ((2810132.372 824698.172, 2810169.540 ... | 337.0 | -0.959260 | -1.018458 |
1 | 0901500010235 | 09015 | 323.0 | 481.0 | 6.0 | 7.0 | POLYGON ((2798881.634 831643.241, 2798843.076 ... | 419.5 | -1.034395 | -0.924539 |
2 | 0900200010097 | 09002 | 448.0 | 780.0 | 5.0 | 37.0 | POLYGON ((2792415.239 836846.390, 2792356.808 ... | 484.5 | -0.892094 | -0.850542 |
3 | 0900200011184 | 09002 | 240.0 | 389.0 | 0.0 | 11.0 | POLYGON ((2792260.139 836768.777, 2792333.695 ... | 409.0 | -1.128883 | -0.936492 |
4 | 0900300011285 | 09003 | 1017.0 | 1291.0 | 0.0 | 23.0 | POLYGON ((2802121.600 817466.682, 2802124.157 ... | 500.0 | -0.244339 | -0.832897 |
La Autocorrelación Espacial Global se relaciona con el patrón geográfico general presente en un conjunto de datos. Los estadísticos diseñados para medir ésta caracterizan a los datos espaciales en términos de su grado de agrupamiento (clustering) y lo resumen. Este resumen puede ser visual o numérico; en esta sección, se estudiará un ejemplo de cada uno de éstos: el Gráfico de Moran, y el Estadístico de I de Moran de Autocorrelación Espacial
El Gráfico de Moran es una forma de visualizar un conjunto de datos espaciales para explorar la naturaleza y robustez de la Autocorrelación Espacial. En escencia se trata de un Diagrama de Dispersión tradicional en el que la variable de interés se contrasta contra su Rezago Espacial. Para poder interpretar los valores en términos de la media, y sus cantidades en términos de desviaciones estándar, la variable de interés usualmente se estandariza.
Desde el punto de vista técnico, crear un Gráfico de Moran en Python sigue un proceso muy similar al de calcular cualquier otro Diagrama de Dispersión, dado que se tenga la variable estandarizada y su Rezago Espacial de antemano:
# Preparación de la Gráfica
fig, filas = plt.subplots(1, figsize = (10,10))
# Generación del Diagrama de Dispersión
sns.regplot(x = 'p_sderech_std' , y = 'rez_espacial_std', data = agebs)
# Adición de lineas horizontal y vertical
plt.axvline(0, c='k', alpha=0.5)
plt.axhline(0, c='k', alpha=0.5)
# Muestra del resultado
plt.show()