Caso 11. Regresión Logística datos de ingresos adultos USA

Objetivo

Realizar e interpretar regresión logísitica con datos de personas e ingresos de USA

Descripción

Construir un modelo de regresión logísitca aplicado a datos de personas y sus ingresos en USA

La variable dependiente es los ingresos identificado por 0 y 1, los ganan por debajo o igual a 50 Mil y los que ganan por encima de 50 Mil.

Proceso

• 1. Cargar librerías
• 2. Cargar datos
• 3. Identificar variables
• 4. Crear datos de entrenamiento y validación
• 5. Crear modelo de regresión logística
• 6. Analizar y/o describir el modelo
• 7. Evaluar el modelo
• 8. Realizar predicciones con el conjunto de datos de validación
• 9. Interpretación

1. Cargar librerías

library(dplyr)
library(ggplot2)
library(knitr)
library(caret)
library(readr)
library(scales)
library(DT)

2. Cargar datos

datos <- read.csv("https://raw.githubusercontent.com/rpizarrog/FundamentosMachineLearning/master/datos/adultos_clean.csv", encoding = "UTF-8")

datatable(datos, caption = "Los datos", options = list(pageLength = 10)) 

3. Identificar variables

• age la edad de la persona
• workclass es un tipo o clase de trabajo de la persona, privado, gobierno, por su cuenta,
• education indica el nivel educativo de la persona
• educational es el valor numérico de education
• marital.status es su estado civil
• race es el tipo de raza de persona
• gender es el género de la persona
• hours.per.week son las horas que trbaja por semana
• income son los ingresos
• age_scale la edad escalada
• hours.per.week escalada
• income10 con valores de 0 gana menos de 50 mil y 1 mas de 50 mil

4. Crear datos de entrenamiento y validación

• 70 % datos de entrenamiento
• 30 % datos de validación

set.seed(2020)
entrena <- createDataPartition(y = datos$income10, p = 0.7, list = FALSE, times = 1)

# Datos entrenamiento
datos.entrenamiento <- datos[entrena, ]  # [renglones, columna]

# Datos validación
datos.validacion <- datos[-entrena, ]

datatable(datos.entrenamiento, caption = "Datos de entrenamiento", options = list(pageLength = 10))

datatable(datos.validacion, caption = "Datos de validación",options = list(pageLength = 10))

5. Crear modelo de regresión logística

Con la regresión logística, dado un conjunto particular de valores de las variables independientes elegidas, se estima la probabilidad de los ingresos de una persona ‘<=50’ o ‘>50’ o lo que es lo mismo ingresos con valores de 0 y 1 en la variable income10.

• Por medio de la función gml() se contruye un modelo de regresión logística
• Variable dependiente o predictiva es ‘income10’, ya que depende de todas las demás variables
• Variables independientes o predictoras, todas las demás: “age”, “workclass”, “education”, “educational.num”, “marital.status”, “race”, “gender”, “hours.per.week”
• Se utiliza el conjunto de datos de entrenamiento
• La finalidad de consruir el modelo de regresión logística es entre otros cosas, para conocer los coeficienes y el nivel de significación de cada variable independiente o predictora así como las pruebas t y F
• La ecuación \[income10=age.scale+workclass+education+marital.status+race+gender+hours.per.week.scale\] que se asigna a una variable y se utiliza para construir el modelo, significa que la variable ingresos * ‘income10’ depende o es dependiente de todas las demás variables
• La fórmula dada por la ecuacuación:

formula = income10 ~ age.scale + workclass + education + marital.status + race + gender + hours.per.week.scale
modelo <- glm(formula, data = datos.entrenamiento, family =  'binomial')

6. Analizar y/o describir el modelo

summary(modelo)

## 
## Call:
## glm(formula = formula, family = "binomial", data = datos.entrenamiento)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.7548  -0.5778  -0.2603  -0.0684   3.3249  
## 
## Coefficients:
##                           Estimate Std. Error z value Pr(>|z|)    
## (Intercept)               -2.76730    0.23517 -11.767  < 2e-16 ***
## age.scale                  2.04720    0.10474  19.546  < 2e-16 ***
## workclassFederal-gov       1.39346    0.12545  11.107  < 2e-16 ***
## workclassLocal-gov         0.76117    0.11093   6.862 6.79e-12 ***
## workclassNever-worked     -7.98118   98.02680  -0.081  0.93511    
## workclassPrivate           0.88906    0.09705   9.161  < 2e-16 ***
## workclassSelf-emp-inc      1.32433    0.11899  11.130  < 2e-16 ***
## workclassSelf-emp-not-inc  0.31448    0.10762   2.922  0.00348 ** 
## workclassState-gov         0.48978    0.12316   3.977 6.99e-05 ***
## workclassWithout-pay       0.15627    0.85256   0.183  0.85457    
## educationCommunity        -1.00153    0.04434 -22.589  < 2e-16 ***
## educationDropout          -2.69532    0.07543 -35.731  < 2e-16 ***
## educationHighGrad         -1.58616    0.04511 -35.164  < 2e-16 ***
## educationMaster            0.59823    0.06103   9.803  < 2e-16 ***
## educationPhD               1.22350    0.13962   8.763  < 2e-16 ***
## marital.statusNot_married -2.53858    0.05407 -46.948  < 2e-16 ***
## marital.statusSeparated   -2.10873    0.05619 -37.531  < 2e-16 ***
## marital.statusWidow       -2.10698    0.12348 -17.064  < 2e-16 ***
## raceAsian-Pac-Islander     0.32012    0.21854   1.465  0.14298    
## raceBlack                  0.34440    0.20941   1.645  0.10005    
## raceOther                  0.26768    0.28964   0.924  0.35540    
## raceWhite                  0.58693    0.20059   2.926  0.00343 ** 
## genderMale                 0.10151    0.04450   2.281  0.02252 *  
## hours.per.week.scale       3.05128    0.13834  22.056  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 37584  on 34189  degrees of freedom
## Residual deviance: 25090  on 34166  degrees of freedom
## AIC: 25138
## 
## Number of Fisher Scoring iterations: 11

• La mayoria de las variables tiene un impacto estadistico significativo dentro del modelo como lo podemos ver gracias a los ’***’

• Acerca del AIC información obtenida de https://support.minitab.com/es-mx/minitab/18/help-and-how-to/modeling-statistics/regression/how-to/fit-binary-logistic-model/interpret-the-results/all-statistics-and-graphs/model-summary-statistics/

  • Nos menciona que el AIC es una medida de la calidad relativa de un modelo que explica el ajuste y el número de términos en el modelo.

  • Mientras más pequeño sea el AIC, mejor se ajustará el modelo a los datos.

  • Sin embargo, el modelo con el criterio AIC más pequeño para un conjunto de predictores no necesariamente ajusta los datos adecuadamente.

  • Lo que quiere decir que en este caso el modelo para este conjunto de datos si se ajusta de muy buena manera, y nos puede servir bastante bien

7. Evaluar el modelo con matriz de confusión

Comparar los valores de income con valores ajustados

• Utilizando los datos de entrenamiento
• Tres variables
• income10 original con valores 0 y 1 s
• valores ajustados
• valores ajustados codificados 0 y 1 s aquellos cuya robabilidad sea > 0.5 o al 50%
• Con las columnas1 y 3 se puede generar la matriz de confusión

comparar <- data.frame(datos.entrenamiento$income10, as.vector(modelo$fitted.values) )

comparar <- comparar %>%
    mutate(income10ajustados = if_else (modelo$fitted.values > 0.5, 1, 0))

colnames(comparar) <- c("income10", "ajuste", 'income10ajustados')


datatable(comparar, caption = "Comparar valores", options = list(pageLength = 10))

Visualizando los valores ajustados

ggplot(data = comparar, aes(x =row.names(comparar), y = ajuste, )) + 
geom_point(aes(colour = factor(income10))) + 
  geom_hline(yintercept = 0.50)

Matriz de confusión

• Con los datos generados en la variable comparar se genera la matriz de confusión

matriz_confusion <- table(comparar$income10, comparar$income10ajustados, dnn = c("income10", "income10ajustados para predicciones"))
kable(matriz_confusion, caption = "Matriz de confusión")

Matriz de confusión

	0	1
0	24135	1892
1	4056	4107

Precisión del modelo

n = nrow(datos.entrenamiento)
exactidud <- (matriz_confusion[1,1] + matriz_confusion[2,2]) / n

exactidud

## [1] 0.826031

8. Realizar predicciones con el conjunto de datos de validación

• Información obtenida de : https://stat.ethz.ch/R-manual/R-devel/library/stats/html/predict.lm.html

• Nos menciona que se.fit es un indicador que muestra si es que los errores Standard requieren ser calculados

predicciones <- predict(modelo, datos.validacion, se.fit = TRUE)

predicciones.ajustadas <- predicciones$fit

kable(head(predicciones.ajustadas, 10))

	x
3	-0.7964486
4	-0.4623941
6	-5.0440275
12	1.0617147
19	-3.8320026
25	0.2488330
26	0.2057424
29	-0.6797358
41	2.2801743
42	0.3710815

Convertir las predicciones en probabilidades

predicciones_prob <- exp(predicciones$fit) / (1 + exp(predicciones$fit))

kable(head(predicciones_prob))

	x
3	0.3107857
4	0.3864180
6	0.0064064
12	0.7430181
19	0.0212067
25	0.5618892

Agregar una columna de la predicción al final del conjunto de datos de validación

• Agregar una columna de los datos de validación con los valores probabilísticos de las predicciones con cbind() y con un uevo conjunto de datos llamado las .predicciones
• Agregar columna con valor 1 cuando la predicción es mayor que 0.5 y 0 cuando la predicción es menor o igual a 0.5
• Verificar las columnas income10 e income10.prediccion

las.predicciones <- cbind(datos.validacion, predicciones_prob)

las.predicciones <- las.predicciones %>%
  mutate(income10.prediccion =  if_else(predicciones_prob > 0.5, 1, 0))
  
kable(head(las.predicciones))

X	age	workclass	education	educational.num	marital.status	race	gender	hours.per.week	income	age.scale	educational.num.scale	hours.per.week.scale	income10	predicciones_prob	income10.prediccion
3	28	Local-gov	Community	12	Married	White	Male	40	>50K	0.1506849	0.7333333	0.3979592	1	0.3107857	0
4	44	Private	Community	10	Married	Black	Male	40	>50K	0.3698630	0.6000000	0.3979592	1	0.3864180	0
6	34	Private	Dropout	6	Not_married	White	Male	30	<=50K	0.2328767	0.3333333	0.2959184	0	0.0064064	0
12	36	Federal-gov	Bachelors	13	Married	White	Male	40	<=50K	0.2602740	0.8000000	0.3979592	0	0.7430181	1
19	37	Private	HighGrad	9	Widow	White	Female	20	<=50K	0.2739726	0.5333333	0.1938776	0	0.0212067	0
25	25	Private	Bachelors	13	Married	White	Male	40	<=50K	0.1095890	0.8000000	0.3979592	0	0.5618892	1

Matriz de confusión de predicciones

matriz_confusion_predicciones <- table(las.predicciones$income10, las.predicciones$income10.prediccion, dnn = c("income10", "predicciones"))
matriz_confusion_predicciones

##         predicciones
## income10     0     1
##        0 10361   767
##        1  1752  1772

Precisión del modelo

n <- sum(matriz_confusion)
exactitud <- round(sum(matriz_confusion[1,1], matriz_confusion[2,2]) / n * 100,2)


paste("Las predicciones la atinan al ", exactitud, "%","de los casos")

## [1] "Las predicciones la atinan al  82.6 % de los casos"

Predicción con tres nuevos registros

Identificar los nuevos datos

• 3 nuevos registros

edad <- c(70,45,60)

clase.empleo <- c('State-gov','Never-worked','Self-emp-inc')
nivel.educacion <- c('HighGrad', 'Bachelors', 'Bachelors')
edo.civil <- c('Married', 'Separated', 'Widow')
raza <- c('White', 'Asian-Pac-Islander', 'Black')
genero <- c('Female', 'Male', 'Female')

horas <- c(58,60,65)

a.predecir <- data.frame(rbind(cbind(edad, clase.empleo, nivel.educacion, edo.civil,raza, genero, horas)))

colnames(a.predecir) <- c('age.scale', 'workclass', 'education', 'marital.status', 'race', 'gender', 'hours.per.week.scale')

a.predecir

##   age.scale    workclass education marital.status               race gender
## 1        70    State-gov  HighGrad        Married              White Female
## 2        45 Never-worked Bachelors      Separated Asian-Pac-Islander   Male
## 3        60 Self-emp-inc Bachelors          Widow              Black Female
##   hours.per.week.scale
## 1                   58
## 2                   60
## 3                   65

Escalar los valores númericos

• Escalar significa centrar conforme con los valores mínimos y máximo de datos originales
• Escalar el valor numérico de la edad igualando y con todos los valores de la edad de los datos originale
• Escalar el valor de horas por semana ‘hours.per.week.scale’ de acuerdo a la columna hours.per.week de todos los datos originales
• Modificar mutate() las columnas de age y hours.per.week
• Se toman para este caso el primero registro [1] de los valores escalados: edad.escalada[1] y horas.escalada[1]

edad.escalada <- rescale(c(edad, min(datos$age), max(datos$age)))
edad.escalada <- edad.escalada[1:3]

# Escalando las horas por semana
horas.escalada<- rescale(c(horas, min(datos$hours.per.week), max(datos$hours.per.week)))
horas.escalada <- horas.escalada[1:3]

a.predecir <- a.predecir %>%
  mutate(age.scale = edad.escalada,
         hours.per.week.scale = horas.escalada)
 
a.predecir

##   age.scale    workclass education marital.status               race gender
## 1 0.7260274    State-gov  HighGrad        Married              White Female
## 2 0.3835616 Never-worked Bachelors      Separated Asian-Pac-Islander   Male
## 3 0.5890411 Self-emp-inc Bachelors          Widow              Black Female
##   hours.per.week.scale
## 1            0.5816327
## 2            0.6020408
## 3            0.6530612

Realizar las predicciones con los nuevos registros

predicciones.nuevos <- predict(modelo, a.predecir, se.fit = TRUE)
predicciones.nuevos

## $fit
##            1            2            3 
## -0.015713200 -9.813362419 -0.007000296 
## 
## $se.fit
##          1          2          3 
##  0.1066299 98.0268280  0.1568007 
## 
## $residual.scale
## [1] 1

prediccion_prob.nuevos <- exp(predicciones.nuevos$fit) / (1 + exp(predicciones.nuevos$fit))
prediccion_prob.nuevos

##            1            2            3 
## 4.960718e-01 5.471257e-05 4.982499e-01

las.predicciones <- if_else(prediccion_prob.nuevos > 0.5, 1, 0)

cat("Son las predicciones para las personas con esas características ")

## Son las predicciones para las personas con esas características

las.predicciones

## [1] 0 0 0

print("Las predicciones en la columna final")

## [1] "Las predicciones en la columna final"

cbind(a.predecir ,las.predicciones)

##   age.scale    workclass education marital.status               race gender
## 1 0.7260274    State-gov  HighGrad        Married              White Female
## 2 0.3835616 Never-worked Bachelors      Separated Asian-Pac-Islander   Male
## 3 0.5890411 Self-emp-inc Bachelors          Widow              Black Female
##   hours.per.week.scale las.predicciones
## 1            0.5816327                0
## 2            0.6020408                0
## 3            0.6530612                0

9. Interpretación

En este conjunto de datos que habla sobre los ingresos de personas que viven en USA se desea saber como es que se determina si una persona recibe más de 50000 dólares o menos o igual a 50000 dólares de ingresos. Para eso la variable objetivo va a ser el income que esta dado como ‘>50k’ o ‘<=50k’ siendo las variables independientes las siguientes:

• age la edad de la persona

• workclass es un tipo o clase de trabajo de la persona, privado, gobierno, por su cuenta

• education indica el nivel educativo de la persona

• educational es el valor numérico de education

• marital.status es su estado civil

• race es el tipo de raza de persona

• gender es el género de la persona

• hours.per.week son las horas que trbaja por semana

• income son los ingresos

• age_scale la edad escalada

• hours.per.week escalada

Para hacer este análisis se decidió usar el modelo de regresión logística en donde se obtuvieron los resultados que dicen que la mayoría de las variables tienen un significado estadístico importante y el valor AIC llega a ser bastante bajo con un valor de 25138 comparado con la cantidad de datos, esto quiere decir que el conjunto de datos se adapta muy bien al modelo. En cuanto a la matriz de confusión nos podemos dar cuenta que:

• El modelo dice que 24135 ganan por debajo de 50 Mil de los 34190 registros totales. De tal forma que le atinó a 70.59 a los VP verdaderos positivos

• El modelo dice que 4107 ganan por encima de 50 Mil de los 34190 registros totales. De tal forma que le atinó a 12.01 a los FP verdaderos negativos

Y en cuanto a la precisión del modelo nos podemos dar cuenta que el modelo predice muy bien los casos en donde se obtiene un ingreso de menos de 50000 dólares y en general el modelo tiene una precisión del 82% con un error del 18% y en cuanto a las predicciones igualmente tiene una certeza del 82% con un error del 18%, lo que quiere decir que el modelo es bastante adecuado para este caso y predice en la mayoría de los casos para recibir un ingreso de menos de 50000 muy bien, aún necesita refinamiento con talvez otro modelo.