library(tidyverse)
library(readxl)
library(openxlsx)
library(plotly)
library(flextable)
require(CGPfunctions)
library(knitr)
library(caret)
datos = read_excel("~/Datos_Rotación.xlsx")

Punto 1

  1. Seleccionar 3 variables categóricas (distintas de rotación) y 3 variables cuantitativas, que consideren estén relacionadas con la rotación. Nota: Justificar por que estas variables están relacionadas y que tipo de relación se espera (Hipótesis). Ejemplo: Se espera que las horas extra se relacionen con la rotación ya que las personas podrían desgastarse mas al trabajar horas extra y descuidan aspectos personales. La hipótesis es que las personas que trabajan horas extra tienen mayor posibilidad de rotar que las que no trabajan extra. (serian 6, una por variable).
glimpse(datos)
## Rows: 1,470
## Columns: 24
## $ Rotación                    <chr> "Si", "No", "Si", "No", "No", "No", "No", …
## $ Edad                        <dbl> 41, 49, 37, 33, 27, 32, 59, 30, 38, 36, 35…
## $ `Viaje de Negocios`         <chr> "Raramente", "Frecuentemente", "Raramente"…
## $ Departamento                <chr> "Ventas", "IyD", "IyD", "IyD", "IyD", "IyD…
## $ Distancia_Casa              <dbl> 1, 8, 2, 3, 2, 2, 3, 24, 23, 27, 16, 15, 2…
## $ Educación                   <dbl> 2, 1, 2, 4, 1, 2, 3, 1, 3, 3, 3, 2, 1, 2, …
## $ Campo_Educación             <chr> "Ciencias", "Ciencias", "Otra", "Ciencias"…
## $ Satisfacción_Ambiental      <dbl> 2, 3, 4, 4, 1, 4, 3, 4, 4, 3, 1, 4, 1, 2, …
## $ Genero                      <chr> "F", "M", "M", "F", "M", "M", "F", "M", "M…
## $ Cargo                       <chr> "Ejecutivo_Ventas", "Investigador_Cientifi…
## $ Satisfación_Laboral         <dbl> 4, 2, 3, 3, 2, 4, 1, 3, 3, 3, 2, 3, 3, 4, …
## $ Estado_Civil                <chr> "Soltero", "Casado", "Soltero", "Casado", …
## $ Ingreso_Mensual             <dbl> 5993, 5130, 2090, 2909, 3468, 3068, 2670, …
## $ Trabajos_Anteriores         <dbl> 8, 1, 6, 1, 9, 0, 4, 1, 0, 6, 0, 0, 1, 0, …
## $ Horas_Extra                 <chr> "Si", "No", "Si", "Si", "No", "No", "Si", …
## $ Porcentaje_aumento_salarial <dbl> 11, 23, 15, 11, 12, 13, 20, 22, 21, 13, 13…
## $ Rendimiento_Laboral         <dbl> 3, 4, 3, 3, 3, 3, 4, 4, 4, 3, 3, 3, 3, 3, …
## $ Años_Experiencia            <dbl> 8, 10, 7, 8, 6, 8, 12, 1, 10, 17, 6, 10, 5…
## $ Capacitaciones              <dbl> 0, 3, 3, 3, 3, 2, 3, 2, 2, 3, 5, 3, 1, 2, …
## $ Equilibrio_Trabajo_Vida     <dbl> 1, 3, 3, 3, 3, 2, 2, 3, 3, 2, 3, 3, 2, 3, …
## $ Antigüedad                  <dbl> 6, 10, 0, 8, 2, 7, 1, 1, 9, 7, 5, 9, 5, 2,…
## $ Antigüedad_Cargo            <dbl> 4, 7, 0, 7, 2, 7, 0, 0, 7, 7, 4, 5, 2, 2, …
## $ Años_ultima_promoción       <dbl> 0, 1, 0, 3, 2, 3, 0, 0, 1, 7, 0, 0, 4, 1, …
## $ Años_acargo_con_mismo_jefe  <dbl> 5, 7, 0, 0, 2, 6, 0, 0, 8, 7, 3, 8, 3, 2, …
  1. Variables Categóricas:
  1. Variables cuantitativas:

Planteamiento de hipótesis

  1. Se espera que la obligación de realizar viajes de negocios y la frecuencia se relacionen con la rotación, puesto los viajes pueden implicar menos tiempo en familia y más desgaste. La hipótesis es que las personas que tienen mayor frecuencia de viajes de negocios tienen mayor probabilidad de rotar que las que tienen menor frecuencia de viaje.

  2. Se espera que el estado civil se relacione con la rotación pues el estar casado, soltero o divorciado puede alterar el riesgo que una persona puede asumir frente al cambio de trabajo debido a la responsabilidad familiar. La hipótesis es que las personas que están solteras tienen mayor probabilidad de rotar que las que tienen un estado civil diferente.

  3. Se espera que el nivel de rendimiento laboral esté relacionado con la rotación pues este puede ser reflejo de la competencia y conformidad del empleado con su trabajo, lo que puede llevar a decidir retirarse o permanecer en su puesto. La hipótesis es que las personas que tienen mayor rendimiento laboral tienen menos probabilidad de rotación que aquellos con menor rendimiento.

  4. Se espera que la edad esté relacionada con la rotación puesto que a medida que las personas envejecen tienen más dificultad para encontrar otros empleos por lo que puede afectar su intención de rotar. La hipótesis es que las personas a medida que aumentan en edad, su probabilidad de rotar disminuye.

  5. Se espera que el Ingreso Mensual esté relacionado con la rotación pues el tener un cierto nivel de ingreso puede ser incentivo para permanecer o querer irse de un puesto de trabajo. La hipótesis es que las personas con mayor nivel de ingresos tienen menos probabilidad de rotar que aquellas con menos nivel de ingresos.

  6. Se espera que la Antiguedad en el Cargo esté relacionada con la rotación pues esta puede ser reflejo de un estancamiento en la posición llevando a inconformidad por parte del empleado. La hipótesis es que las personas con más antiguedad en el cargo tienen mayor probabilidad de rotación que las de menos antiguedad.

Punto 2

  1. Realizar un análisis univariado (caracterización). Nota: Los indicadores o gráficos se usan dependiendo del tipo de variable (cuanti o cuali). Incluir interpretaciones de la rotación.

Empezamos observando la distribución proporcional de cada una de las variables escogidas.

library(patchwork)
p1 <- datos %>% 
  ggplot(aes("Viaje de Negocios")) +
  geom_bar(aes(y = (..count..)/sum(..count..))) +
  theme(axis.text.x = element_text(angle = 15, hjust = 1)) +
  labs(x = "Viaje de Negocios", y = "Proporción")
p2 <- datos %>% 
  ggplot(aes(Estado_Civil)) +
  geom_bar(aes(y = (..count..)/sum(..count..))) +
  labs(x = "Estado Civil", y = "")
p3 <- datos %>% 
  ggplot(aes(Rendimiento_Laboral)) +
  geom_bar(aes(y = (..count..)/sum(..count..))) +
  labs(y = "")
p4 <- datos %>% 
  ggplot(aes(Edad)) +
  geom_histogram(aes(y = ..density..), binwidth = 5) +
  labs(y = "Densidad")
p5 <- datos %>% 
  ggplot(aes(Ingreso_Mensual)) +
  geom_histogram(aes(y = ..density..), binwidth = 2000) +
  labs(y = "")
p6 <- datos %>% 
  ggplot(aes(Antigüedad_Cargo)) +
  geom_histogram(aes(y = ..density..), binwidth = 1) +
  labs(y = "")

p1 + p2 + p3 + p4 + p5 + p6

Se observa que la mayoría de los empleados (70% aprox.) viajan raramente; aproximadamente el 45% de los empleados están casados y el 30% están solteros; la mayoría de empleados tiene un rendimiento de 3 (85% aprox.); la mayoría de los empleados se encuentra en el rango entre 30 y 40 años; la mayor parte de los empleados tiene un ingreso inferior a las 7,000 unidades monetarias; el 25% del personal tiene 2 años de antiguedad en su cargo, el 15% tiene 7 años de antiguedad en su cargo y el 17% tiene menos de un año de antiguedad en el cargo.

datos %>% 
  ggplot(aes(Rotación)) +
  geom_bar(aes(y = (..count..)/sum(..count..))) +
  labs(y = "Proporción")

Se observa que hay una relativamente baja rotación de empleados, correspondiente a menos del 20% del personal.

Punto 3

p_viajes <- datos %>% 
  ggplot(aes("Viaje de Negocios", fill = Rotación)) +
  geom_bar(position = "dodge") +
  labs(y = "Cantidad") +
  theme(axis.text.x = element_text(angle = 15, hjust = 1))
p_viajes_prop <- datos %>% 
  ggplot(aes("Viaje de Negocios", fill = Rotación)) +
  geom_bar(aes(y = ..count..), position = "fill") +
  labs(y = "Proporcion") +
  theme(axis.text.x = element_text(angle = 15, hjust = 1))
  

p_viajes + p_viajes_prop

Se observa que los empleados que la mayor parte de empleados que rotan raramente viajan, no obstante la mayor parte de los empleados que no rotan también pertenecen a los que raramente viajan. Ahora bien, cuando analizamos las categorías proporcionalmente observamos que el grupo que más rotación tiene es el de los que viajan frecuentemente con 25% de rotación, seguido por los que raramente viajan y por último los que no viajan. Por tanto, se apoya la hipótesis de que los que más viajan tienen más probabilidad de rotar.

p_civil <- datos %>% 
  ggplot(aes(Estado_Civil, fill = Rotación)) +
  geom_bar(position = "dodge") +
  labs(y = "Cantidad") +
  theme(axis.text.x = element_text(angle = 15, hjust = 1))

p_civil_prop <- datos %>% 
  ggplot(aes(Estado_Civil, fill = Rotación)) +
  geom_bar(aes(y = ..count..), position = "fill") +
  labs(y = "Proporcion") +
  theme(axis.text.x = element_text(angle = 15, hjust = 1))
  
p_civil + p_civil_prop

Se observa que a pesar de que los casados representan la mayor cantidad de empleados, los solteros representan la mayor cantidad de rotaciones así como también cuentan con la mayor proporción de rotación como categoría (25% aproximadamente). Lo cual apoya la hipótesis planteada de que las personas que están solteras tienen mayor probabilidad de rotar que las que tienen un estado civil diferente

p_rendimiento <- datos %>% 
  ggplot(aes(Rendimiento_Laboral, fill = Rotación)) +
  geom_bar(position = "dodge") +
  labs(y = "Cantidad") +
  theme(axis.text.x = element_text(angle = 15, hjust = 1))

p_rendimiento_prop <- datos %>% 
  ggplot(aes(Rendimiento_Laboral, fill = Rotación)) +
  geom_bar(aes(y = ..count..), position = "fill") +
  labs(y = "Proporcion") +
  theme(axis.text.x = element_text(angle = 15, hjust = 1))
  
p_rendimiento + p_rendimiento_prop

A pesar de que el mayor número de empleados que rotan tienen un rendimiento laboral de 3, al hacer el análisis proporcional no se encuentran diferencias en la proporción de rotación entre los distintos niveles de rendimiento, lo cual no apoya la hipótesis planteada de que las personas que tienen mayor rendimiento laboral tienen menos probabilidad de rotación que aquellos con menor rendimiento.

p_edad <- datos %>% 
  ggplot(aes(Rotación, Edad, fill = Rotación)) +
  geom_boxplot() +
  geom_jitter(alpha = 0.1) +
  guides(fill = "none")
  
ggplotly(p_edad)

Se observa que los individuos que no rotan tienen una edad mediana de 36 años, mientras que los que sí rotan cuentan con una edad mediana de 32 años. Lo cual no apoya la hipótesis planteada de que las personas a medida que aumentan en edad, su probabilidad de rotar disminuye.

p_ingreso <- datos %>% 
  ggplot(aes(Rotación, Ingreso_Mensual, fill = Rotación)) +
  geom_boxplot() +
  geom_jitter(alpha = 0.1) +
  guides(fill = "none")
  
ggplotly(p_ingreso)

Se observa que los empleados que no rotan tienen un ingreso mediano de 5,204 unidades monetarias mesuales, superior al ingreso mediano de los que rotan el cual se encuentra en 3,202 unidades mensuales. Lo cual apoya la hipótesis planteada de que las personas con mayor nivel de ingresos tienen menos probabilidad de rotar que aquellas con menos nivel de ingresos.

p_ant_cargo <- datos %>% 
  ggplot(aes(Rotación, Antigüedad_Cargo, fill = Rotación)) +
  geom_boxplot() +
  geom_jitter(alpha = 0.1) +
  guides(fill = "none")
  
ggplotly(p_ant_cargo)

Se observa que los empleados que no rotan tienen una antiguedad en el cargo mediana de 3 años, mientras que los que rotan tienen una antiguedad mediana de 2 años. Lo cual no apoya la hipótesis planteada de que las personas con más antiguedad en el cargo tienen mayor probabilidad de rotación que las de menos antiguedad.

Punto 4

  1. Realizar la estimación de un modelo de regresión logistico en el cual la variable respuesta es rotación (y=1 es si rotación, y=0 es no rotación) y las covariables las 6 seleccionadas. Interprete los coeficientes del modelo y la significancia de los parametros.
rotacion_dummy <- datos %>% 
  mutate(Rotación = ifelse(Rotación == "Si", 1, 0)) 
set.seed(52)

rotacion_dummy_aleatorizada <- rotacion_dummy %>% 
  sample_frac(size = 1, replace = FALSE)
# Obtenemos el numero de observaciones
n_obs <- nrow(rotacion_dummy_aleatorizada)

# Identificamos la fila donde se partira en dos grupos: particion
particion <- round(n_obs * 0.8) # En este caso decidi utilizar la relacion 80/20

# Creamos un grupo de Entrenamiento (60% de los datos)
entrenamiento <- rotacion_dummy_aleatorizada %>% 
  slice(1:particion)

# Creamos un grupo de Testeo (40% de los datos)
testeo <- rotacion_dummy_aleatorizada %>% 
  slice(particion + 1:n_obs)

nrow(entrenamiento) / nrow(rotacion_dummy_aleatorizada)
## [1] 0.8
modelo_logit <- glm(Rotación ~ factor(`Viaje de Negocios`) + factor(Estado_Civil)  + factor(Rendimiento_Laboral) + Ingreso_Mensual + Antigüedad_Cargo +  Edad, family = binomial, data = entrenamiento)

as_flextable(modelo_logit)

Probamos el modelo con nuestra base de testeo

rotacion_predicho <- predict(modelo_logit, testeo, type = "response")

broom::augment(modelo_logit, newdata = testeo, type.predict = "response")
## # A tibble: 294 × 25
##    Rotación  Edad Viaje d…¹ Depar…² Dista…³ Educa…⁴ Campo…⁵ Satis…⁶ Genero Cargo
##       <dbl> <dbl> <chr>     <chr>     <dbl>   <dbl> <chr>     <dbl> <chr>  <chr>
##  1        0    57 Raramente IyD          24       2 Cienci…       3 M      Repr…
##  2        0    33 Raramente IyD           7       3 Salud         3 M      Dire…
##  3        0    39 Raramente IyD          10       5 Salud         2 M      Dire…
##  4        0    42 Raramente IyD          29       3 Cienci…       2 M      Dire…
##  5        0    34 Frecuent… RH           11       3 Cienci…       3 M      Recu…
##  6        0    38 Raramente IyD          15       2 Cienci…       3 M      Dire…
##  7        1    35 Frecuent… Ventas       18       4 Mercad…       4 F      Ejec…
##  8        0    54 Raramente IyD          19       4 Salud         4 F      Dire…
##  9        0    41 Raramente IyD           1       3 Cienci…       4 M      Inve…
## 10        0    56 No_Viaja  IyD           1       4 Cienci…       3 M      Repr…
## # … with 284 more rows, 15 more variables: Satisfación_Laboral <dbl>,
## #   Estado_Civil <chr>, Ingreso_Mensual <dbl>, Trabajos_Anteriores <dbl>,
## #   Horas_Extra <chr>, Porcentaje_aumento_salarial <dbl>,
## #   Rendimiento_Laboral <dbl>, Años_Experiencia <dbl>, Capacitaciones <dbl>,
## #   Equilibrio_Trabajo_Vida <dbl>, Antigüedad <dbl>, Antigüedad_Cargo <dbl>,
## #   Años_ultima_promoción <dbl>, Años_acargo_con_mismo_jefe <dbl>,
## #   .fitted <dbl>, and abbreviated variable names ¹​`Viaje de Negocios`, …
## # ℹ Use `print(n = ...)` to see more rows, and `colnames()` to see all variable names

Clasificacion binaria (1)

Podemos comenzar a clasificar las observaciones de nuestra base de testeo, la cual hace las veces de conjunto de observaciones (out-of-sample). Pasamos de tener meras probabilidades a asignar una clasificacion definida partiendo de las probabilidades propuestas por nuestro modelo.

modelo_logit_clasificado <- broom::augment(modelo_logit, newdata = testeo, type.predict = "response") %>% 
  mutate(rot_predic = ifelse(.fitted > 0.5, 1, 0)) 

knitr::kable(slice_sample(modelo_logit_clasificado %>% 
  select(Rotación, `Viaje de Negocios`, Estado_Civil, Antigüedad_Cargo, Edad, Ingreso_Mensual, Rendimiento_Laboral, rot_predic), n = 10), align = 'ccccc')
Rotación Viaje de Negocios Estado_Civil Antigüedad_Cargo Edad Ingreso_Mensual Rendimiento_Laboral rot_predic
0 Frecuentemente Divorciado 5 26 6397 4 0
0 Raramente Soltero 0 38 2858 3 0
0 No_Viaja Divorciado 5 33 5368 4 0
0 No_Viaja Divorciado 2 28 2706 3 0
0 Raramente Casado 2 33 2911 3 0
1 Raramente Casado 2 52 4941 3 0
1 Raramente Soltero 2 39 3904 3 0
0 Raramente Soltero 4 25 4487 3 0
0 Raramente Casado 2 23 2819 3 0
0 Frecuentemente Casado 7 36 9738 3 0

Matriz de Confusion: Que tan acertado fue mi modelo para predecir los resultados de rotación en esta base de datos?

Existe una herramienta muy util para observar que tan acertado fue nuestro modelo para predecir una variable categorica (out-of-sample), la Matriz de Confusion.

Esta es una tabla que nos permite observar el desempeno de un algoritmo (tipicamente de aprendizaje supervisado) para predecir una variable categorica.

prediccion_categoria <- ifelse(rotacion_predicho > 0.5, 1, 0)

table(prediccion_categoria, testeo[["Rotación"]])
##                     
## prediccion_categoria   0   1
##                    0 249  40
##                    1   1   4

También podemos usar tambien la funcion confusionMatrix() del paquete caret

caret::confusionMatrix(as.factor(prediccion_categoria), as.factor(testeo[["Rotación"]]))
## Confusion Matrix and Statistics
## 
##           Reference
## Prediction   0   1
##          0 249  40
##          1   1   4
##                                          
##                Accuracy : 0.8605         
##                  95% CI : (0.8156, 0.898)
##     No Information Rate : 0.8503         
##     P-Value [Acc > NIR] : 0.3473         
##                                          
##                   Kappa : 0.1369         
##                                          
##  Mcnemar's Test P-Value : 2.946e-09      
##                                          
##             Sensitivity : 0.99600        
##             Specificity : 0.09091        
##          Pos Pred Value : 0.86159        
##          Neg Pred Value : 0.80000        
##              Prevalence : 0.85034        
##          Detection Rate : 0.84694        
##    Detection Prevalence : 0.98299        
##       Balanced Accuracy : 0.54345        
##                                          
##        'Positive' Class : 0              
##

Curva ROC: Cual debe ser mi matriz de confusion ideal?

Tambien podemos calcular la curva ROC, en la cual cada uno de los puntos representa una matriz de confusion con un umbral de probabilidad distinto. En el eje Y tiene la Sensibilidad de las predicciones y en el eje X tiene la probabilidad de falsos positivos (1 - Especificidad)

caTools::colAUC(rotacion_predicho, testeo[["Rotación"]], plotROC = TRUE)

##              [,1]
## 0 vs. 1 0.7489091

Cada punto sobre la curva es un posible umbral de prediccion, pues, cada uno tiene una sensibilidad y una especificidad distinta. En otras palabras, cada uno de los puntos representa una matriz de confusion que no hemos tenido que evaluar a mano.

Area Bajo la Curva ROC

La curva ROC nos da una buena impresion acerca de la calidad de nuestro modelo predictivo. Es importante notar que los modelos con predicciones aleatorias, tienden a tener curvas ROC cercanas a la diagonal de la grafica. Por otro lado, los modelos de clasificacion que permiten una separacion perfecta de las categorias, tienen curvas que producen una “caja” con un unico punto en (1, 0) la esquina superior izquierda del grafico. Pues representan un modelo en el que es posible alcanzar una tasa de 100% verdaderos positivos (sensibilidad) y 0% falsos negativos (1 - Especificidad).

Asi, cuando se calcula el area bajo la curva ROC de un modelo perfecto, el resultado es exactamente igual a 1. Esto nos da una manera de formalizar una medida de precision (o exactitud) de nuestro modelo utilizando el Area bajo la Curva.

  • La medida va de 0 a 1
    • 0.5 es la medida de un modelo de prediccion aleatorio
    • 1.0 es la medida de un modelo perfecto
    • 0.0 es la medida de un modelo que siempre esta errado
# Usamos el paquete `caret` que facilita el trabajo: Debemos crear un objeto de tipo trainControl 
miControl <- trainControl(
  method = "cv", # Esta funcion hace la validacion cruzada automaticamente a partir de la base de datos original completa.
  number = 10,
  summaryFunction = twoClassSummary,
  classProbs = TRUE, # IMPORTANT!
  verboseIter = TRUE
)

modelo_area_bajo_curva <- train(
    Rotación ~ factor(`Viaje de Negocios`) + factor(Estado_Civil)  + factor(Rendimiento_Laboral) + Ingreso_Mensual + Antigüedad_Cargo +  Edad, 
    datos,
    method = "glm",
    trControl = miControl
    )
## Warning in train.default(x, y, weights = w, ...): The metric "Accuracy" was not
## in the result set. ROC will be used instead.
## + Fold01: parameter=none 
## - Fold01: parameter=none 
## + Fold02: parameter=none 
## - Fold02: parameter=none 
## + Fold03: parameter=none 
## - Fold03: parameter=none 
## + Fold04: parameter=none 
## - Fold04: parameter=none 
## + Fold05: parameter=none 
## - Fold05: parameter=none 
## + Fold06: parameter=none 
## - Fold06: parameter=none 
## + Fold07: parameter=none 
## - Fold07: parameter=none 
## + Fold08: parameter=none 
## - Fold08: parameter=none 
## + Fold09: parameter=none 
## - Fold09: parameter=none 
## + Fold10: parameter=none 
## - Fold10: parameter=none 
## Aggregating results
## Fitting final model on full training set
# Imprimimos nuestro modelo
modelo_area_bajo_curva
## Generalized Linear Model 
## 
## 1470 samples
##    6 predictor
##    2 classes: 'No', 'Si' 
## 
## No pre-processing
## Resampling: Cross-Validated (10 fold) 
## Summary of sample sizes: 1324, 1323, 1323, 1323, 1323, 1323, ... 
## Resampling results:
## 
##   ROC        Sens       Spec      
##   0.7082853  0.9983805  0.07192029

Nuestra medida de Area Bajo la Curva nos resulto en 0.70, lo que implica que nuestro modelo bueno, con una sensibilidad de 99% y una especificidad de 7.1%.