A continuación se presenta todo el código R usado para este trabajo.
Introducción y definición del usuario.
El presente trabajo se realiza basado en la encuesta de calidad de vida publicada por el DANE en el año 2019. Esta encuesta se encuentra dirigida a hogares y a partir de la información obtenida de ellos se calculan indicadores para la medición de varios aspectos económicos y sociales. Además, facilitan el conocimiento y explicación de los determinantes o factores causales del comportamiento de dichos aspectos, lo cual es de gran importancia para el diseño, monitoreo y medición de resultados de las políticas públicas. La Encuesta de Calidad de Vida (ECV) es una investigación que el DANE realiza con el objeto de recoger información sobre diferentes aspectos y dimensiones del bienestar y las condiciones de vida de los hogares, incluyendo temas como: el acceso a bienes y servicios públicos, privados o comunales, salud, educación, atención integral de niños y niñas menores de 5 años, entre otros. (DIMPE,2019).
Para darle solución al problema se define al jefe del hogar como la fuente principal de información con respecto a las características que tiene, valga la redundancia, el hogar. Así el número de hijos hará referencia a los hijos que tenga dicho jefe. Las variables explicativas a usar serán escogidas minuciosamente en términos de correlaciones con la variable respuesta de tal forma que permita mejorar la efectividad del modelo sin caer en sobre ajuste. Finalmente se creará una aplicación web con el modelo predictivo que tendrá como usuario final a cualquier persona del territorio nacional colombiano que desee identificar el número de hijos con base en ciertas variables que describen el hogar; además, es importante aclarar que el usuario previamente debe leer este reporte técnico. Un ejemplo de un usuario podría ser un trabajador comunitario que quiera establecer cuántos hijos hay en ciertos hogares donde realiza su trabajo social, así esta persona puede ingresar a la aplicación y con un número reducido de variables, definir el total de hijos por hogar y con esto establecer estrategias efectivas.
De acuerdo con la estructura de la base de datos, en una vivienda puede haber varios hogares, por tanto, es preciso distinguir con un identificador a los hogares que habitan una vivienda. Para ello se procedió así:
Dado que cada vivienda se identifica con la variable Directorio, y cada hogar de dicha vivienda se identifica con la variable Secuencia_p, se procedió a hacer las respectivas agrupaciones, de tal manera que un hogar sería identificado por el Directorio de la vivienda al cual pertenece, un guión y después la Secuencia_p del respectivo hogar.
II.b. Definición de la variable respuesta. Para la creación de la variable respuesta se define la cantidad de hijos del hogar como el número de hijos que tiene el jefe del hogar.
Debido a que la variable respuesta no estaba explícita en los datos, se procedió al cálculo de ésta de la siguiente manera: se tomó la variable P6051, que corresponde al parentesco del encuestado con el jefe del hogar, y a partir de esto se identificaron 2 posibles casos:
En el hogar hay uno (1) o más hijos: este caso ocurre cuando entre las respuestas a la pregunta P6051 para cada hogar hay una o más ocurrencias de la respuesta número 3, la cual corresponde a que el parentesco del encuestado con el jefe de hogar es hijo(a) o hijastro(a). Por lo tanto, el número de hijos de ese hogar es igual al número de ocurrencias de la respuesta número 3.
En el hogar no hay hijos: ocurre cuando entre las respuestas a la pregunta P6051 no está la número 3.
Basado en lo anteriormente expuesto la variable respuesta corresponderá a: num hijos= num hijos del jefe del hogar
III. Descripción de las variables explicativas.III.a. Definición del procedimiento usado. La selección de las variables explicativas se hizo por medio de un proceso de varios pasos así:
Se revisó cuáles variables tenían definición en la documentación suministrada por el DANE, con el fin de descartar aquellas que carecían de este significado.
Al hacer una revisión de cada variable dentro de las diferentes tablas proporcionadas por el DANE en la encuesta de calidad de vida, se encontró que algunas variables presentes eran creación del DANE formuladas a partir de ciertos indicadores preestablecidos y en otros casos no se encontró descripción alguna que permitiera identificar la pregunta en cuestión. Basado en la descripción del usuario final que hará uso de la aplicación se decide que las variables explicativas usadas deben tener sentido y permitir una respuesta fácil por parte del jefe del hogar, a su vez que estas tengan interpretación y tengan sentido para explicar el número de hijos del hogar.
Se hace imputación de datos, pues para los modelos que se propondrán no pueden existir valores nulos.
Esta imputación de datos se hace creando una nueva categoría ‘000’ en variables categóricas y usando la media de los datos para el caso donde las variables son de tipo numérico.
Se calcularon las correlaciones de las variables predictoras con la variable de salida, con el fin de encontrar una relación significativa entre estas.
Se hace un barrido inicial de la base de datos y se calculan dos correlaciones: Pearson y Spearman. Con el fin de encontrar dependencia lineal o no lineal. Se define un conjunto de cincuenta variables iniciales que cumplen con dos requisitos: Correlación significativa e interpretables. Se propone un primer modelo lineal con resultados no significativos, debido a que, aunque las correlaciones son diferentes de cero, estas no son suficientes para aumentar un factor de determinación R2 altamente significativo y problemas graves en términos de residuales; es por ello que no se define este como un modelo final.
A pesar de que todas las variables resultan significativas, el resultado solo explica 15% de la variabilidad de los hijos del hogar. Se hace un análisis de residuales y no se encuentra validez en los supuestos de normalidad en los errores y varianza constante. Es por ello que se descarta esta alternativa.
Con base en los resultados del punto anterior, se realizó una discusión entre los miembros del equipo para determinar, según el contexto del problema, cuáles variables eran mejores candidatas para la predicción del número de hijos.
Debido a los deficientes resultados del primer modelo lineal realizado, se busca hacer una reducción de dimensionalidad poniendo en contexto con el equipo aquellas variables con correlaciones altas y significativas en el primer modelo para, en conjunto, leer detenidamente cuales tenían sentido para explicar el número de hijos.
Posteriormente, se ajustó un modelo con la metodología XGBOOST buscando hacer una clasificación para múltiples clases, con el fin de encontrar las variables más significativas de las ya preseleccionadas en los pasos previos. Esto nos dio como resultado 40 variables significativas, y a partir de estas se realizó una nueva discusión, en la cual finalmente se escogieron las 18 variables predictoras definitivas.
#datostae <- read.csv2("C:/Users/HP ENVY/Desktop/prediccion hijos tae/Resultado.csv")
#porcentaje de NAs en cada variable
#colMeans(is.na(datostae)) #p8586 es la ultima sin NA
#Eliminando variables que tengan más del 1% de NA
#columnasNA <- which(colMeans(is.na(datostae))>0.01)
#length(columnasNA)
#datostae[,columnasNA] <- NULL
#Eliminando observaciones con NA
#datostae <- na.omit(datostae)
#Eliminando variables que tengan más del 5% de valores faltantes
#columnasfal <- which(colMeans(datostae=="")>0.05)
#length(columnasfal)
#datostae[,columnasfal] <- NULL
#Eliminando observaciones con datos faltantes
#obs_fal <- which(rowMeans(datostae=="")>0)
#length(obs_fal)
#datostae <- datostae[-obs_fal,]
#Eliminando variables que tengan más del 1% de valores guion
#columnasguion <- which(colMeans(datostae=="-")>0.01)
#length(columnasguion)
#datostae[,columnasguion] <- NULL
#Eliminando observaciones con datos guion
#obs_guion <- which(rowMeans(datostae=="-")>0)
#length(obs_guion)
#datostae <- datostae[-obs_guion,]
#for (i in 3:length(datostae)){
# datostae[,i] <- as.numeric(datostae[,i])
#}
#d <- as.matrix(datostae[,3:121])
#sort(cor(d)[1,])
#which(na.omit(cor(d)[1,]),>0.15)
#seleccionadas1 <- data.frame(Hijos=datostae$Hijos,
# P6040=datostae$P6040,
# P6083=as.factor(datostae$P6083),
# P6081=as.factor(datostae$P6081),
# P5502=as.factor(datostae$P5502),
# P6080=as.factor(datostae$P6080))
#P784S1=as.factor(datostae$P784S1),
#P3007S2=as.factor(datostae$P3007S2),
#P3007S1=as.factor(datostae$P3007S1),
#P3007S3=as.factor(datostae$P3007S3),
#P1079=as.factor(datostae$P1079),
#P784S2=as.factor(datostae$P784S2),
#P3007S4=as.factor(datostae$P3007S4),
#P1072S1=as.factor(datostae$P1072S1))
#m1 <- lm(Hijos~.,data = seleccionadas1)
#knitr::kable(summary(m1)$coefficients)
#s <- summary(m1)
#knitr::kable(data.frame(Ajuste="R2 Ajustado",Resultado=s$adj.r.squared))III.b. Presentación de las variables escogidas.
Es importante aclarar que para las variables predictoras escogidas las respuestas tomadas fueron las del jefe del hogar.
Luego del exhaustivo proceso de selección de variables realizado, se escogió el siguiente conjunto de variables:
P6390S2: ¿A qué actividad se dedica principalmente la empresa o negocio en la que … realiza su trabajo?
La respuesta a esta pregunta es un código que se debe de validar con el DANE, por medio del código CIIU que se puede verificar en el siguiente link de la cámara de comercio de bogotá donde se puede buscar la actividad económica. https://linea.ccb.org.co/descripcionciiu/
P1_DEPARTAMENTO: Departamento de nacimiento.
La respuesta a esta pregunta es un código que se debe de validar con el DANE por medio del siguiente enlace: https://fopep.gov.co/wp-content/uploads/2019/02/Tabla-C%C3%B3digos-Dane.pdf
P8587S1: Grado o año aprobado.
Los grados posibles de cursar en cada nivel educativo corresponden con los grados que conforman el nivel. Cinco (5) representa Básica primaria, nueve (9) representa básica secundaria y once (11) o doce (12) representa educación media.
P415: ¿Cuántas horas a la semana trabaja normalmente ….. en ese trabajo ?
Se utilizan para la medición de horas de trabajo, información que hace posible caracterizar las condiciones de empleo diferenciando entre empleo de tiempo completo y de tiempo parcial. Para diligenciar esta pregunta se indaga por el número de horas que la persona trabaja normalmente a la semana. Se aclara que no se debe tener en cuenta las horas que trabajó de más o de menos por cualquier motivo en la semana, sino aquellas horas que trabaja normalmente en su empleo principal.
P6750: ¿cuál fue la ganancia neta o los honorarios netos de… En esa actividad, negocio, profesión o finca, el mes pasado ?
En esta pregunta su respuesta es la Ganancia neta Valor ($ COP) un número positivo.
P756S2: Municipio.
La respuesta a esta pregunta es un código que se debe de validar con el DANE.
P8587: ¿cuál es el nivel educativo más alto alcanzado por … Y el último año o grado aprobado en este nivel?
1 Ninguno. 2 Preescolar. 3 Básica Primaria (1º - 5º). 4 Básica secundaria (6º - 9º). 5 Media (10º - 13º). 6 Técnico sin título . 7 Técnico con título . 8 Tecnológico sin título. 9 Tecnológico con título. 10 Universitario sin título. 11 Universitario con título. 12 Postgrado sin título. 13 Postgrado con título.
P8520S1A1: Estrato para tarifa.
¿Con cuáles de los siguientes servicios públicos, privados o comunales cuenta la vivienda? 1. Energía eléctrica Estrato para tarifa: 1 Bajo - Bajo 2 Bajo 3 Medio - Bajo 4 Medio 5 Medio - Alto 6 Alto 8 Planta eléctrica 9 No conoce el estrato o no cuenta con recibo de pago. 0 Recibos sin estrato o el servicio es pirata.
P1896: En general, qué tan satisfecho(a) se siente … con su ingreso actualmente?
En general, qué tan satisfecho(a) se siente … con su ingreso actualmente? 10 Totalmente satisfecho(a) 9 8 7 6 5 4 3 2 1 0 Totalmente insatisfecho(a).
P1897: En general, qué tan satisfecho(a) se siente … con su salud actualmente?
En general, qué tan satisfecho(a) se siente … con su salud actualmente? 10 Totalmente satisfecho(a) 9 8 7 6 5 4 3 2 1 0 Totalmente insatisfecho(a).
P1899: En general, qué tan satisfecho(a) se siente … con su trabajo/actividad actualmente?
En general, qué tan satisfecho(a) se siente … con su trabajo/actividad actualmente? 10 Totalmente satisfecho(a) 9 8 7 6 5 4 3 2 1 0 Totalmente insatisfecho(a).
P1927: ¿En cuál escalón diría usted que se encuentra parado(a) en este momento?
Imagine una escalera con escalones numerados de 0 a 10, donde 0 es el escalón más bajo y 10 el escalón más alto. El más alto representa la mejor vida que usted podría tener y el más bajo, la peor.
P4015: Material predominante de los pisos.
Material predominante de los pisos 1. Alfombra o tapete de pared a pared 2. Madera pulida y lacada, parqué 3. Mármol 4. Baldosa, vinilo, tableta, ladrillo, laminado 5. Madera burda, tabla, tablón, otro vegetal 6. Cemento, gravilla 7. Tierra, arena.
P4567: ¿Cuál es el material predominante del techo o cubierta?
P5502: Actualmente…
Actualmente…: 1. No está casado(a) y vive en pareja hace menos de dos años 2. No está casado(a) y vive en pareja hace dos años o más 3. Está viudo(a) 4. Está separado(a) o divorciado(a) 5. Está soltero(a) 6. Está casado(a).
P6087: ¿Cuál es o fue el nivel de educación más alto alcanzado por el padre de…..?
Cuál es o fue el nivel de educación más alto alcanzado por el padre de …? 1. Algunos años de primaria 2. Toda la primaria 3. Algunos años de secundaria 4. Toda la secundaria 5. Uno o más años de técnica o tecnológica 6. Técnica o tecnológica completa 7. Uno o más años de universidad 8. Universitaria completa 9. Ninguno 10. No sabe.
P6088: ¿cuál es o fue el nivel de educación más alto alcanzado por la madre de……?
Cuál es o fue el nivel de educación más alto alcanzado por la madre de …? 1. Algunos años de primaria 2. Toda la primaria 3. Algunos años de secundaria 4. Toda la secundaria 5. Uno o más años de técnica o tecnológica 6. Técnica o tecnológica completa 7. Uno o más años de universidad 8. Universitaria completa 9. Ninguno 10. No sabe.
P6040: ¿cuántos años cumplidos tiene…?
Valor numérico de la edad.
#library(xgboost)
#library(shiny)
#library(shinyWidgets)
#library(shinythemes)
#library(V8)
#datos <- read.csv('Resultado1.csv',sep = ';', encoding = 'UTF-8')
#modelopython<-xgb.load('modelo1.bin')
# Define UI for application that draws a histogram
#ui <-(
# navbarPage(
# title="TAE",
# theme=shinytheme("spacelab"),
# inverse=TRUE,
# position = "static-top",
# tabPanel("Clasificacion Hijos",
# fluidPage(
# # titlePanel("------------------------------------------------------------------------------- "),
# flowLayout(
# numericInput("S1",
# label = h4("Actividad Economica"),
# value=1, min = min(datos$ACTIVIDAD_ECONOMICA),
# max=max(datos$ACTIVIDAD_ECONOMICA)),
# numericInput("S2",
# label = h4("Grado Aprobado"),
# value=1, min = min(datos$GRADO_APROBADO),
# max=max(datos$GRADO_APROBADO)),
# numericInput("S3",
# label = h4("Horas Trabaja Semana"),
# value=1, min = min(datos$HORAS_SEM_TRABAJA), max=max(datos$HORAS_SEM_TRABAJA)),
# numericInput("S4",
# label = h4("Ganancia Neta"),
# value=1, min = min(datos$GANANCIA_NETA),
# max=max(datos$GANANCIA_NETA)),
# numericInput("S5",
# label = h4("Municipio Nacimiento"),
# value = 1, min = min(datos$MUNICIPIO_NACIMIENTO),
# max=max(datos$MUNICIPIO_NACIMIENTO)),
# numericInput("S6",
# label = h4("Nivel Educativo"),
# value=1, min = min(datos$NIVEL_EDUCATIVO),
# max=max(datos$NIVEL_EDUCATIVO)),
# numericInput("S7",
# label = h4("Estrato Tarifa"),
# value=1, min = min(datos$ESTRATO_TARIFA),
# max=max(datos$ESTRATO_TARIFA)),
# numericInput("S8",
# label = h4("Departamento"),
# value=1, min = min(datos$DEPARTAMENTO),
# max=max(datos$DEPARTAMENTO)),
# numericInput("S9",
# label = h4("Satisfaccion Ingreso"),
# value=1, min = min(datos$SACTISFACCION_INGRESO),
# max=max(datos$SACTISFACCION_INGRESO)),
# numericInput("S10",
# label = h4("Satisfaccion Salud"),
# value=1, min = min(datos$SACTISFACCION_SALUD),
# max=max(datos$SACTISFACCION_SALUD)),
# numericInput("S11",
# label = h4("Satisfaccion Trabajo"),
# value=1, min = min(datos$SACTISFACCION_TRABAJO),
# max=max(datos$SACTISFACCION_TRABAJO)),
# numericInput("S12",
# label = h4("Escalon Vida"),
# value=1, min = min(datos$ESCALON_VIDA),
# max=max(datos$ESCALON_VIDA)),
# numericInput("S13",
# label = h4("Material Pisos"),
# value=1, min = min(datos$MATERIAL_PISOS),
# max=max(datos$MATERIAL_PISOS)),
# numericInput("S14",
# label = h4("Material Techo"),
# value=1, min = min(datos$MATERIAL_TECHO),
# max=max(datos$MATERIAL_TECHO)),
# numericInput("S15",
# label = h4("Estado Civil"),
# value=1, min = min(datos$ESTADO_CIVIL),
# max=max(datos$ESTADO_CIVIL)),
# numericInput("S16",
# label = h4("Nivel Educativo Padre"),
# value=1, min = min(datos$NIVEL_EDUCACION_PADRE),
# max=max(datos$NIVEL_EDUCACION_PADRE)),
# numericInput("S17",
# label = h4("Nivel Educativo Madre"),
# value=1, min = min(datos$NIVEL_EDUCACION_MADRE),
# max=max(datos$NIVEL_EDUCACION_MADRE)),
# numericInput("S18",
# label = h4("Edad"),
# value=1, min = min(datos$EDAD),
# max=max(datos$EDAD)),
# verbatimTextOutput("text",placeholder = F),
# verbatimTextOutput(""),
# verbatimTextOutput(""),
# verbatimTextOutput(""),
# verbatimTextOutput(""),
# verbatimTextOutput(""),
# img(src="1.png", height="150", width="250")
# )
#)
#)))
#server <- function(input, output,session) {
# output$text <- renderText({
# {
# Data1 = reactive({
# df<-data.frame(
# ACTIVIDAD_ECONOMICA = input$S1,
# GRADO_APROBADO = input$S2,
# HORAS_SEM_TRABAJA = input$S3,
# GANANCIA_NETA = input$S4,
# MUNICIPIO_NACIMIENTO = input$S5,
# NIVEL_EDUCATIVO = input$S6,
# ESTRATO_TARIFA = input$S7,
# DEPARTAMENTO = input$S8,
# SACTISFACCION_INGRESO = input$S9,
# SACTISFACCION_SALUD = input$S10,
# SACTISFACCION_TRABAJO = input$S11,
# ESCALON_VIDA = input$S12,
# MATERIAL_PISOS = input$S13,
# MATERIAL_TECHO = input$S14,
# ESTADO_CIVIL = input$S15,
# NIVEL_EDUCACION_PADRE = input$S16,
# NIVEL_EDUCACION_MADRE = input$S17,
# EDAD = input$S18)
# })
# prediccion <- predict(modelopython, newdata = as.matrix(Data1()))
# #print(prediccion)
# paste(c("Hijos Estimados:" ,prediccion))
# }
# })
#}
# Run the application
#shinyApp(ui = ui, server = server) Definición del Modelo.V.a. Modelo escogido. Luego de realizar algunas pruebas pilotos con enfoques de regresión lineal, hacer sus comparaciones de métricas y verificaciones de supuestos se decide implementar un modelo Gradient Boosting el cual está formado por un conjunto de árboles de decisión individuales, entrenados de forma secuencial, de forma que cada nuevo árbol trata de mejorar los errores de los árboles anteriores. La predicción de una nueva observación se obtiene agregando las predicciones de todos los árboles individuales que forman el modelo.
Un modelo Gradient Boosting Trees está formado por un conjunto (ensemble) de árboles de decisión individuales, entrenados de forma secuencial. Cada nuevo árbol emplea información del árbol anterior para aprender de sus errores, mejorando iteración a iteración. En cada árbol individual, las observaciones se van distribuyendo por bifurcaciones (nodos) generando la estructura del árbol hasta alcanzar un nodo terminal. La predicción de una nueva observación se obtiene agregando las predicciones de todos los árboles individuales que forman el modelo. Para entender cómo funcionan los modelos Gradient Boosting Trees es necesario conocer primero los conceptos de ensemble y boosting.
El algoritmo es implementado en python debido a que tiene una función de regularización mejor que la de R, es decir que a la hora de determinar la probabilidad de escoger una clasificación u otra este presenta mejor desempeño.
# Importando Librerias
#import pandas as pd
#import numpy as np
#from scipy import stats
#import seaborn as sns
#import matplotlib.pyplot as plt
#import sweetviz as sv
#import treeplot
#from numpy import array
#from numpy import sort
#from sklearn.impute import SimpleImputer
#from sklearn.preprocessing import StandardScaler, PolynomialFeatures,LabelEncoder
#from sklearn.model_selection import train_test_split,GridSearchCV,RandomizedSearchCV
#from sklearn.utils import resample
#from sklearn.feature_selection import SelectFromModel,SelectKBest,f_classif
#from sklearn.metrics import confusion_matrix
#from sklearn.linear_model import LogisticRegression
#from sklearn.ensemble import RandomForestClassifier
#from xgboost import XGBClassifier
#from sklearn.metrics import accuracy_score, #precision_score,recall_score,classification_report,roc_curve,roc_auc_score
#import pandas as pd
#import xgboost as xgb
#from xgboost import plot_tree
#from urllib.request import urlretrieve
#datos = pd.read_csv('Resultado1.csv', sep=';',encoding='UTF-8')
#datos.info()
#datos_train, datos_test = train_test_split(datos, test_size=.20, random_state=1)
#datos_train_mat = xgb.DMatrix(datos_train.drop("HIJOS", 1), label=datos_train["HIJOS"])
#datos_test_mat = xgb.DMatrix(datos_test.drop("HIJOS", 1), label=datos_test["HIJOS"])
#rondas =1600
#evaluacion = [(datos_test_mat, "eval"), (datos_train_mat, "train")]
#modelo = xgb.train(parametros, datos_test_mat, rondas, evaluacion)
#prediccion = modelo.predict(datos_test_mat)
#print(classification_report(datos_test["HIJOS"], prediccion))
#cm = confusion_matrix(datos_test["HIJOS"], prediccion)
#cm
#ax = treeplot.xgboost(modelo)
#modelo.save_model('modelo1.bin')Una de las ventajas de este algoritmo es que permite la entrada de datos nulos, adicional a esto permite la exportación en formato binario para ser usado por ejemplo en R. Se define como el número de nodos siete debido a que si es menor a esto la varianza se hace pequeña pero el sesgo crece, es por eso que se define un punto de equilibrio de sesgo y varianza y se definió a partir de ensayo y error, para llegar a este modelo final se hicieron varios con el que se probaba efectividad del ajuste. Se define el learning rate a un valor de 0.01 este es una media que ayuda a que el modelo no se deje influenciar por el paso anterior (o árbol anterior) es por ello que se garantiza que el modelo aprenda verificando siempre que no se cayera en un problema de sobre parametrización.