Logistic Regression Model

¿Qué es logistic regression ?

La regresión logística es un modelo estadístico que se utiliza para determinar la PROBABILIDAD de que ocurra un evento.

Loading libraries

library(foreign)
library(dplyr)        # data manipulation 
library(ggplot2)      # data visualization 
library(psych)        # functions for multivariate analysis 
library(corrplot)     # correlation plots
library(jtools)       # presentation of regression analysis 
library(lmtest)       # diagnostic checks - linear regression analysis 
library(car)          # diagnostic checks - linear regression analysis
library(factoextra)   # provides functions to extract and visualize the output of exploratory multivariate data analyses
library(ggfortify)    # data visualization tools for statistical analysis results
library(tidyverse)    # data manipulation and visualization
library(janitor)      # examining and cleaning datasets 
library(viridis)      # allows you to generate color palettes  
library(scales)       # scaling is a technique for comparing data that isn't measured in the same way
#install.packages("caret")
#library(caret)        # Classification And REgression Training. Set of functions that attempt to streamline the process 
                      # for creating predictive models. 
library(lubridate)    # makes it easier to work with dates and times

Preparación de los datos

Mediante excel, se prepararon los datos: se unieron las bases de bajas y de colaboradores, se mantuvieron las columnas que tuvieran en común (edad, salario, estado civil, genero, puesto), y se agregó la variable de “BAJA”.

rh_logistic<-read.csv("/Users/elenavela/Downloads/logistica.csv")
str(rh_logistic)

## 'data.frame':    336 obs. of  7 variables:
##  $ edad          : int  32 36 23 21 29 46 29 31 50 19 ...
##  $ salario_diario: num  500 152 152 152 152 ...
##  $ estado        : chr  "Nuevo Leon" "Nuevo Leon" "Nuevo Leon" "Nuevo Leon" ...
##  $ civil         : chr  "Soltero" "Union libre" "Matrimonio" "Soltero" ...
##  $ genero        : chr  "MASCULINO" "FEMENINO" "FEMENINO" "FEMENINO" ...
##  $ puesto        : chr  "DISENO" "AYUDANTE GENERAL" "AYUDANTE GENERAL" "AYUDANTE GENERAL" ...
##  $ baja          : int  1 1 1 1 1 1 1 1 1 1 ...

Comprendemos la base de datos nueva

summary(rh_logistic)

##       edad       salario_diario     estado             civil          
##  Min.   :18.00   Min.   :144.4   Length:336         Length:336        
##  1st Qu.:23.00   1st Qu.:176.7   Class :character   Class :character  
##  Median :30.00   Median :180.7   Mode  :character   Mode  :character  
##  Mean   :32.23   Mean   :178.4                                        
##  3rd Qu.:39.00   3rd Qu.:180.7                                        
##  Max.   :67.00   Max.   :500.0                                        
##     genero             puesto               baja       
##  Length:336         Length:336         Min.   :0.0000  
##  Class :character   Class :character   1st Qu.:0.0000  
##  Mode  :character   Mode  :character   Median :1.0000  
##                                        Mean   :0.6905  
##                                        3rd Qu.:1.0000  
##                                        Max.   :1.0000

Convertimos las variables a factor para poder mostrar los estadísticos descriptivos.

rh_logistic$genero<-as.factor(rh_logistic$genero)
rh_logistic$puesto<-as.factor(rh_logistic$puesto)
rh_logistic$civil<-as.factor(rh_logistic$civil)
rh_logistic$estado<-as.factor(rh_logistic$estado)

Se verifica la base de datos

str(rh_logistic)

## 'data.frame':    336 obs. of  7 variables:
##  $ edad          : int  32 36 23 21 29 46 29 31 50 19 ...
##  $ salario_diario: num  500 152 152 152 152 ...
##  $ estado        : Factor w/ 2 levels "Coahuila","Nuevo Leon": 2 2 2 2 2 2 2 2 2 2 ...
##  $ civil         : Factor w/ 4 levels "Divorcio","Matrimonio",..: 3 4 2 3 3 3 4 4 3 3 ...
##  $ genero        : Factor w/ 2 levels "FEMENINO","MASCULINO": 2 1 1 1 1 1 1 2 2 2 ...
##  $ puesto        : Factor w/ 43 levels "ANALISTA DE NOMINAS /AUX DE R.H.",..: 17 9 9 9 9 9 9 9 9 9 ...
##  $ baja          : int  1 1 1 1 1 1 1 1 1 1 ...

summary(rh_logistic)

##       edad       salario_diario         estado            civil    
##  Min.   :18.00   Min.   :144.4   Coahuila  : 17   Divorcio   :  6  
##  1st Qu.:23.00   1st Qu.:176.7   Nuevo Leon:319   Matrimonio :102  
##  Median :30.00   Median :180.7                    Soltero    :149  
##  Mean   :32.23   Mean   :178.4                    Union libre: 79  
##  3rd Qu.:39.00   3rd Qu.:180.7                                     
##  Max.   :67.00   Max.   :500.0                                     
##                                                                    
##        genero                 puesto         baja       
##  FEMENINO :195   AYUDANTE GENERAL:237   Min.   :0.0000  
##  MASCULINO:141   COSTURERA       : 20   1st Qu.:0.0000  
##                  SOLDADOR        : 15   Median :1.0000  
##                  MONTACARGUISTA  :  6   Mean   :0.6905  
##                  RESIDENTE       :  6   3rd Qu.:1.0000  
##                  CHOFER          :  4   Max.   :1.0000  
##                  (Other)         : 48

median(rh_logistic$salario_diario)

## [1] 180.68

# La mediana de salario es 180.7
# Existe un valor que se aleja: 500

Se reemplaza el valor que se sale de lo “normal” en salario diario con la mediana

rh_logistic$salario_diario<-replace(rh_logistic$salario_diario,rh_logistic$salario_diario>499,180.7)
rh_logistic$salario_diario[is.na(rh_logistic$salario_diario)]<-median(rh_logistic$salario_diario,na.rm=TRUE)

summary(rh_logistic)

##       edad       salario_diario         estado            civil    
##  Min.   :18.00   Min.   :144.4   Coahuila  : 17   Divorcio   :  6  
##  1st Qu.:23.00   1st Qu.:176.7   Nuevo Leon:319   Matrimonio :102  
##  Median :30.00   Median :180.7                    Soltero    :149  
##  Mean   :32.23   Mean   :177.4                    Union libre: 79  
##  3rd Qu.:39.00   3rd Qu.:180.7                                     
##  Max.   :67.00   Max.   :337.1                                     
##                                                                    
##        genero                 puesto         baja       
##  FEMENINO :195   AYUDANTE GENERAL:237   Min.   :0.0000  
##  MASCULINO:141   COSTURERA       : 20   1st Qu.:0.0000  
##                  SOLDADOR        : 15   Median :1.0000  
##                  MONTACARGUISTA  :  6   Mean   :0.6905  
##                  RESIDENTE       :  6   3rd Qu.:1.0000  
##                  CHOFER          :  4   Max.   :1.0000  
##                  (Other)         : 48

Se crea una categoría de referencia praa la variaable dependiente: “baja”

rh_logistic$baja<-as.factor(rh_logistic$baja)
rh_logistic$dv_baja<-fct_recode(rh_logistic$baja, "BAJA"="1","NO BAJA"="0")

Observamos una visualización de los datos para veer insights del dataset rh_logistic.

tapply(rh_logistic$salario_diario,
       list(rh_logistic$genero,rh_logistic$civil), mean)

##           Divorcio Matrimonio  Soltero Union libre
## FEMENINO    180.68   179.1338 175.3135    177.1121
## MASCULINO   178.70   181.1262 176.2136    177.3359

K-Means Clustering

Se exploran clusters relacionados con la edad en años.

rh_edad<-rh_logistic
rh_edad<-subset(rh_edad,select = -c(puesto,baja,estado,dv_baja))
summary(rh_edad)

##       edad       salario_diario          civil           genero   
##  Min.   :18.00   Min.   :144.4   Divorcio   :  6   FEMENINO :195  
##  1st Qu.:23.00   1st Qu.:176.7   Matrimonio :102   MASCULINO:141  
##  Median :30.00   Median :180.7   Soltero    :149                  
##  Mean   :32.23   Mean   :177.4   Union libre: 79                  
##  3rd Qu.:39.00   3rd Qu.:180.7                                    
##  Max.   :67.00   Max.   :337.1

Se normaliza

rh_edad_norm<-scale(rh_edad[1:2]) 

Vemos el número óptimo de los clusters

fviz_nbclust(rh_edad_norm, kmeans, method="wss")+ 
  geom_vline(xintercept=5, linetype=2)+           
  labs(subtitle = "Elbow method")  

Se tiene que el número óptimo de clusters es de 5

edad_cluster1<-kmeans(rh_edad_norm,5)

fviz_cluster(edad_cluster1,data=rh_edad_norm)

Se agregan los clusters a la base de datos original

rh_logistic3<-rh_logistic
rh_logistic3$Clusters<-edad_cluster1$cluster

rh_logistic4<-rh_logistic3 %>% group_by(Clusters) %>% summarise(edad=max(edad)) %>% arrange(desc(edad))

Se agrupan por nombre

rh_logistic3$Cluster_Names<-factor(rh_logistic3$Clusters,levels = c(1,2,3,4,5), 
                              labels=c("Mediana Edad","Joven-Mediana","Joven","Jubilación","Alta"))

Mediana el promedio de salario y edad

rh_logistic5<-rh_logistic3 %>% group_by(Cluster_Names) %>% dplyr::summarize(edad_años=mean(edad), 
                                                             salario_imss=mean(salario_diario),
                                                             Count=n())

use table format to display the clusters information

clusters<-as.data.frame(rh_logistic5)
clusters

##   Cluster_Names edad_años salario_imss Count
## 1  Mediana Edad  34.56897     188.3157    58
## 2 Joven-Mediana  27.62500     150.9000    40
## 3         Joven  24.08442     180.1727   154
## 4    Jubilación  53.76923     169.8323    26
## 5          Alta  45.03448     180.8395    58

Regresión logistica

#install.packages('DEoptimR') 
#install.packages("caret", dependencies = TRUE)

library(caret)

Se dividen los datos en training datasets y test datasets

set.seed(123) 
training<-rh_logistic3$Cluster_Names %>%     
  createDataPartition(p=0.75,list=FALSE)
train.data<-rh_logistic3[training, ]
test.data<-rh_logistic3[-training, ]         

Regresión logística simple

model specification

logit_model1<-glm(dv_baja~salario_diario, data=train.data, family=binomial(link='logit'))
summary(logit_model1)

## 
## Call:
## glm(formula = dv_baja ~ salario_diario, family = binomial(link = "logit"), 
##     data = train.data)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.7198  -1.4712   0.9097   0.9097   0.9268  
## 
## Coefficients:
##                 Estimate Std. Error z value Pr(>|z|)  
## (Intercept)     3.419576   1.486788   2.300   0.0214 *
## salario_diario -0.015226   0.008302  -1.834   0.0667 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 322.42  on 253  degrees of freedom
## Residual deviance: 318.51  on 252  degrees of freedom
## AIC: 322.51
## 
## Number of Fisher Scoring iterations: 4

Aunque no presenta mucha significancia, encontramos que el valor p de salario diario no es muy alto, de hecho es bastante cercano a 0.05. Por lo tanto, se considera que el modelo puede dar un buen nivel de predicción.

Predicciones a partir del modelo de regresión

# ¿A que  grupos o categorías hacen referencia estas probabilidades?
contrasts(test.data$dv_baja) 

##         BAJA
## NO BAJA    0
## BAJA       1

probabilities<-logit_model1 %>% predict(test.data,type="response")
predicted.classes<-ifelse(probabilities > 0.5, "bajas", "no_bajas")

# vemos las probabilidades del dataset relacionado con bajas y no bajas
head(probabilities) 

##         1         2        13        14        16        18 
## 0.6611056 0.7523437 0.6745473 0.6745473 0.7523437 0.7523437

# vemos las probabilidades del dataset relacionado con bajas
head(predicted.classes) 

##       1       2      13      14      16      18 
## "bajas" "bajas" "bajas" "bajas" "bajas" "bajas"

Vemos la presición del modelo 1 de regresión logistica

test.data$logit_model1_prob <- predict(logit_model1, test.data, type="response")

test.data <- test.data  %>% mutate(logit_model1_pred = 1*(logit_model1_prob > .50) + 0,
                         bajas_binary_1 = 1*(dv_baja == "BAJA") + 0)

test.data <- test.data %>% mutate(accurate_1=1*(logit_model1_pred == bajas_binary_1))
sum(test.data$accurate_1)/nrow(test.data) 

## [1] 0.7560976

A partir del R2 calculado, nos indica que aproximadamente el 76% de la variación en la variable dependiente es explicada por el logit model specification

El gráfico de la regresión logística demuestra una curva con forma de S

ggplot(rh_logistic3,aes(x=salario_diario, y=as.numeric(dv_baja)-1)) + 
  geom_point(alpha=.5) +
  stat_smooth(method="glm", se=FALSE, fullrange=TRUE, method.args=list(family=binomial)) + 
  ylab("Probabilidad") + xlim(100,500)+
  labs(
    title = "Modelo de Regresión Logística", 
    x = "Salario Diario",
    y = "Probabilidad de Bajas"
  )

## `geom_smooth()` using formula 'y ~ x'

Mediante la presente gráfica observamos que hay una probabilidad más alta de que las personas se den de baja cuando ganan entre $100 y $200 diarios, que aquellas que ganan más de $200. De otra manera, podríamos ver que: las personas que más ganan son las que tienen menos probabilidad de renunciar o ser dados de bajas, entre menos gane una persoona, más alta la probabilidad de que se salga de FORM (relación negativa).

Regresión logística múltiple

logit_model2<-glm(dv_baja~salario_diario+edad, data=train.data, family=binomial(link='logit'))
summary(logit_model2)

## 
## Call:
## glm(formula = dv_baja ~ salario_diario + edad, family = binomial(link = "logit"), 
##     data = train.data)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.9180  -1.2113   0.7491   0.8774   1.3140  
## 
## Coefficients:
##                 Estimate Std. Error z value Pr(>|z|)   
## (Intercept)     4.516971   1.514512   2.982  0.00286 **
## salario_diario -0.013742   0.008205  -1.675  0.09398 . 
## edad           -0.041222   0.012693  -3.248  0.00116 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 322.42  on 253  degrees of freedom
## Residual deviance: 307.67  on 251  degrees of freedom
## AIC: 313.67
## 
## Number of Fisher Scoring iterations: 4

En este modelo de logística, observamos como la edad tiene una mayor significancia y mayor impacto en la bajas que el salario diario. La edad tiene un valor de p de 0.0012, por lo que se considera que tiene un buen nivel de predicción y que puede ser usado en un modelo, como lo vemos a continuación.

Vemos la predicción del modelo 1 de regresión logística

test.data$logit_model2_prob <- predict(logit_model2, test.data, type="response")

test.data <- test.data  %>% mutate(logit_model2_pred = 1*(logit_model2_prob > .50) + 0,
                                   bajas_binary_2 = 1*(dv_baja == "BAJA") + 0)

test.data <- test.data %>% mutate(accurate_2=1*(logit_model2_pred == bajas_binary_2))
sum(test.data$accurate_2)/nrow(test.data) 

## [1] 0.7560976

A partir del R2 calculado, nos indica que aproximadamente el 76% de la variación en la variable dependiente es explicada por el logit model specification

El gráfico de la regresión logística demuestra una curva con forma de S

ggplot(rh_logistic3,aes(x=edad, y=as.numeric(dv_baja)-1)) + 
  geom_point(alpha=.5) +
  stat_smooth(method="glm", se=FALSE, fullrange=TRUE, method.args = list(family=binomial)) + 
  ylab("Probabilidad") + xlim(18,75)+
  labs(
    title = "Modelo de Regresión Logística", 
    x = "Edad",
    y = "Probabilidad de Bajas"
  )

## `geom_smooth()` using formula 'y ~ x'

Con el modelo de regresión múltiple, notando que la edad tiene mayor significancia, es la que se utiliza para verlo en gráfica. En la superior, vemos como existe una mayor probabilidad de que una persona salga de la empresa si es joven. Observamos que las personas después de los 50s, presentan una probabilidad menor a 50% de salir de la empresa.

Consideremos los resultados de K-Means Clustering

logit_model3<-glm(dv_baja~salario_diario+Cluster_Names, data=train.data, family=binomial(link='logit'))
summary(logit_model3)

## 
## Call:
## glm(formula = dv_baja ~ salario_diario + Cluster_Names, family = binomial(link = "logit"), 
##     data = train.data)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.7832  -1.3318   0.7971   0.7971   1.8558  
## 
## Coefficients:
##                            Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                10.88129    4.02143   2.706 0.006813 ** 
## salario_diario             -0.05387    0.02192  -2.457 0.013994 *  
## Cluster_NamesJoven-Mediana -2.05716    0.87721  -2.345 0.019021 *  
## Cluster_NamesJoven         -0.16490    0.42842  -0.385 0.700305    
## Cluster_NamesJubilación    -2.67380    0.72895  -3.668 0.000244 ***
## Cluster_NamesAlta          -0.47422    0.49071  -0.966 0.333847    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 322.42  on 253  degrees of freedom
## Residual deviance: 295.44  on 248  degrees of freedom
## AIC: 307.44
## 
## Number of Fisher Scoring iterations: 5

En este modelo, queremos ver que tanto impacto tiene el pertenecer a un cluster y el salaio diario en la probabilidad de que un empleado salga de laa empresa. Vemos una importance significancia, con valor de 0.0002 de p, en el cluster de Jubilación; las características de este cluster es que el promedio de edad es de 53 y el promedio de salario diario es de $169.

Igualmente, vemos como el cluster de Joven-Mediana tiene un significancia de un p value de 0.019, que aunque no es tan singificativa como el de Jubilación, resulta interesante. Las características que vemos es que tienen un promedio d e edad de 27 y un promedio de salario diario de $150 (el menor de todos los clusters).

Vemos, igualmente, que la variable que menos significancia tiene es el de los jóvenes, seguido por el de la edad “alta”.

test.data$logit_model3_prob <- predict(logit_model3, test.data, type="response")

test.data <- test.data  %>% mutate(logit_model3_pred = 1*(logit_model3_prob > .50) + 0,
                                   bajas_binary_3 = 1*(dv_baja == "BAJA") + 0)

test.data <- test.data %>% mutate(accurate_3=1*(logit_model3_pred == bajas_binary_3))
sum(test.data$accurate_3)/nrow(test.data) 

## [1] 0.7560976

A partir del R2 calculado, nos indica que aproximadamente el 76% de la variación en la varaiable dependiente es explicada por el logit model specification