«

Ideas para modelos de inferencia ecológica

Las tres ideas que describo a continuación buscan resolver el problema que presenta la heterogeneidad espacial en inferencia ecológica. En el caso de las transiciones entre votos en primera vuelta y en ballotage (ver mi post anterior sobre el tema), esto se expresa en que en lugares distintos el comportamiento electoral difiere. El modelo tiene que incorporar esta heterogeneidad. El trabajo de Calvo y Escolar (2003) sirve de inspiración para el planteo del problema y para una de las ideas.$\def\N{\mathbb{N}}\def\R{\mathbb{R}}\def\Err{\text{Err}}\def\E{\mathbb{E}}\def\sen{\text{sen}}$

Idea 1: suavidad espacial

La idea es partir del modelo multinomial-Dirichlet de King, Rosen y Tanner pero hacer un cambio. King et al buscan modelar las diferencias entre mesas de votación, y para eso ponen una fuente de aleatoriedad en la variable $\beta$ de cada mesa (que estima la matriz de transferencia): la hacen tomar valores según una distribución Dirichlet, en cuyos parámetros se puede incluir una covariable que asocia un número a cada mesa, y por lo tanto las distingue. La idea es eliminar esa aleatoriedad, pero incorporar una restricción sobre las $\beta$: si están asociados a mesas que pertenecen a una misma unidad geográfica predeterminada (un mismo departamento, por ejemplo) son iguales; y la variación entre las $\beta$ de unidades cercanas está acotada, de manera que la función que, dada una unidad, asigna su $\beta$, resulta suave en un sentido que cuantificamos. Doy los detalles.

La idea es hacer maximum likelihood sobre un modelo multinomial. En la mesa $p$ los resultados de la primera vuelta son $X^p\in\N^7$ y los de la segunda, $Y^p\in\N^3$; tenemos $N^p=\sum_{i=1}^7X^p_i=\sum_{j=1}^3 Y^p_j$ (la cantidad de gente en esa mesa). Asigno una variable $\beta^p\in\R^{7\times 3}$ a la mesa, que cumple $\beta^p_{ij}\geqq 0$ y $\sum_{j=1}^3\beta^p_{ij}=1$ para $i=1,\ldots,7$. La idea es que $\beta$ no sea exactamente la matriz de transición, pero que el valor esperado de $Y^p$ sea $X^p\beta^p$. Modelamos $Y^p$ como dada por una distribución multinomial de $N^p$ tiradas con probabilidades de que salga $j=1,2,3$ dadas por $p_j=\frac{1}{N^p}\sum_{i=1}^7 X_i^p\beta^p_{ij}$. Entonces tenemos $$P(Y^p=y\mid X^p=x,\beta^p)=\frac{N^p!}{Y_1^p!Y_2^p!Y_3^p!}p_1^{Y_1^p}p_2^{Y_2^p}p_3^{Y_3^p}.$$

Ahora ¿cómo tomo los $\beta^p$? La idea de King-Rosen-Tanner es samplearlos de una distribución Dirichlet con ciertos parámetros. Nuestra idea es tomarlos como parámetros a estimar bajo una restricción de suavidad espacial. Es decir, que los $\beta^p$, como función de una variable espacial (geográfica), se asume que varían de forma “suave”. O, si pensamos que las mesas forman un grafo, de manera que si $p$ y $p’$ son escuelas “vecinas” está definida su distancia $d(p,p’)$, entonces restringimos el espacio de búsqueda a los $\beta^p$ tales que los $\|\beta^p-\beta^{p'}\|$ están controlados de alguna manera por los $d(p,p’)$. Una manera posible: poner una constante $C>0$ y restricciones $\|\beta^p-\beta^{p'}\|\leqq Cd(p,p')$ estilo Lipschitz, donde $\|{\cdot}\|$ es alguna norma. Otra: restringir $\sum_{(p,p')\in E}\frac{1}{d(p,p')}\|\beta^p-\beta^{p'}\|\leqq C$ y $\beta^p=\beta^{p’}$ si $d(p,p’)=0$; o $\sum_{(p,p')\in E}e^{-\lambda d(p,p')}\|\beta^p-\beta^{p'}\|\leqq C$, algo por el estilo. Una pregunta importante va a ser cómo elegir $C$. Creemos que se puede hacer cross-validation sobre la log-likelihood media para obtener el mejor valor de este parámetro.

Entonces, dada $\beta$, uno puede medir cuánto se ajusta a la realidad tomando como medida del “error” lo siguiente: $$\Err(\beta)=-\sum_p\log(P(Y^p=y\mid X^p=x,\beta^p))=$$ $$=-\sum_p\left(\log(N^p!)-\sum_{j=1}^3\log(Y_j^p!)+\sum_{j=1}^3Y_j^p\log\left(\frac{1}{N^p}\sum_{i=1}^7 X_i^p\beta^p_{ij}\right)\right).$$ Maximum likelihood es lo mismo que minimizar esa función. Lo que uno busca es, pues, buscar $\beta$ que minimice $\Err(\beta)$ sujeto a las restricciones de suavidad espacial. Quitando los términos constantes, se busca minimizar $$\Err(\beta)=-\sum_p\sum_{j=1}^3Y_j^p\log\left(\frac{1}{N^p}\sum_{i=1}^7 X_i^p\beta^p_{ij}\right)_{}.$$ Si los $Y_j^p$ sumaran 1 lo que estaríamos minimizando es la distancia de Kullback-Leibler, o cross entropy, entre la distribución del voto en segunda vuelta que da $X^p$ y $\beta^p$ contra la que efectivamente se dio.

Para asegurar la suavidad espacial podemos incorporar el término regularizador a la función a minimizar: $$\Err(\beta)=-\sum_p\sum_{j=1}^3Y_j^p\log\left(\frac{1}{N^p}\sum_{i=1}^7 X_i^p\beta^p_{ij}\right)+\phi\left(\sum_{(p,p')\in E}\frac{1}{d(p,p')}\|\beta^p-\beta^{p'}\|\right),$$ donde $\phi$ es una función penalizadora para imponer la restricción, por ejemplo $\phi(x)=k\max\{x-C,0\}$ o $\phi(x)=\frac1k(e^{k(x-C)}-1)$, con $k\to+\infty$. Otra posibilidad es tomar $\phi(x)=\lambda x$ y buscar (con búsqueda binaria) el $\lambda$ más chico tal que $\text{argmín}_\beta\Err(\beta)$ cumple la restricción.

Para olvidarnos de la restricción $\beta^p_{ij}\geqq 0$, $\sum_{j=1}^3\beta^p_{ij}=1$ para $i=1,\ldots,7$, podemos tomar $\beta^p=\sigma(\gamma^p)$, donde $\sigma:\R^{7\times 3}\to\R^{7\times 3}$ está dada por $$\sigma(\gamma)_{ij}=\frac{e^{\gamma_{ij}}}{\sum_{k=1}^3e^{\gamma_{ik}}}$$ y $\gamma^p$ puede tomar cualquier valor en $\R^{7\times 3}$ (es una biyección entre ambos dominios si fijamos $\gamma^p_{i3}=0$ para $i=1,\ldots,7$ y reforzamos las desigualdades a $\beta^p_{ij}>0$).

Podemos tomar como aristas del grafo los $K$ vecinos más cercanos de cada unidad que tomamos (escuela, circuito electoral o departamento), y $d(p, p’)$ puede ser simplemente la distancia en kilómetros entre los centroides de las unidades. Otra posibilidad: un minimum spanning tree.

Metodología

Tomamos unidades geográficas. Voy a tomar los departamentos como unidades. Son el siguiente nivel de desagregación después de provincias; en CABA son las comunas y en el resto del país son los municipios. El modelo va a asumir que en cada unidad las mesas van a seguir la misma distribución. (El resultado de la segunda vuelta en la unidad $p$ y en la mesa $k$ sigue una distribución multinomial con probabilidades $q_j=\frac{1}{N^p}\sum_{i=1}^7 X_i^k\beta^p_{ij}$, $j=1,2,3$; notar que el $\beta^p$ está dado por la unidad y es el mismo para cada mesa; el $X^k$ -el resultado de la primera vuelta- es propio de la mesa.) ¿Cómo se eligen las unidades, en general? Es un metaparámetro del modelo. Lo que dice el modelo es que en cada unidad las mesas se comportan exactamente igual. Entonces hay que tomar las unidades más grandes posibles en las que pase eso. En la práctica va a depender de la cantidad de datos.

Voy a comparar modelos, por lo que voy a separar un train set de un test set antes de entrenarlos. El entrenamiento de cada $\beta^p$, donde $p$ es una unidad, se hace sobre una muestra de mesas en esa unidad. Por lo tanto la separación de train set y test set va a darse en cada unidad: para cada una voy a tomar un $30\%$ de mesas al azar y separarlas; todas esas mesas van a conformar el test set. Si no hacemos esto podríamos terminar haciendo overfitting en algunas unidades y “underfitting” en otras. (Se podría hacer la prueba para ver si esto tiene un impacto o no.)

Dentro del train set podemos hacer $K$-cross validation. Y vamos a hacer la partición de la misma manera: para cada unidad partimos la muestra de entrenamiento en $K$ partes proporcionales, etc.

Generalización

La suavidad espacial no tiene que ser necesariamente geográfica. Si uno tiene indicadores $Z_1,\ldots,Z_K$ sobre las mesas ($Z_k\in\R$) podemos definir una distancia entre mesas $r$ y $s$ dada por $d(r,s)=\sum_{k=1}^K\eta_k(Z_k^r-Z_k^s)^2$, con $\eta_k\geqq 0$ fijos (o algo por el estilo). Usando algún método de clustering podemos obtener nuestras unidades (la generalización del concepto de “municipio”). Y luego podemos usar como distancia entre unidades ya sea $d(p,p')=\sum_{k=1}^K\eta_k(\E(Z_k^r)-\E(Z_k^s))^2=\sum_{k=1}^K\eta_k\E(Z_k^r-Z_k^s)^2$ o $d(p,p')=\E\left[d(r,s)\right]=\sum_{k=1}^K\eta_k\E\left[(Z_k^r-Z_k^s)^2\right]$, donde $r$ y $s$ son las mesas de las unidades $p$ y $p’$, respectivamente, con la distribución empírica.

Idea 2: regularización

Partimos de la misma base: en la mesa $p$, $X^p\in\N^7$ son los votos en primera vuelta, $Y^p\in\N^p$ son los votos en la segunda vuelta, con $\sum_{i=1}^7X_i^p=\sum_{j=1}^3Y_j^p=N^p$, y $\beta^p\in\R^{7\times 3}$ es una estimación de la matriz de transición, de manera que $Y^p$ sigue una distribución multinomial de $N^p$ tiradas con probabilidades $q_j = \frac{1}{N^p}\sum_{i=1}^7 \beta_{ij}^pX_i^p$. Tenemos $$\beta^p_{ij}=\frac{e^{\gamma_{ij}^p}}{\sum_{k=1}^3e^{\gamma_{ik}^p}},$$ con $\gamma^p\in\R^{7\times 3}$, $\gamma^p_{i3}=0$. El modelo es: $\gamma_{ij}^p = \sum_{k=1}^K \delta^k_{ij} Z_k^p$, donde $Z_1,\ldots,Z_K$ son indicadores o covariables que toman valores en $\R$ para cada mesa.

Si lo pensamos como un problema de estimación, buscamos los números $\delta_{ij}^k$ (con $i=1,\ldots,7$, $j=1,2,3$, $k=1,\ldots,K$) que maximizan la verosimilitud, es decir, buscamos minimizar $$\Err(\delta)=-\sum_p\sum_{j=1}^3Y_j^p\log\left(\frac{1}{N^p}\sum_{i=1}^7 X_i^p\frac{e^{\sum_{k=1}^K \delta^k_{ij} Z_k^p}}{\sum_{j'=1}^3e^{\sum_{k=1}^K \delta^k_{ij'} Z_k^p}}\right)_{}.$$

Si esto lo pensamos como un problema de predicción (no es necesario), buscamos minimizar el error cuadrático medio, es decir, $$\Err(\delta)=\sum_p\sum_{j=1}^3\left(Y_j^p-\sum_{i=1}^7 X_i^p\frac{e^{\sum_{k=1}^K \delta^k_{ij} Z_k^p}}{\sum_{j'=1}^3e^{\sum_{k=1}^K \delta^k_{ij'} Z_k^p}}\right)^2.$$ Imitando a LASSO podemos hacer regularización $L^1$, lo que sería minimizar $$\Err(\delta)=\sum_p\sum_{j=1}^3\left(Y_j^p-\sum_{i=1}^7 X_i^p\frac{e^{\sum_{k=1}^K \delta^k_{ij} Z_k^p}}{\sum_{j'=1}^3e^{\sum_{k=1}^K \delta^k_{ij'} Z_k^p}}\right)^2+\lambda\sum_{i=1}^7\sum_{j=1}^3\sum_{k=1}^k|\delta_{ij}^k|,$$ donde $\lambda$ se determina por cross-validation.

¿Esto basta para capturar el impacto no lineal de la variable geográfica? Si la variable es $U,V$ (latitud y longitud, por ejemplo) se puede probar con indicadores de forma monomial $U^nV^m$ o trigonométrica $\sen(nU)\sen(mV)$, $\sen(nU)\cos(mV)$, etc, o splines.

Idea 3: eliminar la autocorrelación espacial

El modelo de King, Rosen y Tanner dice que $\E(Y_j^p) = \frac{1}{N^p}\sum_{i=1}^7 \beta_{ij}^pX_i^p$. Es decir, $$Y_j^p = \frac{1}{N^p}\sum_{i=1}^7 \beta_{ij}^pX_i^p + \epsilon_j^p,$$ con $\E(\epsilon_j^p)=0$. Supone que $\epsilon_j^p$ y $\epsilon_j^{p'}$ son independientes si $p$ y $p’$ son mesas distintas. Ahora bien, si no incorporamos la variable geográfica al modelo, esto va a fallar: mesas cercanas no van a ser independientes entre sí. Podemos modelar esto así, basándonos en Anselin y Ray (1991): $$\epsilon_j = \lambda W \epsilon_j + \mu_j,$$ donde $\lambda$ es un parámetro a estimar (suficientemente chico), $W\in\R^{P\times P}$ con $W_{p,p'}=d(p,p'),$ y $\mu_j^p$ son errores independientes y con esperanza nula. Es decir, $(I-\lambda W)\epsilon_j = \mu_j$, o $\epsilon_j = (I-\lambda W)^{-1}\mu_j$. Reemplazando y multiplicando por $I-\lambda W$ obtenemos $(I-\lambda W)Y_j = (I-\lambda W)\frac{1}{N}\sum_{i=1}^7 \beta_{ij}X_i+ \mu_j$. Esto da el modelo auto-regresivo $$Y_j^p = \lambda \sum_{p'}d(p,p')Y_j^{p'} + \frac{1}{N^{p}}\sum_{i=1}^7 \beta_{ij}^{p}X_i^{p} - \lambda\sum_{p'}d(p,p')\frac{1}{N^{p'}}\sum_{i=1}^7 \beta_{ij}^{p'}X_i^{p'}+\mu_j^p.$$ Usando mínimos cuadrados podemos estimar $\lambda$ y los parámetros $\gamma_{ij}$ que determinan la distribución de los $\beta_{ij}^p$.