Cuando los datos no tienen ubicación: un método para asignar viviendas censales

Los microdatos del censo y su falta de resolución espacial #

El INE (Instituto Nacional de Estadísticas) publicó los datos del Censo 2024 con microdatos a nivel comunal y 60 indicadores agregados por manzana. Es decir, si quiero ver datos de personas, solo los puedo visualizar a nivel de comuna, si quiero algo más específico como una manzana, solo puedo ver valores grupales. Esta no es la situación ideal en la que podemos saber en cuál manzana vive cada persona. Es entendible esa restricción impuesta por el INE, puesto que la densidad de los datos permitiría identificar personas. Cruzar el censo con otros datos permitiría usos que aporten a las políticas públicas y a la calidad de vida de las personas, pero también puede haber agentes maliciosos que hagan mal uso de los datos. Balancear esa necesidad de privacidad no es fácil. De hecho, investigadores de CEDEUS han cuestionado la decisión del INE porque sin desagregación espacial es difícil medir distintos fenómenos sociales y urbanos:

Sin este detalle será imposible evaluar problemas asociados al deterioro de las viviendas a nivel de distritos, zonas o manzanas censales y así focalizar políticas, planes, programas o proyectos de rehabilitación o regeneración urbana. Tampoco podremos calcular indicadores socioterritoriales para orientar programas de apoyo social. Esto solo por poner algunos ejemplos de las limitaciones que genera esta inédita decisión.

Quizás la mejor manera de que INE disponibilice los datos es utilizando técnicas de privacidad diferencial, lo que permitiría entregar datos que mantengan las propiedades estadísticas para preguntas específicas, a pesar de estar perturbados o randomizados. Hay países que lo han entregado así. Eso no implica que sea fácil de hacer. De hecho, hay que definir de manera previa cuáles preguntas se van a responder. Y, como indican en CEDEUS, esas preguntas son muchas.

¿Qué hacer, entonces? Podemos buscar una solución intermedia: un método de asignación que, dados los datos individuales a nivel comunal y los indicadores agregados a nivel de manzana, pueda distribuir las viviendas de las personas en distintas áreas geográficas, minimizando el error de asignación (medible porque se comparan los indicadores agregados con los indicadores que se pueden calcular después de asignar). En la siguiente imagen vemos comunas (datos individuales), manzanas (datos agregados), y una grilla regular hexagonal (sistema H3, resolución 8) a la que asignaremos viviendas del censo:

Entonces, podemos definir el problema así. Disponemos de:

• Microdatos: cada persona \(p\) pertenece a una vivienda \(v\), y cada vivienda pertenece a una comuna \(c\). No conocemos la posición de \(v\) dentro de la comuna.
• Manzanas: cada manzana \(m\) tiene indicadores agregados (totales de personas por categoría). Una manzana pertenece a una única comuna.
• Grilla \(G\): cada celda \(g\) del sistema H3 cubre un área fija. Asignamos cada manzana a la celda con mayor intersección: \(h(m) = \arg \max_{g} \text{Área}(m \cap g))\).

Y nuestro objetivo es encontrar una asignación \(\phi: V \rightarrow G\) tal que cada vivienda \(v\) de la comuna \(c\) se asigne a una celda \(g\) que intersecte a \(c\), minimizando la discrepancia entre los indicadores agregados en cada celda (calculados a partir de las manzanas) y los indicadores resultantes de la asignación (calculados a partir de las viviendas asignadas):

Donde:

• La variable \( \text{objetivo}_{g,a} \) es el valor esperado del atributo \(a\) en la celda \(g\), derivado de las manzanas asignadas a esa celda.
• \(x_{v,a}\) es el valor del atributo \(a\) en la vivienda \(v\) (por ejemplo, número de personas mayores a 65 años).
• \(\mathcal{A}\) es el conjunto de atributos a optimizar.

Para ilustrar qué tipo de análisis permite este método, veamos el caso de la inmigración en conjunto con resultados del modelo de asignación que describiré más abajo. En los datos por manzana, la variable n_inmigrantes indica cuántos habitantes son migrantes, y en conjunto con la variable n_per, se puede calcular la fracción migrante de la población por manzana. Pero no conocemos otras variables asociadas: ¿Cuál es su nacionalidad? ¿Cuál es el período de llegada al país? Esa información sí está en los microdatos por comuna. Al no haber más atributos por manzana asociados a esa población, puesto que los otros indicadores están agregados considerando población local también, no se puede profundizar en el estudio de la inmigración si se quiere hacer a un nivel más específico que comuna. Un método de asignación permite distribuir espacialmente los atributos que están en los microdatos, como el período de llegada:

Al ver los mapas nos damos cuenta de que la población migrante tiene una distribución espacial diferente si consideramos su período de llegada.

Hay otras variables similares en los microdatos, como p27_nacionalidad, que indica la nacionalidad de las personas. Esta es su distribución:

Noten que en efecto, el valor origen_otra_nacionalidad (mapa a la derecha) tiene una distribución similar a la distribución de migrantes por manzana.

Existen métodos exactos de optimización combinatoria que resuelven este problema. Sin embargo, su costo computacional escala mal con el número de viviendas y celdas. En ocasiones se necesita software especializado, en otras, poder de cómputo y mucho tiempo. Yo les presentaré una solución basada en métodos computacionales que es rápida y donde podemos tener control de la calidad de la solución.

El método de asignación #

Hoy, hablar de datos nos hace pensar enseguida en métodos de Machine Learning o incluso de Deep Learning. Pero este no es el caso. Un modelo de ML o DL necesita datos de entrenamiento. Si quisiéramos hacer algo en esa dirección, tendríamos que tomar un censo pasado en el que existan microdatos por manzana, de modo de implementar un predictor. Pero el Censo 2017 tampoco tuvo microdatos públicos a nivel de manzana, y en 2017 la población era muy diferente (y qué decir en comparación al censo anterior a ese, de 2002).

La respuesta está en un método que aprendí cuando hacía mi doctorado. Trabajaba con la biblioteca d3.js, que se utiliza (incluso hoy) para implementar visualizaciones de datos. Una visualización que ya implementa esta biblioteca es el gráfico de redes:

En este gráfico, los nodos de la red son pelotitas, y los vínculos entre ellas, líneas rectas que las conectan. Una red no tiene necesariamente posiciones, por ejemplo, una red social indica las relaciones sociales o de interacción entre personas. No es necesario conocer las coordenadas geográficas de las personas para visualizar la red. Pero, para poder verla con claridad, debemos asignarle una posición a los nodos. Una posición tal que los nodos conectados entre sí estén cerca y al mismo tiempo, no colisionen entre sí.

¿Se dieron cuenta de que las posiciones finales de cada nodo se calculan de manera iterativa? Al comienzo, la posición de los nodos es aleatoria, y paso a paso se van alejando entre sí. ¿En qué dirección se pueden mover? La idea es desplazarse en una dirección que mejore la calidad del gráfico, que puede medirse como la distancia entre nodos y una métrica de sobreposición de ellos. Esto, nuevamente, es un problema de optimización con múltiples restricciones. Y como tal es un problema costoso. ¿Cómo es que se puede hacer en tiempo real, entonces?

La respuesta está en la documentación de d3.js: utilizan una implementación basada en el método Simulated Annealing. En castellano, significa «enfriamiento simulado», y significa lo siguiente (cito a Wikipedia):

El nombre e inspiración de SA viene del proceso de recocido del acero y cerámicas, una técnica que consiste en calentar y luego enfriar lentamente el material para variar sus propiedades físicas. El calor causa que los átomos aumenten su energía y que puedan así desplazarse de sus posiciones iniciales (un mínimo local de energía); el enfriamiento lento les da mayores probabilidades de recristalizar en configuraciones con menor energía que la inicial (mínimo global).

En el caso de la visualización, en cada iteración se aplica este enfoque: primero, los nodos (átomos) se desplazan de manera aleatoria en cualquier dirección, pero a medida que el sistema se enfría, solo se pueden desplazar en una dirección al azar si es que su nueva posición es a lo sumo, levemente peor con una tolerancia (temperatura) que se va reduciendo con el tiempo. Este método permite implementar una visualización como la que les mostré o incluso resolver el problema del vendedor viajero:

Lo interesante del método es que permite escapar de óptimos locales (soluciones que parecen buenas si solo miras su vecindad, pero que pueden ser mucho peores que otras), aumentando la probabilidad de encontrar soluciones cercanas al óptimo global. En cambio, los métodos que usualmente se utilizan para optimizar funciones (como máxima verosimilitud) tienen el problema de atascarse en óptimos locales.

Como suele suceder en la ciencia, el método fue definido (¿o descubierto?) por distintas personas de manera independiente entre 1970 y 1983. Su nombre proviene del paper publicado en Science (1983) por Kirkpatrick et al¹.

Entonces, ¿cómo aplicarlo al censo? Resumámoslo primero: obtenemos una solución inicial basada en asignación ingenua (lo que efectivamente nos daría un óptimo local), luego, en cada paso, intercambiamos de lugar viviendas de manera aleatoria, y aceptamos ese cambio dependiendo de la tolerancia o temperatura actual del sistema. Veamos ahora cómo construir una asignación inicial razonable y cómo refinarla iterativamente usando SA.

La implementación del método #

Primero necesitamos preprocesar los datos de una comuna: se asigna cada manzana a la celda con mayor intersección, se suman los indicadores de manzanas por celda, y se construyen perfiles de vivienda desde microdatos de personas, puesto que asignamos viviendas directamente, no personas. Estos perfiles deben tener los atributos \(\mathcal{A}\).

Luego viene la etapa de escalado. Como los indicadores totales de cada celda (obtenidos desde sus manzanas) y de los microdatos pueden diferir², para cada atributo \(a \in \mathcal{A}\), usamos la distribución espacial de la grilla pero los totales reales:

Donde:

• \(X_{g,a}^{\text{grilla}}\) es el valor del atributo \(a\) en la celda \(g\), sumado desde las manzanas.
• \(\sum_v x_{v,a}\) es el total real del atributo en los microdatos, calculado a partir de los perfiles de viviendas.

Esto nos permite avanzar a la etapa de asignación inicial. Aquí debemos calcular una distancia entre vivienda y celda. Una distancia pequeña implica que la vivienda está bien asignada, porque sus indicadores son similares a los de la celda. Se define así:

Donde:

• \(w_a\) es el peso del atributo \(a\) (prioridad alta/media/baja). Esto permite decirle al método que un atributo es más relevante que otro para el estudio que quieras hacer.
• \(p_{v,a}\) es la proporción del atributo \(a\) en la vivienda \(v\).
• \(p_{g,a}\) es la proporción del atributo \(a\) en la celda \(g\).

Esta distancia nos permite definir la probabilidad de que una vivienda \(v\) pertenezca a una celda \(g\), utilizando la función softmax:

Donde:

• \(\text{capacidad}_g\) es la cantidad de viviendas esperadas de la celda \(g\). Esto proviene de los datos de manzanas, puesto que indican cuantas viviendas hay, por tipo (casa, departamento o móvil).
• \(\tau\) es un parámetro de suavidad. Mientras más alto es, las probabilidades son uniformes (la vivienda podría estar en muchas celdas); un valor bajo hace que se concentren en las celdas de menor distancia³.

Con la función, asignamos cada vivienda a la celda con mayor probabilidad de contenerla. Lo hacemos estratificado, es decir, cada tipo de vivienda (casa, departamento, etc.) se procesa por separado. Si la celda con mayor probabilidad ya vio agotada su capacidad, y todavía hay viviendas cuya probabilidad máxima las lleva a esa celda, entonces le asignamos la segunda celda, o la tercera, y así sucesivamente. Esto garantiza consistencia de la asignación con la distribución de tipos de viviendas observada.

Después viene la etapa de simulated annealing. Necesitamos definir su energía como el error de asignación total, es decir, la sumatoria del error en cada celda \(g\):

El movimiento aleatorio que haremos es tomar dos celdas al azar: una con probabilidad de elección proporcional a su error; y otra, con proporción inversa. Luego, tomamos una vivienda del mismo tipo en cada una de ellas, y las cambiamos de celda. Medimos la nueva energía del sistema, obteniendo la diferencia \(\Delta E = E_{\text{con cambio}} - E_{\text{sin cambio}}\). Si \(\Delta E\) es negativo, perdemos energía, lo que es bueno porque se reduce el error. Nos preguntamos si aceptamos ese cambio o no utilizando el criterio de Metropolis:

Este criterio contiene el parámetro de temperatura \(T\). Inicialmente, la temperatura es muy alta, por lo que la probabilidad de aceptar malos cambios también lo es. Pero a medida que iteramos, la temperatura va bajando con una tasa de enfriamiento \(\alpha\), y por tanto se vuelve más difícil aceptar un cambio que aumente la energía. El enfriamiento del paso \(k+1)\) se define como:

Estos intercambios se hacen hasta que la temperatura alcance el umbral mínimo. Al terminar, la asignación de viviendas permite generar mapas como los que vimos al comienzo de este post. Ahora bien, ejecutar el algoritmo no es suficiente, también es necesario entender cómo se comporta. Para ello, tenemos que evaluar sus resultados. Lo podemos hacer calculando un error relativo por atributo, definido como:

En la fórmula, el denominador es el total real del atributo \(a\) en la comuna. Esto indica qué fracción del total de viviendas está «mal ubicada». Por ejemplo, si hay 1000 migrantes y el numerador suma 50, el error relativo es 0.05 (o 5%). Conocer los errores por atributo permite determinar si es necesario incorporar más variables para asignación, cambiar sus importancias, o incluso determinar si el método no funciona en una comuna específica. ¡El proceso de definición, implementación y aceptación de los resultados del método también es un método iterativo!

El siguiente diagrama esquematiza la metodología:

Aplicamos esto comuna por comuna. Así, incluso cuando haya errores, las viviendas siguen asignadas a la comuna a la que pertenecen. En ciudades como Santiago esto es importante. A mi juicio, es preferible obtener una distribución aproximada a trabajar con una distribución uniforme.

Finalmente, si deseas que la asignación en sí misma sea probabilística, puedes aplicar bootstrapping con múltiples semillas para la generación de números aleatorios.

Discusión #

Este método heurístico tiene limitaciones importantes. Por un lado, no garantiza encontrar la asignación óptima. Pero eso está bien, probablemente un método más complejo (o completo) tampoco lo encontrará dadas las restricciones de privacidad o la capacidad de cómputo disponible. Este método, en cambio, permite obtener resultados en tiempos razonables (asignar todo el Gran Santiago toma una hora en mi computador de escritorio) e incorporar tus propias variables o datos si lo ves necesario. El método en sí mismo es directo de implementar en Python utilizando operaciones de numpy.

Por otro lado, hay variables de baja frecuencia que son difíciles de optimizar. Además, el conjunto de atributos \(\mathcal{A}\) no necesariamente se correlaciona con todo lo que quisiéramos poder medir. Ejemplifiquemos esto: si bien dentro de indicadores se incluye cuántas personas utilizan cada modo de transporte para ir al trabajo, no hay ningún indicador que hable explícitamente dónde está el lugar de trabajo, algo que sí está disponible en los microdatos (la variable p44_lug_trab). Sí hay indicadores respecto al tipo de trabajo o nivel educacional de las personas, pero lo cierto es que estos no están tan correlacionados con el destino, a sabiendas de que el eje Alameda-Providencia-Apoquindo acapara gran cantidad de viajes de múltiples tipos de trabajo (servicios, profesionales, académicos, etc.). Como resultado, el siguiente mapa muestra la distribución espacial asignada de la variable:

La distribución relativa del lugar de trabajo muestra patrones reconocibles: en el sector oriente, la población tiende a trabajar en la misma comuna, o incluso para otro país. Sin embargo, la variación intracomuna no es fuerte (y tenderíamos a pensar que sí), y eso se debe a que la asignación no considera esta variable ni se utilizó una variable que estuviese fuertemente correlacionada. Por tanto, la distribución se vuelve aleatoria dentro de la comuna.

Esto último no es realmente un problema. Es una limitación que debe tenerse clara. Una manera de resolverla es incorporando datos externos (como datos origen-destino de transporte), que también existen.

Me encantó implementar este algoritmo. Me hace pensar que no todo es Machine Learning, sino que también hay belleza y valor en soluciones computacionales clásicas. Un algoritmo como Simulated Annealing es interpretable, no requiere entrenamiento ni grandes cantidades de datos, y nos permite controlar (o al menos entender) los errores que comete. Es una maravilla cómo una metáfora física (el enfriamiento de metales) se convirtió en una herramienta computacional que puede usarse en disciplinas como la visualización o el análisis de datos sociodemográficos.