Distribución de los ingresos en Colombia

1 Introducción

La base de datos personas.csv disponible en la página web del curso, corresponde a una muestra aleatoria del módulo de Personas (para las que el ingreso total está disponible y es mayor que cero) de la encuesta Medición de Pobreza Monetaria y Desigualdad 2021 llevada a cabo en Colombia por el DANE, la cual se encuentra disponible en este enlace.

El Universo de la encuesta está conformado por la población civil no institucional, residente en todo el territorio nacional; va dirigida a todos los hogares encontrados en la vivienda. La encuesta utiliza informante directo para las personas de 18 años y más, y para aquellas de 10 a 17 años que trabajen o estén buscando trabajo. Para los demás se acepta informante idóneo (persona del hogar mayor de 18 años, que a falta del informante directo pueda responder correctamente las preguntas). No se acepta información de empleados del servicio doméstico, pensionistas, vecinos o menores, excepto cuando el menor de edad es el jefe del hogar o cónyuge.

El objetivo de este trabajo es construir un modelo multinivel completamente Bayesiano, tomando como datos de entrenamiento el ingreso total (ingtot; ingreso total por persona que resulta de sumar cada una de las fuentes de ingresos tanto observadas como imputadas), con el fin de modelar los ingresos por dominio (dominio; cada una de las 24 a.M., otras cabeceras y resto), y establecer un ranking junto con una segmentación de los mismos. Por lo tanto, se toma como variable de agrupamiento el dominio, y como variable respuesta el ingreso total.

2 Modelos

\(\textsf{M}_1\): Modelo Normal

Distribución muestral: \[ y_{i,j}\mid\theta,\sigma^2 \,{\stackrel{\text{i.i.d.}}{\sim}}\,\textsf{N}(\theta,\sigma^2)\,, \] para \(i = 1,\ldots,n_j\) y \(j = 1,\ldots,m\), donde \(y_{i,j}\) es la variable respuesta del individuo \(i\) en el grupo \(j\) y \(\textsf{N}(\theta,\sigma^2)\) denota la distribución Normal con media \(\theta\) y varianza \(\sigma^2\).

Distribución previa: \[ \theta\sim\textsf{N}(\mu_0,\gamma_0^2)\,,\quad \sigma^2\sim\textsf{GI}\left(\frac{\nu_0}{2},\frac{\nu_0\sigma^2_0}{2}\right)\,, \] donde \(\mu_0,\gamma_0^2,\nu_0,\sigma^2_0\) son los hiperparámetros del modelo y \(\textsf{GI}(\alpha,\beta)\) denota la distribución Gamma-Inversa con media \(\frac{\beta}{\alpha-1}\), para \(\alpha > 1\), y varianza \(\frac{\beta^2}{(\alpha-1)^2(\alpha-2)}\), para \(\alpha > 2\).

\(\textsf{M}_2\): Modelo Normal con medias específicas

Distribución muestral: \[ y_{i,j}\mid\theta_j,\sigma^2 \,{\stackrel{\text{ind.}}{\sim}}\,\textsf{N}(\theta_j,\sigma^2)\,. \]

Distribución previa: \[ \theta_j\mid\mu,\tau^2 \,{\stackrel{\text{i.i.d.}}{\sim}}\,\textsf{N}(\mu,\tau^2)\,,\quad \mu\sim\textsf{N}(\mu_0,\gamma_0^2)\,,\quad \tau^2\sim\textsf{GI}\left(\frac{\eta_0}{2},\frac{\eta_0\tau^2_0}{2}\right)\,,\quad \sigma^2\sim\textsf{GI}\left(\frac{\nu_0}{2},\frac{\nu_0\sigma^2_0}{2}\right)\,, \] donde \(\mu_0,\gamma_0^2,\eta_0,\tau^2_0,\nu_0,\sigma^2_0\) son los hiperparámetros del modelo.

\(\textsf{M}_3\): Modelo Normal con medias y varianzas específicas

Distribución muestral: \[ y_{i,j}\mid\theta_j,\sigma_j^2 \,{\stackrel{\text{ind.}}{\sim}}\,\textsf{N}(\theta_j,\sigma^2_j)\,. \]

Distribución previa: \[\begin{align*} &\theta_j\mid\mu,\tau^2 \,{\stackrel{\text{i.i.d.}}{\sim}}\,\textsf{N}(\mu,\tau^2)\,, & &\mu\sim\textsf{N}(\mu_0,\gamma_0^2)\,, & &\tau^2\sim\textsf{GI}\left(\frac{\eta_0}{2},\frac{\eta_0\tau^2_0}{2}\right)\,, &\\ &\sigma^2_j\sim\textsf{GI}\left(\frac{\nu}{2},\frac{\nu\sigma^2}{2}\right)\,, & &\nu\sim\text{Constante}\,,& &\sigma^2\sim\textsf{G}\left(\frac{\alpha_0}{2},\frac{\beta_0}{2}\right)\,,& \end{align*}\] donde \(\mu_0,\gamma_0^2,\eta_0,\tau^2_0,\nu,\alpha_0,\beta_0\) son los hiperparámetros del modelo y \(\textsf{G}(\alpha,\beta)\) denota la distribución Gamma con media \(\frac{\alpha}{\beta}\) y varianza \(\frac{\alpha}{\beta^2}\).

\(\textsf{M}_4\): Modelo t

Distribución muestral: \[ y_{i,j}\mid\theta,\sigma^2 \,{\stackrel{\text{i.i.d.}}{\sim}}\,\textsf{t}_\kappa(\theta,\sigma^2)\,, \] donde \(\textsf{t}_\kappa(\theta,\sigma^2)\) denota la distribución t con \(\kappa\) grados de libertad con media \(\theta\), para \(\kappa > 1\), y varianza \(\frac{\kappa}{\kappa - 2}\,\sigma^2\), para \(\kappa > 2\).

La variable aleatoria \(X\) tiene distribución t con parámetros \(\kappa \in \mathbb{N}\), \(-\infty < \theta < \infty\), \(\sigma^2 > 0\), i.e., \(X\mid\kappa,\theta,\sigma^2\sim\textsf{t}_\kappa(\theta,\sigma^2)\), si su función de densidad de probabilidad es \[ p(x\mid\kappa,\theta,\sigma^2) = \frac{1}{\sqrt{\pi\kappa\sigma^2}} \, \frac{\Gamma\left( (\kappa+1)/2 \right)}{\Gamma\left( \kappa/2 \right)}\,\left( 1 + \frac{( x - \theta )^2}{\kappa\sigma^2} \right)^{-(\kappa+1)/2} \,,\quad -\infty < x < \infty\,. \] Si \(X\mid\kappa,\theta,\sigma^2\sim\textsf{t}_\kappa(\theta,\sigma^2)\), entonces \(\textsf{E}(X) = \theta\), para \(\kappa > 1\), y \(\textsf{Var}(X) = \frac{\kappa}{\kappa - 2}\,\sigma^2\), para \(\kappa > 2\).

Esta distribución es útil para modelar y se encuentra implementada en el paquete de que se encuentra disponible en este enlace.

Para ajustar este modelo de manera directa utilizando el muestreador de Gibbs, se debe tener en cuenta que la distribución muestral \(y_{i,j}\mid\theta,\sigma^2 \,{\stackrel{\text{i.i.d.}}{\sim}}\,\textsf{t}_\kappa(\theta,\sigma^2)\) es equivalente a la distribución jerárquica dada por \[ y_{i,j}\mid\theta,\varsigma^2_{i,j} \,{\stackrel{\text{ind.}}{\sim}}\,\textsf{N}(\theta,\varsigma^2_{i,j})\,,\qquad \varsigma^2_{i,j}\mid\sigma^2\,{\stackrel{\text{i.i.d.}}{\sim}}\,\textsf{GI}\left( \frac{\kappa}{2}, \frac{\kappa\sigma^2}{2}\right)\,, \] donde las variables \(\varsigma^2_{i,j}\) son cantidades auxiliares (desconocidas) cuyo objetivo es facilitar la implementación del muestreador de Gibbs.

La inclusión de las variables \(\varsigma^2_{i,j}\) en el modelo permite que todas las distribuciones condicionales completas de las cantidades desconocidas (incluyendo las mismas variables auxiliares) tengan forma probabilística conocida. En Gelman et al. (2013, pp. 293-294) hay una discusión detallada al respecto. Si no se consideran las variables \(\varsigma^2_{i,j}\) en el modelo, la implementación del muestreador de Gibbs requeriría de otros métodos numéricos más sofisticados como el algoritmo de Metropolis-Hastings o el algoritmo de Monte Carlo Hamiltoniano, dado que la distribuciones condicionales completas tanto de \(\theta\) como \(\sigma^2\) no tendrían forma probabilística conocida.

Esta misma consideración acerca de las variables auxiliares se debe tener en cuenta para la implementación computacional de los modelos 5 y 6.

Distribución previa: \[ \theta\sim\textsf{N}(\mu_0,\gamma_0^2)\,,\quad \sigma^2\sim\textsf{G}\left(\frac{\alpha_0}{2},\frac{\beta_0}{2}\right)\,,\quad \kappa\sim\text{Constante}\,, \] donde \(\mu_0,\gamma_0^2,\alpha_0,\beta_0,\kappa\) son los hiperparámetros del modelo.

\(\textsf{M}_5\): Modelo t con medias específicas

Distribución muestral: \[ y_{i,j}\mid\theta_j,\sigma^2 \,{\stackrel{\text{ind.}}{\sim}}\,\textsf{t}_\kappa(\theta_j,\sigma^2)\,. \]

Distribución previa: \[\begin{align*} &\theta_j\mid\mu,\tau^2 \,{\stackrel{\text{i.i.d.}}{\sim}}\,\textsf{N}(\mu,\tau^2)\,, & &\mu\sim\textsf{N}(\mu_0,\gamma_0^2)\,, & &\tau^2\sim\textsf{GI}\left(\frac{\eta_0}{2},\frac{\eta_0\tau^2_0}{2}\right)\,, &\\ &\sigma^2\sim\textsf{G}\left(\frac{\alpha_0}{2},\frac{\beta_0}{2}\right)\,, & &\kappa\sim\text{Constante}\,, & & & \end{align*}\] donde \(\mu_0,\gamma_0^2,\eta_0,\tau^2_0,\alpha_0,\beta_0,\kappa\) son los hiperparámetros del modelo.

\(\textsf{M}_6\): Modelo t con medias y varianzas específicas

Distribución muestral: \[ y_{i,j}\mid\theta_j,\sigma^2_j \,{\stackrel{\text{ind.}}{\sim}}\,\textsf{t}_\kappa(\theta_j,\sigma_j^2)\,. \]

Distribución previa: \[\begin{align*} &\theta_j\mid\mu,\tau^2 \,{\stackrel{\text{i.i.d.}}{\sim}}\,\textsf{N}(\mu,\tau^2)\,, & &\mu\sim\textsf{N}(\mu_0,\gamma_0^2)\,, & &\tau^2\sim\textsf{GI}\left(\frac{\eta_0}{2},\frac{\eta_0\tau^2_0}{2}\right)\,, &\\ &\sigma_j^2\,{\stackrel{\text{i.i.d.}}{\sim}}\,\textsf{G}\left(\frac{\alpha}{2},\frac{\beta}{2}\right)\,, & &\alpha\sim\text{Constante}\,, & &\beta\sim\textsf{G}\left(\frac{a_\beta}{2},\frac{b_\beta}{2}\right)\,, &\\ & & &\kappa\sim\text{Constante}\,, & & & \end{align*}\] donde \(\mu_0,\gamma_0^2,\eta_0,\tau^2_0,\alpha,a_\beta,b_\beta,\kappa\) son los hiperparámetros del modelo.

3 Grafo acíclico dirigido de \(\textsf{M}_6\)

A continuación se presenta el DAG de \(\textsf{M}_6\) incluyendo las variables auxiliares.

• \(\textsf{N}\): Normal.
• \(\textsf{G}\): Gamma.
• \(\textsf{GI}\): Gamma-Inversa.

4 Ajuste de los modelos

En todos los casos, se utiliza la siguiente convención:

• \(n\): número total de personas en la muestra.
• \(m\): número de dominios.
• \(n_j\): número de personas en el dominio \(j\).
• \(y_{i,j}\): ingreso total en escala logarítmica de la persona \(i\) en el dominio \(j\).

Los modelos presentados anteriormente se ajustan por medio de muestreadores de Gibbs con \(B = 11000\) iteraciones. Las primeras 1000 iteraciones del algoritmo constituyen el periodo de calentamiento del algoritmo, de manera que no se tienen en cuenta para realizar las inferencias. Para tal fin se emplean distribuciones previas empíricas difusas definidas por los siguientes hiperparámetros:

• \(\textsf{M}_1\): \(\mu_0 = 13.495\), \(\gamma_0^2 = 11.382\), \(\nu_0 = 1\), \(\sigma^2_0 = 1.182\).
• \(\textsf{M}_2\): \(\mu_0 = 13.495\), \(\gamma_0^2 = 11.382\), \(\eta_0 = 1\), \(\tau^2_0 = 1.182\), \(\nu_0 = 1\), \(\sigma^2_0 = 1.182\).
• \(\textsf{M}_3\): \(\mu_0 = 13.495\), \(\gamma_0^2 = 11.382\), \(\eta_0 = 1\), \(\tau^2_0 = 1.182\), \(\nu = 1\), \(\alpha_0 = 1\), \(\beta_0 = 0.846\).
• \(\textsf{M}_4\): \(\mu_0 = 13.495\), \(\gamma_0^2 = 11.382\), \(\alpha_0 = 1\), \(\beta_0 = 0.846\), \(\kappa = 3\).
• \(\textsf{M}_5\): \(\mu_0 = 13.495\), \(\gamma_0^2 = 11.382\), \(\eta_0 = 1\), \(\tau^2_0 = 1.182\), \(\alpha_0 = 1\), \(\beta_0 = 0.846\), \(\kappa = 3\).
• \(\textsf{M}_6\): \(\mu_0 = 13.495\), \(\gamma_0^2 = 11.382\), \(\eta_0 = 1\), \(\tau^2_0 = 1.182\), \(\alpha = 1\), \(a_\beta = 1\), \(b_\beta = 1.182\), \(\kappa = 3\).

Los hiperparámetros se establecieron teniendo en cuenta que la media muestral es \(\bar{y} = 13.495\), la varianza muestral es \(s^2_y = 1.182\) y el inverso de la varianza muestral es \(1/s^2_y = 0.846\).

En este repositorio de GitHub se encuentra la código de R para ajustar todos los modelos. A continuación se presenta la cadena de la log-verosimilitud de cada \(\textsf{M}_k\), para \(k=1,\ldots,6\).

Se observa que en todos los casos, basados en la log-verosimilitud, no hay evidencia de falta de convergencia.

5 Comparación de modelos

Con el fin de comparar predictivamente los modelos, se calcula el Criterio de Información de la Devianza (DIC, Deviance Information Criterion) de cada \(\textsf{M}_k\), para \(k=1,\ldots,6\). El DIC es una versión Bayesiana del Criterio de Información de Akaike (AIC, Akaike Infomation Criterion) y se define como: \[ \text{DIC} = -2\,\text{log}\,p(\boldsymbol{y}\mid\hat{\boldsymbol{\theta}}_{\text{Bayes}}) + 2p_{\text{DIC}} \] donde \[ \hat{\boldsymbol{\theta}}_{\text{Bayes}} = \textsf{E}(\boldsymbol{\theta}\mid\boldsymbol{y})\approx\frac{1}{B}\sum_{b=1}^B \boldsymbol{\theta}^{(b)} \] es la media posterior de \(\boldsymbol{\theta}\) y \[ p_{\text{DIC}} = 2\left( \text{log}\,p(\boldsymbol{y}\mid\hat{\boldsymbol{\theta}}_{\text{Bayes}}) - \textsf{E}( \text{log}\,p(\boldsymbol{y}\mid\boldsymbol{\theta}) \mid \boldsymbol{y} ) \right) \approx 2\left( \text{log}\,p(\boldsymbol{y}\mid\hat{\boldsymbol{\theta}}_{\text{Bayes}}) - \frac{1}{B}\sum_{b=1}^B \text{log}\,p\left(\boldsymbol{y}\mid \boldsymbol{\theta}^{(b)}\right) \right) \] es el número efectivo de parámetros.

DIC y número efectivo de parámetros de cada modelo.

	DIC	pDIC
M1	11480.84	1.99
M2	11276.78	24.26
M3	11276.25	47.46
M4	11428.46	1.96
M5	11168.09	24.66
M6	11165.88	47.74

Se observa que el mejor modelo en términos predictivos (penalizado por el número efectivo de parámetros) es el Modelo 6, dado que este modelo minimiza el DIC. También se observa que el número efectivo de parámetros corresponde es aproximadamente 48.

Referencias:

• Spiegelhalter, D. J., Best, N. G., Carlin, B. P., & Van Der Linde, A. (2002). Bayesian measures of model complexity and fit. Journal of the royal statistical society: Series b (statistical methodology), 64(4), 583-639.
• Spiegelhalter, D. J., Best, N. G., Carlin, B. P., & Van der Linde, A. (2014). The deviance information criterion: 12 years on. Journal of the Royal Statistical Society: Series B: Statistical Methodology, 485-493.

6 Bondad de ajuste de \(\textsf{M}_6\) en Bogotá D.C.

Es posible chequear la bondad de ajuste interna del modelo por medio de estadísticos de prueba calculados a partir de la distribución predictiva posterior del estadístico.

Si el estadístico es un valor típico de la distribución predictiva posterior, entonces se dice que el modelo captura adecuadamente la característica de interés que representa el estadístico de prueba. Esta característica se puede cuantificar por medio del valor \(p\) predictivo posterior (ppp), esto es, la probabilidad posterior de que el estadístico sea mayor que el valor observado correspondiente.

Los resultados correspondientes para Bogotá asociados con la media, la mediana, la desviación estándar, el coeficiente de variación, el rango y el rango intercuartílico usando \(\textsf{M}_6\) se presentan a continuación.

Valores \(p\) predictivos posteriores (ppp). DE: desviación estándar. CV: coeficiente de variación. R: rango. RI: rango intercuartílico.

	Media	Mediana	DE	CV	R	RI
M_6	0.861	0.521	0.7	0.677	0.904	0.525

El Modelo 6 permite caracterizar apropiadamente la media, la mediana, la desviación estándar, el coeficiente de variación, el rango y el rango intercuartílico, dado que ningún valor ppp es extremo (menor a 5% o superior a 95%).

7 Ranking

A continuación se presenta un ranking Bayesiano de los dominios basado en los efectos promedio \(\theta_1,\ldots,\theta_m\) de \(\textsf{M}_6\). La siguiente visualización incluye simultáneamente las estimaciones puntuales y los intervalos de credibilidad al 95% de cada \(\theta_j\), para \(j=1,\ldots,m\).

Ranking de los dominios basado en los efectos promedio del Modelo 4. En rojo: efectos promedio significativamente inferiores a 13.830. En negro: efectos promedio que no difieren significativamente de 13.830. En verde: efectos promedio significativamente superiores a 13.830. Se observa que 13.830 corresponde a un SMLMV de 2022 en escala logarítmica (línea vertical en color gris).

El Top 5 del ranking está conformado por 1. Medellin, 2. Manizales, 3. Tunja, 4. Bogotá, 5. Bucaramanga.Solamente Medellin y Manizales tiene efectos promedio significativamente superiores a 1 SMLMV de 2022 (valor de referencia). Aunque los efectos promedio de Tunja y Bogotá no difieren significativamente del valor de referencia, las estimaciones puntuales correspondientes son claramente superiores. Similarmente, el efecto promedio de Bucaramanga no difiere significativamente del valor de referencia, pero su estimación puntual es aproximadamente igual a este valor. Las últimas posiciones del ranking las ocupan 21. Quibdo, 22. Rioacha, 23. Resto Urbano, 24. Sincelejo, 25. Rural. Este último dominio (Rural) se encuentra muy por debajo de los demás dominios. Finalmente, la incertidumbre asociada con las estimaciones de los efectos promedio, la cual se ve reflejada en la amplitud de los intervalos, varia considerablemente.

8 Top 5

A continuación se estima puntualmente la media, la desviación estándar y el coeficiente de variación de los ingresos para el Top 5 del ranking usando \(\textsf{M}_6\). Para lograr tal fin en escala real se emplea el método delta dado que todos los modelos se ajustaron en escala log.

Estimaciones de la media, la desviación estándar (DE) y el coeficiente de variación (CV) de los ingresos para el Top 5 del ranking usando el Modelo 6. Las estimaciones de la media y la desviación estándar estánd dadas en la escala real en pesos y la del CV en puntos porcentuales.

	Ranking	Media	DE	CV (%)
MEDELLIN	1	1227401	1425492	116.1
MANIZALES	2	1158390	1111272	95.9
TUNJA	3	1137319	1342489	118.0
BOGOTA	4	1131058	1570543	138.9
BUCARAMANGA	5	1023315	1325420	129.5

Se observa que la volatilidad (variabilidad) de los ingresos en estos dominios es muy alta a pesar de ocupar las primeras posiciones del ranking.

9 Segmentación

Con el fin de segmentar los 25 dominios de acuerdo con los ingresos totales de las personas, se conforma un arreglo que contenga las estimaciones puntuales de los \(\theta_j\) y los \(\sigma_j\) de todos los dominios usando \(\textsf{M}_6\). A continuación se usa este arreglo como insumo, para segmentar los dominios por medio de agrupación jerárquica con cuatro grupos.

Repitiendo este mismo protocolo en cada iteración del muestreador de Gibbs se obtiene que los mismos dominios tienen una probabilidad posterior alta de pertenecer al mismo grupo (colores más oscuros), pero se distinguen dos grandes grupos, uno del dominio 1 al 14 y otro del dominio 15 al 25 (usando el ordenamiento del ranking).