Informe para el SEREMI de desarrollo social

P1

Cargue los paquetes “data.table”, “ggplot2” y “caret”, junto con la base de datos. (1 punto)

library(data.table)
library(ggplot2)
library(caret)

esi2019<-fread("C:/Users/Constanza/Downloads/esi-2019---personas.csv")

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)

esi2019[,ss_t:=as.numeric(ss_t)]
str(esi2019$ss_t)
##  num [1:96240] 0 229822 0 0 1100000 ...
esi2019[,region:=as.character(region)]
str(esi2019$region)
##  chr [1:96240] "13" "13" "11" "11" "11" "11" "11" "11" "8" "8" "8" "11" ...
ggplot(data=esi2019[ss_t>0 & ss_t<2000000], aes(x=ss_t)) + geom_histogram() + 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.

esi2019[,ss_t:=as.numeric(ss_t)]
esi2019[,region:=as.character(region)]

ggplot(data=esi2019[ss_t>0 & ss_t<2000000], aes(x=ss_t, fill=region)) + 
  geom_histogram(color="black") + 
  facet_wrap(~region, scale="free_y") + theme(axis.text.x = element_text(angle=75, vjust=0.6)) + labs(x="Ingresos", y="", title = "Ingresos por región", caption = "Fuente: INE" )

P4

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

esi2019<-esi2019[!is.na(esi2019$ss_t)]

ggplot(esi2019,  aes(x=edad, y=ss_t)) + geom_point() + labs(title="Relación Edad-Ingreso", x="Edad", y="Ingresos", caption="Fuente: INE") + scale_y_continuous(labels=function(n){format(n, scientific = FALSE)})

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)

reg_1 <- lm(ss_t~ sexo + cine , data = esi2019)
get_regression_table(reg_1) %>% kable()
term estimate std_error statistic p_value lower_ci upper_ci
intercept 170227.110 3498.074 48.663 0 163370.906 177083.31
sexo -54407.560 2163.093 -25.153 0 -58647.209 -50167.91
cine 202.736 26.849 7.551 0 150.112 255.36

P4.1

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

esi2019[,prediccion_1:=predict(reg_1)]

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)

esi2019_2<-esi2019[cine==7]
esi2019_2[,mean(prediccion_1), by= "sexo"]
##    sexo        V1
## 1:    1 117238.70
## 2:    2  62831.14

La predicción de salario para una mujer con educación universitaria será de 62831,14 y para un hombre con educación universitaria será de 117238,7.

P5

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. (7 puntos)

reg_2 <- lm(ss_t~ sexo+cine+edad+region , data = esi2019, na.action(esi2019))
get_regression_table(reg_2) %>% kable()
term estimate std_error statistic p_value lower_ci upper_ci
intercept 129401.481 6986.795 18.521 0.000 115707.406 143095.556
sexo -57576.095 2153.554 -26.735 0.000 -61797.048 -53355.142
cine 166.899 26.688 6.254 0.000 114.591 219.207
edad 956.663 43.742 21.871 0.000 870.929 1042.398
region10 347.659 7153.907 0.049 0.961 -13673.956 14369.273
region11 18748.520 8819.281 2.126 0.034 1462.782 36034.258
region12 57918.561 9621.285 6.020 0.000 39060.899 76776.223
region13 53681.685 6490.757 8.270 0.000 40959.839 66403.530
region14 3595.879 8245.925 0.436 0.663 -12566.084 19757.842
region15 -19177.859 8270.124 -2.319 0.020 -35387.253 -2968.466
region16 -22530.338 8771.861 -2.568 0.010 -39723.134 -5337.542
region2 57365.309 8414.961 6.817 0.000 40872.036 73858.581
region3 10311.325 8752.041 1.178 0.239 -6842.623 27465.272
region4 -15232.933 7442.489 -2.047 0.041 -29820.167 -645.699
region5 8488.458 6748.328 1.258 0.208 -4738.225 21715.140
region6 -10040.376 7510.290 -1.337 0.181 -24760.499 4679.747
region7 -18659.604 7322.094 -2.548 0.011 -33010.864 -4308.343
region8 -13940.893 6926.813 -2.013 0.044 -27517.404 -364.382
region9 -19130.623 7573.619 -2.526 0.012 -33974.871 -4286.374

P5.1

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

esi2019[,prediccion_2:=predict(reg_2)]

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)

prediccion_1<-predict(reg_1) 
esi2019[,prediccion_1:=predict(reg_1)]

predicciones1<-data.table(RMSE=RMSE(prediccion_1,esi2019$ss_t), MAE=MAE(prediccion_1,esi2019$ss_t))

prediccion_2<-predict(reg_2) 

predicciones2<-data.table(RMSE=RMSE(prediccion_2,esi2019$ss_t), MAE=MAE(prediccion_2,esi2019$ss_t))

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)

La raíz del error cuadrático medio de la predicción del modelo 1 (de nuestro colega) es 303322.1 mientras que la de la predicción del modelo 2 es 301121.3, siendo la diferencia entre modelo 1 y modelo 2 de 2200,8. Por otro lado, la media absoluta del error de la predicción del modelo 1 es de 147367 y la de la predicción del modelo 2 es de 144524.1, siendo la diferencia entre el modelo 1 y el modelo 2 de 2842,9. Estos datos nos muestran que los indicadores son menores en el segundo modelo y, por ende, este hace una mejor predicción de la muestra.

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)
Importante: Utilice set.seed(12345). Si no sabe qué hacer con los NA, lea la documentación de la fución para saber como omitirlos.

set.seed(12345) 

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

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

predkfolds2<-train(ss_t ~ sexo+cine+edad+region, data=esi2019, method="lm", trControl= setup)

print(predkfolds1)
## Linear Regression 
## 
## 79088 samples
##     2 predictor
## 
## No pre-processing
## Resampling: Cross-Validated (5 fold) 
## Summary of sample sizes: 63270, 63270, 63271, 63271, 63270 
## Resampling results:
## 
##   RMSE      Rsquared     MAE     
##   302710.7  0.008607358  147381.8
## 
## Tuning parameter 'intercept' was held constant at a value of TRUE
print(predkfolds2)
## Linear Regression 
## 
## 79088 samples
##     4 predictor
## 
## No pre-processing
## Resampling: Cross-Validated (5 fold) 
## Summary of sample sizes: 63270, 63270, 63270, 63271, 63271 
## Resampling results:
## 
##   RMSE      Rsquared    MAE     
##   301062.8  0.02257009  144564.2
## 
## Tuning parameter 'intercept' was held constant at a value of TRUE

Dado lo anterior, podemos concluir que se mantienen las conclusiones obtenidas dentro de la muestra ya que siguen siendo menores los errores del modelo 2, sugiriendo que es un mejor predictor.