Descripción de las personas encuestadas según la pertenencia a programas sociales en Colombia

Análisis de la Encuesta de Calidad de Vida 2023

Mostrar código 
rm(list=ls())

Ruta1 <- c("D:\\OneDrive - Universidad Nacional de Colombia\\Maestr_Est_Ap\\ECV\\2023\\1BD")
Ruta2 <-c("C:\\Users\\jaarc\\OneDrive - Universidad Nacional de Colombia\\Maestr_Est_Ap\\ECV\\2023\\1BD")
setwd(Ruta2)

library(ggplot2)
library(knitr)
library(kableExtra)
library(dplyr)
library(survey)
library(srvyr)

options(scipen = 999)

bd <- readRDS("bd1.rds")
bd0 <- bd|>
  dplyr::filter(!is.na(PROG_DUMMY))

bd01 <- bd0|>
  dplyr::filter(PERCAPITA>0 & !is.na(PERCAPITA))

bd02 <- bd0|>
  dplyr::filter(nivel_edu4>0 & !is.na(nivel_edu4))

dim(bd)
[1] 240212    101

Introducción

En Colombia operan diversos programas sociales dirigidos a poblaciones específicas, con objetivos y estrategias diferenciadas. Algunos están enfocados en personas mayores, otros en niños, jóvenes o en hogares en situación de vulnerabilidad económica, entre otros grupos. El análisis de las características de estas poblaciones suele basarse en registros administrativos particulares de cada programa, lo que limita las posibilidades de comparación con la población general.

La Encuesta Nacional de Calidad de Vida (ECV) 2023 permite identificar a personas que reportan ser beneficiarias de distintos programas sociales. No obstante, debido a la forma en que se aborda esta temática en el instrumento, no es posible conocer con precisión el nombre de muchos de estos programas (DANE, 2024).

Dado lo anterior, este ejercicio propone caracterizar a los beneficiarios y no beneficiarios de programas sociales en Colombia, utilizando la ECV 2023 y considerando su diseño muestral complejo. Esta aproximación busca aportar evidencia que complemente las fuentes administrativas, permitiendo una visión más integral de la población beneficiaria a nivel nacional.

Metodología

El presente ejercicio de naturaleza retrospectiva sigue un diseño de corte-transversal a partir de la ECV para el año 2023. De acuerdo con el DANE (2024), el universo de este instrumento refiere a la población civil no institucional residente en todo el territorio nacional, exceptuando la población de Providencia y el centro poblado y rural disperso de San Andrés.

En este contexto, las estimaciones se realizarán teniendo en cuenta el diseño del diseño muestral complejo de la ECV para 2023, siguiendo la estrategia planteada por Gutiérrez et al. (2024). De acuerdo a lo descrito en la ficha metodológica de la encuesta, se encuentran que como estratos grupos de municipios, en donde adicionalmente se encuentra la ciudad de Bogotá por inclusión forzosa; así mismo, se describe que hay una estratificación a nivel de área (cabecera, centro poblado y rural disperso) (DANE, 2024).

A continuación, se ilustra la declaración de dicho diseño muestral en la base de datos de la encuesta:

Mostrar código 
bd_dis <- bd %>% # Base de datos.
  as_survey_design(
    strata = c(ESTRATO2020, CLASE), # Id de los estratos.
    ids = c(DIRECTORIO, SECUENCIA_P, ORDEN), # Id para las observaciones. #
    weights = FEX_C, # Factores de expansión.
    pps="brewer",
    nest = TRUE # Anidado dentro del estrato
  )

bd_dis2 <- bd0 %>% # Base de datos.
  as_survey_design(
    strata = c(ESTRATO2020, CLASE), # Id de los estratos.
    ids = c(DIRECTORIO, SECUENCIA_P, ORDEN), # Id para las observaciones. #
    weights = FEX_C, # Factores de expansión.
    pps="brewer",
    nest = TRUE # Anidado dentro del estrato
  )

bd_dis3 <- bd01 %>% # Base de datos.
  as_survey_design(
    strata = c(ESTRATO2020, CLASE), # Id de los estratos.
    ids = c(DIRECTORIO, SECUENCIA_P, ORDEN), # Id para las observaciones. #
    weights = FEX_C, # Factores de expansión.
    pps="brewer",
    nest = TRUE # Anidado dentro del estrato
  )


bd_dis4 <- bd02 %>% # Base de datos.
  as_survey_design(
    strata = c(ESTRATO2020, CLASE), # Id de los estratos.
    ids = c(DIRECTORIO, SECUENCIA_P, ORDEN), # Id para las observaciones. #
    weights = FEX_C, # Factores de expansión.
    pps="brewer",
    nest = TRUE # Anidado dentro del estrato
  )

Resultados

Total población

La población representada en la ECV para el año 2023 asciende a los 52.3 millones de personas, consistentemente con lo descrito por el DANE (DANE, 2024). Puntualmente, se encuentra que un poco más de 5 millones de personas han sido beneficiarias de alguno de los programas sociales en el país (coeficiente de variación [CV]: 0.5%), tal como se ilustra en la siguiente tabla (Table 1). Las personas sobre quienes no se indagó corresponden a más de 30 millones, las cuales serán excluidas de los análisis posteriores.

En este contexto, las estimaciones sobre las personas beneficiarias o no de los programas sociales en Colombia corresponden a más de 21 millones, de las cuales el equivalente al 23.7% mencionó haber recibido algún apoyo gubernamental (Table 1).

Mostrar código 
bd_dis %>%
  mutate(PROG_DUMMY = as.character(PROG_DUMMY))%>%
group_by(PROG_DUMMY) %>%
  cascade(
    N_d = survey_total(vartype = c("cv", "se","ci")),
    .fill = "Total")%>%
  kable(format.args = list(decimal.mark = ',', big.mark = "."), digits=2)%>%
  kable_styling(bootstrap_options = "striped", font_size = 12, full_width = F)%>%
  kable_classic(full_width = F, html_font = "Arial Narrow")%>%
  row_spec(0, color = "white", background = "navy", align = "center", bold= T)
bd_dis %>%
  filter(!is.na(PROG_DUMMY))|>
  mutate(PROG_DUMMY = as.character(PROG_DUMMY))%>%
group_by(PROG_DUMMY) %>%
  cascade(
    Prop_d = survey_prop(vartype = c("cv", "se","ci")),
    .fill = "Total")%>%
  kable(format.args = list(decimal.mark = ',', big.mark = "."), digits=3)%>%
  kable_styling(bootstrap_options = "striped", font_size = 12, full_width = F)%>%
  kable_classic(full_width = F, html_font = "Arial Narrow")%>%
  row_spec(0, color = "white", background = "navy", align = "center", bold= T)
Table 1: Personas en Colombia según recepción de programa social
(a) Total
PROG_DUMMY N_d N_d_cv N_d_se N_d_low N_d_upp
0 16.177.741 0,01 100.758,57 15.980.255 16.375.227
1 5.022.700 0,01 58.196,23 4.908.636 5.136.764
Total 52.314.391 0,01 287.228,61 51.751.425 52.877.357
NA 31.113.951 0,01 215.112,07 30.692.333 31.535.569
(b) Distribución porcentual
PROG_DUMMY Prop_d Prop_d_cv Prop_d_se Prop_d_low Prop_d_upp
0 0,763 0,003 0,002 0,758 0,768
1 0,237 0,010 0,002 0,232 0,242
Total 1,000 0,000 0,000 1,000 1,000
NA 0,000 0,000 0,000 0,000 0,000

Población estudiada

Específicamente, más de 4 millones de personas se han beneficiado de 1 programa social (CV 1.2%), y alrededor de 600 personas han recibido apoyos de 4 programas. No obstante, frente a estos últimos las estimaciones asociadas presentan alta incertidumbre, considerando que su CV está alrededor del 50% (Table 2).

En este sentido, las personas beneficiarias de 1 programa social representan el 21.5% (CV 1.2%) de la población estudiada.

Mostrar código 
bd_dis2 %>%
  mutate(PROG_NUM = as.character(PROG_NUM))%>%
group_by(PROG_NUM) %>%
  cascade(
    N_d = survey_total(vartype = c("cv", "se","ci")),
    .fill = "Total")%>%
  kable(format.args = list(decimal.mark = ',', big.mark = "."), digits=2)%>%
  kable_styling(bootstrap_options = "striped", font_size = 12, full_width = F)%>%
  kable_classic(full_width = F, html_font = "Arial Narrow")%>%
  row_spec(0, color = "white", background = "navy", align = "center", bold= T)
bd_dis2 %>%
  mutate(PROG_NUM = as.character(PROG_NUM))%>%
group_by(PROG_NUM) %>%
  cascade(
    Prop_d = survey_prop(vartype = c("cv", "se","ci")),
    .fill = "Total")%>%
  kable(format.args = list(decimal.mark = ',', big.mark = "."), digits=3)%>%
  kable_styling(bootstrap_options = "striped", font_size = 12, full_width = F)%>%
  kable_classic(full_width = F, html_font = "Arial Narrow")%>%
  row_spec(0, color = "white", background = "navy", align = "center", bold= T)
Table 2: Número de programas a los que han sido beneficiarios
(a) Total
PROG_NUM N_d N_d_cv N_d_se N_d_low N_d_upp
0 16.177.740,56 0,01 100.758,57 15.980.254,60 16.375.226,52
1 4.561.736,00 0,01 54.427,01 4.455.059,52 4.668.412,48
2 434.725,92 0,04 17.072,38 401.264,20 468.187,65
3 25.631,65 0,14 3.590,07 18.595,15 32.668,15
4 606,09 0,50 302,67 12,86 1.199,32
Total 21.200.440,22 0,01 111.178,30 20.982.531,68 21.418.348,76
(b) Distribución porcentual
PROG_NUM Prop_d Prop_d_cv Prop_d_se Prop_d_low Prop_d_upp
0 0,763 0,003 0,002 0,758 0,768
1 0,215 0,011 0,002 0,211 0,220
2 0,021 0,039 0,001 0,019 0,022
3 0,001 0,140 0,000 0,001 0,002
4 0,000 0,499 0,000 0,000 0,000
Total 1,000 0,000 0,000 1,000 1,000

Características demográficas

Sexo

Según sexo se observan ligeras diferencias, se evidencia que una mayor proporción de personas no beneficiarias son hombres (55.8%, CV <1%), mientras que en el caso de las personas beneficiarias tan solo el 48% eran de sexo masculino (CV 1.1%) (Table 3).

Mostrar código 
bd_dis2 %>%
  mutate(PROG_DUMMY = as.character(PROG_DUMMY),
         Sexo=as.character(Sexo))%>%
group_by(PROG_DUMMY, Sexo) %>%
  cascade(
    N_d = survey_total(vartype = c("cv", "se","ci"), na.rm=TRUE),
    .fill = "Total")%>%
  kable(format.args = list(decimal.mark = ',', big.mark = "."), digits=2)%>%
  kable_styling(bootstrap_options = "striped", font_size = 12, full_width = F)%>%
  kable_classic(full_width = F, html_font = "Arial Narrow")%>%
  row_spec(0, color = "white", background = "navy", align = "center", bold= T)
bd_dis2 %>%
  mutate(PROG_DUMMY = as.character(PROG_DUMMY),
         Sexo=as.character(Sexo))%>%
  filter(PROG_DUMMY=="0")%>%
group_by(Sexo) %>%
  cascade(
    Prop_d = survey_prop(vartype = c("cv", "se","ci"), na.rm=TRUE),
    .fill = "Total")%>%
  kable(format.args = list(decimal.mark = ',', big.mark = "."), digits=3)%>%
  kable_styling(bootstrap_options = "striped", font_size = 12, full_width = F)%>%
  kable_classic(full_width = F, html_font = "Arial Narrow")%>%
  row_spec(0, color = "white", background = "navy", align = "center", bold= T)
bd_dis2 %>%
  mutate(PROG_DUMMY = as.character(PROG_DUMMY),
         Sexo=as.character(Sexo))%>%
  filter(PROG_DUMMY=="1")%>%
group_by(Sexo) %>%
  cascade(
    Prop_d = survey_prop(vartype = c("cv", "se","ci"), na.rm=TRUE),
    .fill = "Total")%>%
  kable(format.args = list(decimal.mark = ',', big.mark = "."), digits=3)%>%
  kable_styling(bootstrap_options = "striped", font_size = 12, full_width = F)%>%
  kable_classic(full_width = F, html_font = "Arial Narrow")%>%
  row_spec(0, color = "white", background = "navy", align = "center", bold= T)
Table 3: Sexo según pertenencia a programas sociales
(a) Total
PROG_DUMMY Sexo N_d N_d_cv N_d_se N_d_low N_d_upp
0 Hombre 9.026.435 0,01 83.568,07 8.862.642 9.190.227
0 Mujer 7.151.306 0,01 72.044,82 7.010.099 7.292.513
0 Total 16.177.741 0,01 100.758,57 15.980.255 16.375.227
1 Hombre 2.406.388 0,02 37.159,37 2.333.556 2.479.220
1 Mujer 2.616.311 0,02 41.881,87 2.534.223 2.698.400
1 Total 5.022.700 0,01 58.196,23 4.908.636 5.136.764
Total Total 21.200.440 0,01 111.178,30 20.982.532 21.418.349
(b) Distribución porcentual: No beneficiarios
Sexo Prop_d Prop_d_cv Prop_d_se Prop_d_low Prop_d_upp
Hombre 0,558 0,007 0,004 0,551 0,565
Mujer 0,442 0,008 0,004 0,435 0,449
Total 1,000 0,000 0,000 1,000 1,000
(c) Distribución porcentual: Beneficiarios
Sexo Prop_d Prop_d_cv Prop_d_se Prop_d_low Prop_d_upp
Hombre 0,479 0,011 0,005 0,469 0,490
Mujer 0,521 0,010 0,005 0,510 0,531
Total 1,000 0,000 0,000 1,000 1,000

Así mismo, se cuenta con evidencia para afirmar que la pertenencia entre la pertenencia a un programa social y el sexo no son independientes como lo ilustra la siguiente Table 4. Es decir, que el sexo se perfila como un factor determinante para ser incluido en los programas sociales en el país. Esto puede estar explicado en gran medida porque muchos programas están dirigidos a personas cabeza de hogar -típicamente mujeres-, como familias en acción (cf. https://dps2018.prosperidadsocial.gov.co/inf/not/Paginas/Mujeres-beneficiarias-de-Prosperidad-Social,-protagonistas-en-reducci%C3%B3n-de-pobreza.aspx)

Mostrar código 
svychisq(~Sexo + PROG_DUMMY, bd_dis2, statistic="F")
Table 4: Prueba de independencia: Sexo y pertenencia a programas sociales 

    Pearson's X^2: Rao & Scott adjustment

data:  NextMethod()
F = 146.43, ndf = 1, ddf = 85707, p-value < 0.00000000000000022

Edad

La edad promedio fue de 42.7 años (CV <1%), en donde las personas beneficiarias tendieron a ser menores que las no beneficiarias, con una edad promedio de 37.9 y 44.2 años, respectivamente. Consistentemente, la mediana tendió a situarse cerca del promedio, tal como se ilustra en la siguiente tabla (Table 5).

Por su parte, la desviación estándar permite evidenciar una mayor dispersión de la edad en el caso de las personas beneficiarias (25.7 años en los beneficiarios vs. 18.9 en los no beneficiarios).

Mostrar código 
bd_dis2 %>%
  mutate(PROG_DUMMY = as.character(PROG_DUMMY))%>%
group_by(PROG_DUMMY) %>%
  cascade(
    Mean_d = survey_mean(Edad, vartype = c("cv", "se","ci"), na.rm=TRUE),
    Median_d=survey_median(Edad, vartype = c("cv", "se"), na.rm=TRUE),
    Quantile_d=survey_quantile(Edad, quantiles=c(0.25, 0.75), vartype = c("cv"), na.rm=TRUE),
    Sd = sqrt(survey_var(Edad, level = 0.95, vartype = c("se"))),
    .fill = "Total")%>%
  kable(format.args = list(decimal.mark = ',', big.mark = "."), digits=2)%>%
  kable_styling(bootstrap_options = "striped", font_size = 12, full_width = F)%>%
  kable_classic(full_width = F, html_font = "Arial Narrow")%>%
  row_spec(0, color = "white", background = "navy", align = "center", bold= T)
Table 5: Edad según pertenencia a programas sociales
PROG_DUMMY Mean_d Mean_d_cv Mean_d_se Mean_d_low Mean_d_upp Median_d Median_d_cv Median_d_se Quantile_d_q25 Quantile_d_q75 Quantile_d_q25_cv Quantile_d_q75_cv Sd Sd_se
0 44,20 0,00 0,16 43,89 44,52 44 0,01 0,26 31 58 0,01 0,00 18,86 1,85
1 37,85 0,01 0,29 37,28 38,42 37 0,01 0,26 15 59 0,09 0,01 25,69 2,56
Total 42,70 0,00 0,14 42,42 42,98 43 0,01 0,26 29 58 0,01 0,00 20,86 1,73

Visualmente, el histograma y el boxplot permiten ilustrar las tendencias descritas anteriormente; además evidencia una ligera concentración en niños en la población estudiada (Figure 1).

Mostrar código 
svyhist(~ Edad , bd_dis2,
main = "", col = "grey80",
xlab = "Edad (años)", probability = FALSE )

svyboxplot( formula= Edad~as.character(PROG_DUMMY) , design=bd_dis2,
col = c("goldenrod", "forestgreen"),
ylab = "Edad", xlab = "Beneficiarios o no" )
(a) Histograma
(b) Boxplot
Figure 1: Gráficos de la edad

Una prueba t, señala que las diferencias en la edad entre las personas beneficiarias y las no beneficiarias son estadísticamente significativas (Table 6). Lo anterior implica que los programas sociales en el país han tendido a concentrarse en personas más jóvenes que las no beneficiarias.

Mostrar código 
svyttest(Edad ~ as.character(PROG_DUMMY), bd_dis2)
Table 6: T-test: Edad según pertenencia a programas sociales 

    Design-based t-test

data:  Edad ~ as.character(PROG_DUMMY)
t = -19.585, df = 85706, p-value < 0.00000000000000022
alternative hypothesis: true difference in mean is not equal to 0
95 percent confidence interval:
 -6.988953 -5.717383
sample estimates:
difference in mean 
         -6.353168

Años de escolaridad

Frente a la escolaridad de la población en estudio, se encuentra que el promedio para el país en 2023 fue de 10.6 años (CV < 1%) (Table 7). Particularmente, las personas beneficiarias de programas sociales presentaron un promedio de 8 años de escolaridad (CV 1%), mientras que las no beneficiarias alcanzaron un promedio de 11.2 años (CV <1%). Con relación a la mediana, se observa que las personas beneficiarias de programas sociales esta fue de 6 años (CV 17%), mientras que las no beneficiarias alcanzaron una mediana de 12 años (CV 2%) (Table 7).

Mostrar código 
bd_dis2 %>%
  mutate(PROG_DUMMY = as.character(PROG_DUMMY))%>%
group_by(PROG_DUMMY) %>%
  cascade(
    Mean_d = survey_mean(nivel_edu4, vartype = c("cv", "se","ci"), na.rm=TRUE),
    Median_d=survey_median(nivel_edu4, vartype = c("cv", "se"), na.rm=TRUE),
    Quantile_d=survey_quantile(nivel_edu4, quantiles=c(0.25, 0.75), vartype = c("cv"), na.rm=TRUE),
    Sd = sqrt(survey_var(nivel_edu4, level = 0.95, vartype = c("se"))),
    .fill = "Total")%>%
  kable(format.args = list(decimal.mark = ',', big.mark = "."), digits=2)%>%
  kable_styling(bootstrap_options = "striped", font_size = 12, full_width = F)%>%
  kable_classic(full_width = F, html_font = "Arial Narrow")%>%
  row_spec(0, color = "white", background = "navy", align = "center", bold= T)
Table 7: Años de escolaridad según pertenencia a programas sociales
PROG_DUMMY Mean_d Mean_d_cv Mean_d_se Mean_d_low Mean_d_upp Median_d Median_d_cv Median_d_se Quantile_d_q25 Quantile_d_q75 Quantile_d_q25_cv Quantile_d_q75_cv Sd Sd_se
0 11,21 0,00 0,03 11,15 11,28 12 0,02 0,26 6 14 0,17 0,02 NA NaN
1 8,04 0,01 0,05 7,94 8,14 6 0,17 1,02 6 12 0,17 0,02 NA NaN
Total 10,61 0,00 0,03 10,55 10,67 12 0,02 0,26 6 13 0,17 0,02 NA NaN

Este comportamiento se ilustra en la siguiente figura (Figure 2), el panel (a) retrata el histograma de los años de escolaridad para toda la población, mientras que el panel (b) ilustra el boxplot según pertenencia a algún programa social. Puntualmente en toda la población de estudio, se destaca una concentración no solo cercana a los 10 años de escolaridad sino cercana a los 6 años, lo que puede estar correlacionado con el número de años asociados al logro de completar la educación primaria (Figure 2).

Mostrar código 
svyhist(~ nivel_edu4 , bd_dis2,
main = "", col = "grey80",
xlab = "Años de escolaridad", probability = FALSE )

svyboxplot( formula= nivel_edu4~as.character(PROG_DUMMY) , design=bd_dis2,
col = c("goldenrod", "forestgreen"),
ylab = "Años de escolaridad", xlab = "Beneficiarios o no" )
(a) Histograma
(b) Boxplot
Figure 2: Gráficos de los años de escolaridad

Una prueba t, señala que las diferencias en los años de escolaridad entre las personas beneficiarias y las no beneficiarias son estadísticamente significativas (Table 8). Lo anterior implica que los programas sociales en el país han tendido a concentrarse en personas con menos años de escolaridad que las no beneficiarias.

Mostrar código 
svyttest(nivel_edu4 ~ as.character(PROG_DUMMY), bd_dis2)
Table 8: T-test: Años de escolaridad según pertenencia a programas sociales 

    Design-based t-test

data:  nivel_edu4 ~ as.character(PROG_DUMMY)
t = -51.099, df = 83787, p-value < 0.00000000000000022
alternative hypothesis: true difference in mean is not equal to 0
95 percent confidence interval:
 -3.294944 -3.051515
sample estimates:
difference in mean 
         -3.173229

Número de dificultades físicas

Con relación al número de dificultades físicas, en donde se encuentran dificultades relacionadas con los sentidos, funciones ejecutivas y habiliades motoras como la visión, audición, habla, movilidad, aprendizaje e interactuar con otros. Se observa que promedio en la población de estudio fue de 0.4 dificultades (CV 1%),en donde las personas beneficiarias tuvieron un promedio ligeramente superior (0.51, CV 2%) y las personas no beneficiarias levemente inferior al promedio de la población de estudio (0.36, CV 2%) (Table 9).

Mostrar código 
bd_dis2 %>%
  mutate(PROG_DUMMY = as.character(PROG_DUMMY))%>%
group_by(PROG_DUMMY) %>%
  cascade(
    Mean_d = survey_mean(num_dif_fisica, vartype = c("cv", "se","ci"), na.rm=TRUE),
    Median_d=survey_median(num_dif_fisica, vartype = c("cv", "se"), na.rm=TRUE),
    Quantile_d=survey_quantile(num_dif_fisica, quantiles=c(0.25, 0.75), vartype = c("cv"), na.rm=TRUE),
    Sd = sqrt(survey_var(num_dif_fisica, level = 0.95, vartype = c("se"))),
    .fill = "Total")%>%
  kable(format.args = list(decimal.mark = ',', big.mark = "."), digits=2)%>%
  kable_styling(bootstrap_options = "striped", font_size = 12, full_width = F)%>%
  kable_classic(full_width = F, html_font = "Arial Narrow")%>%
  row_spec(0, color = "white", background = "navy", align = "center", bold= T)
Table 9: Número de dificultades físicas según pertenencia a programas sociales
PROG_DUMMY Mean_d Mean_d_cv Mean_d_se Mean_d_low Mean_d_upp Median_d Median_d_cv Median_d_se Quantile_d_q25 Quantile_d_q75 Quantile_d_q25_cv Quantile_d_q75_cv Sd Sd_se
0 0,36 0,02 0,01 0,35 0,37 0 Inf 0,26 0 0 Inf Inf 0,91 0,16
1 0,51 0,02 0,01 0,49 0,53 0 Inf 0,26 0 1 Inf 0,26 1,16 0,21
Total 0,40 0,01 0,01 0,38 0,41 0 Inf 0,26 0 0 Inf Inf 0,97 0,15

Análogamente, la figura 3 (Figure 3) presenta el comportamiento de esta variable, en donde se destaca una alta concentración de personas con 0 dificultades físicas.

Mostrar código 
svyhist(~ num_dif_fisica , bd_dis2,
main = "", col = "grey80",
xlab = "Número de dificultades físicas", probability = FALSE )

svyboxplot( formula= num_dif_fisica~as.character(PROG_DUMMY) , design=bd_dis2,
col = c("goldenrod", "forestgreen"),
ylab = "Número de dificultades físicas", xlab = "Beneficiarios o no" )
(a) Histograma
(b) Boxplot
Figure 3: Gráficos del Número de dificultades físicas

Una prueba t, para evaluar si las diferencias entre los beneficiarios y no beneficiarios son estadísticamente significativas, señala que efectivamente las diferencias en el número de dificultades físicas entre las personas beneficiarias y las no beneficiarias son estadísticamente significativas (Table 10). Por lo que se argumenta que los programas sociales en el país han tendido a concentrarse en personas con más dificultades físicas que las no beneficiarias.

Mostrar código 
svyttest(num_dif_fisica ~ as.character(PROG_DUMMY), bd_dis2)
Table 10: T-test: Años de escolaridad según pertenencia a programas sociales 

    Design-based t-test

data:  num_dif_fisica ~ as.character(PROG_DUMMY)
t = 11.929, df = 85706, p-value < 0.00000000000000022
alternative hypothesis: true difference in mean is not equal to 0
95 percent confidence interval:
 0.1240047 0.1727669
sample estimates:
difference in mean 
         0.1483858

Características económicas

Ingresos

Respecto a los ingresos, se evidencia que en promedio las personas de la población de estudio reportan ingresos per cápita de $1.237.974 (CV 2%), en donde puntualmente, las personas no beneficiarias reportan un ingreso per cápita superior ($1.396.141, CV 2%) al de las personas beneficiarias ($573.176, CV 3%) (Table 11). Es decir, las personas no beneficiarias presentan en promedio un ingreso per cápita 2.4 veces superior al de las personas beneficiarias de programas sociales en el país.

Mostrar código 
bd_dis2 %>%
  mutate(PROG_DUMMY = as.character(PROG_DUMMY))%>%
group_by(PROG_DUMMY) %>%
  cascade(
    Mean_d = survey_mean(PERCAPITA, vartype = c("cv", "se","ci"), na.rm=TRUE),
    Median_d=survey_median(PERCAPITA, vartype = c("cv", "se"), na.rm=TRUE),
    Quantile_d=survey_quantile(PERCAPITA, quantiles=c(0.25, 0.75), vartype = c("cv"), na.rm=TRUE),
    Sd = sqrt(survey_var(PERCAPITA, level = 0.95, vartype = c("se"))),
    .fill = "Total")%>%
  kable(format.args = list(decimal.mark = ',', big.mark = "."), digits=2)%>%
  kable_styling(bootstrap_options = "striped", font_size = 12, full_width = F)%>%
  kable_classic(full_width = F, html_font = "Arial Narrow")%>%
  row_spec(0, color = "white", background = "navy", align = "center", bold= T)
Table 11: Ingresos según pertenencia a programas sociales
PROG_DUMMY Mean_d Mean_d_cv Mean_d_se Mean_d_low Mean_d_upp Median_d Median_d_cv Median_d_se Quantile_d_q25 Quantile_d_q75 Quantile_d_q25_cv Quantile_d_q75_cv Sd Sd_se
0 1.396.141,0 0,02 28.971,16 1.339.357,8 1.452.924,2 755.000,0 0,01 6.356,29 400.000,0 1.485.000,0 0,01 0,01 NA NaN
1 573.176,1 0,03 15.809,62 542.189,4 604.162,8 376.666,7 0,01 3.826,38 232.666,7 594.166,7 0,01 0,01 NA NaN
Total 1.237.974,4 0,02 23.639,95 1.191.640,3 1.284.308,5 648.888,9 0,01 4.145,43 337.500,0 1.300.000,0 0,01 0,01 NA NaN

Gráficamente se puede evidenciar una asimetría positiva de la variable ingresos, estándo altamente concentrados alrededor del cero y que en este contexto las medianas según pertenencia a algún programa social son mucho menores al promedio (Figure 4).

Mostrar código 
svyhist(~ PERCAPITA , bd_dis2,
main = "", col = "grey80",
xlab = "Ingreso per cápita", probability = FALSE )

svyboxplot( formula= PERCAPITA~as.character(PROG_DUMMY) , design=bd_dis2,
col = c("goldenrod", "forestgreen"),
ylab = "Ingreso per cápita", xlab = "Beneficiarios o no" )
(a) Histograma
(b) Boxplot
Figure 4: Gráficos del Número de dificultades físicas

Complementariamente, una prueba para evaluar si las diferencias de medias entre los grupos son estadísticamente significativas, indica que esta brecha en los ingresos entre las personas beneficiarias y las no beneficiarias son estadísticamente significativas (Table 12). En este contexto es posible señalar que los programas sociales en el país han tendido a concentrarse en personas con menos ingresos que las no beneficiarias.

Mostrar código 
svyttest(PERCAPITA ~ as.character(PROG_DUMMY), bd_dis2)
Table 12: T-test: Ingresos según pertenencia a programas sociales 

    Design-based t-test

data:  PERCAPITA ~ as.character(PROG_DUMMY)
t = -24.832, df = 85706, p-value < 0.00000000000000022
alternative hypothesis: true difference in mean is not equal to 0
95 percent confidence interval:
 -887920.7 -758009.0
sample estimates:
difference in mean 
         -822964.9

Índice de hacinamiento

Frente al índice de hacinamiento (número de personas por cuarto), se evidencia que el promedio en la población de estudio fue de 1.65 personas por cuarto (CV < 1%), en donde las personas beneficiarias presentaron un promedio superior (1.91, CV 1%) y las personas no beneficiarias un valor inferior al promedio de la población de estudio (1.58, CV < 1%) (Table 13).

Mostrar código 
bd_dis2 %>%
  mutate(PROG_DUMMY = as.character(PROG_DUMMY))%>%
group_by(PROG_DUMMY) %>%
  cascade(
    Mean_d = survey_mean(i_hacinamiento, vartype = c("cv", "se","ci"), na.rm=TRUE),
    Median_d=survey_median(i_hacinamiento, vartype = c("cv", "se"), na.rm=TRUE),
    Quantile_d=survey_quantile(i_hacinamiento, quantiles=c(0.25, 0.75), vartype = c("cv"), na.rm=TRUE),
    Sd = sqrt(survey_var(i_hacinamiento, level = 0.95, vartype = c("se"))),
    .fill = "Total")%>%
  kable(format.args = list(decimal.mark = ',', big.mark = "."), digits=2)%>%
  kable_styling(bootstrap_options = "striped", font_size = 12, full_width = F)%>%
  kable_classic(full_width = F, html_font = "Arial Narrow")%>%
  row_spec(0, color = "white", background = "navy", align = "center", bold= T)
Table 13: Ingresos según pertenencia a programas sociales
PROG_DUMMY Mean_d Mean_d_cv Mean_d_se Mean_d_low Mean_d_upp Median_d Median_d_cv Median_d_se Quantile_d_q25 Quantile_d_q75 Quantile_d_q25_cv Quantile_d_q75_cv Sd Sd_se
0 1,58 0,00 0,01 1,57 1,59 1,50 0,03 0,04 1,00 2 0,06 0,06 NA NaN
1 1,91 0,01 0,01 1,89 1,93 1,75 0,05 0,09 1,25 2 0,07 0,04 NA NaN
Total 1,65 0,00 0,00 1,64 1,66 1,50 0,03 0,04 1,00 2 0,06 0,04 NA NaN

En la siguiente figura (Figure 5) se ilustra el comportamiento de esta variable, en donde se destaca una alta concentración de personas con 1 persona por cuarto. De su lado, el cuartil 75 en todos los casos se encuentra en 2 personas por cuarto (Table 13).

Mostrar código 
svyhist(~ i_hacinamiento , bd_dis2,
main = "", col = "grey80",
xlab = "Índice de hacinamiento", probability = FALSE )

svyboxplot( formula= i_hacinamiento~as.character(PROG_DUMMY) , design=bd_dis2,
col = c("goldenrod", "forestgreen"),
ylab = "Índice de hacinamiento", xlab = "Beneficiarios o no" )
(a) Histograma
(b) Boxplot
Figure 5: Gráficos del Índice de hacinamiento

Las diferencias entre los beneficiarios y no beneficiarios en el índice de hacinamiento son estadísticamente significativas (Table 14).

Mostrar código 
svyttest(i_hacinamiento ~ as.character(PROG_DUMMY), bd_dis2)
Table 14: T-test: Índice de hacinamiento según pertenencia a programas sociales 

    Design-based t-test

data:  i_hacinamiento ~ as.character(PROG_DUMMY)
t = 27.581, df = 85706, p-value < 0.00000000000000022
alternative hypothesis: true difference in mean is not equal to 0
95 percent confidence interval:
 0.3046295 0.3512379
sample estimates:
difference in mean 
         0.3279337

Número de bienes

Consistentemente, con las características económicas presentadas anteriormente, las personas no beneficiarias tienen en promedio más bienes (6,2, CV < 1) que las beneficiarias de programas sociales en el país (4,8, CV 1%) (Table 15).

Mostrar código 
bd_dis2 %>%
  mutate(PROG_DUMMY = as.character(PROG_DUMMY))%>%
group_by(PROG_DUMMY) %>%
  cascade(
    Mean_d = survey_mean(num_bienes, vartype = c("cv", "se","ci"), na.rm=TRUE),
    Median_d=survey_median(num_bienes, vartype = c("cv", "se"), na.rm=TRUE),
    Quantile_d=survey_quantile(num_bienes, quantiles=c(0.25, 0.75), vartype = c("cv"), na.rm=TRUE),
    Sd = sqrt(survey_var(num_bienes, level = 0.95, vartype = c("se"))),
    .fill = "Total")%>%
  kable(format.args = list(decimal.mark = ',', big.mark = "."), digits=2)%>%
  kable_styling(bootstrap_options = "striped", font_size = 12, full_width = F)%>%
  kable_classic(full_width = F, html_font = "Arial Narrow")%>%
  row_spec(0, color = "white", background = "navy", align = "center", bold= T)
Table 15: Número de bienes según pertenencia a programas sociales
PROG_DUMMY Mean_d Mean_d_cv Mean_d_se Mean_d_low Mean_d_upp Median_d Median_d_cv Median_d_se Quantile_d_q25 Quantile_d_q75 Quantile_d_q25_cv Quantile_d_q75_cv Sd Sd_se
0 6,21 0,00 0,02 6,17 6,25 6 0,04 0,26 4 8 0,06 0,03 NA NaN
1 4,78 0,01 0,03 4,72 4,83 5 0,05 0,26 3 6 0,09 0,04 NA NaN
Total 5,94 0,00 0,02 5,90 5,97 6 0,04 0,26 4 8 0,06 0,03 NA NaN

Estas diferencias son estadísticamente significativas, tal como lo ilustra la siguiente prueba t (Table 16).

Mostrar código 
svyttest(num_bienes ~ as.character(PROG_DUMMY), bd_dis2)
Table 16: T-test: Número de bienes según pertenencia a programas sociales 

    Design-based t-test

data:  num_bienes ~ as.character(PROG_DUMMY)
t = -38.679, df = 85706, p-value < 0.00000000000000022
alternative hypothesis: true difference in mean is not equal to 0
95 percent confidence interval:
 -1.510107 -1.364446
sample estimates:
difference in mean 
         -1.437276

Análisis exploratorios

Correlaciones

A manera exploratoria, se presentan las correlaciones entre las variables de edad y el número de bienes según la pertenencia a programas sociales. Específicamente, se observa que a nivel general la correlación es de 0.03 positiva y estadísticamente positiva (Intevalo de confianza [IC] 95% 0.02; 0.04) aunque de baja intensidad (Table 17). No obstante al desagregar según la pertenencia a programas sociales, se observa que la correlación es negativa y estadísticamente significativa para el caso de las personas beneficiarias (-0.12; IC 95% -0.15; -0.10) y positiva y estadísticamente significativa para el caso de las personas no beneficiarias (0.09, IC 95% 0.08-0.1) (Table 17).

Mostrar código 
bd_dis2 %>%
  mutate(PROG_DUMMY = as.character(PROG_DUMMY))%>%
group_by(PROG_DUMMY) %>%
  cascade(
  Correlation=survey_corr(x=Edad, y=num_bienes, vartype = c("ci", "se"), na.rm=TRUE),
  .fill = "Total"
  )%>%
kable(format.args = list(decimal.mark = ',', big.mark = "."), digits=2)%>%
  kable_styling(bootstrap_options = "striped", font_size = 12, full_width = F)%>%
  kable_classic(full_width = F, html_font = "Arial Narrow")%>%
  row_spec(0, color = "white", background = "navy", align = "center", bold= T)
Table 17: Correlación: Número de bienes y edad
PROG_DUMMY Correlation Correlation_low Correlation_upp Correlation_se
0 0,09 0,08 0,10 0,01
1 -0,12 -0,15 -0,10 0,01
Total 0,03 0,02 0,04 0,01

Regresión

Determinantes del ingreso y la educación

Con el ánimo de identificar los determinantes del ingreso como de la educación, se plantean 3 modelos de regresión. Dos de naturaleza gamma y uno lineal para el caso de la educación. En todos los casos se incluyen como covariables el sexo, la edad, el número de bienes, el número de dificultades físicas y la pertenencia a un programa social.

Para este fin y teniendo en cuenta lo observado previamente se plantea un modelo gamma para el ingreso y se plantean 2 modelos para los años de escolaridad uno de naturaleza gamma 1 y otro lineal. En este sentido, se espera que el modelo gamma para la educación sea más robusto que el modelo lineal, dado que la variable dependiente presenta una distribución asimétrica positiva.

En la siguiente tabla se presentan los resultados del modelo para los ingresos ingresos, bajo una especificación exploratoria, mencionada anteriormente (Table 18) 2. Allí se evidencia que el ingreso promedio ajustado por las diferentes covariables es de $123.264 pesos, un valor sustancialmente inferior al descrito anteriormente en la Table 11. De su parte, los coeficientes exponenciados de las covariables 3 indican que las personas de sexo femenino tienen un ingreso un 11% inferior respecto al de los hombres, que por cada año adicional de vida el ingreso crece un 2%, que pertenecer a un programa social se asocia con un ingreso 36% inferior al de las personas no beneficiarias, que por cada aumento en el número de dificultades físicas el ingreso disminuye un 2% y que por cada bien adicional se asocia con un aumento del 10%. Por último, un año de escolaridad adicional se asocia con un aumento del ingreso del 9%.

Mostrar código 
mod <- svyglm(formula = PERCAPITA ~ Sexo + Edad +PROG_DUMMY+num_dif_fisica+nivel_edu4+num_bienes,
                 design = bd_dis3,
                 family = Gamma(link = "log")) #log <-  exp; inverse <- 1/coef

mod %>%
  broom::tidy() %>%
  mutate(exponentiated_estimate = exp(estimate))%>%
  kable(format.args = list(decimal.mark = ',', big.mark = "."), digits=2)%>%
  kable_styling(bootstrap_options = "striped", font_size = 12, full_width = F)%>%
  kable_classic(full_width = F, html_font = "Arial Narrow")%>%
  row_spec(0, color = "white", background = "navy", align = "center", bold= T)
Table 18: Modelo de regresión: Ingreso (Gamma)
term estimate std.error statistic p.value exponentiated_estimate
(Intercept) 11,72 0,06 195,24 0,00 123.264,72
SexoMujer -0,11 0,02 -6,10 0,00 0,89
Edad 0,02 0,00 18,91 0,00 1,02
PROG_DUMMY -0,45 0,02 -24,62 0,00 0,64
num_dif_fisica -0,02 0,01 -2,62 0,01 0,98
nivel_edu4 0,08 0,00 17,61 0,00 1,09
num_bienes 0,10 0,01 18,89 0,00 1,10

De su parte, los modelos de regresión asociados a identificar los determinantes del número de años de educación, señalan que un mayor número de bienes y ser de sexo femenino (respecto a los de sexo masculino), inciden positivamente sobre dicha variable dependiente; mientras que la edad, el número de dificultades físicas y pertenecer a un programa social se asocian negativamente con el número de años de escolaridad. (Table 19).

Mostrar código 
mod1 <- svyglm(formula = nivel_edu4 ~ Sexo + Edad +PROG_DUMMY+num_dif_fisica+PERCAPITA+num_bienes,
                 design = bd_dis4) 
mod1 %>%
  broom::tidy() %>%
  kable(format.args = list(decimal.mark = ',', big.mark = "."), digits=2)%>%
  kable_styling(bootstrap_options = "striped", font_size = 12, full_width = F)%>%
  kable_classic(full_width = F, html_font = "Arial Narrow")%>%
  row_spec(0, color = "white", background = "navy", align = "center", bold= T)
mod2 <- svyglm(formula = nivel_edu4 ~ Sexo + Edad +PROG_DUMMY+num_dif_fisica+PERCAPITA+num_bienes,
                 design = bd_dis4,
                 family = Gamma(link = "log")) #log <-  exp; inverse <- 1/coef

mod2 %>%
  broom::tidy() %>%
  mutate(exponentiated_estimate = exp(estimate))%>%
  kable(format.args = list(decimal.mark = ',', big.mark = "."), digits=2)%>%
  kable_styling(bootstrap_options = "striped", font_size = 12, full_width = F)%>%
  kable_classic(full_width = F, html_font = "Arial Narrow")%>%
  row_spec(0, color = "white", background = "navy", align = "center", bold= T)
Table 19: Modelos de regresión: Educación
(a) Educación (Lineal)
term estimate std.error statistic p.value
(Intercept) 10,73 0,10 103,74 0,00
SexoMujer 0,29 0,05 6,30 0,00
Edad -0,08 0,00 -44,95 0,00
PROG_DUMMY -1,22 0,07 -17,27 0,00
num_dif_fisica -0,05 0,03 -2,09 0,04
PERCAPITA 0,00 0,00 2,29 0,02
num_bienes 0,66 0,03 24,57 0,00
(b) Educación (Gamma)
term estimate std.error statistic p.value exponentiated_estimate
(Intercept) 2,40 0,01 288,60 0 11,01
SexoMujer 0,03 0,00 7,53 0 1,03
Edad -0,01 0,00 -51,90 0 0,99
PROG_DUMMY -0,12 0,01 -22,68 0 0,89
num_dif_fisica -0,01 0,00 -3,09 0 0,99
PERCAPITA 0,00 0,00 15,47 0 1,00
num_bienes 0,05 0,00 57,19 0 1,06

En la siguiente tabla se presentan diferentes métricas de ajuste de los modelos. Considerando que los únicos modelos a comparar tienen que ver con los relacionados a los que emplean a la variable de años de escolaridad como variable dependiente, se observa que desde la óptica del AIC, AICc y BIC el modelo lineal presenta un mejor ajuste que el modelo gamma; no obstante, desde la eprspectiva del RMSE, el modelo gamma presenta un menor error (Table 20).

Mostrar código 
performance::model_performance(mod)|>
  broom::tidy() %>%
  dplyr::select(column, mean)|>
  rename(Value=mean) %>%
  kable(format.args = list(decimal.mark = ',', big.mark = "."), digits=2)%>%
  kable_styling(bootstrap_options = "striped", font_size = 12, full_width = F)%>%
  kable_classic(full_width = F, html_font = "Arial Narrow")%>%
  row_spec(0, color = "white", background = "navy", align = "center", bold= T)
performance::model_performance(mod1)|>
  broom::tidy() %>%
  dplyr::select(column, mean)|>
  rename(Value=mean) %>%
  kable(format.args = list(decimal.mark = ',', big.mark = "."), digits=2)%>%
  kable_styling(bootstrap_options = "striped", font_size = 12, full_width = F)%>%
  kable_classic(full_width = F, html_font = "Arial Narrow")%>%
  row_spec(0, color = "white", background = "navy", align = "center", bold= T)
performance::model_performance(mod2)|>
  broom::tidy() %>%
  dplyr::select(column, mean)|>
  rename(Value=mean) %>%
  kable(format.args = list(decimal.mark = ',', big.mark = "."), digits=2)%>%
  kable_styling(bootstrap_options = "striped", font_size = 12, full_width = F)%>%
  kable_classic(full_width = F, html_font = "Arial Narrow")%>%
  row_spec(0, color = "white", background = "navy", align = "center", bold= T)
Table 20: Métricas de ajuste de los modelos
(a) Ingreso (Gamma)
column Value
AIC 60.730,79
AICc 2.434.306,24
BIC 2.434.380,88
R2 0,40
R2_adjusted 0,39
RMSE 1.919.841,77
Sigma 2,88
(b) Educación (Gamma)
column Value
AIC 409.806,05
AICc 609.338,89
BIC 609.413,02
R2 0,36
R2_adjusted 0,35
RMSE 3,22
Sigma 10,94
(c) Educación (Lineal)
column Value
AIC 7.655,29
AICc 406.308,66
BIC 406.382,79
R2 0,35
R2_adjusted 0,34
RMSE 315,58
Sigma 0,09

Verificación de supuestos

A continuación se presentan diferentes pruebas para verificar los supuestos de los modelos de regresión planteados. En este sentido, se presentan pruebas gráficas y formales para verificar la normalidad de los residuos, la homocedasticidad, la independencia de los residuos y presencia de observaciones influyentes.

1. Distribución de los residuos

Pruebas de normalidad - Modelo lineal (Educación)

A continuación se presenta una prueba gráfica para identificar si los residuos del modelo lineal de la educación siguen una distribución normal. En este caso, la línea azul representa los residuos del modelo y la línea roja la distribución normal. Se observa que si bien hay zonas en las que los residuos se ajustan a una distribución normal, esto no ocurre en toda la distribución (Figure 6).

Mostrar código 
library(svydiags)
library(magrittr)

stdresids <- as.numeric(svystdres(mod1)$stdresids)

bd_dis4$variables %<>%
  mutate(stdresids = stdresids)

# Histograma de los residuales

p1_hist <- ggplot(data = bd_dis4$variables,
                  aes(x = stdresids)) +
  geom_histogram(
    aes(y = ..density..),
    colour = "black",
    bins = 30,
    fill = "blue",
    alpha = 0.3) + 
  geom_density(linewidth = 2, colour = "blue") +
  geom_function(fun = dnorm, colour = "red",
                linewidth = 2) +
  theme_bw() + labs(y = "")

p1_hist
Figure 6: Prueba gráfica: Normalidad de los residuos del modelo lineal (Educación)

Lo anterior es posible verificarlo con pruebas formales como el test de Anderson Darling y el test de Lilliefors. En ambos casos, se observa que los residuos no siguen una distribución normal, lo que sugiere que el modelo no se ajusta adecuadamente a los datos (Table 21).

Mostrar código 
nortest::ad.test(stdresids)|>
  broom::tidy() %>%
  kable(format.args = list(decimal.mark = ',', big.mark = "."), digits=5)%>%
  kable_styling(bootstrap_options = "striped", font_size = 12, full_width = F)%>%
  kable_classic(full_width = F, html_font = "Arial Narrow")%>%
  row_spec(0, color = "white", background = "navy", align = "center", bold= T)
nortest::lillie.test(stdresids)|>
  broom::tidy() %>%
  kable(format.args = list(decimal.mark = ',', big.mark = "."), digits=5)%>%
  kable_styling(bootstrap_options = "striped", font_size = 12, full_width = F)%>%
  kable_classic(full_width = F, html_font = "Arial Narrow")%>%
  row_spec(0, color = "white", background = "navy", align = "center", bold= T)
Table 21: Test de normalidad de los residuos (modelo lineal- Educación)
(a) Anderson Darling
statistic p.value method
56,98471 0 Anderson-Darling normality test
(b) Lilliefors
statistic p.value method
0,01615 0 Lilliefors (Kolmogorov-Smirnov) normality test

Pruebas de distribución - Modelo gamma (Ingreso)

Para revisar si los residuos del modelo del ingreso siguen el patrón esperado para la distribución Gamma, se presenta una prueba gráfica que similar a la anterior busca identificar si los residuos del modelo (línea azul) siguen una distribución gamma (línea verde). En este caso, se observa que los residuos no siguen una distribución gamma, lo que sugiere que el modelo no se ajusta adecuadamente a los datos (Figure 7).

Mostrar código 
library(svydiags)
library(magrittr)
# Extraer residuos
residuos <- as.numeric(svystdres(mod)$stdresids)
residuos2 <- scales::rescale(residuos,  to=c(1e-8,max(residuos)), from = range(residuos, na.rm = TRUE, finite = TRUE))

ParGamma0 <- fitdistrplus::fitdist(residuos2, "gamma", method="mge", optim.method="BFGS")
#ParGamma0$estimate

shape_gamma <- ParGamma0$estimate[1]
rate_gamma <- ParGamma0$estimate[2]

# shape_gamma/rate_gamma
# mean(residuos2)
# var(residuos2)

ParGamma <- MASS::fitdistr(residuos2, densfun = "gamma", start=list(shape=1, scale=1))
#ParGamma$estimate
shape_gamma <- ParGamma$estimate[1]
scale_gamma <- ParGamma$estimate[2]



# Visualizar residuos
ggplot(data = tibble(residuos), aes(x = residuos)) +
  geom_histogram(
    aes(y = ..density..),
    colour = "black",
    bins = 30,
    fill = "royalblue",
    alpha = 0.3
  ) +
  geom_density(linewidth = 1, colour = "royalblue") +
  geom_function(fun = dnorm, colour = "orange",linewidth = 2)+
  stat_function(fun = dgamma, args = list(shape = shape_gamma, scale = scale_gamma), color = "limegreen", size = 1)+
  labs(title = "", x = "Residuos", y = "Densidad")
Figure 7: Prueba gráfica: Distribución Gamma de los residuos del modelo del ingreso

2.a. Homocedasticidad (modelo lineal)

Respecto al supuesto de varianza constante u homocedasticidad, se presenta una prueba gráfica que busca identificar si los residuos del modelo siguen un patrón de varianza constante según cada una de las covariables (Figure 8). En este caso, se observa que los residuos no siguen un patrón de varianza constante, lo que sugiere que los intervalos de confianza y las pruebas de hipótesis pueden ser incorrectos, por lo que no es apropiado generar conclusiones acerca de la significancia de los coeficientes.

Mostrar código 
library(patchwork)

bd_dis4$variables %<>%
  mutate(pred = predict(mod1))

#summary(mod1)

g2 <- ggplot(data = bd_dis4$variables,
             aes(x = Sexo, y = stdresids))+
  geom_point() +
  geom_hline(yintercept = 0) + theme_bw()

g3 <- ggplot(data = bd_dis4$variables,
             aes(x = Edad , y = stdresids))+
  geom_point() +
  geom_hline(yintercept = 0) + 
  theme_bw()

g4 <- ggplot(data = bd_dis4$variables,
             aes(x = PROG_DUMMY, y = stdresids))+
  geom_point() +
  geom_hline(yintercept = 0) + 
  theme_bw()

g5 <- ggplot(data = bd_dis4$variables,
             aes(x = num_dif_fisica, y = stdresids))+
  geom_point() + geom_hline(yintercept = 0) +
  theme_bw()

g6 <- ggplot(data = bd_dis4$variables,
             aes(x = PERCAPITA, y = stdresids))+
  geom_point() + geom_hline(yintercept = 0) +
  theme_bw()

g7 <- ggplot(data = bd_dis4$variables,
             aes(x = num_bienes, y = stdresids))+
  geom_point() + geom_hline(yintercept = 0) +
  theme_bw()

(g2|g3)/(g4|g5)/(g6|g7)
Figure 8: Prueba gráfica: Homocedasticidad de los residuos del modelo lineal (Educación)

2.b. Relación entre la media y la varianza (modelo gamma)

Teniendo en cuenta que la media de un modelo gamma corresponde a \[\alpha \times \beta\] y la varianza \[\alpha \times \beta^2\], se espera que la varianza de los datos debería crecer exponencialmente respecto de la media. De tal manera la siguiente figura presenta los valores estimados de la variable dependiente y su error, en donde se espera que a mayor valor estimado, se evidencien mayores residuos (Figure 9). Sin embargo, esto no ocurre por lo que se argumenta que este supuesto no se cumple en el modelo del ingreso.

Mostrar código 
predicciones <- predict(mod, type = "response")

ggplot(data = tibble(predicciones, residuos), aes(x = predicciones, y = residuos)) +
  geom_point() +
  labs(title = "Residuos vs Predicciones", x = "Predicciones", y = "Residuos")
Figure 9: Prueba gráfica: Relación entre la media y la varianza del modelo gamma (Ingreso)

3. Multicolinealidad (VIF)

La especificación para el modelo del ingreso no presenta problemas de multicolinealidad, como lo ilustran factores de inflación de varianza (VIF) menores a 10, sin embargo, el relacionado con la educación presenta problemas de multicolinealidad con la variable PERCAPITA, como se ilustra al exihibir un VIF superior a 10 (Table 22). Especificaciones futuras debe considerar la transformación o supresión de dicha variable.

Mostrar código 
performance::check_collinearity(mod)|>
  data.frame()|>
  #broom::tidy() %>%
  kable(format.args = list(decimal.mark = ',', big.mark = "."), digits=2)%>%
  kable_styling(bootstrap_options = "striped", font_size = 12, full_width = F)%>%
  kable_classic(full_width = F, html_font = "Arial Narrow")%>%
  row_spec(0, color = "white", background = "navy", align = "center", bold= T)
performance::check_collinearity(mod1)|>
  data.frame()|>
  #broom::tidy() %>%
  kable(format.args = list(decimal.mark = ',', big.mark = "."), digits=2)%>%
  kable_styling(bootstrap_options = "striped", font_size = 12, full_width = F)%>%
  kable_classic(full_width = F, html_font = "Arial Narrow")%>%
  row_spec(0, color = "white", background = "navy", align = "center", bold= T)
Table 22: Multicolinealidad (VIF)
Term VIF VIF_CI_low VIF_CI_high SE_factor Tolerance Tolerance_CI_low Tolerance_CI_high
Sexo 1,29 1,28 1,30 1,14 0,78 0,77 0,78
Edad 2,95 2,91 2,98 1,72 0,34 0,34 0,34
PROG_DUMMY 1,11 1,10 1,11 1,05 0,90 0,90 0,91
num_dif_fisica 1,49 1,47 1,50 1,22 0,67 0,67 0,68
nivel_edu4 2,93 2,90 2,97 1,71 0,34 0,34 0,34
num_bienes 2,20 2,17 2,22 1,48 0,46 0,45 0,46
Term VIF VIF_CI_low VIF_CI_high SE_factor Tolerance Tolerance_CI_low Tolerance_CI_high
Sexo 1,02 1,01 1,03 1,01 0,98 0,97 0,99
Edad 1,64 1,62 1,65 1,28 0,61 0,60 0,62
PROG_DUMMY 2,26 2,24 2,29 1,50 0,44 0,44 0,45
num_dif_fisica 1,28 1,27 1,30 1,13 0,78 0,77 0,78
PERCAPITA 11,12 10,97 11,27 3,33 0,09 0,09 0,09
num_bienes 8,92 8,80 9,04 2,99 0,11 0,11 0,11

4. Observaciones influyentes

Frente a la presencia de observaciones influyentes -puntos de datos que tienen un impacto sustancial en los coeficientes del modelo-, se presentan diferentes métricas para identificarlas. En este caso, se presentan las distancias de Cook, los DFBETAS, los DFFITS y la presencia de puntos de apalancamiento.

4.a. Distancias de Cook

Las distancias de Cook permiten evidenciar la presencia de múltiples observaciones influyentes en el modelo del ingreso, que bien pueden deberse a problemas de medición de la variable, una especificación insuficiente o que sencillamente el modelo gamma no representa adecuadamente la estructura de la relación de los datos para todas las observaciones y una forma funcional diferente es requerida (Figure 10).

Mostrar código 
d_cook = data.frame(
  cook = svyCooksD(mod),
  id = 1:length(svyCooksD(mod)))

ggplot(d_cook, aes(y = cook, x = id)) +
  geom_point() +
  theme_bw(20)
Figure 10: Observaciones influyentes: Distancias de Cook

4.b. Detección de observaciones influyentes (Df Betas(𝑖)𝑗)

Análogamente, una aproximación con la técnica de DfBetas que mide cuánto cambia un coeficiente específico cuando se elimina una observación, permite identificar observaciones influyentes. Consistentemente con lo evidenciado en la anterior prueba, se tiene evidencia de la presencia de observaciones influyentes (Figure 11). Por lo anterior se argumenta que el modelo planteado para el ingreso es frágil y no modela adecuadamente la relación entre la variable dependiente y las covariables en todas sus observaciones.

Mostrar código 
d_dfbetas = data.frame(t(svydfbetas(mod)$Dfbetas))
colnames(d_dfbetas) <- paste0("Beta_", 1:7)
#d_dfbetas %>% slice(1:5L)

d_dfbetas$id <- 1:nrow(d_dfbetas)

d_dfbetas <- reshape2::melt(d_dfbetas,
                            id.vars = "id")
cutoff <- svydfbetas(mod)$cutoff

d_dfbetas %<>%
  mutate(
    Criterio = ifelse(
      abs(value) > cutoff, "Si", "No"))

tex_label <- d_dfbetas %>%
  filter(Criterio == "Si") %>%
  arrange(desc(abs(value))) %>%
  slice(1:10L)

ggplot(d_dfbetas, aes(y = abs(value), x = id)) +
  geom_point(aes(col = Criterio)) +
  geom_text(data = tex_label,
            angle = 45,
            vjust = -1,
            aes(label = id)) +
  geom_hline(aes(yintercept = cutoff)) +
  facet_wrap(. ~ variable, nrow = 2) +
  scale_color_manual(
    values = c("Si" = "red", "No" = "black")) +
  theme_bw()
Figure 11: Observaciones influyentes: Df Betas(𝑖)𝑗

4.c. Detección de observaciones influyentes (𝐷𝑓 𝐹𝑖𝑡𝑠(𝑖))

Resultados similares se presen tan al emplear la técnica de DFFITS o diferencia en los valores estimados al excluir un valor y estimar el modelo, consistentemente con las gráficas anteriores, evidencia la presencia de múltiples observaciones influeyentes en el modelo (Figure 12).

Mostrar código 
d_dffits = data.frame(
  dffits = svydffits(mod)$Dffits,
  id = 1:length(svydffits(mod)$Dffits))

cutoff <- svydffits(mod)$cutoff

d_dffits %<>% mutate(
  C_cutoff = ifelse(abs(dffits) > cutoff, "Si", "No"))

ggplot(d_dffits, aes(y = abs(dffits), x = id)) +
  geom_point(aes(col = C_cutoff)) +
  geom_hline(yintercept = cutoff) +
  scale_color_manual(
    values = c("Si" = "red", "No" = "black"))+
  theme_bw()
Figure 12: Observaciones influyentes: 𝐷𝑓 𝐹𝑖𝑡𝑠(𝑖)

4.d. Detección de observaciones influyentes: Apalancamiento (ℎ𝑖𝑖)

En esta misma línea, la detección de puntos de apalancamiento resalta la presencia de observaciones influyentes en el modelo del ingreso, tal como se ilustra a continuación (Figure 13).

Mostrar código 
vec_hat <- svyhat(mod, doplot = FALSE)

d_hat = data.frame(hat = vec_hat,
                   id = 1:length(vec_hat))
d_hat %<>% mutate(
  C_cutoff = ifelse(hat > (3 * mean(hat)),
                    "Si", "No"))

ggplot(d_hat, aes(y = hat, x = id)) +
  geom_point(aes(col = C_cutoff)) +
  geom_hline(yintercept = (3 * mean(d_hat$hat))) +
  scale_color_manual(
    values = c("Si" = "red", "No" = "black"))+
  theme_bw()
Figure 13: Observaciones influyentes: Apalancamiento

Inferencia

Considerando que los modelos previamente descritos no cumplen con uno o más supuestos, se argumenta que hacer inferencia sobre la significancia de los paramétros o del modelo en su conjunto sería inapropiado. Por este motivo no se realiza una lectura inferencial de los resultados obtenidos.

Conclusiones

Las personas beneficiarias de los programas sociales en Colombia a corte de 2023, tienden a ser de sexo femenino, más jóvenes que los no beneficiarios, con bajos niveles de escolaridad, un mayor número de dificultades físicas, menores ingresos y un mayor hacinamiento. Estos hallazgos por un lado implican que la focalización de los programas sociales está centrada en personas con mayor vulnerabilidad socioeconómica.

Complementariamente, asociados a los modelos de regresión, estrategias orientadas a mejorar el ingreso y la escolaridad deben considerar centrarse en personas de sexo femenino y en intervenciones dirigidas a personas con un mayor número de dificultades físicas. Sin embargo, no es posible hacer inferencia respecto a dichos hallazgos, considerando que los modelos planteados no cumplen con los supuestos requeridos para este propósito.

Referencias

DANE. (2024). Encuesta Nacional de Calidad de Vida 2023. https://www.dane.gov.co/index.php/estadisticas-por-tema/salud/calidad-de-vida-ecv/encuesta-nacional-de-calidad-de-vida-ecv-2023
Gutiérrez, A., Téllez, C., & Guerrero, S. (2024). Análisis de encuestas de hogares con R.

Footnotes

  1. Se empleó como función de enlace el logaritmo, por lo que los coeficientes se exponenciaran para facilitar su interpretación↩︎

  2. Si bien a continuación se hace una lectura de los coeficientes, esta no debe tomarse como definitiva, puesto que aún no se ha validado si los supuestos subyacentes a este modelo se cumplen, sin embargo se hace por efectos ilustrativos.↩︎

  3. Que se pueden leer como un ratio del ingreso.↩︎