Informe para la SEREMI de desarrollo social

Suponga que a usted le han contratado desde la SEREMI de desarrollo social para encontrar algún modelo que haga una predicción del salario de las personas. Para ello, le han proporcionado algunos datos sobre el salario para distintos individuos, con los cuales usted deberá encontrar algún vínculo entre alguna variable interesante y el salario. Para esto:

rm(list=ls())

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

Importante:Borre todas las observaciones que tienen datos ´NA´ en alguna de las variables. El diccionario de variables se encuentra aquí

library(data.table)
library(ggplot2)
library(moderndive)
library(knitr)
library(lattice)
library(caret)
seremi <- fread("/Users/kiwixx/Desktop/Data Science/datos_tarea_5.csv")
seremi <- as.data.table(seremi)
seremi <- seremi[,ss_t:=as.integer(seremi$ss_t)]
## Warning in eval(jsub, SDenv, parent.frame()): NAs introduced by coercion
library(tidyr)
seremi <- seremi[!is.na(seremi$ss_t)]
seremi <- seremi[!is.na(seremi$sexo)]
seremi <- seremi[!is.na(seremi$nacionalidad)]
seremi <- seremi[!is.na(seremi$edad)]
seremi <- seremi[!is.na(seremi$r_p_c)]
seremi <- seremi[!is.na(seremi$region)]
seremi <- seremi[!is.na(seremi$cise)]

2. 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(seremi$region)
## [1] "integer"
ingreso_por_region <- seremi[ss_t < 2000000 & ss_t > 0, ss_t, by = "region"]
#ippr <- ingreso_por_region[, (mean(ss_t)), by = "region"] ## areglar
ggplot(ingreso_por_region, aes(x=ss_t)) +
  geom_histogram(bins = 10) +  
  facet_wrap(facets="region")

3. 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.

ggplot(ingreso_por_region, aes(x=ss_t, fill=region)) +
  geom_histogram() +
  theme(axis.text.x = element_text(angle=75, vjust=0.5)) +
  facet_wrap(facets= "region", scales="free_y")
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

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

Grafico con filtro en los ingresos (solicitado en preguntas anteriores)

ggplot(seremi[ss_t > 0 & ss_t < 2000000 & !is.na(ss_t)], aes(edad, ss_t)) +
  geom_point(color = "Blue", size = 0.8)

Grafico sin filtro, y mostrando la completa disperción de la edad~ingreso

ggplot(seremi[!is.na(ss_t)], aes(edad, ss_t)) +
  geom_point(color = "Blue", size = 0.8) +
  labs(x = "Edad", y = "Ingreso Promedio")

4. 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)

modelo1 <- lm(ss_t ~ sexo + nivel, data = seremi)
get_regression_table(modelo1) %>% kable()
term estimate std_error statistic p_value lower_ci upper_ci
intercept 169925.514 3498.895 48.565 0 163067.701 176783.326
sexo -54460.931 2162.967 -25.179 0 -58700.333 -50221.529
nivel 298.479 36.080 8.273 0 227.762 369.197

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

seremi <- seremi[,predicciones_salario:=predict(modelo1)]

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

HOMBRE

regressionHombre <- lm(ss_t ~ (sexo==1)+ (nivel==9) ,data = seremi)
get_regression_table(regressionHombre) %>% kable()
term estimate std_error statistic p_value lower_ci upper_ci
intercept 36690.96 1487.700 24.663 0 33775.08 39606.84
sexo == 1TRUE 53098.04 2107.393 25.196 0 48967.56 57228.52
nivel == 9TRUE 205527.95 3133.791 65.584 0 199385.73 211670.16 
seremi <- seremi[,predicciones_salariohombre_universitario:=predict(regressionHombre)]

MUJER

regressionMujer <- lm(ss_t ~  (sexo==2) + (nivel==9) ,data = seremi)
get_regression_table(regressionMujer) %>% kable()
term estimate std_error statistic p_value lower_ci upper_ci
intercept 89789.00 1598.816 56.160 0 86655.33 92922.67
sexo == 2TRUE -53098.04 2107.393 -25.196 0 -57228.52 -48967.56
nivel == 9TRUE 205527.95 3133.791 65.584 0 199385.73 211670.16 
seremi[,predicciones_salario_mujer_universitario:=predict(regressionMujer)]

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

propuesta <- lm(ss_t ~ sexo + 
                  nivel + 
                  nacionalidad + 
                  edad + 
                  r_p_c+region + 
                  cise  ,data = seremi)

get_regression_table(propuesta) %>% kable()
term estimate std_error statistic p_value lower_ci upper_ci
intercept 36646.039 4252.314 8.618 0 28311.529 44980.549
sexo -23329.579 1901.089 -12.272 0 -27055.701 -19603.456
nivel 125.410 31.477 3.984 0 63.715 187.105
nacionalidad -1.014 0.280 -3.615 0 -1.564 -0.464
edad -254.345 39.041 -6.515 0 -330.864 -177.825
r_p_c -35.683 5.782 -6.171 0 -47.016 -24.350
region 36757.988 5738.331 6.406 0 25510.893 48005.082
cise 107120.533 684.223 156.558 0 105779.460 108461.606

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

seremi[,Predic_5:=predict(propuesta)]
seremi[1:10,.(Predic_5)] %>% kable()
Predic_5
-18307.858
340877.452
-23676.547
6387.883
330798.095
307844.747
323160.700
-21021.787
302167.665
10487.218

6. 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)

p1 <- predict(modelo1)
pr1 <-data.table(RMSE=RMSE(p1,seremi$ss_t), MAE=MAE(p1,seremi$ss_t))
pr1
##        RMSE      MAE
## 1: 303300.2 147359.3

Lo anterior los explica la regression en relación a sexo~y estudios. ahora volvemos a predecir con lo anterior y lo propuesto anteriormente

p2 <- predict(propuesta)
pr2 <-data.table(RMSE=RMSE(p2,seremi$ss_t),
                        MAE=MAE(p1,seremi$ss_t))
pr2
##        RMSE      MAE
## 1: 264053.7 147359.3

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

Dentro de los dos modelos anteriores, podemos ver que el Error Absoluto Medio de la primera prediccion era mayor a el de la segunda. Esto quiere decir que la prediccion de el segundo modelo es “mejor” que el primero, dado los datos y el nivel predictivo obtenido.

8. 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)
setupKCV <- trainControl(method = "cv" , number = 5)
predkfolds1 <-train(ss_t ~ sexo + 
                      nivel ,data = seremi, method="lm", trControl= setupKCV)

predkfolds2<-train(ss_t ~ sexo + 
                     nivel + 
                     nacionalidad + 
                     edad + 
                     r_p_c + 
                     region + 
                     cise  ,data = seremi, method="lm", trControl= setupKCV)

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     
##   302679.1  0.008774494  147373.4
## 
## Tuning parameter 'intercept' was held constant at a value of TRUE

Dado la siguente prediccion, y con el cross validation, determinamos que dejamos todo tal como fue dicho en la pregunta 7, ya que el cross validation nos da un valor de Error Cuadratico Medio menor a los anteriores.

Con esto mismo, nos damos cuetna que nuestro RMSE en ambos casos sigue un comportamiento similar, por lo que nuevamente podemos concluir que las concluciones anteriores siguen iguales.