En la primera parte de la práctica vimos cómo la escala de análisis tiene influencia sobre los resultados del mismo. En esta segunda parte vamos a ver como diferentes formas de agregar los datos, en una misma escala, también pueden tener una influencia sobre los resultados de un análisis.
Para esta segunda parte vamos a volver a trabajar con los datos de mezcla de usos de suelo en la Ciudad de México y vamos a tratar de estimar su efecto en la generación de viajes. Los datos de uso de suelo están calculados a partir del DENUE, mientras que la información sobre los viajes es de la Encuesta Origen Destino del 2007, por lo que la escala de análisis serán los distritos de tráfico de dicha encuesta.
Como siempre, el primer paso es leer los datos:
In [3]:
from geopandas import GeoDataFrame
datos = GeoDataFrame.from_file('data/distritos_variables.shp')
datos.head()
Out[3]:
El shapefile que leimos tiene columnas para cada uno de los tipos de uso de suelo que nos interesan (vivienda, comercio, servicios y ocio), adicionalmente tiene tres columnas con información de los viajes:
Como pueden recordar de la práctica anterior en que usamos estos datos, la mezcla de usos de suelo se mide utilizando el índice de entropía, en esta ocasión vmos a calcularlo de acuerdo a la siguiente fórmula:
\begin{equation} E = \sum\limits_{j}{\frac{p_{j}*ln(p_{j})}{ln(J)}} \end{equation}Donde $p_{j}$ representa la proporción del $j-ésimo$ uso de suelo con respecto al total y $J$ es el número de usos de suelo.
In [4]:
import numpy as np
intensidad = datos['comercio'] + datos['viv'] + datos['ocio'] + datos['servicios']
prop_comercio = datos['comercio'] / intensidad
prop_viv = datos['viv'] / intensidad
prop_ocio = datos['ocio'] / intensidad
prop_servicios = datos['servicios'] / intensidad
entropia = (prop_comercio*np.log(prop_comercio) + prop_viv*np.log(prop_viv) + prop_ocio*np.log(prop_ocio)
+ prop_servicios*np.log(prop_servicios))/np.log(4)
entropia.head()
Out[4]:
Lo que hicimos para calcular la entropía es relativamente fácil. Primero creamos la variable intensidad, que es la suma de todos los usos de suelo y luego fuimos calculando cada una de las proporciones, finalmente las sumamos y las dividimos por el logaritmo natural del número de usus de suelo. Lo que tenemos ahora es una Serie que contiene los datos de entropía, ahora lo que necesitamos es unir esa serie a nuestros datos originales:
In [5]:
datos['entropia'] = entropia
datos.head()
Out[5]:
Ahora que ya tenemos todos los datos, probemos un modelo que nos trate de explicar, por ejemplo, la cantidad de viajes que terminan en cada distrito. Lo primero que vamos a hacer, es explorar la colinearidad de las variables que tenemos:
In [6]:
corr = datos[['ocio','comercio','servicios','viv','entropia']].corr()
w, v = np.linalg.eig(corr)
print w
Parece ser que, si utilizamos las variables de arriba, vamos a tener problemas de colinearidad, entonces, para seleccionar las variables, observemos la matriz de correlación:
In [7]:
print corr
Como puden ver, la entropía está bastante correlacionada con las demás variables, sobre todo con comercio y ocio, seleccionesmos entonces, por lo pronto, entropía y vivienda como variables explicativas.
Ahora, antes de correr nuestro modelo, vamos a normalizar las variables porque tiene escalas de variación muy diferentes y ya sabemos que eso nos puede traer problemas:
In [8]:
variables = datos[['entrada','viv','entropia']]
vars_norm = (variables - variables.mean())/variables.std()
print vars_norm.mean()
print vars_norm.std()
Ahora que ya seleccionamos las varibles y las normalizamos, vamos a correr un primer modelo:
In [9]:
import statsmodels.formula.api as sm
model = sm.ols(formula="entrada ~ viv + entropia",data=vars_norm).fit()
print model.summary()
In [27]:
import clusterpy
datos_cp = clusterpy.importArcData('data/distritos_variables')
datos_cp.cluster('random',[datos_cp.fieldNames[0]],50,dissolve=0)
Recuerden que para ver las regiones resultantes pueden exportar los resultados como shapefile usando datos_cp.exportArcData('ruta/para/guardar'). Por lo pronto omitiremos este paso y más bien vamos a unir el resultado de los clusters a los datos originales para poder hacer algunas operaciones.
El primer paso es saber en qué columna guardo el algoritmo los resultados de la regionalización, para esto, vamos a ver qué columnas tiene nuestro objeto:
In [28]:
datos_cp.fieldNames
Out[28]:
Como pueden ver, el algoritmo le pone un nombre al azar a la columna donde va a guardar los resultados, random_20160405130236, en este caso.
Ahora lo que vamos a hacer es extraer sólo el identificador de distrito y el identificador de región, que es lo que necesitamos para hacer un merge con nuestros datros originales:
In [29]:
resultado = datos_cp.getVars(['cve_dist','random_20160406101421'])
print resultado
El resultado es un diccionario, es decir, un conjunto de parejas llave : valor, en este caso las llaves son el id de los polígonos y los valores son la clave del distrito y el identificador de la región. Como queremos unir estos resultados con nuestros datos originales, necesitamos convertirlos en un DataFrame, afortunadamete esto es relativamente fácil:
In [30]:
import pandas as pd
df = pd.DataFrame(resultado)
df.head()
Out[30]:
Casi lo logramos, el único problema es que los datos están "al revés", es decir renglones y columnas están intercambiados. Lo que necesitamos es trasponerlos:º
In [31]:
df = df.transpose()
df.head()
Out[31]:
Ahora pongámosle nombre a las columnas, la primera es el identificador de distrito y la segunda el de región:
In [32]:
df.columns=['cve_dist','id_region']
df.head()
Out[32]:
Ahora sí podemos hacer un merge con los datos originales:
In [33]:
region_1 = datos.merge(df,how='inner',on='cve_dist')
region_1.head()
Out[33]:
Ahora ya tenemos casi todo lo que necesitamos para correr un nuevo modelo, esta vez sobre las variables agregadas en regiones, lo único que necesitamos es, justamente, calcular las variables agregadas, para esto, vamos a hacer un "gropup by":
In [34]:
agregados = region_1.groupby(by='id_region').sum()
agregados.head()
Out[34]:
Aquí simplemente agrupamos los datos por su id de región y calculamos la suma de las variables sobre cada grupo.
Sólo hay un problema, la entropía quedo calculada como la suma de las entropías individuales y eso no nos sirve, necesitamos volver a calcularla:
In [35]:
intensidad = agregados['comercio'] + agregados['viv'] + agregados['ocio']
prop_comercio = agregados['comercio'] / intensidad
prop_viv = agregados['viv'] / intensidad
prop_ocio = agregados['ocio'] / intensidad
entropia = (prop_comercio*np.log(prop_comercio) + prop_viv*np.log(prop_viv) + prop_ocio*np.log(prop_ocio))/np.log(3)
agregados.ix[:,'entropia']= entropia
agregados.head()
Out[35]:
Ahora sí, podemos correr el mismo modelo de antes pero sobre nuestros datos agregados (recordemos que antes hay que normalizarlos):
In [43]:
variables = agregados[['entrada','viv','entropia']]
vars_norm = (variables - variables.mean())/variables.std()
model = sm.ols(formula="entrada ~ viv + entropia",data=vars_norm).fit()
print model.summary()
Ahora, el objetivo es comparar los resultados usando diferentes regionalizaciones aleatorias de los mismos datos. El proceso de pasar los datos a regiones es un poco largo, pero no se preocupen, para aliviar ese trabajo escribi una función que se encarga de encapsular todo el proceso:
In [40]:
def regionaliza_agrega(dataframe,shp_path='data/distritos_variables',n_regiones=50):
datos = clusterpy.importArcData(shp_path)
datos.cluster('random',[datos.fieldNames[0]],50,dissolve=0)
resultado = datos.getVars(['cve_dist',datos.fieldNames[-1]])
df = pd.DataFrame(resultado)
df = df.transpose()
df.columns=['cve_dist','id_region']
region = dataframe.merge(df,how='inner',on='cve_dist')
agregados = region.groupby(by='id_region').sum()
intensidad = agregados['comercio'] + agregados['viv'] + agregados['ocio']
prop_comercio = agregados['comercio'] / intensidad
prop_viv = agregados['viv'] / intensidad
prop_ocio = agregados['ocio'] / intensidad
entropia = (prop_comercio*np.log(prop_comercio) + prop_viv*np.log(prop_viv) + prop_ocio*np.log(prop_ocio))/np.log(3)
agregados.ix[:,'entropia']= entropia
return agregados
La función nos regresa los datos agregados en regiones aleatorias, lo único que necesitamos es selecionar las varables que vaos a usar, normalizarlas y correr un modelo:
In [45]:
agregados = regionaliza_agrega(datos)
variables = agregados[['entrada','viv','entropia']]
vars_norm = (variables - variables.mean())/variables.std()
model = sm.ols(formula="entrada ~ viv + entropia",data=vars_norm).fit()
print model.summary()
Como pueden ver, cada vez que corran la celda anterior, van a tener un resultado un poco diferente, es decir, cada modelo arroja resultados diferentes en la misma escala. En este caso, podemos pensar que no es muy grave, nosotros sabemos que las agregaciones son aleatorias y tenemos alguna certeza de que las desviaciones observadas en los resultados son aleatorias, pero ¿Qué pasa cuando no tenemos control sobre las agregaciones?
Antes de examinar esa pregunta, hay un...
Encuentren un modelo que explique razonablemente bien los viajes entrantes a la escala agregada (para esto usen una sóla regionalización aleatoria, es decir, no usen la función que les acabo de explicar). A aprtir de ese modelo, ahora sí, usando la función regionaliza_agrega, haz diferentes experimentos y registra los parámetros más importantes: la $R^2$ y los coeficientes de las variables. Con los resultados de varias repeticiones, trata de mostrar que las diferencias son aleatorias (es decir, que siguen una distribución normal).
Ahora vamos a trabajar sobre tres agregaciones diferentes de los datos, en la misma escala que las que hicimos, pero sobre las que no tenemos control. Simplemente son las unidades que tenemos (más o menos como las AGEBS, por ejemplo)
Lo primero que vamos a hacer es leer los nuevos datos:
In [47]:
datos_nuevos = GeoDataFrame.from_file('data/variables_regiones.shp')
datos_nuevos.head()
Out[47]:
Como pueden ver, son exáctamente igual a los originales, salvo porque tienen tres columnas extra: region_1, region_2 y region_3. Esta son, jústamente, las nuevas agregaciones de los datos, por lo demás, todo es exáctamente igiual que antes.
Entonces, construyamos tres agregaciones, una para cada regionalización, y ajustemos nuestro modelo lineal para ver qué sucede. Empecemos con la primera regionalización:
In [49]:
r_1 = datos_nuevos.groupby(by='region_1').sum()
intensidad = r_1['comercio'] + r_1['viv'] + r_1['ocio']
prop_comercio = r_1['comercio'] / intensidad
prop_viv = r_1['viv'] / intensidad
prop_ocio = r_1['ocio'] / intensidad
entropia = (prop_comercio*np.log(prop_comercio) + prop_viv*np.log(prop_viv) + prop_ocio*np.log(prop_ocio))/np.log(3)
r_1.ix[:,'entropia']= entropia
variables = r_1[['entrada','viv','entropia']]
vars_norm = (variables - variables.mean())/variables.std()
model_1 = sm.ols(formula="entrada ~ viv + entropia",data=vars_norm).fit()
print model_1.summary()
Ahora con la segunda:
In [50]:
r_2 = datos_nuevos.groupby(by='region_2').sum()
intensidad = r_2['comercio'] + r_2['viv'] + r_2['ocio']
prop_comercio = r_2['comercio'] / intensidad
prop_viv = r_2['viv'] / intensidad
prop_ocio = r_2['ocio'] / intensidad
entropia = (prop_comercio*np.log(prop_comercio) + prop_viv*np.log(prop_viv) + prop_ocio*np.log(prop_ocio))/np.log(3)
r_2.ix[:,'entropia']= entropia
variables = r_2[['entrada','viv','entropia']]
vars_norm = (variables - variables.mean())/variables.std()
model_2 = sm.ols(formula="entrada ~ viv + entropia",data=vars_norm).fit()
print model_2.summary()
Ahora ustedes pueden hacer solitos el tercer modelo.
¿Qué es lo que está pasando? ¿Por qué estamos obteniendo ahora valores tan diferentes a lo que teníamos antes?
La parte final de su tarea es tratar de explicar esto.
In [ ]: