Estudio de clasificación según los ingresos

0. - Parametros

Desactivamos la notación científica para facilitar el visualizado de los numeros.

options(scipen=999)#Desactiva la notacion cientifica

#por defecto no muestro ni warnings ni mensajes de la consola en la salida
knitr::opts_chunk$set(warning = FALSE, message = FALSE)

1. - Preparacion del entorno

1.1 - Cargamos las librerias que vamos a utilizar

Comrpobamos si todos los paquetes que necesitamos están instalados en el entorno de R, y si no lo están los instalamos. Después los cargamos todos.

1.2 - Cargamos los datos

Vamos a utilizar unos datos obtenidos de UCI Machine Learning Repository. En ellos se tienen datos demograficos y personales, con los que predecir quienes tienen ingresos mayores a 50.000 dolares o menores a 50.000 dolares.

Usamos fread de data.table para una lectura mucho mas rapida del fichero con los datos. Con esta instrucción tendremos los datos en un dataframe para su manejo y gestión.

Este set de datos ya viene en dos partes, una para training y otra para test. Cargamos ambas y las juntaremos, después borramos las variables temporales. Más adelante ya separaremos entre training y test nosotros.

temp1 <- fread('adult_data.csv')
temp <- fread('adult_test.csv')
df <- rbind(temp1,temp)
rm(temp)
rm(temp1)
knitr::kable(head(df))
age workclass fnlwgt education educacion-num marital-status occupation relationship race sex capital-gain capital-loss hours-per-week native-country income
39 State-gov 77516 Bachelors 13 Never-married Adm-clerical Not-in-family White Male 2174 0 40 United-States <=50K
50 Self-emp-not-inc 83311 Bachelors 13 Married-civ-spouse Exec-managerial Husband White Male 0 0 13 United-States <=50K
38 Private 215646 HS-grad 9 Divorced Handlers-cleaners Not-in-family White Male 0 0 40 United-States <=50K
53 Private 234721 11th 7 Married-civ-spouse Handlers-cleaners Husband Black Male 0 0 40 United-States <=50K
28 Private 338409 Bachelors 13 Married-civ-spouse Prof-specialty Wife Black Female 0 0 40 Cuba <=50K
37 Private 284582 Masters 14 Married-civ-spouse Exec-managerial Wife White Female 0 0 40 United-States <=50K

1.3. - Muestreo

No es necesario el muestreo, ya que tenemos una cantidad de registros manejable.

2 - Analisis exploratorio

2.1 - Analisis exploratorio general y tipo de datos

as.data.frame(sort(names(df)))
##    sort(names(df))
## 1              age
## 2     capital-gain
## 3     capital-loss
## 4    educacion-num
## 5        education
## 6           fnlwgt
## 7   hours-per-week
## 8           income
## 9   marital-status
## 10  native-country
## 11      occupation
## 12            race
## 13    relationship
## 14             sex
## 15       workclass
str(df)
## Classes 'data.table' and 'data.frame':   48842 obs. of  15 variables:
##  $ age           : int  39 50 38 53 28 37 49 52 31 42 ...
##  $ workclass     : chr  "State-gov" "Self-emp-not-inc" "Private" "Private" ...
##  $ fnlwgt        : int  77516 83311 215646 234721 338409 284582 160187 209642 45781 159449 ...
##  $ education     : chr  "Bachelors" "Bachelors" "HS-grad" "11th" ...
##  $ educacion-num : int  13 13 9 7 13 14 5 9 14 13 ...
##  $ marital-status: chr  "Never-married" "Married-civ-spouse" "Divorced" "Married-civ-spouse" ...
##  $ occupation    : chr  "Adm-clerical" "Exec-managerial" "Handlers-cleaners" "Handlers-cleaners" ...
##  $ relationship  : chr  "Not-in-family" "Husband" "Not-in-family" "Husband" ...
##  $ race          : chr  "White" "White" "White" "Black" ...
##  $ sex           : chr  "Male" "Male" "Male" "Male" ...
##  $ capital-gain  : int  2174 0 0 0 0 0 0 0 14084 5178 ...
##  $ capital-loss  : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ hours-per-week: int  40 13 40 40 40 40 16 45 50 40 ...
##  $ native-country: chr  "United-States" "United-States" "United-States" "United-States" ...
##  $ income        : chr  "<=50K" "<=50K" "<=50K" "<=50K" ...
##  - attr(*, ".internal.selfref")=<externalptr>
glimpse(df)
## Observations: 48,842
## Variables: 15
## $ age              <int> 39, 50, 38, 53, 28, 37, 49, 52, 31, 42, 37, 3...
## $ workclass        <chr> "State-gov", "Self-emp-not-inc", "Private", "...
## $ fnlwgt           <int> 77516, 83311, 215646, 234721, 338409, 284582,...
## $ education        <chr> "Bachelors", "Bachelors", "HS-grad", "11th", ...
## $ `educacion-num`  <int> 13, 13, 9, 7, 13, 14, 5, 9, 14, 13, 10, 13, 1...
## $ `marital-status` <chr> "Never-married", "Married-civ-spouse", "Divor...
## $ occupation       <chr> "Adm-clerical", "Exec-managerial", "Handlers-...
## $ relationship     <chr> "Not-in-family", "Husband", "Not-in-family", ...
## $ race             <chr> "White", "White", "White", "Black", "Black", ...
## $ sex              <chr> "Male", "Male", "Male", "Male", "Female", "Fe...
## $ `capital-gain`   <int> 2174, 0, 0, 0, 0, 0, 0, 0, 14084, 5178, 0, 0,...
## $ `capital-loss`   <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...
## $ `hours-per-week` <int> 40, 13, 40, 40, 40, 40, 16, 45, 50, 40, 80, 4...
## $ `native-country` <chr> "United-States", "United-States", "United-Sta...
## $ income           <chr> "<=50K", "<=50K", "<=50K", "<=50K", "<=50K", ...

Variables que tenemos en el data frame:

  • income: variable categorica. >50K, <=50K. candidata de pasar a factor y a 0/1 para facilitar su manejo.
  • age: variable integer. continuous.
  • workclass: variable categorica. Private, Self-emp-not-inc, Self-emp-inc, Federal-gov, Local-gov, State-gov, Without-pay, Never-worked. candidata de pasar a factor.
  • fnlwgt: integer. continuous.
  • education: variable categorica. Categorías: Bachelors, Some-college, 11th, HS-grad, Prof-school, Assoc-acdm, Assoc-voc, 9th, 7th-8th, 12th, Masters, 1st-4th, 10th, Doctorate, 5th-6th, Preschool. Candidata de pasar a factor.
  • education-num: variable entera. continuous.
  • marital-status: variable categorica. Categorías: Married-civ-spouse, Divorced, Never-married, Separated, Widowed, Married-spouse-absent, Married-AF-spouse. Candidata a pasar a factor.
  • occupation: Variable categorica. Categorías: Tech-support, Craft-repair, Other-service, Sales, Exec-managerial, Prof-specialty, Handlers-cleaners, Machine-op-inspct, Adm-clerical, Farming-fishing, Transport-moving, Priv-house-serv, Protective-serv, Armed-Forces. Candidata de pasar a factor.
  • relationship: variable categorica. Wife, Own-child, Husband, Not-in-family, Other-relative, Unmarried. Candidata de pasar a factor.
  • race: variable categorica. Categorías: White, Asian-Pac-Islander, Amer-Indian-Eskimo, Other, Black. Candidata de pasar a factor.
  • sex: variable categorica. Categorías: Female, Male. Candidata de pasar a factor.
  • capital-gain: variable integer. continuous. Aparecen algunos ceros.
  • capital-loss: variable integer. continuous. Aparecen algunos ceros.
  • hours-per-week: variable integer. Horas de trabajo a la semana. continuous.
  • native-country: variable categorica. United-States, Cambodia, England, Puerto-Rico, Canada, Germany, Outlying-US(Guam-USVI-etc), India, Japan, Greece, South, China, Cuba, Iran, Honduras, Philippines, Italy, Poland, Jamaica, Vietnam, Mexico, Portugal, Ireland, France, Dominican-Republic, Laos, Ecuador, Taiwan, Haiti, Columbia, Hungary, Guatemala, Nicaragua, Scotland, Thailand, Yugoslavia, El-Salvador, Trinadad&Tobago, Peru, Hong, Holand-Netherlands. Candidata de pasar a factor.

Hacemos el primer análisis (str, glimpse). Conclusiones:

  • Hay bastantes variables categoricas, que debemos pasar a factores.
  • Se han de investigar los ceros de capital-gain y capital-loss
  • Deberiamos normalizar la variable income como <=50K (0) y >50K (1).
  • En los primeros valores de capital-loss solo aparecen ceros. Podemos ver que valores distintos tiene esta variable y también capital-gain.
Capitalgain <- table(df$'capital-gain')
Capitalloss <- table(df$'capital-loss')
print(Capitalgain)
## 
##     0   114   401   594   914   991  1055  1086  1111  1151  1173  1264 
## 44807     8     5    52    10     6    37     8     1    13     5     2 
##  1409  1424  1455  1471  1506  1639  1731  1797  1831  1848  2009  2036 
##    10     4     4     9    24     1     1    10     9     9     3     5 
##  2050  2062  2105  2174  2176  2202  2228  2290  2329  2346  2354  2387 
##     5     3    15    74    31    28     5    10     7     8    21     1 
##  2407  2414  2463  2538  2580  2597  2635  2653  2829  2885  2907  2936 
##    25    10    15     5    20    31    14    11    42    30    18     4 
##  2961  2964  2977  2993  3103  3137  3273  3325  3411  3418  3432  3456 
##     4    14    11     3   152    51     7    81    34     8     4     6 
##  3464  3471  3674  3781  3818  3887  3908  3942  4064  4101  4386  4416 
##    33    11    22    16    11     8    42    18    54    29   108    24 
##  4508  4650  4687  4787  4865  4931  4934  5013  5060  5178  5455  5556 
##    23    63     4    35    25     4    10   117     2   146    18     6 
##  5721  6097  6360  6418  6497  6514  6612  6723  6767  6849  7262  7298 
##     7     2     3    16    15    10     1     5     6    42     1   364 
##  7430  7443  7688  7896  7978  8614  9386  9562 10520 10566 10605 11678 
##    15     7   410     4     2    82    31     5    64     8    19     4 
## 13550 14084 14344 15020 15024 15831 18481 20051 22040 25124 25236 27828 
##    42    49    34    10   513     8     2    49     1     6    14    58 
## 34095 41310 99999 
##     6     3   244
print(Capitalloss)
## 
##     0   155   213   323   419   625   653   810   880   974  1092  1138 
## 46560     1     5     5     3    17     4     2     6     2    11     4 
##  1258  1340  1380  1408  1411  1421  1429  1485  1504  1510  1539  1564 
##     6    11    10    35     4     1     3    71    26     3     1    43 
##  1573  1579  1590  1594  1602  1617  1628  1648  1651  1668  1669  1672 
##    12    30    62     9    62    11    24     3    11     9    35    50 
##  1719  1721  1726  1735  1740  1741  1755  1762  1816  1825  1844  1848 
##    38    28     9     3    58    44     2    20     4     5     3    67 
##  1870  1876  1887  1902  1911  1944  1974  1977  1980  2001  2002  2042 
##     1    59   233   304     1     3    28   253    36    35    33    12 
##  2051  2057  2080  2129  2149  2163  2174  2179  2201  2205  2206  2231 
##    29    16     1     7     5     2    10    20     1    19     6     7 
##  2238  2246  2258  2267  2282  2339  2352  2377  2392  2415  2444  2457 
##     4     8    39     3     2    27     2    25    11    72    20     4 
##  2465  2467  2472  2489  2547  2559  2603  2754  2824  3004  3175  3683 
##     1     2     4     1     5    17     7     2    14     5     2     2 
##  3770  3900  4356 
##     4     2     3

Vemos que la mayoría de los registros de estas dos variables son cero. Hay que investigarlos más a fondo.

Reservamos las variables a pasar a factores en el vector a_factores:

a_factores <- c("income","workclass","education","marital-status","occupation",
                "relationship","race","sex","native-country")

2.2 - Calidad de datos: Estadísticos básicos

Hacemos un summary, con lapply que sale en formato de lista y se lee mejor. También lo guardo en una variable y lo muestro en pantalla para mejorar su lectura cuando renderizamos el notebook a HTML, y además lo podemos acceder individualmente como vector.

est_basicos <- lapply(df,summary)
print(est_basicos)
## $age
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   17.00   28.00   37.00   38.64   48.00   90.00 
## 
## $workclass
##    Length     Class      Mode 
##     48842 character character 
## 
## $fnlwgt
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   12285  117551  178145  189664  237642 1490400 
## 
## $education
##    Length     Class      Mode 
##     48842 character character 
## 
## $`educacion-num`
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    1.00    9.00   10.00   10.08   12.00   16.00 
## 
## $`marital-status`
##    Length     Class      Mode 
##     48842 character character 
## 
## $occupation
##    Length     Class      Mode 
##     48842 character character 
## 
## $relationship
##    Length     Class      Mode 
##     48842 character character 
## 
## $race
##    Length     Class      Mode 
##     48842 character character 
## 
## $sex
##    Length     Class      Mode 
##     48842 character character 
## 
## $`capital-gain`
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##       0       0       0    1079       0   99999 
## 
## $`capital-loss`
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##     0.0     0.0     0.0    87.5     0.0  4356.0 
## 
## $`hours-per-week`
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    1.00   40.00   40.00   40.42   45.00   99.00 
## 
## $`native-country`
##    Length     Class      Mode 
##     48842 character character 
## 
## $income
##    Length     Class      Mode 
##     48842 character character

Conclusiones: - La edad (age)va de 17 a 90 años, aunque con la media en 38 y el 3rd Qu en 48, parece que habrá valores atípicos por la parte de arriba. Esta variable la discretizaremos. - La variable fnlwgt tiene un rango muy amplio, con una gran diferenciaentre la media y el máximo. Habrá que hacer una análisis de atípicos para ver como es la distribución de la variable, y considerar su discretización. - La variable education-num va entre 1 y 16, con media y mediana en 10. Es posible que haya valores atípicos por la parte de abajo. Analizaremos la distribución de la variable, para considerar su discretización. - Las variables capital-gain y capital-loss tienen más de la mitad de los valores a cero, con una media muy baja y un máximo muy alto. Hay que hacer un análisis de atípicos y ceros, para ver que decisión se toma con estas variables. - La variable hours_per_week tiene una media de 40, con un 1st Qu de 40 y un 3rd Qu de 45, con lo que está muy agrupada alrededor de la media. El maximo es 99, por lo que parece que hay atípicos en la parte alta. Hay que hacer un análisis de atipicos y considerar la discretización.

2.3 - Calidad de datos: Análisis de nulos

Lo pongo en una variable para tener más facilidad su inclusión en un informe automático.

num_nulos <- data.frame(colSums(is.na(df)))
print(num_nulos)
##                colSums.is.na.df..
## age                             0
## workclass                       0
## fnlwgt                          0
## education                       0
## educacion-num                   0
## marital-status                  0
## occupation                      0
## relationship                    0
## race                            0
## sex                             0
## capital-gain                    0
## capital-loss                    0
## hours-per-week                  0
## native-country                  0
## income                          0

Conclusion: no hay valores nulos en ninguna de las variables.

2.4 - Calidad de datos: Análisis de ceros

Hemos detectado muchos ceros en las variables capital-gain y capital-loss, vamos a hacer un analisis de ceros en esas dos variables.

contar_ceros <- function(variable) {
    temp <- transmute(df,if_else(variable==0,1,0))
    sum(temp)
}

num_ceros <- sapply(df,contar_ceros) #aplico la función de contar ceros a cada variable de nuestro                                            conjunto de datos
#transformamos el vector en un data frame para mejorar su visualización
num_ceros <- data.frame(VARIABLE=names(num_ceros),CEROS = as.numeric(num_ceros),stringsAsFactors = F) #el as.numeric es para sacar solo el valor de num_ceros, sin el nombre
num_ceros <- num_ceros %>%
  arrange(desc(CEROS)) %>%
  mutate(PORCENTAJE = CEROS / nrow(df) * 100) #añado a los valores absolutos el porcentajes del total de                                                 observaciones
num_ceros
##          VARIABLE CEROS PORCENTAJE
## 1    capital-loss 46560   95.32779
## 2    capital-gain 44807   91.73867
## 3             age     0    0.00000
## 4       workclass     0    0.00000
## 5          fnlwgt     0    0.00000
## 6       education     0    0.00000
## 7   educacion-num     0    0.00000
## 8  marital-status     0    0.00000
## 9      occupation     0    0.00000
## 10   relationship     0    0.00000
## 11           race     0    0.00000
## 12            sex     0    0.00000
## 13 hours-per-week     0    0.00000
## 14 native-country     0    0.00000
## 15         income     0    0.00000

Conclusion: Estas dos variables, capital-gain y capital-loss contienen datos de ganancias y pérdidas por productos de capital. PArece ser que gran parte de la muestra no tiene productos de capital, como acciones, y por eso hay un valor elevado de ceros. Pero puede ser un punto importante para precedir nivel de ingresos. Como no lo podemos saber, incluiremos ambas variables, pero quizás haya que combinarlas de alguna manera para sintetizar mejor la información para la predicción.

2.5 - Calidad de datos - atipicos

Dividimos el análisis entre las variables tipo numerico y tipo integer. Aquí solo tenemos variables tipo integer, así que solo hacemos el análisis integer.

2.5.1 - Analizamos las que son de tipo numerico

out <- function(variable){
  t(t(head(sort(variable,decreasing = T),20))) #la doble traspuesta es un truco para que se visualice la salida, si no lo que crearia es una coleccion de dataframes que no se ven bien
}
atipicos_variables_numeric <- lapply(df,function(x){
  if(is.double(x)) out(x)
})
print(atipicos_variables_numeric)
## $age
## NULL
## 
## $workclass
## NULL
## 
## $fnlwgt
## NULL
## 
## $education
## NULL
## 
## $`educacion-num`
## NULL
## 
## $`marital-status`
## NULL
## 
## $occupation
## NULL
## 
## $relationship
## NULL
## 
## $race
## NULL
## 
## $sex
## NULL
## 
## $`capital-gain`
## NULL
## 
## $`capital-loss`
## NULL
## 
## $`hours-per-week`
## NULL
## 
## $`native-country`
## NULL
## 
## $income
## NULL

2.5.2 - Analizamos las que son de tipo integer

out_int <- function(variable){
  t(t(table(variable))) #la doble traspuesta es un truco para que se visualice la salida, si no lo que crearia es una coleccion de dataframes que no se ven bien
}
atipicos_variables_integer <- lapply(df,function(x){
  if(is.integer(x)) out_int(x)
})
#hacemos algo más facil la visualización mostrando solo algunos datos de cada variable
for(i in 1:length(atipicos_variables_integer)){
  print(names(atipicos_variables_integer[i]))
  print(head(atipicos_variables_integer[[i]]))
  cat("\n","----------------","\n","\n")
}
## [1] "age"
##         
## variable [,1]
##       17  595
##       18  862
##       19 1053
##       20 1113
##       21 1096
##       22 1178
## 
##  ---------------- 
##  
## [1] "workclass"
## NULL
## 
##  ---------------- 
##  
## [1] "fnlwgt"
##         
## variable [,1]
##    12285    1
##    13492    1
##    13769    3
##    13862    1
##    14878    1
##    18827    1
## 
##  ---------------- 
##  
## [1] "education"
## NULL
## 
##  ---------------- 
##  
## [1] "educacion-num"
##         
## variable [,1]
##        1   83
##        2  247
##        3  509
##        4  955
##        5  756
##        6 1389
## 
##  ---------------- 
##  
## [1] "marital-status"
## NULL
## 
##  ---------------- 
##  
## [1] "occupation"
## NULL
## 
##  ---------------- 
##  
## [1] "relationship"
## NULL
## 
##  ---------------- 
##  
## [1] "race"
## NULL
## 
##  ---------------- 
##  
## [1] "sex"
## NULL
## 
##  ---------------- 
##  
## [1] "capital-gain"
##         
## variable [,1] 
##      0   44807
##      114     8
##      401     5
##      594    52
##      914    10
##      991     6
## 
##  ---------------- 
##  
## [1] "capital-loss"
##         
## variable [,1] 
##      0   46560
##      155     1
##      213     5
##      323     5
##      419     3
##      625    17
## 
##  ---------------- 
##  
## [1] "hours-per-week"
##         
## variable [,1]
##        1   27
##        2   53
##        3   59
##        4   84
##        5   95
##        6   92
## 
##  ---------------- 
##  
## [1] "native-country"
## NULL
## 
##  ---------------- 
##  
## [1] "income"
## NULL
## 
##  ---------------- 
##

Conclusiones:

  • Vamos a discretizar la edad, el fnlwgt, hours_per_week.
  • En la variable education_num, el valor 1 es un valor atípico. Tiene muy pocos valores para discretizar, asi que la dejaremos tal y como está, eliminando los registros con valor 1.
  • Las variables capital-gain y capital-loss las combinaremos para calcular las ganancias netas. Esta nueva variable la discretizaremos.

2.6 - Analisis longitudinal

Analisis de coherencia de las variables dentro de si mismas.

longi <- df %>% 
  summarise_all(mean) %>% #calcular la media de cada variable
  t() %>% #transponerlo para tenerlo en una sola columna y leerlo mejor
  as.data.frame() #reconvertirlo a dataframe porque t() lo pasa a matriz
data.frame(variable = rownames(longi), media = longi$V1) %>% #crear un nuevo dataframe para poder ordenar por el nombre
  arrange(desc(variable)) #ordenar por el nombre para tener la vision longitudinal
##          variable        media
## 1       workclass           NA
## 2             sex           NA
## 3    relationship           NA
## 4            race           NA
## 5      occupation           NA
## 6  native-country           NA
## 7  marital-status           NA
## 8          income           NA
## 9  hours-per-week     40.42238
## 10         fnlwgt 189664.13460
## 11      education           NA
## 12  educacion-num     10.07809
## 13   capital-loss     87.50231
## 14   capital-gain   1079.06763
## 15            age     38.64359

Conclusiones: no se ven anomalias en las medias, a simple vista. La media de jornada laboral es 40h/semana, la edad media son 38 años. En cuanto a ganancias y salario, la media puede ser normal, aunque lo miraremos más en detalle.

Para ver si es por unos pocos atipicos sumamos los 10 valores más altos de cada variable

df %>% 
  select(fnlwgt,'capital-gain','capital-loss') %>% 
  map(~ sum(head(sort(.x,decreasing = T),10)))
## $fnlwgt
## [1] 12973684
## 
## $`capital-gain`
## [1] 999990
## 
## $`capital-loss`
## [1] 39631

Conclusion: efectivamente es debido a algunos atipicos, como vamos a discretizar no nos preocupa y podemos seguir adelante

2.7 - Analisis de coherencia

En este caso no tenemos variables que podamos contrastar unas contra otras para ver si presentan incoherencias.

2.8 Otros

2.8.1 Analisis de la fnlwgt, capital-gain y capital-loss

Hemos visto que en la variables fnlwgt, capital-gain, capital-loss hay bastantes atipicos. Los analizaremos con un histograma.

hist(df$fnlwgt,breaks = 50)

hist(df$`capital-gain`,breaks = 50)

hist(df$`capital-loss`, breaks = 50)

La distribución de fnlwgt es bastante normal, disminuyendo hacia valores altos. Es un comportamiento normal, y la variable no tiene valores atipicos de gran frecuencia. La variable capital-gain indica que hay muchas personas que acumulan ganancias pequeñas, y que las ganancias grandes, aunque existen, tienen una muy baja frecuencia. Lo mismo pasa con el histograma de capital-loss. Ambos son comportamientos normales, y por tanto, seguiremos con el plan de combinar ambas variables y discretizar.

2.9 - Acciones resultado del analisis de calidad de datos y exploratorio

Aquí ejecutamos todas las acciones que hemos ido anotando del análisis de la calidad de los datos.

Vamos a hacer lo siguiente: - eliminar los registros con valores 1 en la variable education-num. - transformar a factor las variables de ‘a_factores’ - debemos cambiar la variable income, donde <=50K (0) y >50K (1).

df <- filter(df,df$`educacion-num`>1)

df <- df %>%
  mutate(income = ifelse(income == "<=50K" | income =="<=50K.",0,1)) %>%
  mutate_at(a_factores,funs(factor))
print(head(df))
##   age        workclass fnlwgt education educacion-num     marital-status
## 1  39        State-gov  77516 Bachelors            13      Never-married
## 2  50 Self-emp-not-inc  83311 Bachelors            13 Married-civ-spouse
## 3  38          Private 215646   HS-grad             9           Divorced
## 4  53          Private 234721      11th             7 Married-civ-spouse
## 5  28          Private 338409 Bachelors            13 Married-civ-spouse
## 6  37          Private 284582   Masters            14 Married-civ-spouse
##          occupation  relationship  race    sex capital-gain capital-loss
## 1      Adm-clerical Not-in-family White   Male         2174            0
## 2   Exec-managerial       Husband White   Male            0            0
## 3 Handlers-cleaners Not-in-family White   Male            0            0
## 4 Handlers-cleaners       Husband Black   Male            0            0
## 5    Prof-specialty          Wife Black Female            0            0
## 6   Exec-managerial          Wife White Female            0            0
##   hours-per-week native-country income
## 1             40  United-States      0
## 2             13  United-States      0
## 3             40  United-States      0
## 4             40  United-States      0
## 5             40           Cuba      0
## 6             40  United-States      0

Ahora que ya tenemos los factores, vemos que niveles tenemos en cada factor, para asegurar que no haya niveles raros en los factores.

levels(df$workclass)
## [1] "?"                "Federal-gov"      "Local-gov"       
## [4] "Never-worked"     "Private"          "Self-emp-inc"    
## [7] "Self-emp-not-inc" "State-gov"        "Without-pay"
levels(df$education)
##  [1] "10th"         "11th"         "12th"         "1st-4th"     
##  [5] "5th-6th"      "7th-8th"      "9th"          "Assoc-acdm"  
##  [9] "Assoc-voc"    "Bachelors"    "Doctorate"    "HS-grad"     
## [13] "Masters"      "Prof-school"  "Some-college"
levels(df$`marital-status`)
## [1] "Divorced"              "Married-AF-spouse"     "Married-civ-spouse"   
## [4] "Married-spouse-absent" "Never-married"         "Separated"            
## [7] "Widowed"
levels(df$occupation)
##  [1] "?"                 "Adm-clerical"      "Armed-Forces"     
##  [4] "Craft-repair"      "Exec-managerial"   "Farming-fishing"  
##  [7] "Handlers-cleaners" "Machine-op-inspct" "Other-service"    
## [10] "Priv-house-serv"   "Prof-specialty"    "Protective-serv"  
## [13] "Sales"             "Tech-support"      "Transport-moving"
levels(df$relationship)
## [1] "Husband"        "Not-in-family"  "Other-relative" "Own-child"     
## [5] "Unmarried"      "Wife"
levels(df$race)
## [1] "Amer-Indian-Eskimo" "Asian-Pac-Islander" "Black"             
## [4] "Other"              "White"
levels(df$sex)
## [1] "Female" "Male"
levels(df$`native-country`)
##  [1] "?"                          "Cambodia"                  
##  [3] "Canada"                     "China"                     
##  [5] "Columbia"                   "Cuba"                      
##  [7] "Dominican-Republic"         "Ecuador"                   
##  [9] "El-Salvador"                "England"                   
## [11] "France"                     "Germany"                   
## [13] "Greece"                     "Guatemala"                 
## [15] "Haiti"                      "Holand-Netherlands"        
## [17] "Honduras"                   "Hong"                      
## [19] "Hungary"                    "India"                     
## [21] "Iran"                       "Ireland"                   
## [23] "Italy"                      "Jamaica"                   
## [25] "Japan"                      "Laos"                      
## [27] "Mexico"                     "Nicaragua"                 
## [29] "Outlying-US(Guam-USVI-etc)" "Peru"                      
## [31] "Philippines"                "Poland"                    
## [33] "Portugal"                   "Puerto-Rico"               
## [35] "Scotland"                   "South"                     
## [37] "Taiwan"                     "Thailand"                  
## [39] "Trinadad&Tobago"            "United-States"             
## [41] "Vietnam"                    "Yugoslavia"

Vemos que las variables workclass, occupation y native-country contienen el caracter “?”, y eso es porque el dato no está informado. Vamos a ver la ocurrencia en cada variable del nivel “?”

count(df,df$workclass=="?")
## # A tibble: 2 x 2
##   `df$workclass == "?"`     n
##   <lgl>                 <int>
## 1 FALSE                 45970
## 2 TRUE                   2789
count(df,df$occupation=="?")
## # A tibble: 2 x 2
##   `df$occupation == "?"`     n
##   <lgl>                  <int>
## 1 FALSE                  45960
## 2 TRUE                    2799
count(df,df$`native-country`=="?")
## # A tibble: 2 x 2
##   `df$\`native-country\` == "?"`     n
##   <lgl>                          <int>
## 1 FALSE                          47903
## 2 TRUE                             856

Son pocos registros (workclass son 2789, en occupation con 2799 y en native-country son 856). Los eliminamos de los datos.

df <- filter(df,df$workclass!="?")
df$workclass <- factor(df$workclass)
df <- filter(df,df$occupation!="?")
df$occupation <- factor(df$occupation)
df <- filter(df,df$`native-country`!="?")
df$`native-country` <- factor(df$`native-country`)

Vamos a cambiar los guiones de los nombres de las variables, ya que pueden confundir a las funcionar que usaemos en el futuro, ya que pueden considerarlo una resta, dependiendo de como estén programadas. O sustituimos el guion medio por un guion bajo, o simplificamos el nombre de la variable, dependiendo del caso.

names(df)[5] <- "education_num"
names(df)[6] <- "marital_status"
names(df)[11] <- "gain"
names(df)[12] <- "loss"
names(df)[13] <- "hours_per_week"
names(df)[14] <- "country"
names(df)
##  [1] "age"            "workclass"      "fnlwgt"         "education"     
##  [5] "education_num"  "marital_status" "occupation"     "relationship"  
##  [9] "race"           "sex"            "gain"           "loss"          
## [13] "hours_per_week" "country"        "income"

Ahora ya tenemos en df el dataframe con los datos listos después de la parte de calidad de datos. Eliminamos algunas variables que ya no utilizaremos.

rm(list = c("test","Capitalgain","Capitalloss"))

3 - Trasformación de datos

3.1 - Creacion de la variable target

En este caso, la variable ya la tenemos en formato dicotómico, y nos han marcado el valor del umbral (50K dolares). Antes ya hemos traducido a formato 0/1, por lo que no necesitamos hacer nada en este punto.

3.2 - Preparacion de las variables independientes

3.2.1 - Preseleccion de variables independientes

Creamos una lista larga con todas las variables independientes. En este punto, no descartaremos ninguna variable. Al no tener variables históricas ni tampoco un código númerico para identificar cada observación, no tenemos por que deseleccionar ninguna. Solo eliminamos la variable target.

ind_larga<-names(df) #lista con todas las variables
no_usar <- c('income') #identificamos las que no queremos usar
ind_larga <- setdiff(ind_larga,no_usar)

En ind_larga tenemos la lista de variables final que utilizaremos para la preselección.

Creamos una muestra m menor para que los calculos sean mas rapidos. Reducir la muestra no restará eficacia a los algoritmos de preselección. El tamaño de la muestra se determina para que, estadisticamente, sea relevante. Seleccionaremos 20000 registros, que son casi la mitad de los registros y por tanto, será estadísticamente relevante con un error muy bajo.

set.seed(12345)
m <- sample_n(df,20000)

3.2.1.1 - Preseleccion con RandomForest

Vamos ha hacer la primera preselección con RandomForest.

pre_rf <- randomForest(formula = reformulate(ind_larga,"income"),data= m,mtry=2,ntree=50, importance = T)
imp_rf <- importance(pre_rf)[,4] #como importance devuelve una matriz con varias metricas, tenemos que extraer asi el decrecimiento en Gini (metrica 4) que es el que mas nos interesa, porque la experiencia demuestra que funciona bien, pero también se pueden probar otras metricas de las que devuelve en la matriz de importancia, o una combinación de ambas
imp_rf <- data.frame(VARIABLE = names(imp_rf), IMP_RF = imp_rf) #lo transformamos a dataframe
imp_rf <- imp_rf %>% arrange(desc(IMP_RF)) %>% mutate(RANKING_RF = 1:nrow(imp_rf)) #creamos el ranking
imp_rf
##          VARIABLE    IMP_RF RANKING_RF
## 1    relationship 654.76379          1
## 2            gain 620.53506          2
## 3             age 611.51912          3
## 4  marital_status 595.57102          4
## 5      occupation 532.07896          5
## 6       education 431.54808          6
## 7          fnlwgt 406.11714          7
## 8   education_num 400.17564          8
## 9  hours_per_week 347.33342          9
## 10      workclass 190.36565         10
## 11           loss 184.13443         11
## 12        country 143.24987         12
## 13            sex 100.97497         13
## 14           race  70.41341         14

3.2.1.2 - Preseleccion con Information Value

La función de IV necesita que la variable target sea numerica. Como nosotros la tenemos como factor, creamos un dataset diferente y le cambiamos la valiable target de factor a numeric.

m2 <- mutate(m,income = as.numeric(as.character(income))) #transformo la target a numerico temporalmente porque este algoritmo necesita que este en numerico, y el as.character es para que lo convierta a 0 y 1, y no a 1 y 2
imp_iv <- smbinning.sumiv(m2[c(ind_larga,'income')],y="income")
##  
## 
  |                                                        
  |                                                  |   0%
  |                                                        
  |---                                               |   7%
  |                                                        
  |-------                                           |  13%
  |                                                        
  |----------                                        |  20%
  |                                                        
  |-------------                                     |  27%
  |                                                        
  |-----------------                                 |  33%
  |                                                        
  |--------------------                              |  40%
  |                                                        
  |-----------------------                           |  47%
  |                                                        
  |---------------------------                       |  53%
  |                                                        
  |------------------------------                    |  60%
  |                                                        
  |---------------------------------                 |  67%
  |                                                        
  |-------------------------------------             |  73%
  |                                                        
  |----------------------------------------          |  80%
  |                                                        
  |-------------------------------------------       |  87%
  |                                                        
  |-----------------------------------------------   |  93%
  |                                                        
  |--------------------------------------------------| 100%
## 
imp_iv <- imp_iv %>% mutate(Ranking = 1:nrow(imp_iv)) %>% select(-Process)
names(imp_iv) <- c('VARIABLE','IMP_IV','RANKING_IV')
imp_iv
##          VARIABLE IMP_IV RANKING_IV
## 1    relationship 1.4508          1
## 2  marital_status 1.2841          2
## 3             age 1.0972          3
## 4   education_num 0.6880          4
## 5            gain 0.6409          5
## 6  hours_per_week 0.4603          6
## 7             sex 0.3059          7
## 8       workclass 0.1244          8
## 9            race 0.0712          9
## 10         fnlwgt     NA         10
## 11      education     NA         11
## 12     occupation     NA         12
## 13           loss     NA         13
## 14        country     NA         14

3.2.1.3 - Preseleccion final

En este punto ya tenemos dos metodologías distintas con las que hemos identificado a las variables que van a ser más predictivas. Lo que vamos a hacer es ver cómo las podemos unificar y obtener un criterio unificado para estar muy seguros de que estamos seleccionando las variables correctas. Y hecho eso vamos a dejar ya identificada la lista corta de variables que pasará a la siguiente fase.

imp_final <- inner_join(imp_rf,imp_iv,by='VARIABLE') %>% 
  select(VARIABLE,IMP_RF,IMP_IV,RANKING_RF,RANKING_IV) %>% #ponerlos en orden mas legible
  mutate(RANKING_TOT = RANKING_RF + RANKING_IV) %>% 
  arrange(RANKING_TOT)
imp_final
##          VARIABLE    IMP_RF IMP_IV RANKING_RF RANKING_IV RANKING_TOT
## 1    relationship 654.76379 1.4508          1          1           2
## 2             age 611.51912 1.0972          3          3           6
## 3  marital_status 595.57102 1.2841          4          2           6
## 4            gain 620.53506 0.6409          2          5           7
## 5   education_num 400.17564 0.6880          8          4          12
## 6  hours_per_week 347.33342 0.4603          9          6          15
## 7      occupation 532.07896     NA          5         12          17
## 8       education 431.54808     NA          6         11          17
## 9          fnlwgt 406.11714     NA          7         10          17
## 10      workclass 190.36565 0.1244         10          8          18
## 11            sex 100.97497 0.3059         13          7          20
## 12           race  70.41341 0.0712         14          9          23
## 13           loss 184.13443     NA         11         13          24
## 14        country 143.24987     NA         12         14          26

¿Son los metodos fiables? Si los resultados fueran muy diferentes, podrían que ambos metodos se anularan y la clasificación final no se corresponde con el nivel de predicción de la variable target. Vamos a hacer una correlacion entre ellos a ver si dan cosas similares.

cor(imp_final$IMP_RF,imp_final$IMP_IV,use = 'complete.obs')
## [1] 0.8879993

Si nos dan fiabilidad los resultados ya que la correlacion es normal y por tanto los resultados de ambos valores se refuerzan.

Decision: vamos a descartar aquellas variables que no hayan salido entre las 10 mas importantes en ninguno de los dos sistemas de seleccion de variables. Con esto solo descartamos la variable capital-loss y la variable country, por lo que abandonamos la idea de combinar la primera con la variable capital-gain y nos quedamos con esta ultima, que el metodo de RF dice que es bastante predictiva.

ind_corta <- imp_final %>%
  filter(RANKING_RF <= 10 | RANKING_IV <= 10) %>% 
  select(VARIABLE) #nos quedamos solo con el nombre
ind_corta <- as.character(ind_corta$VARIABLE) #lo pasamos a un vector en vez de un dataframe

Estas son las variables predictoras con las que vamos a trabajar finalmente. Son 13 variables.

ind_corta
##  [1] "relationship"   "age"            "marital_status" "gain"          
##  [5] "education_num"  "hours_per_week" "occupation"     "education"     
##  [9] "fnlwgt"         "workclass"      "sex"            "race"

3.2.2 - Seleccionar la lista de variables finales del proyecto

En este caso, no hemos quitado variables originalmente, así que no hay que meter ninguna variable.

Solo tenemos que unir la variable target, que quitamos al inicio.

finales <- union(ind_corta,c('income'))#creamos la lista de variables finales, incluyendo las que no tienen historia y la target
finales
##  [1] "relationship"   "age"            "marital_status" "gain"          
##  [5] "education_num"  "hours_per_week" "occupation"     "education"     
##  [9] "fnlwgt"         "workclass"      "sex"            "race"          
## [13] "income"

3.3 - Fichero final y limpieza del entorno

3.3.1 - Fichero final

Generamos el dataset final con las variables seleccionadas.

df <- df %>% 
  select(one_of(finales))

3.3.2 - Limpieza del entorno

Durante todo el proceso anterior hemos creado muchas variables y ficheros temporales, vamos a aprovechar para limpiarlo todo y dejarlo organizado antes de pasar a la siguiente fase

ls() #Vemos todo lo que tenemos cargado en el entorno
##  [1] "a_factores"                 "atipicos_variables_integer"
##  [3] "atipicos_variables_numeric" "contar_ceros"              
##  [5] "df"                         "est_basicos"               
##  [7] "finales"                    "i"                         
##  [9] "imp_final"                  "imp_iv"                    
## [11] "imp_rf"                     "ind_corta"                 
## [13] "ind_larga"                  "instalados"                
## [15] "longi"                      "m"                         
## [17] "m2"                         "no_usar"                   
## [19] "num_ceros"                  "num_nulos"                 
## [21] "out"                        "out_int"                   
## [23] "paquetes"                   "pre_rf"
rm(list=setdiff(ls(),c('df','contar_ceros','out_int'))) #borramos todo excepto el nuevo df
# y vamos a dejar preparadas unas variables que nos van a facilitar cosas en el futuro:
target <- 'income'
indep <- setdiff(names(df),c(target))

Vamos a guardar un cache temporal de datos, de forma que cuando queramos seguir trabajando desde aqui no tengamos que volver a ejecutar todo

saveRDS(df,'cache1.rds')

4 - Creacion de variables sinteticas

Cargamos el cache que hemos grabado de la fase anterior, y volvemos a crear la variable target y el conjunto de variables independientes:

df <- readRDS(file = 'cache1.rds')
target <- 'income'
indep <- setdiff(names(df),c(target))
no_borrar <- c("indep","target","df")
no_borrar <- append(no_borrar,"no_borrar")

4.1 - Creacion de variables sinteticas

No tenemos datos históricos, asi que no tiene sentido crear variables de tenencia, contratación, cancelación, medias o tendencias. Solo tenemos una foto estática con la que tenemos que predecir la variable target income.

Aparte de la discretización de las variables continuas, que haremos más adelante, podemos intentar combinar algunas de las variables categoricas que tenemos, y que tenga sentido que estén relacionadas. Las variables que tenemos y sus niveles son los siguientes:

for(i in 1:length(indep)) {
  if(is.factor(df[[i]])) {
      cat("Variable -",indep[i],"\n")
      cat("Clase -",class(df[[i]]),"\n")
      niveles <- levels(df[[i]])
      cat("Niveles -",niveles,"\n")
      print("-------------------------------")
  }
}
## Variable - relationship 
## Clase - factor 
## Niveles - Husband Not-in-family Other-relative Own-child Unmarried Wife 
## [1] "-------------------------------"
## Variable - marital_status 
## Clase - factor 
## Niveles - Divorced Married-AF-spouse Married-civ-spouse Married-spouse-absent Never-married Separated Widowed 
## [1] "-------------------------------"
## Variable - occupation 
## Clase - factor 
## Niveles - Adm-clerical Armed-Forces Craft-repair Exec-managerial Farming-fishing Handlers-cleaners Machine-op-inspct Other-service Priv-house-serv Prof-specialty Protective-serv Sales Tech-support Transport-moving 
## [1] "-------------------------------"
## Variable - education 
## Clase - factor 
## Niveles - 10th 11th 12th 1st-4th 5th-6th 7th-8th 9th Assoc-acdm Assoc-voc Bachelors Doctorate HS-grad Masters Prof-school Some-college 
## [1] "-------------------------------"
## Variable - workclass 
## Clase - factor 
## Niveles - Federal-gov Local-gov Never-worked Private Self-emp-inc Self-emp-not-inc State-gov Without-pay 
## [1] "-------------------------------"
## Variable - sex 
## Clase - factor 
## Niveles - Female Male 
## [1] "-------------------------------"
## Variable - race 
## Clase - factor 
## Niveles - Amer-Indian-Eskimo Asian-Pac-Islander Black Other White 
## [1] "-------------------------------"

Por la naturaleza de las variables, podríamos hacer las siguientes transformaciones:

  • Combinar relationship con marital_status de alguna manera.
  • Reducir el número de niveles de occupation, de tal manera que agrupen por tipo de trabajo u otra agrupación, y que mantenga la distribución original.
  • Reducir los niveles de la variable education, agrupando en niveles de primaria, secundaria, graduados, doctorados, master, formación profesional.
  • Reducir workclass a empleo publico, privado, sin remuneración, autonomo.

Vamos con la combinación de relantionship con marital_status. Vemos como se distribuyen una con respecto de la otra:

table(df$relationship,df$marital_status)
##                 
##                  Divorced Married-AF-spouse Married-civ-spouse
##   Husband               0                11              18636
##   Not-in-family      3434                 0                 18
##   Other-relative      166                 1                182
##   Own-child           429                 1                125
##   Unmarried          2266                 0                  0
##   Wife                  0                19               2070
##                 
##                  Married-spouse-absent Never-married Separated Widowed
##   Husband                            0             0         0       0
##   Not-in-family                    280          6663       587     686
##   Other-relative                    44           819        75      58
##   Own-child                         56          5858       130      20
##   Unmarried                        168          1221       617     510
##   Wife                               0             0         0       0
ggplot(df,aes(relationship,marital_status)) + geom_count()

Vamos a agrupar de la siguiente manera, para reducir los niveles a cuatro: with_children, married, unmarried, other. Mantenemos más o menos la distribución, y reducimos en 1 la dimensionalidad.:

df <- df %>% mutate(c_status = as.factor(case_when(
      (relationship == "Other-relative") | (relationship =="Not-in-family" & marital_status =="Married-civ-spouse") ~ '01_other', 
      relationship == "Own-child"  ~ '02_with_children',
      (relationship == "Unmarried" | relationship =="Not-in-family") & (marital_status =="Married-spouse-absent" | marital_status =="Divorced" | marital_status =="Separated" | marital_status =="Never-married"  | marital_status =="Widowed") ~ '03_unmarried',
      (relationship == "Wife" | relationship =="Husband") & (marital_status =="Married-AF-spouse" | marital_status =="Married-civ-spouse") ~ '04_married',
      TRUE ~ '00_ERROR'))
)
#vemos la distribución, nos aseguramos que no tenemos ningún caso en error
table(df$c_status) 
## 
##         01_other 02_with_children     03_unmarried       04_married 
##             1363             6619            16432            20736
ggplot(df,aes(c_status)) + geom_bar()

#eliminamos las variables originales
df <- select(df, -relationship)
df <- select(df, -marital_status)

Ahora pasamos a reducir el numero de niveles de occupation y mantener una distribución lo mas homogenea y monotonica posible entre los nuevos niveles.

table(df$occupation)
## 
##      Adm-clerical      Armed-Forces      Craft-repair   Exec-managerial 
##              5537                14              6014              5983 
##   Farming-fishing Handlers-cleaners Machine-op-inspct     Other-service 
##              1463              2041              2958              4787 
##   Priv-house-serv    Prof-specialty   Protective-serv             Sales 
##               230              6007               976              5406 
##      Tech-support  Transport-moving 
##              1420              2314
ggplot(df,aes(occupation)) + geom_bar()

df <- df %>% mutate(worktype = as.factor(case_when(
      occupation == "Adm-clerical"  ~ '01_admin',
      occupation== "Exec-managerial" ~ '02_Executives',
      (occupation== "Sales" | occupation== "Tech-support" | occupation== "Armed-Forces") ~ '03_Technical',
      (occupation== "Other-service" | occupation== "Priv-house-serv" | occupation== "Prof-specialty" | occupation== "Protective-serv") ~ '04_services',
      (occupation== "Craft-repair" | occupation== "Farming-fishing" | occupation== "Handlers-cleaners" | occupation== "Machine-op-inspct" | occupation== "Transport-moving") ~ '05_manual_work',
      TRUE ~ '00_ERROR'))
)
#vemos la distribución, nos aseguramos que no tenemos ningún caso en error
table(df$worktype)
## 
##       01_admin  02_Executives   03_Technical    04_services 05_manual_work 
##           5537           5983           6840          12000          14790
ggplot(df,aes(worktype)) + geom_bar()

#eliminamos la variable original
df <- select(df, -occupation)

La siguiene variable es education, intentaremos agrupar por nivel de estudios, primaria, secundaria, graduados, doctorados, master, formación profesional.

table(df$education)
## 
##         10th         11th         12th      1st-4th      5th-6th 
##         1223         1619          577          222          449 
##      7th-8th          9th   Assoc-acdm    Assoc-voc    Bachelors 
##          823          676         1507         1959         7570 
##    Doctorate      HS-grad      Masters  Prof-school Some-college 
##          544        14783         2514          785         9899
ggplot(df,aes(education)) + geom_bar()

df <- df %>% mutate(educational_level = as.factor(case_when(
      (education == "1st-4th" | education == "5th-6th" | education == "7th-8th" | education == "9th" | education == "10th" | education == "11th" | education == "12th")  ~ '01_primary',
      education == "HS-grad"  ~ '02_high_school',
      (education == "Prof-school" | education == "Some-college" | education == "Assoc-acdm")  ~ '03_professional',
      education == "Bachelors"  ~ '04_bachelor',
      (education == "Masters" | education == "Assoc-voc")  ~ '05_master',
      education == "Doctorate"  ~ '06_doctor',
      TRUE ~ '00_ERROR'))
)

#vemos la distribución, nos aseguramos que no tenemos ningún caso en error
table(df$educational_level)
## 
##      01_primary  02_high_school 03_professional     04_bachelor 
##            5589           14783           12191            7570 
##       05_master       06_doctor 
##            4473             544
ggplot(df,aes(educational_level)) + geom_bar()

#eliminamos la variable original
df <- select(df, -education)

Ahora vamos a reducir el numero de niveles de la variable workclass, agrupandolo por el origen del dinero entrante.

table(df$workclass)
## 
##      Federal-gov        Local-gov     Never-worked          Private 
##             1406             3096                0            33245 
##     Self-emp-inc Self-emp-not-inc        State-gov      Without-pay 
##             1646             3791             1945               21
ggplot(df,aes(workclass)) + geom_bar()

df <- df %>% mutate(work = as.factor(case_when(
      workclass == "Self-emp-inc"  ~ '01_self_emp',
      (workclass == "Self-emp-not-inc" | workclass == "Without-pay")  ~ '02_not_inc',
      (workclass == "Federal-gov" | workclass == "Local-gov" | workclass == "State-gov")  ~ '03_public',
      workclass == "Private"  ~ '04_private',
      TRUE ~ '00_ERROR'))
)

#vemos la distribución, nos aseguramos que no tenemos ningún caso en error
table(df$work)
## 
## 01_self_emp  02_not_inc   03_public  04_private 
##        1646        3812        6447       33245
ggplot(df,aes(work)) + geom_bar()

#eliminamos la variable original
df <- select(df, -work)

Si vemos la relación entre education_num y la nueva variable sintética educational_level, vemos que la relación es practicaemnte lineal y que a cada grupo de numeros de education_num se corresponde un nivel de educational_level. Por lo tanto, podemos eliminar education_num.

plot(df$education_num, df$educational_level, xlab="Education Num", ylab="Education Level", pch=19)

df <- select(df,-education_num)

Con esto ya hemos simplificado aquellas variables categoricas que tenía sentido simplificar, para reducir el numero de factores y porque el comportamiento no será muy diferente con respecto a la variable terget, creando unas nuevas variables sintéticas. Ahora ya podemos pasar a la discretización.

rm(list=setdiff(ls(),no_borrar))

4.2. Discretización

4.2.1. Discretización automática

Primero vamos a crear la funcion que va a discretizar de forma automatica maximizando la capacidad predictiva de la nueva variable. Ademas, como vamos a usar en la modelizacion un algoritmo lineal, que es la regresion logistica, vamos a intentar que la discretizacion sea monotonica.

Las variables que tenemos para discretizar en nuestro dataset son age, hours_per_week, fnlwgt y gain.

Creamos la función de discretización automatica:

discretizar <- function(vi,target){
  temp_df <- data.frame(vi = vi, target = target)
  #smbinning necesita que la target sea numerica
  temp_df$target <- as.numeric(as.character(temp_df$target))
  disc <- smbinning(temp_df, y = 'target', x = 'vi')
  return(disc)
}

Ahora hacemos las discretizaciones automáticamente:

#age:

disc_temp_age <- discretizar(df$age,df$income)
df_temp <- select(df,age,income) #creamos este temporal porque smbinning.gen necesita que el df tenga el mismo numero de columnas que la salida de la funcion discretizar
df_temp <- smbinning.gen(df_temp,disc_temp_age,chrname = 'age_disc')
#Metemos en df la nueva variable discretizada y eliminamos la original
df <- cbind(df,df_temp[3]) %>% select(-age)

#hours_per_week

disc_temp_hours <- discretizar(df$hours_per_week,df$income)
df_temp <- select(df,hours_per_week,income) #creamos este temporal porque smbinning.gen necesita que el df tenga el mismo numero de columnas que la salida de la funcion discretizar
df_temp <- smbinning.gen(df_temp,disc_temp_hours,chrname = 'hours_per_week_disc')
#Metemos en df la nueva variable discretizada y eliminamos la original
df <- cbind(df,df_temp[3]) %>% select(-hours_per_week)

#gain

disc_temp_gain <- discretizar(df$gain,df$income)
df_temp <- select(df,gain,income) #creamos este temporal porque smbinning.gen necesita que el df tenga el mismo numero de columnas que la salida de la funcion discretizar
df_temp <- smbinning.gen(df_temp,disc_temp_gain,chrname = 'gain_disc')
#Metemos en df la nueva variable discretizada y eliminamos la original
df <- cbind(df,df_temp[3]) %>% select(-gain)

Como la variable fnlwgt ha dado un error con la función smbinning, procederemos a discretizarla manualmente.

4.2.2. Discretización manual

Lo primero es ver que distribucion tiene la variable fnlwgt. La variable discretizada debería de ser parecida a la distribución de la variable original.

También la distribución de la variable target en cada categoria a realizar sea aproximadamente la misma, para ayudar a los modelos y aproximarnos a tener una variable monotonica.

ggplot(df,aes(fnlwgt)) + geom_density()

#No se ve muy bien la distribucion por culpa de los atipicos. Vamos a limitar el eje x
ggplot(df,aes(fnlwgt)) + geom_density() + scale_x_continuous(limits = c(0, 500000))

#Parece que la gran mayoria esta por debajo de 500000€, vamos a limitar un poco mas
ggplot(df,aes(fnlwgt)) + geom_density() + scale_x_continuous(limits = c(0, 400000))

#Parece que la gran mayoria esta por debajo de 400000€, vamos a limitar un poco mas
ggplot(df,aes(fnlwgt)) + geom_density() + scale_x_continuous(limits = c(0, 300000))

#Parece que la gran mayoria esta por debajo de 300000€, vamos a limitar un poco mas
ggplot(df,aes(fnlwgt)) + geom_density() + scale_x_continuous(limits = c(0, 250000))

#Parece que la gran mayoria esta por debajo de 250000€, vamos a limitar un poco mas
ggplot(df,aes(fnlwgt)) + geom_density() + scale_x_continuous(limits = c(0, 210000))

#Pedimos un histograma para aproximar mejor la forma que queremos conseguir, que es una discretización monotonica
ggplot(df,aes(fnlwgt)) + geom_histogram(bins = 40) + scale_x_continuous(limits = c(0, 210000))

#Ya sabemos que queremos una forma decreciente ahora veamos como se comporta la variable target para ver si podremos generar algo monotonico, y con la distribucion de la variable target similar en cada tramo.  Hay que probar con varios valores para el numero de tramos  y coger el que mejor se aproxime.
ggplot(df,aes(fnlwgt,fill=income)) + geom_histogram(bins = 40,position='fill') + scale_x_continuous(limits = c(0, 210000)) 

#Sabiendo ambas cosas vamos a apoyarnos en los deciles para intuir donde podemos hacer buenos cortes y no dejar a muchos clientes en el mismo tramo y asi distribuir los target entre todos los tramos
as.data.frame(quantile(df$fnlwgt,prob = seq(0, 1, length = 11)))
##      quantile(df$fnlwgt, prob = seq(0, 1, length = 11))
## 0%                                              13492.0
## 10%                                             66118.0
## 20%                                            105936.0
## 30%                                            130404.7
## 40%                                            157950.6
## 50%                                            178250.0
## 60%                                            196342.8
## 70%                                            220237.0
## 80%                                            260335.6
## 90%                                            328466.0
## 100%                                          1490400.0

Procedemos a discretizarla manualmente

df <- df %>% mutate(fnlwgt = as.factor(case_when(
      fnlwgt <= 65000 ~ '01_MENOR_65K',
      fnlwgt > 65000 & fnlwgt <= 120000 ~ '02_DE_65K_A_120K',
      fnlwgt > 120000 & fnlwgt <= 180000 ~ '03_DE_105K_A_180K',
      fnlwgt > 180000 & fnlwgt <= 320000 ~ '04_DE_160K_A_320K',
      fnlwgt > 320000  ~ '05_MAS_320K',
      TRUE ~ '00_ERROR'))
)

Veamos si la distribucion ha quedado similar a la original, que no la hemos cambiado con la discretización.

ggplot(df,aes(fnlwgt)) + geom_bar()

Y ahora vamos a comprobar si la penetracion de la target es similar también a la original.

ggplot(df,aes(fnlwgt,fill=income)) + geom_bar(position='fill')

Cambiamos el nombre de la variable para indicar que es discretizada:

names(df)[1] <- "fnlwgt_disc"

4.2.3. Revisión de la distribución de las variables discretizadas

Vamos a hacer una inspeccion visual de todas las variables a ver si han salido bien. Vemos que todas siguen una distribución coherente con su significado.

df %>% 
  select_if(is.factor) %>% 
  gather() %>% 
  ggplot(aes(value)) +
    geom_bar() +
    facet_wrap(~ key, scales = "free") + #hace los graficos para todas las variables
    theme(axis.text=element_text(size=4))#esto es para cambiar el tamaño del texto del eje y que se lea bien

Ahora vamos a analizar la penetración de la target en cada categoría para ver si las variables han salido monotónicas

a <- function(var1,var2) {
  df_temp <- data.frame(var1 = df[[var1]],var2 = df[[var2]])
  df_temp %>% 
    group_by(var1) %>% 
    summarise(Conteo = n(), Porc = mean(as.numeric(as.character(var2)))) %>% 
  ggplot(aes(var1,Porc)) + geom_bar(stat='identity') + xlab(var1)
}
df2_nombres <- df %>% select_if(is.factor) %>% names()
lapply(df2_nombres,function(x){a(x,'income')})
## [[1]]

## 
## [[2]]

## 
## [[3]]

## 
## [[4]]

## 
## [[5]]

## 
## [[6]]

## 
## [[7]]

## 
## [[8]]

## 
## [[9]]

## 
## [[10]]

## 
## [[11]]

La mayoria han salido monotonicas o con la misma distribución original, coherente con su significado. Podemos dejar el dataset aqui, ya preparado con sus variables sinteticas y discretizaciones listas para la modelización.

4.2.4. Organización del dataframe y limpieza del entorno

Antes de continuar vamos a guardar en un objeto de R las discretizaciones, porque las necesitaremos despues para poner el modelo en produccion

#Vamos a crear un objeto de tipo lista que es lo ideal para guardar objetos complejos como las discretizaciones
discretizaciones <- list(
 disc_temp_age = disc_temp_age,
 disc_temp_gain = disc_temp_gain,
 disc_temp_hours = disc_temp_hours
)
saveRDS(discretizaciones,'02_CortesDiscretizaciones.rds')

Vamos a ver como ha quedado nuestro fichero antes de pasar a la fase de modelizacion

glimpse(df)
## Observations: 45,150
## Variables: 11
## $ fnlwgt_disc         <fct> 02_DE_65K_A_120K, 02_DE_65K_A_120K, 04_DE_...
## $ workclass           <fct> State-gov, Self-emp-not-inc, Private, Priv...
## $ sex                 <fct> Male, Male, Male, Male, Female, Female, Fe...
## $ race                <fct> White, White, White, Black, Black, White, ...
## $ income              <fct> 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 0, 0, ...
## $ c_status            <fct> 03_unmarried, 04_married, 03_unmarried, 04...
## $ worktype            <fct> 01_admin, 02_Executives, 05_manual_work, 0...
## $ educational_level   <fct> 04_bachelor, 04_bachelor, 02_high_school, ...
## $ age_disc            <fct> 08 <= 41, 10 <= 54, 07 <= 38, 10 <= 54, 04...
## $ hours_per_week_disc <fct> 03 <= 43, 01 <= 34, 03 <= 43, 03 <= 43, 03...
## $ gain_disc           <fct> 01 <= 4650, 01 <= 4650, 01 <= 4650, 01 <= ...

Vamos a reordernar las variables simplemente por una cuestion estetica

#creamos un vector con las variables centrales
centrales <- setdiff(names(df),c('income'))
df <- df %>% select(
  one_of(centrales),
  income)

Comprobamos de nuevo

glimpse(df)
## Observations: 45,150
## Variables: 11
## $ fnlwgt_disc         <fct> 02_DE_65K_A_120K, 02_DE_65K_A_120K, 04_DE_...
## $ workclass           <fct> State-gov, Self-emp-not-inc, Private, Priv...
## $ sex                 <fct> Male, Male, Male, Male, Female, Female, Fe...
## $ race                <fct> White, White, White, Black, Black, White, ...
## $ c_status            <fct> 03_unmarried, 04_married, 03_unmarried, 04...
## $ worktype            <fct> 01_admin, 02_Executives, 05_manual_work, 0...
## $ educational_level   <fct> 04_bachelor, 04_bachelor, 02_high_school, ...
## $ age_disc            <fct> 08 <= 41, 10 <= 54, 07 <= 38, 10 <= 54, 04...
## $ hours_per_week_disc <fct> 03 <= 43, 01 <= 34, 03 <= 43, 03 <= 43, 03...
## $ gain_disc           <fct> 01 <= 4650, 01 <= 4650, 01 <= 4650, 01 <= ...
## $ income              <fct> 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 0, 0, ...

Limpiamos el entorno de cualquier cosa que no sea el dataframe. Primero lo salvamos en un fichero.

Guardamos otro cache temporal

saveRDS(df,'cache3.rds')

Limpiamos la memoria de todas las variables.

a_borrar <- setdiff(ls(),no_borrar)
rm(list=c(a_borrar,'a_borrar'))

Cargamos el cache temporal

df <- readRDS('cache3.rds')

5. Modelizacion

5.1 - Preparar las funciones que vamos a necesitar

Funcion para crear una matriz de confusion.

confusion<-function(real,scoring,umbral){ 
  conf<-table(real,scoring>=umbral)
  if(ncol(conf)==2) return(conf) else return(NULL)
}

Funcion para calcular las metricas de los modelos: acierto, precision, cobertura y F1

metricas<-function(matriz_conf){
  acierto <- (matriz_conf[1,1] + matriz_conf[2,2]) / sum(matriz_conf) *100
  precision <- matriz_conf[2,2] / (matriz_conf[2,2] + matriz_conf[1,2]) *100
  cobertura <- matriz_conf[2,2] / (matriz_conf[2,2] + matriz_conf[2,1]) *100
  F1 <- 2*precision*cobertura/(precision+cobertura)
  salida<-c(acierto,precision,cobertura,F1)
  return(salida)
}

Funcion para probar distintos umbrales y ver el efecto sobre precision y cobertura. Genera una tabla con las diferentes metricas para cada umbral de corte, variando un 5% este umbral.

umbrales<-function(real,scoring){
  umbrales<-data.frame(umbral=rep(0,times=19),acierto=rep(0,times=19),precision=rep(0,times=19),cobertura=rep(0,times=19),F1=rep(0,times=19))
  cont <- 1
  for (cada in seq(0.05,0.95,by = 0.05)){
    datos<-metricas(confusion(real,scoring,cada))
    registro<-c(cada,datos)
    umbrales[cont,]<-registro
    cont <- cont + 1
  }
  return(umbrales)
}

Funciones que calculan la curva ROC y el AUC

roc<-function(prediction){
  r<-performance(prediction,'tpr','fpr')
  plot(r)
}

auc<-function(prediction){
  a<-performance(prediction,'auc')
  return(a@y.values[[1]])
}

5.2 - Creamos las particiones de training (70%) y test (30%)

Establecemos una semilla para que nos salgan los mismos resultados en todas las ejecuciones:

set.seed(12345)

Generamos una variable aleatoria con una distribucion 70-30

df$random<-sample(0:1,size = nrow(df),replace = T,prob = c(0.3,0.7))

Creamos los dos dataframes

train<-filter(df,random==1)
test<-filter(df,random==0)
#Eliminamos ya la random para que no moleste
df$random <- NULL

5.3 - Creación del modelo de propensión

Nota: Vamos a probar tres algoritmos diferentes para ver cual funciona mejor y aprender como se comparan

5.3.1 - Identificamos las variables

#Las independientes seran todas menos la target
independientes <- setdiff(names(df),c('income'))
target <- 'income'

5.3.2 - Creamos la formula para usar en el modelo

formula <- reformulate(independientes,target)

5.3.3 - Modelizamos con regresion logistica

Primero vamos a hacer un modelo con todas las variables. Cargamos la formula en una copia para la regresion logistica, ya que iremos haciendo cambios específicos para la regresión logistica. La preselección de variables es básica para la regresión logística, ya que funciona mejor con pocas variables y que no tengan correlación entre ellas.

formula_rl <- formula
rl <- glm(formula_rl,train,family=binomial(link='logit'))
summary(rl)
## 
## Call:
## glm(formula = formula_rl, family = binomial(link = "logit"), 
##     data = train)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -3.2109  -0.5249  -0.2159  -0.0114   3.4592  
## 
## Coefficients:
##                                  Estimate Std. Error z value
## (Intercept)                      -8.45290    0.48807 -17.319
## fnlwgt_disc02_DE_65K_A_120K       0.13404    0.06924   1.936
## fnlwgt_disc03_DE_105K_A_180K      0.34743    0.06512   5.335
## fnlwgt_disc04_DE_160K_A_320K      0.29068    0.06220   4.673
## fnlwgt_disc05_MAS_320K            0.43857    0.07654   5.730
## workclassLocal-gov               -0.57403    0.10761  -5.334
## workclassPrivate                 -0.46292    0.09057  -5.111
## workclassSelf-emp-inc            -0.26490    0.11860  -2.234
## workclassSelf-emp-not-inc        -0.82243    0.10483  -7.846
## workclassState-gov               -0.60667    0.11857  -5.116
## workclassWithout-pay             -0.21001    0.93868  -0.224
## sexMale                           0.17254    0.05118   3.371
## raceAsian-Pac-Islander            0.84568    0.25791   3.279
## raceBlack                         0.49229    0.24796   1.985
## raceOther                         0.67501    0.33634   2.007
## raceWhite                         0.88775    0.23848   3.723
## c_status02_with_children         -0.66259    0.22364  -2.963
## c_status03_unmarried             -0.05222    0.18575  -0.281
## c_status04_married                2.16644    0.18385  11.784
## worktype02_Executives             0.81048    0.07449  10.880
## worktype03_Technical              0.37188    0.07626   4.877
## worktype04_services               0.41506    0.07100   5.846
## worktype05_manual_work           -0.23493    0.07196  -3.265
## educational_level02_high_school   1.00052    0.08115  12.329
## educational_level03_professional  1.58623    0.08296  19.121
## educational_level04_bachelor      2.13927    0.08701  24.588
## educational_level05_master        2.04583    0.09088  22.511
## educational_level06_doctor        3.10922    0.16690  18.629
## age_disc02 <= 24                  0.86264    0.39333   2.193
## age_disc03 <= 27                  1.41853    0.37903   3.743
## age_disc04 <= 29                  1.95258    0.37784   5.168
## age_disc05 <= 33                  2.26103    0.37264   6.068
## age_disc06 <= 35                  2.47013    0.37521   6.583
## age_disc07 <= 38                  2.72706    0.37293   7.313
## age_disc08 <= 41                  2.79977    0.37318   7.502
## age_disc09 <= 45                  2.86792    0.37257   7.698
## age_disc10 <= 54                  3.04759    0.37124   8.209
## age_disc11 <= 61                  2.89946    0.37352   7.762
## age_disc12 > 61                   2.39894    0.37707   6.362
## hours_per_week_disc02 <= 39       0.63538    0.10299   6.169
## hours_per_week_disc03 <= 43       0.82123    0.07608  10.794
## hours_per_week_disc04 > 43        1.25944    0.07747  16.257
## gain_disc02 > 4650                3.19619    0.09350  34.183
##                                              Pr(>|z|)    
## (Intercept)                      < 0.0000000000000002 ***
## fnlwgt_disc02_DE_65K_A_120K                  0.052900 .  
## fnlwgt_disc03_DE_105K_A_180K     0.000000095471838707 ***
## fnlwgt_disc04_DE_160K_A_320K     0.000002963748604231 ***
## fnlwgt_disc05_MAS_320K           0.000000010059974617 ***
## workclassLocal-gov               0.000000095980947105 ***
## workclassPrivate                 0.000000320272262640 ***
## workclassSelf-emp-inc                        0.025513 *  
## workclassSelf-emp-not-inc        0.000000000000004307 ***
## workclassState-gov               0.000000311397433488 ***
## workclassWithout-pay                         0.822965    
## sexMale                                      0.000748 ***
## raceAsian-Pac-Islander                       0.001042 ** 
## raceBlack                                    0.047109 *  
## raceOther                                    0.044753 *  
## raceWhite                                    0.000197 ***
## c_status02_with_children                     0.003049 ** 
## c_status03_unmarried                         0.778598    
## c_status04_married               < 0.0000000000000002 ***
## worktype02_Executives            < 0.0000000000000002 ***
## worktype03_Technical             0.000001079825092472 ***
## worktype04_services              0.000000005049873389 ***
## worktype05_manual_work                       0.001096 ** 
## educational_level02_high_school  < 0.0000000000000002 ***
## educational_level03_professional < 0.0000000000000002 ***
## educational_level04_bachelor     < 0.0000000000000002 ***
## educational_level05_master       < 0.0000000000000002 ***
## educational_level06_doctor       < 0.0000000000000002 ***
## age_disc02 <= 24                             0.028296 *  
## age_disc03 <= 27                             0.000182 ***
## age_disc04 <= 29                 0.000000237009798587 ***
## age_disc05 <= 33                 0.000000001298139906 ***
## age_disc06 <= 35                 0.000000000045996338 ***
## age_disc07 <= 38                 0.000000000000262082 ***
## age_disc08 <= 41                 0.000000000000062633 ***
## age_disc09 <= 45                 0.000000000000013852 ***
## age_disc10 <= 54                 0.000000000000000223 ***
## age_disc11 <= 61                 0.000000000000008332 ***
## age_disc12 > 61                  0.000000000199086000 ***
## hours_per_week_disc02 <= 39      0.000000000686119204 ***
## hours_per_week_disc03 <= 43      < 0.0000000000000002 ***
## hours_per_week_disc04 > 43       < 0.0000000000000002 ***
## gain_disc02 > 4650               < 0.0000000000000002 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 35675  on 31720  degrees of freedom
## Residual deviance: 21256  on 31678  degrees of freedom
## AIC: 21342
## 
## Number of Fisher Scoring iterations: 8

Análisis de la salida:

  • La distribución de los residuos parece que sigue una distribucion más o menos normal, aunque la mediana está un poco desplazada.
  • Coeficientes: todos tienen sentido en los signos como en los valores. Por ejemplo, la workclass hace que conforme vamos a trabajos sin remuneración, más va en contra del income que se genera, la race no tienen ningún impacto, pero si el sex, en el educational-level, a mayor nivel mas ingresos (signo positivo y coeficientes más altos). Y asi todas las variables.

Vamos a seleccionar aquellas variables que tiene el valor máximo de significación en algunos de sus tramos, para volver a generar el modelo solo con estas variables más significativas.

a_mantener <- c(
  'workclass',
  'sex',
  'c_status',
  'worktype',
  'educational_level',
  'age_disc',
  'hours_per_week_disc',
  'gain_disc',
  'fnlwgt_disc'
)

formula_rl <- reformulate(a_mantener,target)
rl<- glm(formula_rl,train,family=binomial(link='logit'))
summary(rl)
## 
## Call:
## glm(formula = formula_rl, family = binomial(link = "logit"), 
##     data = train)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -3.2013  -0.5261  -0.2177  -0.0145   3.4803  
## 
## Coefficients:
##                                  Estimate Std. Error z value
## (Intercept)                      -7.70257    0.43040 -17.896
## workclassLocal-gov               -0.54265    0.10694  -5.075
## workclassPrivate                 -0.41727    0.08977  -4.648
## workclassSelf-emp-inc            -0.20690    0.11792  -1.755
## workclassSelf-emp-not-inc        -0.76888    0.10403  -7.391
## workclassState-gov               -0.58065    0.11802  -4.920
## workclassWithout-pay             -0.14151    0.93976  -0.151
## sexMale                           0.18974    0.05102   3.719
## c_status02_with_children         -0.63565    0.22328  -2.847
## c_status03_unmarried             -0.02185    0.18522  -0.118
## c_status04_married                2.20866    0.18336  12.045
## worktype02_Executives             0.81273    0.07438  10.927
## worktype03_Technical              0.37680    0.07614   4.949
## worktype04_services               0.41222    0.07084   5.819
## worktype05_manual_work           -0.23613    0.07182  -3.288
## educational_level02_high_school   1.01268    0.08111  12.485
## educational_level03_professional  1.60008    0.08292  19.297
## educational_level04_bachelor      2.16497    0.08690  24.912
## educational_level05_master        2.07077    0.09077  22.814
## educational_level06_doctor        3.14619    0.16641  18.906
## age_disc02 <= 24                  0.86292    0.39329   2.194
## age_disc03 <= 27                  1.41475    0.37897   3.733
## age_disc04 <= 29                  1.94248    0.37777   5.142
## age_disc05 <= 33                  2.25395    0.37256   6.050
## age_disc06 <= 35                  2.45794    0.37513   6.552
## age_disc07 <= 38                  2.72303    0.37286   7.303
## age_disc08 <= 41                  2.79022    0.37311   7.478
## age_disc09 <= 45                  2.85111    0.37249   7.654
## age_disc10 <= 54                  3.03981    0.37117   8.190
## age_disc11 <= 61                  2.89874    0.37345   7.762
## age_disc12 > 61                   2.40241    0.37703   6.372
## hours_per_week_disc02 <= 39       0.62441    0.10282   6.073
## hours_per_week_disc03 <= 43       0.81354    0.07594  10.713
## hours_per_week_disc04 > 43        1.26790    0.07733  16.395
## gain_disc02 > 4650                3.18229    0.09330  34.109
## fnlwgt_disc02_DE_65K_A_120K       0.14345    0.06897   2.080
## fnlwgt_disc03_DE_105K_A_180K      0.34933    0.06489   5.383
## fnlwgt_disc04_DE_160K_A_320K      0.28807    0.06193   4.652
## fnlwgt_disc05_MAS_320K            0.42317    0.07607   5.563
##                                              Pr(>|z|)    
## (Intercept)                      < 0.0000000000000002 ***
## workclassLocal-gov               0.000000388416382708 ***
## workclassPrivate                 0.000003344965257847 ***
## workclassSelf-emp-inc                        0.079332 .  
## workclassSelf-emp-not-inc        0.000000000000145370 ***
## workclassState-gov               0.000000866286557233 ***
## workclassWithout-pay                         0.880305    
## sexMale                                      0.000200 ***
## c_status02_with_children                     0.004414 ** 
## c_status03_unmarried                         0.906092    
## c_status04_married               < 0.0000000000000002 ***
## worktype02_Executives            < 0.0000000000000002 ***
## worktype03_Technical             0.000000747580565823 ***
## worktype04_services              0.000000005906751946 ***
## worktype05_manual_work                       0.001010 ** 
## educational_level02_high_school  < 0.0000000000000002 ***
## educational_level03_professional < 0.0000000000000002 ***
## educational_level04_bachelor     < 0.0000000000000002 ***
## educational_level05_master       < 0.0000000000000002 ***
## educational_level06_doctor       < 0.0000000000000002 ***
## age_disc02 <= 24                             0.028227 *  
## age_disc03 <= 27                             0.000189 ***
## age_disc04 <= 29                 0.000000271916541877 ***
## age_disc05 <= 33                 0.000000001448779926 ***
## age_disc06 <= 35                 0.000000000056660876 ***
## age_disc07 <= 38                 0.000000000000281076 ***
## age_disc08 <= 41                 0.000000000000075255 ***
## age_disc09 <= 45                 0.000000000000019451 ***
## age_disc10 <= 54                 0.000000000000000262 ***
## age_disc11 <= 61                 0.000000000000008361 ***
## age_disc12 > 61                  0.000000000186626511 ***
## hours_per_week_disc02 <= 39      0.000000001258484460 ***
## hours_per_week_disc03 <= 43      < 0.0000000000000002 ***
## hours_per_week_disc04 > 43       < 0.0000000000000002 ***
## gain_disc02 > 4650               < 0.0000000000000002 ***
## fnlwgt_disc02_DE_65K_A_120K                  0.037528 *  
## fnlwgt_disc03_DE_105K_A_180K     0.000000073156353588 ***
## fnlwgt_disc04_DE_160K_A_320K     0.000003287975621009 ***
## fnlwgt_disc05_MAS_320K           0.000000026522660906 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 35675  on 31720  degrees of freedom
## Residual deviance: 21302  on 31682  degrees of freedom
## AIC: 21380
## 
## Number of Fisher Scoring iterations: 8

Ahora todas las variables son significativas, y los coeficientes siguen siendo coherentes con el significado de las variables.

Y calculamos el pseudo R cuadrado, a partir de 0.5 el modelo puede considerarse bueno, y nos dice cuantos motivos podemos explicar de porque se contrata el deposito:

pr2_rl <- 1 -(rl$deviance / rl$null.deviance)
pr2_rl
## [1] 0.4029022

Aplicamos el modelo al conjunto de test, generando un vector con las probabilidades

rl_predict <- predict(rl,test,type = 'response')

Vemos que pinta tiene mirando la variable target con la predicción:

plot(rl_predict~test$income)

Aunque las variables están separadas y el modelo predice los 0/1, hay muy poca separación entre las predicciones de 0 y de 1.

Ahora tenemos que transformar la probabilidad en una decision del volumen de ingresos, y esto lo hacemos con los umbrales.

Con la funcion umbrales probamos diferentes cortes

umb_rl<-umbrales(test$income,rl_predict)
umb_rl
##    umbral  acierto precision cobertura       F1
## 1    0.05 59.65448  37.47807  97.77235 54.18569
## 2    0.10 69.34247  44.07115  95.27006 60.26445
## 3    0.15 74.94229  49.27750  91.57766 64.07601
## 4    0.20 78.04006  53.01471  88.00732 66.16955
## 5    0.25 80.61658  56.89726  84.83369 68.11221
## 6    0.30 82.24738  60.20571  80.37839 68.84475
## 7    0.35 83.21543  62.95265  75.86207 68.80709
## 8    0.40 83.66967  65.10591  71.28471 68.05535
## 9    0.45 84.20582  68.09643  66.37168 67.22299
## 10   0.50 84.46645  70.79986  61.85536 66.02606
## 11   0.55 84.58560  73.76920  57.15594 64.40853
## 12   0.60 84.34731  76.78979  51.38847 61.57221
## 13   0.65 84.04200  79.93664  46.20079 58.55734
## 14   0.70 83.20798  82.84062  39.33476 53.34161
## 15   0.75 82.18780  86.00488  32.25511 46.91522
## 16   0.80 81.05592  89.79370  25.23650 39.39971
## 17   0.85 79.95383  92.45283  19.43851 32.12305
## 18   0.90 79.17194  93.63636  15.71559 26.91403
## 19   0.95 78.43473  95.46539  12.20629 21.64502

Seleccionamos el umbral que maximiza la F1

umbral_final_rl<-umb_rl[which.max(umb_rl$F1),1]
umbral_final_rl
## [1] 0.3

Evaluamos la matriz de confusion y las metricas con el umbral optimizado

confusion(test$income,rl_predict,umbral_final_rl)
##     
## real FALSE TRUE
##    0  8411 1741
##    1   643 2634
rl_metricas<-filter(umb_rl,umbral==umbral_final_rl)
rl_metricas
##   umbral  acierto precision cobertura       F1
## 1    0.3 82.24738  60.20571  80.37839 68.84475

Evaluamos la ROC

#creamos el objeto prediction
rl_prediction<-prediction(rl_predict,test$income)
#visualizamos la ROC
roc(rl_prediction)

Sacamos las metricas definitivas incluyendo el AUC, que es un 69%.

rl_metricas<-cbind(rl_metricas,AUC=round(auc(rl_prediction),2)*100)
print(t(rl_metricas))
##               [,1]
## umbral     0.30000
## acierto   82.24738
## precision 60.20571
## cobertura 80.37839
## F1        68.84475
## AUC       90.00000

También dibujamos el gain chart:

plot(performance(rl_prediction,'tpr','rpp'))

5.3.4 - Modelizamos con Arboles de decision

5.3.4.1. Modelización

Los arboles de decisión permiten explicar muy bien las decisiones que toman sobre los datos.

Creamos el primer modelo

formula_ar <- formula
ar <- rpart(formula_ar, train, method = 'class', parms = list(
  split = "information"), 
  control = rpart.control(cp = 0.00001))

Revisamos donde el error de validacion cruzada empieza a crecer. Hay que fijarse en el xerror, donde cambia de decreciente a creciente.

printcp(ar)
## 
## Classification tree:
## rpart(formula = formula_ar, data = train, method = "class", parms = list(split = "information"), 
##     control = rpart.control(cp = 0.00001))
## 
## Variables actually used in tree construction:
##  [1] age_disc            c_status            educational_level  
##  [4] fnlwgt_disc         gain_disc           hours_per_week_disc
##  [7] race                sex                 workclass          
## [10] worktype           
## 
## Root node error: 7930/31721 = 0.24999
## 
## n= 31721 
## 
##             CP nsplit rel error  xerror      xstd
## 1  0.109079445      0   1.00000 1.00000 0.0097252
## 2  0.033228247      2   0.78184 0.78184 0.0089063
## 3  0.029760404      4   0.71538 0.72585 0.0086558
## 4  0.024085750      5   0.68562 0.68562 0.0084641
## 5  0.019672131      6   0.66154 0.66154 0.0083442
## 6  0.002143758      7   0.64187 0.64124 0.0082401
## 7  0.001471206     10   0.63480 0.63758 0.0082211
## 8  0.001134931     14   0.62812 0.63367 0.0082006
## 9  0.000945776     15   0.62699 0.63443 0.0082046
## 10 0.000882724     21   0.62106 0.63380 0.0082012
## 11 0.000798655     22   0.62018 0.63253 0.0081946
## 12 0.000756620     30   0.61248 0.63165 0.0081900
## 13 0.000693569     32   0.61097 0.63052 0.0081840
## 14 0.000630517     35   0.60845 0.63001 0.0081813
## 15 0.000567465     36   0.60782 0.63077 0.0081853
## 16 0.000535939     40   0.60555 0.63039 0.0081833
## 17 0.000504414     44   0.60340 0.63039 0.0081833
## 18 0.000441362     48   0.60126 0.62989 0.0081807
## 19 0.000409836     50   0.60038 0.63102 0.0081866
## 20 0.000378310     61   0.59445 0.63329 0.0081986
## 21 0.000353090     88   0.58398 0.63607 0.0082131
## 22 0.000336276     93   0.58222 0.63783 0.0082224
## 23 0.000315259    100   0.57970 0.63796 0.0082230
## 24 0.000277427    114   0.57503 0.63808 0.0082237
## 25 0.000252207    120   0.57327 0.63834 0.0082250
## 26 0.000236444    148   0.56608 0.64161 0.0082421
## 27 0.000210172    159   0.56305 0.64174 0.0082427
## 28 0.000189155    165   0.56179 0.64275 0.0082480
## 29 0.000176545    192   0.55649 0.64603 0.0082649
## 30 0.000168138    197   0.55561 0.64817 0.0082760
## 31 0.000162133    203   0.55460 0.65006 0.0082857
## 32 0.000151324    210   0.55347 0.65044 0.0082877
## 33 0.000144118    215   0.55271 0.65044 0.0082877
## 34 0.000126103    225   0.55095 0.65233 0.0082974
## 35 0.000112092    319   0.53821 0.65347 0.0083032
## 36 0.000100883    328   0.53720 0.65561 0.0083141
## 37 0.000094578    338   0.53619 0.65675 0.0083199
## 38 0.000084069    358   0.53430 0.66028 0.0083378
## 39 0.000075662    374   0.53291 0.66267 0.0083500
## 40 0.000063052    379   0.53253 0.66772 0.0083753
## 41 0.000054044    421   0.52989 0.66885 0.0083810
## 42 0.000050441    428   0.52951 0.67087 0.0083911
## 43 0.000042034    441   0.52875 0.67503 0.0084118
## 44 0.000036030    453   0.52825 0.67491 0.0084112
## 45 0.000031526    460   0.52799 0.67881 0.0084306
## 46 0.000025221    464   0.52787 0.67881 0.0084306
## 47 0.000010000    484   0.52724 0.68247 0.0084486
plotcp(ar)

Parece que minimiza aprox en 0.0004 de complejidad, fila 18 con un error de 0.62989. Generamos un nuevo arbol con ese parametro. Ademas vamos a incluir un nuevo paramtero para que el arbol no tenga mas de 7 niveles.

ar <- rpart(formula_ar, train, method = 'class', parms = list(
  split = "information"), 
  control = rpart.control(cp = 0.0004,maxdepth = 7))

Revisamos de nuevo la complejidad.

printcp(ar)
## 
## Classification tree:
## rpart(formula = formula_ar, data = train, method = "class", parms = list(split = "information"), 
##     control = rpart.control(cp = 0.0004, maxdepth = 7))
## 
## Variables actually used in tree construction:
## [1] age_disc            c_status            educational_level  
## [4] fnlwgt_disc         gain_disc           hours_per_week_disc
## [7] race                workclass           worktype           
## 
## Root node error: 7930/31721 = 0.24999
## 
## n= 31721 
## 
##            CP nsplit rel error  xerror      xstd
## 1  0.10907945      0   1.00000 1.00000 0.0097252
## 2  0.03322825      2   0.78184 0.78184 0.0089063
## 3  0.02976040      4   0.71538 0.72850 0.0086681
## 4  0.02408575      5   0.68562 0.68562 0.0084641
## 5  0.01967213      6   0.66154 0.66154 0.0083442
## 6  0.00214376      7   0.64187 0.64351 0.0082519
## 7  0.00147121      9   0.63758 0.63909 0.0082290
## 8  0.00094578     12   0.63317 0.63821 0.0082244
## 9  0.00053594     16   0.62938 0.64073 0.0082375
## 10 0.00050441     20   0.62724 0.64035 0.0082355
## 11 0.00040000     21   0.62673 0.64035 0.0082355
plotcp(ar)

Ahora parece bastante estable, ya que el error es solo decreciente y no vuelve a crecer. Habría que variar el parametro cp para ver que luego volvería a crecer el error y comprobar que este es el mejor corte.

5.3.4.2. Visualización del arbol

Vamos a crear el grafico del arbol para analizarlo.

rpart.plot(ar,type=2,extra = 7, under = TRUE,under.cex = 0.7,fallen.leaves=F,gap = 0,cex=0.2,yesno = 2,box.palette = "GnYlRd",branch.lty = 3)

Otra visualizacion alternativa

prp(ar, type=2, extra=104, nn=TRUE, fallen.leaves=F,
faclen=0, varlen=0, shadow.col="grey", branch.lty=3)

5.3.4.3. Extracción de las reglas del modelo

Vamos a sacar las reglas que podrian ser utilizadas por ejemplo para hacer una implantacion del arbol. Normalmente los data scientist generan el modelo y los insights con un lenguaje analitico como R, y luego un equipo de data engineers generan esas reglas en C++, Java o SQL para que se ejecuten sobre una base de datos, donde tendremos los datos de origen.

reglas <- rpart.rules(ar,style = 'tall',cover = T)
#style sirve para que la salida sea mas legible y cover añade el % de casos e los que aplica la regla

Podemos llevarnos el nodo final de cada cliente a un data.frame para poder hacer una explotacion posterior.

#Para ello usarmos el predict especficio de rpart y con el parametro nn
ar_numnodos <- rpart.predict(ar,test,nn = T)
head(ar_numnodos)
##            0         1 nn
## 1 0.94894050 0.0510595  4
## 2 0.31894540 0.6810546 59
## 3 0.94894050 0.0510595  4
## 4 0.89095574 0.1090443 24
## 5 0.17125382 0.8287462 11
## 6 0.05109489 0.9489051 15

5.3.4.4. Evaluación del modelo

Vamos a calcular los scorings y evaluar el modelo

ar_predict <- predict(ar,test,type = 'prob')[,2]

Vemos que pinta tiene

plot(ar_predict~test$income)

Con la funcion umbrales probamos diferentes cortes

umb_ar<-umbrales(test$income,ar_predict)
umb_ar
##    umbral  acierto precision cobertura       F1
## 1    0.05 24.40986  24.40042  99.96948 39.22649
## 2    0.10 72.11259  46.30565  89.50259 61.03423
## 3    0.15 75.76141  50.19298  87.30546 63.74067
## 4    0.20 78.60600  53.97794  83.64358 65.61341
## 5    0.25 79.92404  56.21390  80.19530 66.09658
## 6    0.30 80.22191  56.76323  79.52395 66.24301
## 7    0.35 83.23777  63.80517  72.35276 67.81067
## 8    0.40 83.34202  64.21542  71.68142 67.74333
## 9    0.45 84.35475  70.53073  61.64175 65.78733
## 10   0.50 84.35475  70.53073  61.64175 65.78733
## 11   0.55 84.25795  71.41805  59.16997 64.71963
## 12   0.60 84.11646  72.50197  56.24046 63.34422
## 13   0.65 84.11646  72.50197  56.24046 63.34422
## 14   0.70 80.04319  91.40083  20.10986 32.96648
## 15   0.75 79.92404  91.67862  19.49954 32.15903
## 16   0.80 79.92404  91.67862  19.49954 32.15903
## 17   0.85 78.56877  94.63087  12.90815 22.71751
## 18   0.90 78.56877  94.63087  12.90815 22.71751
## 19   0.95  0.95000   0.95000   0.95000  0.95000

Seleccionamos automaticamente el mejor umbral

umbral_final_ar<-umb_ar[which.max(umb_ar$F1),1]
umbral_final_ar
## [1] 0.35

Evaluamos la matriz de confusion y las metricas con el umbral optimizado

confusion(test$income,ar_predict,umbral_final_ar)
##     
## real FALSE TRUE
##    0  8807 1345
##    1   906 2371
ar_metricas<-filter(umb_ar,umbral==umbral_final_ar)
ar_metricas
##   umbral  acierto precision cobertura       F1
## 1   0.35 83.23777  63.80517  72.35276 67.81067

Evaluamos la ROC

#creamos el objeto prediction
ar_prediction<-prediction(ar_predict,test$income)
#visualizamos la ROC
roc(ar_prediction)

Sacamos las metricas definitivas incluyendo el AUC

ar_metricas<-cbind(ar_metricas,AUC=round(auc(ar_prediction),2)*100)
print(t(ar_metricas))
##               [,1]
## umbral     0.35000
## acierto   83.23777
## precision 63.80517
## cobertura 72.35276
## F1        67.81067
## AUC       86.00000

5.3.5 - Modelizamos con Random Forest

5.3.5.1. Modelización

Creamos el modelo:

formula_rf <- formula
rf <- randomForest(formula_rf,train,importance=T)
rf
## 
## Call:
##  randomForest(formula = formula_rf, data = train, importance = T) 
##                Type of random forest: classification
##                      Number of trees: 500
## No. of variables tried at each split: 3
## 
##         OOB estimate of  error rate: 15.9%
## Confusion matrix:
##       0    1 class.error
## 0 21973 1818  0.07641545
## 1  3225 4705  0.40668348

Visualizamos las variables mas importantes:

varImpPlot(rf)

Como hay dos criterios (decrecimiento de la media Accuracy y decrecimiento de la media de Gini) vamos a crear una unica variable agregada y visualizarla para tener una mejor idea de la importancia de cada variable

importancia <- importance(rf)[,3:4]
#normalizamos para poner las dos variables en la misma scala. los valores negativos son las que menos predicen y los positivos las que mas
importancia_norm <- as.data.frame(scale(importancia))
#creamos una unica variable como suma de las otras
importancia_norm <- importancia_norm %>% mutate(
  Variable = rownames(importancia_norm),
  Imp_tot = MeanDecreaseAccuracy + MeanDecreaseGini) %>%
  mutate(Imp_tot = Imp_tot + abs(min(Imp_tot))) %>% 
  arrange(desc(Imp_tot)) %>% 
  select(Variable,Imp_tot,MeanDecreaseAccuracy,MeanDecreaseGini)
#hacemos un grafico para ver la curva de caida de importancia
ggplot(importancia_norm, aes(reorder(Variable,-Imp_tot),Imp_tot)) + geom_bar(stat = "identity") + theme(axis.text.x = element_text(angle = 90,size = 7))

importancia_norm
##               Variable   Imp_tot MeanDecreaseAccuracy MeanDecreaseGini
## 1             c_status 6.0980267           1.78869971        2.1791233
## 2            gain_disc 4.0722811           1.49699338        0.4450840
## 3    educational_level 3.1403525           0.51665560        0.4934932
## 4             age_disc 2.9101775          -0.02229203        0.8022658
## 5             worktype 1.9560700          -0.03848760       -0.1356462
## 6  hours_per_week_disc 1.1389420          -0.38885042       -0.6024113
## 7            workclass 0.9343396          -0.62743422       -0.5684299
## 8                  sex 0.6236535          -0.49038639       -1.0161638
## 9          fnlwgt_disc 0.4281942          -1.18034253       -0.5216670
## 10                race 0.0000000          -1.05455548       -1.0756482

La caida es bastante gradual, pero vemos que la variable race no influye practicamente y hay un salto grande desd e la anterior. Asi que las eliminamos y volvemos a calcular el random forest con el resto.

a_mantener <- importancia_norm %>% 
  filter(Imp_tot > 0.3) %>% 
  select(Variable)

#Extraemos los nombres como un vector
a_mantener <- as.character((a_mantener$Variable))

Creamos de nuevo el modelo con las nuevas variables

formula_rf <- reformulate(a_mantener,target)
rf <- randomForest(formula_rf,train,importance=T)
rf
## 
## Call:
##  randomForest(formula = formula_rf, data = train, importance = T) 
##                Type of random forest: classification
##                      Number of trees: 500
## No. of variables tried at each split: 3
## 
##         OOB estimate of  error rate: 15.89%
## Confusion matrix:
##       0    1 class.error
## 0 21982 1809  0.07603716
## 1  3232 4698  0.40756620

5.3.5.2. Aplicación al conjunto de test y visualizacion

Aplicamos el modelo al conjunto de test, generando un vector con las probabilidades Notar que por el metodo predict de randomforest hay que poner el type=prob para tener el scoring, lo cual nos dara una matriz que nos tenemos que quedar con la segunda columna

rf_predict <- predict(rf,test,type = 'prob')[,2]

Vemos que pinta tiene

plot(rf_predict~test$income)

5.3.5.3. Evaluación del modelo

Con la funcion umbrales probamos diferentes cortes

umb_rf<-umbrales(test$income,rf_predict)
umb_rf
##    umbral  acierto precision cobertura       F1
## 1    0.05 72.53705  46.72719  89.53311 61.40645
## 2    0.10 77.20605  51.99483  85.90174 64.77966
## 3    0.15 80.05808  56.29070  81.78212 66.68325
## 4    0.20 81.62931  59.31891  78.66951 67.63741
## 5    0.25 82.72396  61.99549  75.46536 68.07047
## 6    0.30 83.34202  64.12276  72.04760 67.85458
## 7    0.35 83.94519  66.39368  69.27067 67.80167
## 8    0.40 84.31008  68.13391  67.07354 67.59957
## 9    0.45 84.57815  70.04654  64.29661 67.04853
## 10   0.50 84.54837  71.47963  61.03143 65.84362
## 11   0.55 84.63028  73.42603  58.01038 64.81418
## 12   0.60 84.61538  74.86653  55.63015 63.83053
## 13   0.65 84.31752  76.15007  52.02930 61.82016
## 14   0.70 84.01966  77.65281  48.45896 59.67681
## 15   0.75 83.46117  79.04290  43.85108 56.40824
## 16   0.80 82.89523  80.66333  39.33476 52.88205
## 17   0.85 82.01653  82.11624  33.62832 47.71596
## 18   0.90 81.23464  84.00719  28.53219 42.59681
## 19   0.95 79.57406  85.98383  19.46903 31.74919

Seleccionamos automaticamente el mejor umbral

umbral_final_rf<-umb_rf[which.max(umb_rf$F1),1]
umbral_final_rf
## [1] 0.25

Evaluamos la matriz de confusion y las metricas con el umbral optimizado

confusion(test$income,rf_predict,umbral_final_rf)
##     
## real FALSE TRUE
##    0  8636 1516
##    1   804 2473
rf_metricas <- filter(umb_rf,umbral==umbral_final_rf)
rf_metricas
##   umbral  acierto precision cobertura       F1
## 1   0.25 82.72396  61.99549  75.46536 68.07047

Evaluamos la ROC

#creamos el objeto prediction
rf_prediction <- prediction(rf_predict,test$income)
#visualizamos la ROC
roc(rf_prediction)

Sacamos las metricas definitivas incluyendo el AUC

rf_metricas<-cbind(rf_metricas,AUC=round(auc(rf_prediction),2)*100)
print(t(rf_metricas))
##               [,1]
## umbral     0.25000
## acierto   82.72396
## precision 61.99549
## cobertura 75.46536
## F1        68.07047
## AUC       88.00000

5.3.6. - Comparamos los 3 metodos

comparativa <- rbind(rl_metricas,ar_metricas,rf_metricas)
rownames(comparativa) <- c('Regresion Logistica','Arbol Decision','Random Forest')
t(comparativa) #t simplemente transpone para leerlo mejor
##           Regresion Logistica Arbol Decision Random Forest
## umbral                0.30000        0.35000       0.25000
## acierto              82.24738       83.23777      82.72396
## precision            60.20571       63.80517      61.99549
## cobertura            80.37839       72.35276      75.46536
## F1                   68.84475       67.81067      68.07047
## AUC                  90.00000       86.00000      88.00000

Conclusion: todos serían practicamente igual predictivos, vamos a quedarnos con la regresion logistica por tener algo mejor al AUC(90).

5.3.7. Realizando una cross-validation de los resultados

5.3.7.1. Haciendo las particiones de los datos

Vamos a particionar los datos en 4 datasets iguales, para ir rotando los datasets de entrenamiento y de validación. Esto nos dará un dataset de entrenamiento en cada iteración del 75% de las observaciones y del 25% para el dataset de test.

inTrain <- createDataPartition(df$income,p=0.5,list=FALSE)
df_temp_1 <- df[inTrain,]
df_temp_2 <- df[-inTrain,]

inTrain <- createDataPartition(df_temp_1$income,p=0.5,list=FALSE)
df1 <- df_temp_1[inTrain,]
df2 <- df_temp_1[-inTrain,]

inTrain <- createDataPartition(df_temp_2$income,p=0.5,list=FALSE)
df3 <- df_temp_2[inTrain,]
df4 <- df_temp_2[-inTrain,]

rm(list=c("df_temp_1","df_temp_2","inTrain"))

Ahora hay que hacer 4 modelos de cada tipo, donde 3 de las particiones serán la muestra de entrenamiento y la particion restante será la muestra de test. Y la muestra de test va rotando.

Vamos a realizar la cross-validation solo con la regresión logistica, que en la comparativa anterior ha salido como el mejor modelo. Si se confirman las metricas, lo dejaremos con la elección de la regresión logistica. Si hay mucha variabilidad, haremos lo mismo con arboles de decisión y con random forest, para elegir el mejor modelo.

5.3.7.2. Realizando los modelos con regresión logistica

Modelo 1: df1, df2 y df3 serán la muestra de entrenamiento, y df4 la muestra de test.

train <- rbind(df1,df2,df3)
test <- df4

formula_rl_x <- formula
rl_x1 <- glm(formula_rl_x,train,family=binomial(link='logit'))
summary(rl_x1)
## 
## Call:
## glm(formula = formula_rl_x, family = binomial(link = "logit"), 
##     data = train)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -3.2277  -0.5297  -0.2195  -0.0200   3.5142  
## 
## Coefficients:
##                                  Estimate Std. Error z value
## (Intercept)                      -7.89676    0.46435 -17.006
## fnlwgt_disc02_DE_65K_A_120K       0.16950    0.06712   2.525
## fnlwgt_disc03_DE_105K_A_180K      0.33365    0.06314   5.284
## fnlwgt_disc04_DE_160K_A_320K      0.35092    0.06032   5.818
## fnlwgt_disc05_MAS_320K            0.44102    0.07439   5.929
## workclassLocal-gov               -0.67612    0.10324  -6.549
## workclassPrivate                 -0.55190    0.08694  -6.348
## workclassSelf-emp-inc            -0.30037    0.11416  -2.631
## workclassSelf-emp-not-inc        -0.97807    0.10080  -9.703
## workclassState-gov               -0.83246    0.11521  -7.225
## workclassWithout-pay             -0.61845    0.94324  -0.656
## sexMale                           0.13931    0.04932   2.824
## raceAsian-Pac-Islander            0.55133    0.22267   2.476
## raceBlack                         0.11477    0.21180   0.542
## raceOther                         0.09936    0.31052   0.320
## raceWhite                         0.47149    0.20136   2.342
## c_status02_with_children         -0.71265    0.20826  -3.422
## c_status03_unmarried             -0.15809    0.17143  -0.922
## c_status04_married                2.06566    0.16956  12.183
## worktype02_Executives             0.71092    0.07146   9.949
## worktype03_Technical              0.30477    0.07321   4.163
## worktype04_services               0.41367    0.06809   6.076
## worktype05_manual_work           -0.25228    0.06912  -3.650
## educational_level02_high_school   0.97666    0.07813  12.500
## educational_level03_professional  1.49104    0.08012  18.609
## educational_level04_bachelor      2.06767    0.08387  24.654
## educational_level05_master        2.01437    0.08733  23.066
## educational_level06_doctor        3.15549    0.16296  19.364
## age_disc02 <= 24                  0.85388    0.39272   2.174
## age_disc03 <= 27                  1.61355    0.37656   4.285
## age_disc04 <= 29                  2.09429    0.37576   5.573
## age_disc05 <= 33                  2.35046    0.37099   6.336
## age_disc06 <= 35                  2.46748    0.37343   6.608
## age_disc07 <= 38                  2.80419    0.37138   7.551
## age_disc08 <= 41                  2.85363    0.37159   7.680
## age_disc09 <= 45                  2.94837    0.37097   7.948
## age_disc10 <= 54                  3.10136    0.36976   8.387
## age_disc11 <= 61                  2.98366    0.37188   8.023
## age_disc12 > 61                   2.53824    0.37540   6.761
## hours_per_week_disc02 <= 39       0.66681    0.10002   6.667
## hours_per_week_disc03 <= 43       0.85427    0.07398  11.547
## hours_per_week_disc04 > 43        1.32864    0.07524  17.660
## gain_disc02 > 4650                3.24051    0.09124  35.516
##                                              Pr(>|z|)    
## (Intercept)                      < 0.0000000000000002 ***
## fnlwgt_disc02_DE_65K_A_120K                  0.011558 *  
## fnlwgt_disc03_DE_105K_A_180K      0.00000012621561252 ***
## fnlwgt_disc04_DE_160K_A_320K      0.00000000596374496 ***
## fnlwgt_disc05_MAS_320K            0.00000000305683123 ***
## workclassLocal-gov                0.00000000005800315 ***
## workclassPrivate                  0.00000000021768957 ***
## workclassSelf-emp-inc                        0.008507 ** 
## workclassSelf-emp-not-inc        < 0.0000000000000002 ***
## workclassState-gov                0.00000000000049963 ***
## workclassWithout-pay                         0.512037    
## sexMale                                      0.004736 ** 
## raceAsian-Pac-Islander                       0.013287 *  
## raceBlack                                    0.587916    
## raceOther                                    0.748983    
## raceWhite                                    0.019204 *  
## c_status02_with_children                     0.000622 ***
## c_status03_unmarried                         0.356409    
## c_status04_married               < 0.0000000000000002 ***
## worktype02_Executives            < 0.0000000000000002 ***
## worktype03_Technical              0.00003142911051192 ***
## worktype04_services               0.00000000123387970 ***
## worktype05_manual_work                       0.000262 ***
## educational_level02_high_school  < 0.0000000000000002 ***
## educational_level03_professional < 0.0000000000000002 ***
## educational_level04_bachelor     < 0.0000000000000002 ***
## educational_level05_master       < 0.0000000000000002 ***
## educational_level06_doctor       < 0.0000000000000002 ***
## age_disc02 <= 24                             0.029686 *  
## age_disc03 <= 27                  0.00001827693792117 ***
## age_disc04 <= 29                  0.00000002497134275 ***
## age_disc05 <= 33                  0.00000000023644700 ***
## age_disc06 <= 35                  0.00000000003907686 ***
## age_disc07 <= 38                  0.00000000000004330 ***
## age_disc08 <= 41                  0.00000000000001596 ***
## age_disc09 <= 45                  0.00000000000000190 ***
## age_disc10 <= 54                 < 0.0000000000000002 ***
## age_disc11 <= 61                  0.00000000000000103 ***
## age_disc12 > 61                   0.00000000001366551 ***
## hours_per_week_disc02 <= 39       0.00000000002616725 ***
## hours_per_week_disc03 <= 43      < 0.0000000000000002 ***
## hours_per_week_disc04 > 43       < 0.0000000000000002 ***
## gain_disc02 > 4650               < 0.0000000000000002 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 37953  on 33863  degrees of freedom
## Residual deviance: 22739  on 33821  degrees of freedom
## AIC: 22825
## 
## Number of Fisher Scoring iterations: 8
#cogemos solo las variables relevantes segun el modelo anterior
a_mantener <- c(
  'workclass',
  'sex',
  'c_status',
  'worktype',
  'educational_level',
  'age_disc',
  'hours_per_week_disc',
  'gain_disc',
  'fnlwgt_disc'
)

formula_rl_x <- reformulate(a_mantener,target)
rl_x1 <- glm(formula_rl_x,train,family=binomial(link='logit'))
summary(rl_x1)
## 
## Call:
## glm(formula = formula_rl_x, family = binomial(link = "logit"), 
##     data = train)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -3.2246  -0.5313  -0.2223  -0.0201   3.5369  
## 
## Coefficients:
##                                  Estimate Std. Error z value
## (Intercept)                      -7.52943    0.42224 -17.832
## workclassLocal-gov               -0.65438    0.10269  -6.372
## workclassPrivate                 -0.51789    0.08630  -6.001
## workclassSelf-emp-inc            -0.25818    0.11366  -2.271
## workclassSelf-emp-not-inc        -0.93560    0.10016  -9.341
## workclassState-gov               -0.81576    0.11479  -7.107
## workclassWithout-pay             -0.55544    0.94475  -0.588
## sexMale                           0.15572    0.04918   3.166
## c_status02_with_children         -0.69423    0.20816  -3.335
## c_status03_unmarried             -0.13991    0.17115  -0.817
## c_status04_married                2.09623    0.16929  12.382
## worktype02_Executives             0.71646    0.07136  10.040
## worktype03_Technical              0.31125    0.07312   4.257
## worktype04_services               0.41174    0.06796   6.059
## worktype05_manual_work           -0.25415    0.06900  -3.683
## educational_level02_high_school   0.98881    0.07804  12.671
## educational_level03_professional  1.50634    0.08004  18.819
## educational_level04_bachelor      2.09510    0.08372  25.025
## educational_level05_master        2.04187    0.08718  23.422
## educational_level06_doctor        3.19721    0.16247  19.679
## age_disc02 <= 24                  0.84797    0.39280   2.159
## age_disc03 <= 27                  1.60631    0.37659   4.265
## age_disc04 <= 29                  2.08127    0.37578   5.539
## age_disc05 <= 33                  2.34202    0.37101   6.313
## age_disc06 <= 35                  2.45834    0.37344   6.583
## age_disc07 <= 38                  2.79751    0.37140   7.532
## age_disc08 <= 41                  2.84270    0.37160   7.650
## age_disc09 <= 45                  2.93388    0.37098   7.908
## age_disc10 <= 54                  3.09074    0.36979   8.358
## age_disc11 <= 61                  2.97866    0.37190   8.009
## age_disc12 > 61                   2.53884    0.37544   6.762
## hours_per_week_disc02 <= 39       0.65858    0.09991   6.592
## hours_per_week_disc03 <= 43       0.84965    0.07390  11.497
## hours_per_week_disc04 > 43        1.33545    0.07516  17.768
## gain_disc02 > 4650                3.23135    0.09114  35.455
## fnlwgt_disc02_DE_65K_A_120K       0.17705    0.06693   2.645
## fnlwgt_disc03_DE_105K_A_180K      0.33186    0.06300   5.268
## fnlwgt_disc04_DE_160K_A_320K      0.34250    0.06013   5.696
## fnlwgt_disc05_MAS_320K            0.42003    0.07399   5.677
##                                              Pr(>|z|)    
## (Intercept)                      < 0.0000000000000002 ***
## workclassLocal-gov                0.00000000018640013 ***
## workclassPrivate                  0.00000000195692929 ***
## workclassSelf-emp-inc                        0.023121 *  
## workclassSelf-emp-not-inc        < 0.0000000000000002 ***
## workclassState-gov                0.00000000000118743 ***
## workclassWithout-pay                         0.556587    
## sexMale                                      0.001544 ** 
## c_status02_with_children                     0.000853 ***
## c_status03_unmarried                         0.413646    
## c_status04_married               < 0.0000000000000002 ***
## worktype02_Executives            < 0.0000000000000002 ***
## worktype03_Technical              0.00002074850239553 ***
## worktype04_services               0.00000000137321607 ***
## worktype05_manual_work                       0.000230 ***
## educational_level02_high_school  < 0.0000000000000002 ***
## educational_level03_professional < 0.0000000000000002 ***
## educational_level04_bachelor     < 0.0000000000000002 ***
## educational_level05_master       < 0.0000000000000002 ***
## educational_level06_doctor       < 0.0000000000000002 ***
## age_disc02 <= 24                             0.030867 *  
## age_disc03 <= 27                  0.00001995368049742 ***
## age_disc04 <= 29                  0.00000003049981545 ***
## age_disc05 <= 33                  0.00000000027431330 ***
## age_disc06 <= 35                  0.00000000004610593 ***
## age_disc07 <= 38                  0.00000000000004984 ***
## age_disc08 <= 41                  0.00000000000002012 ***
## age_disc09 <= 45                  0.00000000000000261 ***
## age_disc10 <= 54                 < 0.0000000000000002 ***
## age_disc11 <= 61                  0.00000000000000115 ***
## age_disc12 > 61                   0.00000000001357418 ***
## hours_per_week_disc02 <= 39       0.00000000004349183 ***
## hours_per_week_disc03 <= 43      < 0.0000000000000002 ***
## hours_per_week_disc04 > 43       < 0.0000000000000002 ***
## gain_disc02 > 4650               < 0.0000000000000002 ***
## fnlwgt_disc02_DE_65K_A_120K                  0.008162 ** 
## fnlwgt_disc03_DE_105K_A_180K      0.00000013813808693 ***
## fnlwgt_disc04_DE_160K_A_320K      0.00000001224328274 ***
## fnlwgt_disc05_MAS_320K            0.00000001374434400 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 37953  on 33863  degrees of freedom
## Residual deviance: 22774  on 33825  degrees of freedom
## AIC: 22852
## 
## Number of Fisher Scoring iterations: 8
#calculamos el pseudo R cuadrado
pr2_rl_x1 <- 1 -(rl_x1$deviance / rl_x1$null.deviance)
pr2_rl_x1
## [1] 0.3999542
#aplicamos el modelo al conjunto de test
rl_predict_x1 <- predict(rl_x1,test,type = 'response')

#hacemos una primera visualización de como es el modelo
plot(rl_predict_x1~test$income)

#Calculamos los umbrales y cual es el optimo segun la F1
umb_rl_x1<-umbrales(test$income,rl_predict_x1)
umb_rl_x1
##    umbral  acierto precision cobertura       F1
## 1    0.05 59.73773  37.97101  98.21492 54.76807
## 2    0.10 69.88304  44.95784  95.18029 61.06975
## 3    0.15 75.58923  50.45133  91.78865 65.11333
## 4    0.20 78.54864  54.15573  88.39700 67.16398
## 5    0.25 80.56885  57.34654  84.71974 68.39602
## 6    0.30 82.64221  61.57229  79.97144 69.57602
## 7    0.35 83.28903  63.85104  75.29454 69.10223
## 8    0.40 83.90927  66.76881  70.01071 68.35134
## 9    0.45 84.56495  70.06442  66.01214 67.97794
## 10   0.50 84.78646  72.79226  61.79936 66.84688
## 11   0.55 84.96367  76.28571  57.19386 65.37441
## 12   0.60 84.42318  78.57534  51.19600 61.99741
## 13   0.65 83.71434  81.20544  44.73402 57.68877
## 14   0.70 82.92575  84.57278  38.16494 52.59533
## 15   0.75 81.58781  86.18619  30.73902 45.31579
## 16   0.80 80.44480  89.49468  24.02713 37.88348
## 17   0.85 79.53216  91.82283  19.24313 31.81818
## 18   0.90 78.77902  93.37607  15.60157 26.73600
## 19   0.95 78.01701  95.97701  11.92431 21.21308
umbral_final_rl_x1<-umb_rl[which.max(umb_rl_x1$F1),1]
umbral_final_rl_x1
## [1] 0.3
#evaluamos la matriz de confusión y las métricas
confusion(test$income,rl_predict_x1,umbral_final_rl_x1)
##     
## real FALSE TRUE
##    0  7087 1398
##    1   561 2240
rl_metricas_x1 <- filter(umb_rl_x1,umbral==umbral_final_rl_x1)
rl_metricas_x1
##   umbral  acierto precision cobertura       F1
## 1    0.3 82.64221  61.57229  79.97144 69.57602
#Calculamos la ROC y la AUC

#creamos el objeto prediction
rl_prediction_x1 <- prediction(rl_predict_x1,test$income)
#visualizamos la ROC
roc(rl_prediction_x1)

rl_metricas_x1 <- cbind(rl_metricas_x1,AUC=round(auc(rl_prediction_x1),2)*100)
print(t(rl_metricas_x1))
##               [,1]
## umbral     0.30000
## acierto   82.64221
## precision 61.57229
## cobertura 79.97144
## F1        69.57602
## AUC       90.00000

El resultado ha sido de una AUC de 90 con este modelo.

Modelo 2: df1, df2 y df4 serán la muestra de entrenamiento, y df3 la muestra de test.

train <- rbind(df1,df2,df4)
test <- df3

formula_rl_x <- formula
rl_x2 <- glm(formula_rl_x,train,family=binomial(link='logit'))
summary(rl_x2)
## 
## Call:
## glm(formula = formula_rl_x, family = binomial(link = "logit"), 
##     data = train)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -3.1815  -0.5273  -0.2219  -0.0159   3.9627  
## 
## Coefficients:
##                                  Estimate Std. Error z value
## (Intercept)                      -8.50919    0.54756 -15.540
## fnlwgt_disc02_DE_65K_A_120K       0.15300    0.06760   2.263
## fnlwgt_disc03_DE_105K_A_180K      0.35090    0.06358   5.519
## fnlwgt_disc04_DE_160K_A_320K      0.32453    0.06081   5.337
## fnlwgt_disc05_MAS_320K            0.44684    0.07457   5.992
## workclassLocal-gov               -0.56288    0.10460  -5.381
## workclassPrivate                 -0.48606    0.08843  -5.496
## workclassSelf-emp-inc            -0.21105    0.11519  -1.832
## workclassSelf-emp-not-inc        -0.84312    0.10219  -8.250
## workclassState-gov               -0.71173    0.11568  -6.152
## workclassWithout-pay             -0.77162    0.84043  -0.918
## sexMale                           0.16151    0.04946   3.265
## raceAsian-Pac-Islander            0.53781    0.23196   2.319
## raceBlack                         0.16497    0.22230   0.742
## raceOther                         0.43186    0.31511   1.371
## raceWhite                         0.58558    0.21216   2.760
## c_status02_with_children         -0.72543    0.21445  -3.383
## c_status03_unmarried             -0.11203    0.17704  -0.633
## c_status04_married                2.04570    0.17551  11.656
## worktype02_Executives             0.78907    0.07190  10.975
## worktype03_Technical              0.34292    0.07338   4.673
## worktype04_services               0.39526    0.06859   5.763
## worktype05_manual_work           -0.25657    0.06962  -3.685
## educational_level02_high_school   0.95448    0.07813  12.216
## educational_level03_professional  1.49880    0.08004  18.725
## educational_level04_bachelor      2.04451    0.08360  24.456
## educational_level05_master        1.98428    0.08733  22.721
## educational_level06_doctor        2.98666    0.15712  19.009
## age_disc02 <= 24                  1.30540    0.48046   2.717
## age_disc03 <= 27                  1.97243    0.46869   4.208
## age_disc04 <= 29                  2.49444    0.46764   5.334
## age_disc05 <= 33                  2.81692    0.46372   6.075
## age_disc06 <= 35                  2.96128    0.46570   6.359
## age_disc07 <= 38                  3.20861    0.46407   6.914
## age_disc08 <= 41                  3.27841    0.46422   7.062
## age_disc09 <= 45                  3.40088    0.46364   7.335
## age_disc10 <= 54                  3.52915    0.46274   7.627
## age_disc11 <= 61                  3.44737    0.46440   7.423
## age_disc12 > 61                   2.94915    0.46714   6.313
## hours_per_week_disc02 <= 39       0.70117    0.10058   6.972
## hours_per_week_disc03 <= 43       0.85636    0.07482  11.446
## hours_per_week_disc04 > 43        1.24784    0.07606  16.405
## gain_disc02 > 4650                3.23561    0.09075  35.652
##                                              Pr(>|z|)    
## (Intercept)                      < 0.0000000000000002 ***
## fnlwgt_disc02_DE_65K_A_120K                  0.023620 *  
## fnlwgt_disc03_DE_105K_A_180K       0.0000000340518368 ***
## fnlwgt_disc04_DE_160K_A_320K       0.0000000944022305 ***
## fnlwgt_disc05_MAS_320K             0.0000000020705806 ***
## workclassLocal-gov                 0.0000000739673568 ***
## workclassPrivate                   0.0000000387769993 ***
## workclassSelf-emp-inc                        0.066921 .  
## workclassSelf-emp-not-inc        < 0.0000000000000002 ***
## workclassState-gov                 0.0000000007628404 ***
## workclassWithout-pay                         0.358555    
## sexMale                                      0.001094 ** 
## raceAsian-Pac-Islander                       0.020418 *  
## raceBlack                                    0.458011    
## raceOther                                    0.170527    
## raceWhite                                    0.005778 ** 
## c_status02_with_children                     0.000717 ***
## c_status03_unmarried                         0.526874    
## c_status04_married               < 0.0000000000000002 ***
## worktype02_Executives            < 0.0000000000000002 ***
## worktype03_Technical               0.0000029653357086 ***
## worktype04_services                0.0000000082667299 ***
## worktype05_manual_work                       0.000228 ***
## educational_level02_high_school  < 0.0000000000000002 ***
## educational_level03_professional < 0.0000000000000002 ***
## educational_level04_bachelor     < 0.0000000000000002 ***
## educational_level05_master       < 0.0000000000000002 ***
## educational_level06_doctor       < 0.0000000000000002 ***
## age_disc02 <= 24                             0.006588 ** 
## age_disc03 <= 27                   0.0000257243176351 ***
## age_disc04 <= 29                   0.0000000960110784 ***
## age_disc05 <= 33                   0.0000000012426978 ***
## age_disc06 <= 35                   0.0000000002034108 ***
## age_disc07 <= 38                   0.0000000000047085 ***
## age_disc08 <= 41                   0.0000000000016393 ***
## age_disc09 <= 45                   0.0000000000002213 ***
## age_disc10 <= 54                   0.0000000000000241 ***
## age_disc11 <= 61                   0.0000000000001142 ***
## age_disc12 > 61                    0.0000000002732633 ***
## hours_per_week_disc02 <= 39        0.0000000000031342 ***
## hours_per_week_disc03 <= 43      < 0.0000000000000002 ***
## hours_per_week_disc04 > 43       < 0.0000000000000002 ***
## gain_disc02 > 4650               < 0.0000000000000002 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 37950  on 33861  degrees of freedom
## Residual deviance: 22755  on 33819  degrees of freedom
## AIC: 22841
## 
## Number of Fisher Scoring iterations: 8
#cogemos solo las variables relevantes segun el modelo anterior
a_mantener <- c(
  'workclass',
  'sex',
  'c_status',
  'worktype',
  'educational_level',
  'age_disc',
  'hours_per_week_disc',
  'gain_disc',
  'fnlwgt_disc'
)

formula_rl_x <- reformulate(a_mantener,target)
rl_x2 <- glm(formula_rl_x,train,family=binomial(link='logit'))
summary(rl_x2)
## 
## Call:
## glm(formula = formula_rl_x, family = binomial(link = "logit"), 
##     data = train)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -3.1773  -0.5302  -0.2246  -0.0161   3.8762  
## 
## Coefficients:
##                                  Estimate Std. Error z value
## (Intercept)                      -8.03775    0.50759 -15.835
## workclassLocal-gov               -0.54275    0.10403  -5.217
## workclassPrivate                 -0.44538    0.08773  -5.077
## workclassSelf-emp-inc            -0.16354    0.11466  -1.426
## workclassSelf-emp-not-inc        -0.79361    0.10149  -7.819
## workclassState-gov               -0.68442    0.11512  -5.945
## workclassWithout-pay             -0.70889    0.83985  -0.844
## sexMale                           0.17774    0.04931   3.605
## c_status02_with_children         -0.70777    0.21433  -3.302
## c_status03_unmarried             -0.08787    0.17682  -0.497
## c_status04_married                2.08331    0.17532  11.883
## worktype02_Executives             0.79161    0.07179  11.027
## worktype03_Technical              0.34946    0.07326   4.770
## worktype04_services               0.39142    0.06844   5.719
## worktype05_manual_work           -0.25865    0.06946  -3.724
## educational_level02_high_school   0.96503    0.07806  12.363
## educational_level03_professional  1.51165    0.07997  18.902
## educational_level04_bachelor      2.06636    0.08349  24.750
## educational_level05_master        2.00794    0.08718  23.032
## educational_level06_doctor        3.02027    0.15688  19.253
## age_disc02 <= 24                  1.30061    0.48031   2.708
## age_disc03 <= 27                  1.96476    0.46853   4.193
## age_disc04 <= 29                  2.47991    0.46747   5.305
## age_disc05 <= 33                  2.80412    0.46354   6.049
## age_disc06 <= 35                  2.94673    0.46552   6.330
## age_disc07 <= 38                  3.19685    0.46390   6.891
## age_disc08 <= 41                  3.26348    0.46405   7.033
## age_disc09 <= 45                  3.38118    0.46346   7.296
## age_disc10 <= 54                  3.51458    0.46256   7.598
## age_disc11 <= 61                  3.43989    0.46422   7.410
## age_disc12 > 61                   2.94678    0.46697   6.310
## hours_per_week_disc02 <= 39       0.68604    0.10042   6.832
## hours_per_week_disc03 <= 43       0.84692    0.07471  11.336
## hours_per_week_disc04 > 43        1.25535    0.07597  16.525
## gain_disc02 > 4650                3.22787    0.09061  35.623
## fnlwgt_disc02_DE_65K_A_120K       0.15912    0.06738   2.361
## fnlwgt_disc03_DE_105K_A_180K      0.35114    0.06338   5.540
## fnlwgt_disc04_DE_160K_A_320K      0.31750    0.06056   5.243
## fnlwgt_disc05_MAS_320K            0.42384    0.07412   5.718
##                                              Pr(>|z|)    
## (Intercept)                      < 0.0000000000000002 ***
## workclassLocal-gov                0.00000018157091337 ***
## workclassPrivate                  0.00000038382497891 ***
## workclassSelf-emp-inc                        0.153788    
## workclassSelf-emp-not-inc         0.00000000000000531 ***
## workclassState-gov                0.00000000276349440 ***
## workclassWithout-pay                         0.398634    
## sexMale                                      0.000312 ***
## c_status02_with_children                     0.000959 ***
## c_status03_unmarried                         0.619223    
## c_status04_married               < 0.0000000000000002 ***
## worktype02_Executives            < 0.0000000000000002 ***
## worktype03_Technical              0.00000184225791073 ***
## worktype04_services               0.00000001070443840 ***
## worktype05_manual_work                       0.000196 ***
## educational_level02_high_school  < 0.0000000000000002 ***
## educational_level03_professional < 0.0000000000000002 ***
## educational_level04_bachelor     < 0.0000000000000002 ***
## educational_level05_master       < 0.0000000000000002 ***
## educational_level06_doctor       < 0.0000000000000002 ***
## age_disc02 <= 24                             0.006772 ** 
## age_disc03 <= 27                  0.00002747319383468 ***
## age_disc04 <= 29                  0.00000011269071988 ***
## age_disc05 <= 33                  0.00000000145446898 ***
## age_disc06 <= 35                  0.00000000024515518 ***
## age_disc07 <= 38                  0.00000000000552813 ***
## age_disc08 <= 41                  0.00000000000202721 ***
## age_disc09 <= 45                  0.00000000000029738 ***
## age_disc10 <= 54                  0.00000000000003006 ***
## age_disc11 <= 61                  0.00000000000012623 ***
## age_disc12 > 61                   0.00000000027837500 ***
## hours_per_week_disc02 <= 39       0.00000000000839329 ***
## hours_per_week_disc03 <= 43      < 0.0000000000000002 ***
## hours_per_week_disc04 > 43       < 0.0000000000000002 ***
## gain_disc02 > 4650               < 0.0000000000000002 ***
## fnlwgt_disc02_DE_65K_A_120K                  0.018202 *  
## fnlwgt_disc03_DE_105K_A_180K      0.00000003016939851 ***
## fnlwgt_disc04_DE_160K_A_320K      0.00000015828553913 ***
## fnlwgt_disc05_MAS_320K            0.00000001076876133 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 37950  on 33861  degrees of freedom
## Residual deviance: 22799  on 33823  degrees of freedom
## AIC: 22877
## 
## Number of Fisher Scoring iterations: 8
#calculamos el pseudo R cuadrado
pr2_rl_x2 <- 1 -(rl_x2$deviance / rl_x2$null.deviance)
pr2_rl_x2
## [1] 0.3992297
#aplicamos el modelo al conjunto de test
rl_predict_x2 <- predict(rl_x2,test,type = 'response')

#hacemos una primera visualización de como es el modelo
plot(rl_predict_x2~test$income)

#Calculamos los umbrales y cual es el optimo segun la F1
umb_rl_x2 <- umbrales(test$income,rl_predict_x2)
umb_rl_x2
##    umbral  acierto precision cobertura       F1
## 1    0.05 60.20553  38.22792  97.93005 54.98998
## 2    0.10 69.79093  44.90104  95.53890 61.09083
## 3    0.15 75.42523  50.27090  92.71949 65.19448
## 4    0.20 79.23459  55.09342  88.40114 67.88161
## 5    0.25 81.44933  58.82793  84.18986 69.26013
## 6    0.30 82.97307  62.45753  78.72948 69.65583
## 7    0.35 83.94755  65.81470  73.51892 69.45381
## 8    0.40 84.45252  68.64980  68.77231 68.71100
## 9    0.45 84.70943  71.62379  63.59743 67.37240
## 10   0.50 84.58540  73.91892  58.56531 65.35245
## 11   0.55 84.41708  76.78480  53.35475 62.96062
## 12   0.60 83.95641  79.83143  47.32334 59.42191
## 13   0.65 83.55776  83.26301  42.25553 56.06061
## 14   0.70 82.53898  85.97403  35.43897 50.18954
## 15   0.75 81.33416  89.17700  28.22984 42.88425
## 16   0.80 80.16478  91.45803  22.16274 35.67940
## 17   0.85 79.18143  93.62934  17.30906 29.21687
## 18   0.90 78.54359  95.23810  14.27552 24.82930
## 19   0.95 77.65769  95.75163  10.45682 18.85457
umbral_final_rl_x2 <- umb_rl[which.max(umb_rl_x2$F1),1]
umbral_final_rl_x2
## [1] 0.3
#evaluamos la matriz de confusión y las métricas
confusion(test$income,rl_predict_x2,umbral_final_rl_x2)
##     
## real FALSE TRUE
##    0  7160 1326
##    1   596 2206
rl_metricas_x2 <- filter(umb_rl_x2,umbral==umbral_final_rl_x2)
rl_metricas_x2
##   umbral  acierto precision cobertura       F1
## 1    0.3 82.97307  62.45753  78.72948 69.65583
#Calculamos la ROC y la AUC

#creamos el objeto prediction
rl_prediction_x2 <- prediction(rl_predict_x2,test$income)
#visualizamos la ROC
roc(rl_prediction_x2)

rl_metricas_x2 <- cbind(rl_metricas_x2,AUC=round(auc(rl_prediction_x2),2)*100)
print(t(rl_metricas_x2))
##               [,1]
## umbral     0.30000
## acierto   82.97307
## precision 62.45753
## cobertura 78.72948
## F1        69.65583
## AUC       90.00000

El resultado ha sido de una AUC de 90 con este modelo también.

Modelo 3: df1, df3 y df4 serán la muestra de entrenamiento, y df2 la muestra de test.

train <- rbind(df1,df3,df4)
test <- df2

formula_rl_x <- formula
rl_x3 <- glm(formula_rl_x,train,family=binomial(link='logit'))
summary(rl_x3)
## 
## Call:
## glm(formula = formula_rl_x, family = binomial(link = "logit"), 
##     data = train)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -3.1848  -0.5277  -0.2176  -0.0160   3.9440  
## 
## Coefficients:
##                                   Estimate Std. Error z value
## (Intercept)                       -8.50431    0.53847 -15.793
## fnlwgt_disc02_DE_65K_A_120K        0.07229    0.06736   1.073
## fnlwgt_disc03_DE_105K_A_180K       0.27923    0.06310   4.425
## fnlwgt_disc04_DE_160K_A_320K       0.20157    0.06042   3.336
## fnlwgt_disc05_MAS_320K             0.35797    0.07512   4.765
## workclassLocal-gov                -0.47444    0.10433  -4.547
## workclassPrivate                  -0.45129    0.08807  -5.125
## workclassSelf-emp-inc             -0.19750    0.11465  -1.723
## workclassSelf-emp-not-inc         -0.82788    0.10178  -8.134
## workclassState-gov                -0.67255    0.11605  -5.795
## workclassWithout-pay             -12.40721  114.61978  -0.108
## sexMale                            0.20235    0.04994   4.052
## raceAsian-Pac-Islander             0.58891    0.22339   2.636
## raceBlack                          0.26536    0.21236   1.250
## raceOther                          0.27380    0.31491   0.869
## raceWhite                          0.61049    0.20186   3.024
## c_status02_with_children          -0.55703    0.21923  -2.541
## c_status03_unmarried               0.01352    0.18489   0.073
## c_status04_married                 2.22644    0.18337  12.142
## worktype02_Executives              0.69896    0.07127   9.808
## worktype03_Technical               0.22682    0.07290   3.111
## worktype04_services                0.30345    0.06782   4.474
## worktype05_manual_work            -0.33053    0.06873  -4.809
## educational_level02_high_school    0.93821    0.07802  12.025
## educational_level03_professional   1.50783    0.07978  18.899
## educational_level04_bachelor       2.13068    0.08364  25.476
## educational_level05_master         1.97457    0.08749  22.569
## educational_level06_doctor         3.04237    0.16244  18.729
## age_disc02 <= 24                   1.25569    0.47194   2.661
## age_disc03 <= 27                   2.04399    0.45954   4.448
## age_disc04 <= 29                   2.43934    0.45955   5.308
## age_disc05 <= 33                   2.65123    0.45561   5.819
## age_disc06 <= 35                   2.89437    0.45759   6.325
## age_disc07 <= 38                   3.12310    0.45585   6.851
## age_disc08 <= 41                   3.28275    0.45590   7.201
## age_disc09 <= 45                   3.27842    0.45542   7.199
## age_disc10 <= 54                   3.48808    0.45443   7.676
## age_disc11 <= 61                   3.33535    0.45616   7.312
## age_disc12 > 61                    2.85582    0.45886   6.224
## hours_per_week_disc02 <= 39        0.60635    0.10112   5.996
## hours_per_week_disc03 <= 43        0.80639    0.07405  10.889
## hours_per_week_disc04 > 43         1.25424    0.07544  16.625
## gain_disc02 > 4650                 3.25193    0.09180  35.424
##                                              Pr(>|z|)    
## (Intercept)                      < 0.0000000000000002 ***
## fnlwgt_disc02_DE_65K_A_120K                  0.283198    
## fnlwgt_disc03_DE_105K_A_180K     0.000009628765569812 ***
## fnlwgt_disc04_DE_160K_A_320K                 0.000849 ***
## fnlwgt_disc05_MAS_320K           0.000001887504773538 ***
## workclassLocal-gov               0.000005431750469946 ***
## workclassPrivate                 0.000000298323426771 ***
## workclassSelf-emp-inc                        0.084970 .  
## workclassSelf-emp-not-inc        0.000000000000000415 ***
## workclassState-gov               0.000000006825239304 ***
## workclassWithout-pay                         0.913800    
## sexMale                          0.000050846205082087 ***
## raceAsian-Pac-Islander                       0.008383 ** 
## raceBlack                                    0.211452    
## raceOther                                    0.384610    
## raceWhite                                    0.002492 ** 
## c_status02_with_children                     0.011056 *  
## c_status03_unmarried                         0.941719    
## c_status04_married               < 0.0000000000000002 ***
## worktype02_Executives            < 0.0000000000000002 ***
## worktype03_Technical                         0.001863 ** 
## worktype04_services              0.000007662622683221 ***
## worktype05_manual_work           0.000001514323201394 ***
## educational_level02_high_school  < 0.0000000000000002 ***
## educational_level03_professional < 0.0000000000000002 ***
## educational_level04_bachelor     < 0.0000000000000002 ***
## educational_level05_master       < 0.0000000000000002 ***
## educational_level06_doctor       < 0.0000000000000002 ***
## age_disc02 <= 24                             0.007798 ** 
## age_disc03 <= 27                 0.000008671162674552 ***
## age_disc04 <= 29                 0.000000110792426849 ***
## age_disc05 <= 33                 0.000000005914738205 ***
## age_disc06 <= 35                 0.000000000252806310 ***
## age_disc07 <= 38                 0.000000000007327078 ***
## age_disc08 <= 41                 0.000000000000599371 ***
## age_disc09 <= 45                 0.000000000000608276 ***
## age_disc10 <= 54                 0.000000000000016454 ***
## age_disc11 <= 61                 0.000000000000263558 ***
## age_disc12 > 61                  0.000000000485627251 ***
## hours_per_week_disc02 <= 39      0.000000002019398814 ***
## hours_per_week_disc03 <= 43      < 0.0000000000000002 ***
## hours_per_week_disc04 > 43       < 0.0000000000000002 ***
## gain_disc02 > 4650               < 0.0000000000000002 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 37950  on 33861  degrees of freedom
## Residual deviance: 22589  on 33819  degrees of freedom
## AIC: 22675
## 
## Number of Fisher Scoring iterations: 12
#cogemos solo las variables relevantes segun el modelo anterior
a_mantener <- c(
  'workclass',
  'sex',
  'c_status',
  'worktype',
  'educational_level',
  'age_disc',
  'hours_per_week_disc',
  'gain_disc',
  'fnlwgt_disc'
)

formula_rl_x <- reformulate(a_mantener,target)
rl_x3 <- glm(formula_rl_x,train,family=binomial(link='logit'))
summary(rl_x3)
## 
## Call:
## glm(formula = formula_rl_x, family = binomial(link = "logit"), 
##     data = train)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -3.1816  -0.5278  -0.2195  -0.0167   3.8729  
## 
## Coefficients:
##                                   Estimate Std. Error z value
## (Intercept)                       -8.00865    0.50236 -15.942
## workclassLocal-gov                -0.45352    0.10379  -4.370
## workclassPrivate                  -0.41411    0.08736  -4.740
## workclassSelf-emp-inc             -0.14896    0.11407  -1.306
## workclassSelf-emp-not-inc         -0.78339    0.10108  -7.750
## workclassState-gov                -0.64704    0.11555  -5.600
## workclassWithout-pay             -12.36036  114.64172  -0.108
## sexMale                            0.21508    0.04980   4.318
## c_status02_with_children          -0.52821    0.21902  -2.412
## c_status03_unmarried               0.03957    0.18460   0.214
## c_status04_married                 2.26637    0.18306  12.381
## worktype02_Executives              0.70140    0.07119   9.853
## worktype03_Technical               0.23175    0.07281   3.183
## worktype04_services                0.30049    0.06769   4.439
## worktype05_manual_work            -0.33202    0.06861  -4.839
## educational_level02_high_school    0.94914    0.07796  12.175
## educational_level03_professional   1.52076    0.07972  19.076
## educational_level04_bachelor       2.15518    0.08351  25.808
## educational_level05_master         1.99842    0.08735  22.879
## educational_level06_doctor         3.08233    0.16229  18.993
## age_disc02 <= 24                   1.25395    0.47205   2.656
## age_disc03 <= 27                   2.03942    0.45962   4.437
## age_disc04 <= 29                   2.42862    0.45963   5.284
## age_disc05 <= 33                   2.64280    0.45567   5.800
## age_disc06 <= 35                   2.88267    0.45766   6.299
## age_disc07 <= 38                   3.11798    0.45592   6.839
## age_disc08 <= 41                   3.27457    0.45597   7.181
## age_disc09 <= 45                   3.26700    0.45550   7.172
## age_disc10 <= 54                   3.48107    0.45451   7.659
## age_disc11 <= 61                   3.33297    0.45624   7.305
## age_disc12 > 61                    2.85918    0.45895   6.230
## hours_per_week_disc02 <= 39        0.59364    0.10098   5.879
## hours_per_week_disc03 <= 43        0.79830    0.07395  10.795
## hours_per_week_disc04 > 43         1.25944    0.07537  16.711
## gain_disc02 > 4650                 3.24259    0.09157  35.411
## fnlwgt_disc02_DE_65K_A_120K        0.07732    0.06717   1.151
## fnlwgt_disc03_DE_105K_A_180K       0.27871    0.06295   4.428
## fnlwgt_disc04_DE_160K_A_320K       0.19636    0.06020   3.262
## fnlwgt_disc05_MAS_320K             0.34273    0.07469   4.589
##                                              Pr(>|z|)    
## (Intercept)                      < 0.0000000000000002 ***
## workclassLocal-gov                0.00001243444023726 ***
## workclassPrivate                  0.00000213411069910 ***
## workclassSelf-emp-inc                         0.19159    
## workclassSelf-emp-not-inc         0.00000000000000916 ***
## workclassState-gov                0.00000002149504638 ***
## workclassWithout-pay                          0.91414    
## sexMale                           0.00001571695807378 ***
## c_status02_with_children                      0.01588 *  
## c_status03_unmarried                          0.83029    
## c_status04_married               < 0.0000000000000002 ***
## worktype02_Executives            < 0.0000000000000002 ***
## worktype03_Technical                          0.00146 ** 
## worktype04_services               0.00000903919540601 ***
## worktype05_manual_work            0.00000130353103872 ***
## educational_level02_high_school  < 0.0000000000000002 ***
## educational_level03_professional < 0.0000000000000002 ***
## educational_level04_bachelor     < 0.0000000000000002 ***
## educational_level05_master       < 0.0000000000000002 ***
## educational_level06_doctor       < 0.0000000000000002 ***
## age_disc02 <= 24                              0.00790 ** 
## age_disc03 <= 27                  0.00000911322238670 ***
## age_disc04 <= 29                  0.00000012645419339 ***
## age_disc05 <= 33                  0.00000000664110970 ***
## age_disc06 <= 35                  0.00000000029997911 ***
## age_disc07 <= 38                  0.00000000000798384 ***
## age_disc08 <= 41                  0.00000000000068958 ***
## age_disc09 <= 45                  0.00000000000073698 ***
## age_disc10 <= 54                  0.00000000000001876 ***
## age_disc11 <= 61                  0.00000000000027662 ***
## age_disc12 > 61                   0.00000000046715678 ***
## hours_per_week_disc02 <= 39       0.00000000412925208 ***
## hours_per_week_disc03 <= 43      < 0.0000000000000002 ***
## hours_per_week_disc04 > 43       < 0.0000000000000002 ***
## gain_disc02 > 4650               < 0.0000000000000002 ***
## fnlwgt_disc02_DE_65K_A_120K                   0.24967    
## fnlwgt_disc03_DE_105K_A_180K      0.00000953240090393 ***
## fnlwgt_disc04_DE_160K_A_320K                  0.00111 ** 
## fnlwgt_disc05_MAS_320K            0.00000445781720487 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 37950  on 33861  degrees of freedom
## Residual deviance: 22625  on 33823  degrees of freedom
## AIC: 22703
## 
## Number of Fisher Scoring iterations: 12
#calculamos el pseudo R cuadrado
pr2_rl_x3 <- 1 -(rl_x3$deviance / rl_x3$null.deviance)
pr2_rl_x3
## [1] 0.4038343
#aplicamos el modelo al conjunto de test
rl_predict_x3 <- predict(rl_x3,test,type = 'response')

#hacemos una primera visualización de como es el modelo
plot(rl_predict_x3~test$income)

#Calculamos los umbrales y cual es el optimo segun la F1
umb_rl_x3 <- umbrales(test$income,rl_predict_x3)
umb_rl_x3
##    umbral  acierto precision cobertura       F1
## 1    0.05 60.56875  38.47337  98.21556 55.28880
## 2    0.10 70.10099  45.15474  95.28908 61.27367
## 3    0.15 75.44295  50.29389  91.61313 64.93802
## 4    0.20 78.46386  54.05642  88.22270 67.03729
## 5    0.25 80.37739  57.16028  83.61884 67.90320
## 6    0.30 82.06060  60.71133  78.58672 68.50210
## 7    0.35 82.91106  63.21526  74.51820 68.40295
## 8    0.40 83.27427  65.56540  68.70093 67.09655
## 9    0.45 83.70836  68.58356  63.41899 65.90024
## 10   0.50 84.08930  71.62511  59.45753 64.97660
## 11   0.55 84.21332  74.85380  54.81799 63.28801
## 12   0.60 83.66407  76.78971  49.00071 59.82571
## 13   0.65 83.00850  79.46667  42.54104 55.41609
## 14   0.70 82.19348  81.78170  36.36688 50.34585
## 15   0.75 81.43161  85.80122  30.19272 44.66737
## 16   0.80 80.52799  90.15957  24.19700 38.15419
## 17   0.85 79.58894  92.63699  19.30764 31.95511
## 18   0.90 78.65875  92.99781  15.16774 26.08162
## 19   0.95 78.02977  95.73864  12.02712 21.36969
umbral_final_rl_x3 <- umb_rl[which.max(umb_rl_x3$F1),1]
umbral_final_rl_x3
## [1] 0.3
#evaluamos la matriz de confusión y las métricas
confusion(test$income,rl_predict_x3,umbral_final_rl_x3)
##     
## real FALSE TRUE
##    0  7061 1425
##    1   600 2202
rl_metricas_x3 <- filter(umb_rl_x3,umbral==umbral_final_rl_x3)
rl_metricas_x3
##   umbral acierto precision cobertura      F1
## 1    0.3 82.0606  60.71133  78.58672 68.5021
#Calculamos la ROC y la AUC

#creamos el objeto prediction
rl_prediction_x3 <- prediction(rl_predict_x3,test$income)
#visualizamos la ROC
roc(rl_prediction_x3)

rl_metricas_x3 <- cbind(rl_metricas_x3,AUC=round(auc(rl_prediction_x3),2)*100)
print(t(rl_metricas_x3))
##               [,1]
## umbral     0.30000
## acierto   82.06060
## precision 60.71133
## cobertura 78.58672
## F1        68.50210
## AUC       90.00000

El resultado ha sido de una AUC de 90 con este modelo también.

Modelo 4: df2, df3 y df4 serán la muestra de entrenamiento, y df1 la muestra de test.

train <- rbind(df2,df3,df4)
test <- df1

formula_rl_x <- formula
rl_x4 <- glm(formula_rl_x,train,family=binomial(link='logit'))
summary(rl_x4)
## 
## Call:
## glm(formula = formula_rl_x, family = binomial(link = "logit"), 
##     data = train)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -3.2225  -0.5169  -0.2128  -0.0211   3.8022  
## 
## Coefficients:
##                                  Estimate Std. Error z value
## (Intercept)                      -8.08860    0.45401 -17.816
## fnlwgt_disc02_DE_65K_A_120K       0.07764    0.06756   1.149
## fnlwgt_disc03_DE_105K_A_180K      0.28689    0.06337   4.527
## fnlwgt_disc04_DE_160K_A_320K      0.28284    0.06047   4.677
## fnlwgt_disc05_MAS_320K            0.48824    0.07462   6.543
## workclassLocal-gov               -0.60836    0.10470  -5.811
## workclassPrivate                 -0.45732    0.08833  -5.177
## workclassSelf-emp-inc            -0.19087    0.11495  -1.660
## workclassSelf-emp-not-inc        -0.84737    0.10230  -8.284
## workclassState-gov               -0.59973    0.11640  -5.152
## workclassWithout-pay             -1.13270    0.86473  -1.310
## sexMale                           0.11417    0.04939   2.311
## raceAsian-Pac-Islander            0.59791    0.22812   2.621
## raceBlack                         0.21225    0.21793   0.974
## raceOther                         0.39833    0.31055   1.283
## raceWhite                         0.61420    0.20763   2.958
## c_status02_with_children         -0.57810    0.21213  -2.725
## c_status03_unmarried             -0.07989    0.17750  -0.450
## c_status04_married                2.17216    0.17575  12.359
## worktype02_Executives             0.81923    0.07231  11.330
## worktype03_Technical              0.36470    0.07423   4.913
## worktype04_services               0.47980    0.06869   6.985
## worktype05_manual_work           -0.23222    0.06995  -3.320
## educational_level02_high_school   1.02534    0.07997  12.821
## educational_level03_professional  1.58551    0.08158  19.436
## educational_level04_bachelor      2.18600    0.08539  25.599
## educational_level05_master        2.02616    0.08930  22.689
## educational_level06_doctor        3.14233    0.16437  19.118
## age_disc02 <= 24                  0.60743    0.37496   1.620
## age_disc03 <= 27                  1.39598    0.35746   3.905
## age_disc04 <= 29                  1.98728    0.35662   5.573
## age_disc05 <= 33                  2.26629    0.35124   6.452
## age_disc06 <= 35                  2.34583    0.35423   6.622
## age_disc07 <= 38                  2.60072    0.35179   7.393
## age_disc08 <= 41                  2.77200    0.35189   7.877
## age_disc09 <= 45                  2.79649    0.35134   7.959
## age_disc10 <= 54                  3.02029    0.34997   8.630
## age_disc11 <= 61                  2.86139    0.35218   8.125
## age_disc12 > 61                   2.34964    0.35577   6.604
## hours_per_week_disc02 <= 39       0.60191    0.10027   6.003
## hours_per_week_disc03 <= 43       0.77616    0.07342  10.572
## hours_per_week_disc04 > 43        1.24334    0.07483  16.614
## gain_disc02 > 4650                3.28039    0.09286  35.326
##                                              Pr(>|z|)    
## (Intercept)                      < 0.0000000000000002 ***
## fnlwgt_disc02_DE_65K_A_120K                   0.25043    
## fnlwgt_disc03_DE_105K_A_180K     0.000005977168616729 ***
## fnlwgt_disc04_DE_160K_A_320K     0.000002906946608564 ***
## fnlwgt_disc05_MAS_320K           0.000000000060403501 ***
## workclassLocal-gov               0.000000006223740020 ***
## workclassPrivate                 0.000000225223782365 ***
## workclassSelf-emp-inc                         0.09683 .  
## workclassSelf-emp-not-inc        < 0.0000000000000002 ***
## workclassState-gov               0.000000257289316632 ***
## workclassWithout-pay                          0.19023    
## sexMale                                       0.02081 *  
## raceAsian-Pac-Islander                        0.00877 ** 
## raceBlack                                     0.33009    
## raceOther                                     0.19961    
## raceWhite                                     0.00310 ** 
## c_status02_with_children                      0.00643 ** 
## c_status03_unmarried                          0.65265    
## c_status04_married               < 0.0000000000000002 ***
## worktype02_Executives            < 0.0000000000000002 ***
## worktype03_Technical             0.000000896521508270 ***
## worktype04_services              0.000000000002848947 ***
## worktype05_manual_work                        0.00090 ***
## educational_level02_high_school  < 0.0000000000000002 ***
## educational_level03_professional < 0.0000000000000002 ***
## educational_level04_bachelor     < 0.0000000000000002 ***
## educational_level05_master       < 0.0000000000000002 ***
## educational_level06_doctor       < 0.0000000000000002 ***
## age_disc02 <= 24                              0.10524    
## age_disc03 <= 27                 0.000094130994624736 ***
## age_disc04 <= 29                 0.000000025107742631 ***
## age_disc05 <= 33                 0.000000000110222391 ***
## age_disc06 <= 35                 0.000000000035368415 ***
## age_disc07 <= 38                 0.000000000000143668 ***
## age_disc08 <= 41                 0.000000000000003340 ***
## age_disc09 <= 45                 0.000000000000001728 ***
## age_disc10 <= 54                 < 0.0000000000000002 ***
## age_disc11 <= 61                 0.000000000000000448 ***
## age_disc12 > 61                  0.000000000039933848 ***
## hours_per_week_disc02 <= 39      0.000000001941594260 ***
## hours_per_week_disc03 <= 43      < 0.0000000000000002 ***
## hours_per_week_disc04 > 43       < 0.0000000000000002 ***
## gain_disc02 > 4650               < 0.0000000000000002 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 37950  on 33861  degrees of freedom
## Residual deviance: 22446  on 33819  degrees of freedom
## AIC: 22532
## 
## Number of Fisher Scoring iterations: 8
#cogemos solo las variables relevantes segun el modelo anterior
a_mantener <- c(
  'workclass',
  'sex',
  'c_status',
  'worktype',
  'educational_level',
  'age_disc',
  'hours_per_week_disc',
  'gain_disc',
  'fnlwgt_disc'
)

formula_rl_x <- reformulate(a_mantener,target)
rl_x4 <- glm(formula_rl_x,train,family=binomial(link='logit'))
summary(rl_x4)
## 
## Call:
## glm(formula = formula_rl_x, family = binomial(link = "logit"), 
##     data = train)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -3.2174  -0.5175  -0.2158  -0.0213   3.7150  
## 
## Coefficients:
##                                  Estimate Std. Error z value
## (Intercept)                      -7.59325    0.40832 -18.596
## workclassLocal-gov               -0.58383    0.10414  -5.606
## workclassPrivate                 -0.41815    0.08768  -4.769
## workclassSelf-emp-inc            -0.14050    0.11441  -1.228
## workclassSelf-emp-not-inc        -0.79954    0.10164  -7.866
## workclassState-gov               -0.57540    0.11592  -4.964
## workclassWithout-pay             -1.07255    0.86473  -1.240
## sexMale                           0.13196    0.04923   2.680
## c_status02_with_children         -0.55151    0.21192  -2.602
## c_status03_unmarried             -0.05008    0.17717  -0.283
## c_status04_married                2.21466    0.17543  12.624
## worktype02_Executives             0.82244    0.07221  11.390
## worktype03_Technical              0.36956    0.07410   4.988
## worktype04_services               0.47465    0.06853   6.927
## worktype05_manual_work           -0.23533    0.06980  -3.372
## educational_level02_high_school   1.03454    0.07991  12.946
## educational_level03_professional  1.59729    0.08152  19.593
## educational_level04_bachelor      2.20873    0.08528  25.899
## educational_level05_master        2.04747    0.08917  22.960
## educational_level06_doctor        3.17983    0.16409  19.379
## age_disc02 <= 24                  0.60640    0.37510   1.617
## age_disc03 <= 27                  1.38975    0.35754   3.887
## age_disc04 <= 29                  1.97498    0.35668   5.537
## age_disc05 <= 33                  2.26050    0.35131   6.435
## age_disc06 <= 35                  2.33565    0.35429   6.592
## age_disc07 <= 38                  2.59428    0.35186   7.373
## age_disc08 <= 41                  2.76252    0.35196   7.849
## age_disc09 <= 45                  2.78488    0.35141   7.925
## age_disc10 <= 54                  3.00853    0.35005   8.595
## age_disc11 <= 61                  2.85608    0.35225   8.108
## age_disc12 > 61                   2.35143    0.35586   6.608
## hours_per_week_disc02 <= 39       0.58992    0.10010   5.893
## hours_per_week_disc03 <= 43       0.76704    0.07327  10.469
## hours_per_week_disc04 > 43        1.24802    0.07471  16.705
## gain_disc02 > 4650                3.27396    0.09262  35.348
## fnlwgt_disc02_DE_65K_A_120K       0.08127    0.06735   1.207
## fnlwgt_disc03_DE_105K_A_180K      0.28498    0.06320   4.509
## fnlwgt_disc04_DE_160K_A_320K      0.27330    0.06025   4.536
## fnlwgt_disc05_MAS_320K            0.46671    0.07419   6.291
##                                              Pr(>|z|)    
## (Intercept)                      < 0.0000000000000002 ***
## workclassLocal-gov               0.000000020700745713 ***
## workclassPrivate                 0.000001853101397855 ***
## workclassSelf-emp-inc                        0.219452    
## workclassSelf-emp-not-inc        0.000000000000003660 ***
## workclassState-gov               0.000000692195263397 ***
## workclassWithout-pay                         0.214854    
## sexMale                                      0.007355 ** 
## c_status02_with_children                     0.009257 ** 
## c_status03_unmarried                         0.777426    
## c_status04_married               < 0.0000000000000002 ***
## worktype02_Executives            < 0.0000000000000002 ***
## worktype03_Technical             0.000000611653079797 ***
## worktype04_services              0.000000000004308174 ***
## worktype05_manual_work                       0.000747 ***
## educational_level02_high_school  < 0.0000000000000002 ***
## educational_level03_professional < 0.0000000000000002 ***
## educational_level04_bachelor     < 0.0000000000000002 ***
## educational_level05_master       < 0.0000000000000002 ***
## educational_level06_doctor       < 0.0000000000000002 ***
## age_disc02 <= 24                             0.105958    
## age_disc03 <= 27                             0.000102 ***
## age_disc04 <= 29                 0.000000030746349734 ***
## age_disc05 <= 33                 0.000000000123836235 ***
## age_disc06 <= 35                 0.000000000043254776 ***
## age_disc07 <= 38                 0.000000000000166637 ***
## age_disc08 <= 41                 0.000000000000004197 ***
## age_disc09 <= 45                 0.000000000000002284 ***
## age_disc10 <= 54                 < 0.0000000000000002 ***
## age_disc11 <= 61                 0.000000000000000515 ***
## age_disc12 > 61                  0.000000000039047041 ***
## hours_per_week_disc02 <= 39      0.000000003784551620 ***
## hours_per_week_disc03 <= 43      < 0.0000000000000002 ***
## hours_per_week_disc04 > 43       < 0.0000000000000002 ***
## gain_disc02 > 4650               < 0.0000000000000002 ***
## fnlwgt_disc02_DE_65K_A_120K                  0.227520    
## fnlwgt_disc03_DE_105K_A_180K     0.000006502918265238 ***
## fnlwgt_disc04_DE_160K_A_320K     0.000005728131186358 ***
## fnlwgt_disc05_MAS_320K           0.000000000315798550 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 37950  on 33861  degrees of freedom
## Residual deviance: 22489  on 33823  degrees of freedom
## AIC: 22567
## 
## Number of Fisher Scoring iterations: 8
#calculamos el pseudo R cuadrado
pr2_rl_x4 <- 1 -(rl_x4$deviance / rl_x4$null.deviance)
pr2_rl_x4
## [1] 0.4074112
#aplicamos el modelo al conjunto de test
rl_predict_x4 <- predict(rl_x4,test,type = 'response')

#hacemos una primera visualización de como es el modelo
plot(rl_predict_x4~test$income)

#Calculamos los umbrales y cual es el optimo segun la F1
umb_rl_x4 <- umbrales(test$income,rl_predict_x4)
umb_rl_x4
##    umbral  acierto precision cobertura       F1
## 1    0.05 59.78030  37.95064  97.68023 54.66347
## 2    0.10 69.66690  44.76959  95.00357 60.85963
## 3    0.15 74.85826  49.64899  90.86367 64.21185
## 4    0.20 77.95889  53.43244  87.22341 66.26898
## 5    0.25 80.09391  56.78318  82.90507 67.40171
## 6    0.30 81.68852  60.02182  78.55103 68.04761
## 7    0.35 82.79589  63.05404  74.12562 68.14304
## 8    0.40 83.38058  65.61025  69.45039 67.47573
## 9    0.45 83.87668  68.65502  64.48965 66.50718
## 10   0.50 84.10702  71.52007  59.77873 65.12442
## 11   0.55 83.97413  73.85872  54.85368 62.95310
## 12   0.60 83.44259  75.87354  48.82227 59.41368
## 13   0.65 82.90220  78.23834  43.11206 55.59135
## 14   0.70 82.36180  81.75411  37.25910 51.18902
## 15   0.75 81.44047  84.62292  30.83512 45.20010
## 16   0.80 80.52799  88.71795  24.69665 38.63763
## 17   0.85 79.34975  91.09948  18.62955 30.93333
## 18   0.90 78.64989  93.36283  15.06067 25.93731
## 19   0.95 77.82601  94.89489  11.27766 20.15949
umbral_final_rl_x4 <- umb_rl[which.max(umb_rl_x4$F1),1]
umbral_final_rl_x4
## [1] 0.35
#evaluamos la matriz de confusión y las métricas
confusion(test$income,rl_predict_x4,umbral_final_rl_x4)
##     
## real FALSE TRUE
##    0  7269 1217
##    1   725 2077
rl_metricas_x4 <- filter(umb_rl_x4,umbral==umbral_final_rl_x4)
rl_metricas_x4
##   umbral  acierto precision cobertura       F1
## 1   0.35 82.79589  63.05404  74.12562 68.14304
#Calculamos la ROC y la AUC

#creamos el objeto prediction
rl_prediction_x4 <- prediction(rl_predict_x4,test$income)
#visualizamos la ROC
roc(rl_prediction_x4)

rl_metricas_x4 <- cbind(rl_metricas_x4,AUC=round(auc(rl_prediction_x4),2)*100)
print(t(rl_metricas_x4))
##               [,1]
## umbral     0.35000
## acierto   82.79589
## precision 63.05404
## cobertura 74.12562
## F1        68.14304
## AUC       89.00000

El resultado ha sido de una AUC de 90 con este modelo también.

Juntamos los 4 set de metricas para hacer una comparativas de los 4 modelos generados:

rl_metricas_x <- rbind(rl_metricas_x1,rl_metricas_x2,rl_metricas_x3,rl_metricas_x4)
rl_metricas_x
##   umbral  acierto precision cobertura       F1 AUC
## 1   0.30 82.64221  61.57229  79.97144 69.57602  90
## 2   0.30 82.97307  62.45753  78.72948 69.65583  90
## 3   0.30 82.06060  60.71133  78.58672 68.50210  90
## 4   0.35 82.79589  63.05404  74.12562 68.14304  89

Salvo pequeñas variaciones, las metricas son similares en los 4 modelos. Por lo tanto, hemos validado la regresión logistica con esta cross-validation, y como tiene la AUC más alta, elegimos este modelo como el que mejor ajusta los datos. No evaluamos con cross-validation los otros dos modelos con arboles de decision y con random forest.

5.3.8. - Escribimos el scoring final en el dataset y guardamos el modelo

df$scoring <- predict(rl,df,type = 'response')
saveRDS(rl,'03_modelo_final.rds') #guardamos el modelo
saveRDS(df,'cache_final.rds') #guardamos el dataset con el scoring

6. Evaluacion y analisis de negocio

6.1. Visualización por tramos

Vamos a visualizar la penetración de la target por tramos de scoring. Este grafico es muy potente para ver que el modelo es consistente, ya que debe presentar una linea descente en la tasa de contratacion conforme se desciende en el scoring.

#Creamos una funcion para visualizar la contratacion real por percentiles de scoring
vis <- function(scoring,real) {
    #Preparar el dataframe de visualizaciÛn
    vis_df <- data.frame(Scoring = scoring, Perc_Scoring = cut_number(scoring, 20), Real = real)
    levels(vis_df$Perc_Scoring) <- seq(from = 100,to = 5,by = -5)
    vis_gr <- vis_df %>% group_by(Perc_Scoring) %>% summarise(Tasa_Contr = mean(as.numeric(as.character(Real)))) %>% arrange(Perc_Scoring)
    #ordenar el factor para el grafico
    vis_gr$Perc_Scoring <- factor(vis_gr$Perc_Scoring, levels = vis_gr$Perc_Scoring[order(vis_gr$Perc_Scoring, decreasing = T)])
    #hacemos el grafico
    ggplot(vis_gr,aes(Perc_Scoring, Tasa_Contr)) + 
      geom_col(fill='grey') + 
      geom_hline(aes(yintercept =      mean(as.numeric(as.character(vis_df$Real)))), col = 'black') +
      labs(title = 'Contratación real por tramo de scoring', x = 'Tramo de Scoring', y = 'Contratación real')
}
vis(df$scoring,df$income)

El gráfico es coherente con un modelo bien ajustado, siendo completamente descencente.

6.2. Cuantas personas cobran más de 50K dolares

Puede haber interes en saber quien puede estar en un nivel de ingresos por encima de 50K dolares, por ejemplo, para hacer promociones de productos financieros o publicidad de productos de alta gama. Solo con este parametros no se puede saber la cantidad de dinero disponible, pero puede ser un primer indicador.

Para ver a partir de que punto están las personas cuyo scoring está por encima de la penetración de la variable target (aprox 25%). Parece ser que son 16500 personas aproximadamente.

penetracion_target <- mean(as.numeric(as.character(df$income)))
df %>% 
  arrange(desc(scoring)) %>% 
  ggplot(aes(y = scoring, x = seq_along(scoring))) +
  geom_line() + 
  geom_vline(xintercept = 16500, col = 'orange') +
  geom_hline(yintercept = penetracion_target,col='blue') +
  labs(x = 'PERSONAS ORDENADOS POR SCORING', y = 'SCORING')

Estos son muchas personas, y solo por el nivel de ingresos no podemos inferir cuanto dinero tienen disponible. Una de las variables es las ganancias financieras, que la discretizamos en la variable gain_disc. Podemos filtrar el rango más alto de esta variable discretizada, y vemos que personas tienen el mayor scoring y además tienen ingresos de productos financieros, y es más probable que tengan más dinero disponible:

df %>% 
  arrange(desc(scoring)) %>% 
  filter(gain_disc == "02 > 4650") %>%
  ggplot(aes(y = scoring, x = seq_along(scoring))) +
  geom_line() + 
  geom_vline(xintercept = 2375, col = 'orange') +
  geom_hline(yintercept = penetracion_target,col='blue') +
  labs(x = 'PERSONAS ORDENADAS POR SCORING CON GANANCIAS FINANCIERAS', y = 'SCORING')

Vemos que todos estos clientes tienen un muy alto scoring y que la curva cae muy rapidamente, como era de esperar. Ahora ya solo tenemos unos 2375 personas que tienen ingresos de más de 50K dolares y además tienen ganancias financieras altas. Con estas dos variables podemos tener una mejor aproximación que estas personas puede tener dinero disponible, y por tanto, ser target de nuestra publicidad y promociones.

A partir de estas personas, podemos filtrar, para ser un poco más conservadores, aquellos clientes con ganancias financieras y con scoring mayor de 0.35, que es el umbral optimo que salio de las metricas de la regresión logistica. Con este grupo, y con las variables que tenemos, podríamos hacer segmentaciones tipicas, como por ejemplo, por edad o por nivel de estudios, dependiendo de la promoción o publicidad que queremos hacer.

Vamos a visualizar la media del scoring para los clientes con mayores ganancias financieras y con scoring mayor de 0.35 en cada una de las franjas de edad:

df %>%
  arrange(desc(scoring)) %>%
  filter(gain_disc == "02 > 4650") %>%
  filter(scoring > 0.35) %>%
  group_by(age_disc) %>%
  summarise(Conteo = n(), porc = mean(as.numeric(as.character(scoring)))) %>%
  ggplot(aes(age_disc,porc)) + geom_bar(stat='identity') + xlab("age")

Podemos ver que los mayores scorings están entre 38 y 61 años, con esto podemos focalizar mucho más nuestra campaña.

Podemos hacer la misma visualización pero por nivel de educación:

df %>%
  arrange(desc(scoring)) %>%
  filter(gain_disc == "02 > 4650") %>%
  filter(scoring > 0.35) %>%
  group_by(educational_level) %>%
  summarise(Conteo = n(), porc = mean(as.numeric(as.character(scoring)))) %>%
  ggplot(aes(educational_level,porc)) + geom_bar(stat='identity') + xlab("Nivel de educación")

Vemos que el mayor scoring es para las personas con educación superior.

Hacemos los filtros anteriores y tendremos el dataset con los clientes foco de la campaña:

df_final <- df %>% 
            arrange(desc(scoring)) %>% 
            filter(gain_disc == "02 > 4650") %>%
            filter(scoring > 0.35) %>%
            filter(age_disc == "07 <= 38" | age_disc == "08 <= 41" | age_disc == "09 <= 45" | age_disc == "10 <= 54" | age_disc == "11 <= 61") %>%
            filter(educational_level == "03_professional" | educational_level == "04_bachelor" | educational_level == "05_master" | educational_level =="06_doctor")

#dibujamos el scoring del set de clientes escogidos
ggplot(df_final,aes(y = scoring, x = seq_along(scoring))) +
  geom_line() + 
  labs(x = 'PERSONAS ORDENADAS POR SCORING CON GANANCIAS FINANCIERAS', y = 'SCORING')

print(head(df_final))
##         fnlwgt_disc    workclass  sex  race   c_status      worktype
## 1       05_MAS_320K      Private Male White 04_married 02_Executives
## 2       05_MAS_320K Self-emp-inc Male White 04_married 02_Executives
## 3  02_DE_65K_A_120K Self-emp-inc Male White 04_married 02_Executives
## 4 04_DE_160K_A_320K  Federal-gov Male White 04_married   04_services
## 5 04_DE_160K_A_320K      Private Male White 04_married 02_Executives
## 6 04_DE_160K_A_320K      Private Male White 04_married 02_Executives
##   educational_level age_disc hours_per_week_disc gain_disc income
## 1         06_doctor 10 <= 54             04 > 43 02 > 4650      1
## 2         06_doctor 08 <= 41             04 > 43 02 > 4650      1
## 3         06_doctor 10 <= 54             04 > 43 02 > 4650      1
## 4         06_doctor 10 <= 54             04 > 43 02 > 4650      1
## 5         06_doctor 10 <= 54             04 > 43 02 > 4650      1
## 6         06_doctor 10 <= 54             04 > 43 02 > 4650      1
##     scoring
## 1 0.9978724
## 2 0.9977875
## 3 0.9977200
## 4 0.9976058
## 5 0.9975655
## 6 0.9975655
cat("Número de clientes objetivo:",nrow(df_final),"\n")
## Número de clientes objetivo: 1322

Grabamos la lista de personas objetivo en un fichero:

fwrite(df_final,"personas_objetivo.csv")

FINAL DEL DOCUMENTO

Andrés López Peñaranda

06/12/2019