rm (list = ls())

Pregunta 1

library(data.table)
library(caret)
## Loading required package: lattice
## Loading required package: ggplot2
library(ggplot2)
library(knitr)
library(moderndive)
library(rpart.plot)
## Loading required package: rpart
library(rpart)
library(readr)
datos <- read_delim("datos_tarea_5.csv", 
    ";", escape_double = FALSE, na = "NA", 
    trim_ws = TRUE)
## 
## -- Column specification --------------------------------------------------------
## cols(
##   .default = col_double(),
##   a6_otro = col_character(),
##   b16_otro = col_character(),
##   e3_7 = col_logical(),
##   e3_8 = col_logical(),
##   e3_11 = col_logical(),
##   e3_12 = col_logical(),
##   e19_otro = col_character(),
##   d1_monto = col_number(),
##   ss_t = col_number(),
##   d2_1_monto = col_number(),
##   d2_1_porcentaje = col_logical(),
##   d2_2_monto = col_number(),
##   d2_2_porcentaje = col_logical(),
##   d2_3_monto = col_number(),
##   d2_3_porcentaje = col_logical(),
##   d2_3_opcionb = col_logical(),
##   d2_4_monto = col_number(),
##   d2_4_porcentaje = col_logical(),
##   d2_5_monto = col_number(),
##   d2_5_porcentaje = col_logical()
##   # ... with 98 more columns
## )
## i Use `spec()` for the full column specifications.

P2

Muestre un histograma de ingreso por region. Limite el histograma a las observaciones con ingresos (ss_t) mayores a 0 y menores a 2 millones de pesos. Recuerde ver la clase de las variables (4 puntos)

class(datos$ss_t)
## [1] "numeric"
datos=as.data.table(datos)
ggplot(datos[0<ss_t&ss_t<2000000], aes(x=ss_t,fill=region))+
  geom_histogram(bins=30)+facet_wrap(facets = "region")

### P3 ### Arregle los gráficos para que cada uno tenga un eje x legible, además arregle el eje y. Agregue color a cada gráfico:(4 puntos) Hint: utilice la función scales = ‘free_y’ dentro de facet_wrap para dejar libre el eje y. También recuerde chequear la clase de la variable region para que se asigne correctamente un color a cada región.

datos$region=as.character(datos$region)
ggplot(datos[0<ss_t&ss_t<2000000], aes(x=ss_t,fill=region))+
  geom_histogram(bins=30)+facet_wrap(facets = "region",scales = 'free_y')+
  theme(axis.text.x = element_text(angle=90,vjust=1))

P4

Realice un gráfico de dispersión que muestre la relación entre el ingreso y la edad. (4 puntos)

ggplot(datos,aes(x=ss_t, y=edad)) + geom_point()+scale_x_continuous(labels=scales::comma)+
  labs(title = "Diagrama de dispersión",
       subtitle = "ScatterPlot",
       caption = "Fuente: INE",
       x = "Edad",
       y = "Ingresos")

P4

Un miembro de su equipo propone un modelo que calcule la incidencia del sexo y la educación en el ingreso, sin considerar variables adicionales. Usted quiere mostrarle que aquel modelo está incompleto. Para lograr esto, haga el modelo de regresión anteriormente mencionado. (5 puntos)

datos= datos[!is.na(datos$ss_t),]
datos= datos[!is.na(datos$sexo),]
datos= datos[!is.na(datos$cine),]
reg_1 <- lm(ss_t~sexo + cine  , data = datos)
get_regression_table(reg_1) %>% kable()
term estimate std_error statistic p_value lower_ci upper_ci
intercept 1.004020e+14 1.768835e+12 56.762 0.000 9.693508e+13 1.038689e+14
sexo -2.326727e+13 1.103550e+12 -21.084 0.000 -2.543021e+13 -2.110432e+13
cine 4.564607e+10 1.359551e+10 3.357 0.001 1.899902e+10 7.229312e+10

Otra forma

summary(reg_1)
## 
## Call:
## lm(formula = ss_t ~ sexo + cine, data = datos)
## 
## Residuals:
##        Min         1Q     Median         3Q        Max 
## -1.227e+14 -7.727e+13 -5.410e+13 -5.400e+13  9.452e+14 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.004e+14  1.769e+12  56.762  < 2e-16 ***
## sexo        -2.327e+13  1.104e+12 -21.084  < 2e-16 ***
## cine         4.565e+10  1.360e+10   3.357 0.000787 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.71e+14 on 96237 degrees of freedom
## Multiple R-squared:  0.004709,   Adjusted R-squared:  0.004688 
## F-statistic: 227.6 on 2 and 96237 DF,  p-value: < 2.2e-16

P4.1

Calcule el salario predicho para cada observación.

pred1<-predict(reg_1) ## calcular los y predichos para cada observación
datos[,pred1:=predict(reg_1)]

En adición para tener una respuesta más completa, agregaremos la siguiente tabla.

predicciones1<-data.table(RMSE=RMSE(pred1,datos$ss_t),
                         MAE=MAE(pred1,datos$ss_t))

predicciones1
##           RMSE          MAE
## 1: 1.71014e+14 1.076885e+14

P4.2

Dado el modelo anterior, calcule la predicción del salario para una mujer y para un hombre con educación Universitario (5 puntos)

datos[cine==7,.N,by=.(sexo,pred1)]
##    sexo        pred1    N
## 1:    1 7.745423e+13 6224
## 2:    2 5.418696e+13 6922

Como grupo nos llamo la atencion que los salarios predichos den tan bajos (pese a tener educacion universitaria), por lo que decidimos comprobar cuantas mujeres y hombres en la base de datos que poseen educacion universitaria poseen un ss_t=0. Esto se puede evidenciar en la siguiente tabla

datos[ss_t==0 & cine==7,.N,by=sexo]
##    sexo    N
## 1:    1 3079
## 2:    2 3508

Podemos comprobar que efectivamente, existen 3079 datos de hombres con educación universitaria que tienen ss_t = 0, y tambien 3508 datos de mujeres con educación universitaria que tienen ss_t = 0.

P.5

Teniendo ya el modelo propuesto por su colega, contrarreste usted con un modelo de predicción en base a los datos con los que ya cuenta. Utilice el argumento na.action por los datos nulos en la función de la regresión.

datos[, factsexo := as.factor(sexo)]
datos[, facttipo := as.factor(tipo)]

reg_2 <- lm(ss_t~factsexo +  cine + facttipo +edad, data = datos,na.action=na.exclude)

summary(reg_2)
## 
## Call:
## lm(formula = ss_t ~ factsexo + cine + facttipo + edad, data = datos, 
##     na.action = na.exclude)
## 
## Residuals:
##        Min         1Q     Median         3Q        Max 
## -1.423e+14 -7.349e+13 -6.005e+13 -3.956e+13  9.652e+14 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  6.116e+13  1.254e+12  48.772  < 2e-16 ***
## factsexo2   -2.497e+13  1.102e+12 -22.668  < 2e-16 ***
## cine         3.162e+10  1.356e+10   2.332  0.01971 *  
## facttipo2   -4.733e+12  1.623e+12  -2.917  0.00353 ** 
## facttipo3   -1.826e+13  1.393e+12 -13.105  < 2e-16 ***
## edad         5.504e+11  2.381e+10  23.115  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.704e+14 on 96234 degrees of freedom
## Multiple R-squared:  0.01175,    Adjusted R-squared:  0.0117 
## F-statistic: 228.9 on 5 and 96234 DF,  p-value: < 2.2e-16

P 5.1

Calcule el salario predicho para cada observación. (4 puntos)

pred2<-predict(reg_2) ## calcular los y predichos para cada observación
datos[,pred1:=predict(reg_2)]

En adición para tener una respuesta más completa, agregaremos la siguiente tabla.

predicciones2<-data.table(RMSE=RMSE(pred2,datos$ss_t),
                         MAE=MAE(pred2,datos$ss_t))

predicciones2
##            RMSE          MAE
## 1: 1.704077e+14 1.067312e+14

P6

Contando ya con su modelo, usted puede realizar una comparación directa entre ambos. Calcule los errores de predicción dentro de muestra para ambos modelos. (8 puntos)

predicciones1<-data.table(RMSE=RMSE(pred1,datos$ss_t),
                         MAE=MAE(pred1,datos$ss_t))

predicciones1
##           RMSE          MAE
## 1: 1.71014e+14 1.076885e+14
predicciones2<-data.table(RMSE=RMSE(pred2,datos$ss_t),
                         MAE=MAE(pred2,datos$ss_t))

predicciones2
##            RMSE          MAE
## 1: 1.704077e+14 1.067312e+14

P7

Interprete la diferencia en los errores de predicción entre el Modelo 1 y el Modelo 2. ¿Qué modelo hace una mejor predicción dentro de muestra? (5 puntos)

Se puede ver que tanto el error cuadrático medio como error absoluto medio de la segunda regresion estan un poco mejor estimados que el de la regresion 1, sin embargo estos errores igual son grandes, por lo que habría que agregar mas variables a la regresion para tener una regresion mejor estimada.El RMSE, podemos analizar la diferencia entre el valor predicho y el valor observado, dando asi, el residuo. Se puede visualizar que efectivamente es menor el RMSE del modelo 2 con respecto al modelo 1, pero de todas formas, creemos que se debe disminuir el area al cuadrado de los residuos, para que asi, la prediccion sea mucho mas precisa.Con respecto al MAE, nos ayuda a cuantificar nuestra prediccion tomando dos variables, la predicha y la observada, entonces nos da un promedio de las diferencias en valor absoluto Cuando todas las diferencias individuales se ponderán por igual en el promedio, entonces al disminuir se puede analizar que existe una mayor precision del modelo 2 con respecto al 1.

P8

Realice validación cruzada (CV) a los modelos de la pregunta anterior por el método K-folds con 5 folds. ¿Se mantienen las conclusiones obtenidas en el análisis dentro de muestra? (9 puntos)

library(jtools)
set.seed(12345) ## setear una semilla

setupKCV <- trainControl(method = "cv" , number = 5)

predkfolds1<-train(ss_t~sexo + cine ,data=datos,method="lm",trControl= setupKCV)

predkfolds2<-train(ss_t~factsexo +  cine + facttipo +edad,data=datos,method="lm",trControl= setupKCV)

print(predkfolds1)
## Linear Regression 
## 
## 96240 samples
##     2 predictor
## 
## No pre-processing
## Resampling: Cross-Validated (5 fold) 
## Summary of sample sizes: 76992, 76992, 76992, 76992, 76992 
## Resampling results:
## 
##   RMSE          Rsquared     MAE        
##   1.710092e+14  0.004790278  1.07691e+14
## 
## Tuning parameter 'intercept' was held constant at a value of TRUE
print(predkfolds2)
## Linear Regression 
## 
## 96240 samples
##     4 predictor
## 
## No pre-processing
## Resampling: Cross-Validated (5 fold) 
## Summary of sample sizes: 76992, 76992, 76992, 76992, 76992 
## Resampling results:
## 
##   RMSE          Rsquared    MAE         
##   1.703979e+14  0.01171366  1.067353e+14
## 
## Tuning parameter 'intercept' was held constant at a value of TRUE

Teniendo en cuenta los calculos realizados de Validacion Cruzada en el modelo de regresion 1 y 2, y podemos interpretar que el modelo 2 tiene un R cuadrado de un 1.9%, esto quiere decir, que el modelo con las variables de (factsexo-edad-fatctipo y cine) explica la variable independiente en un 1.9%. Por otro lado, el modelo 1, con tan solo considerar las variables (sexo-cine) explica la variable independiente en un 0.8%. A partir de lo aprendido en Econometria, podemos agregar que desde la teoria, al agregar mas variables dependientes al modelo, en este caso edad y tipo, aumenta tu r cuadrado, debido a que el modelo tiene mas variables a considerar siendo menos sesgado y con mayor precision.

En el lado de la validación cruzada, tiene el objetivo de ver que tan ajustado esta el modelo, entonces se debe ir entrenando el modelo con respecto a la muestra y cuando ya este entrenado, aplicarlo a otro grupo de muestra. Para lograr esto se debe realizar un proceso iterativo y su propósito es verificar si nuestro modelo puede ser aplicable para otro tipo de muestras, es decir, si aparte de ser aplicable para la muestra de la base de datos otorgado por el INE, puede ser igualmente aplicable a una base de datos diferente a la del INE y tener resultados similares, por lo que la eficacia del modelo sería aplicable para otras muestras.

En conclusion, podemos decir que si se mantienen las conclusiones obtenidas en el analisis debido a que el RMSE y el MAE disminuyeron del modelo 1 al modelo 2, no de forma significativa, pero disminuyeron, y eso tal como se analizo anteriormente,nos dice que el modelo 2 es mas preciso que el modelo 1. Ademas, podemos ver que el modelo 2 sigue teniendo un R-cuadrado fuera de muestra más alto que el R-cuadrado fuera de muestra del modelo original. Por lo que por un lado predice mejor el modelo 2.