Grupo conformado por: Biel Arenas Bosch, Benjamín Gálvez Megía, Jose Antonio Doyague Hernández, Alberto Fernández Hernández

Nota: para la realización de la práctica, y a lo largo del concurso, se ha tomado como base el score obtenido en el modelo 2 realizado en clase:

  1. Selección de variables numéricas y categóricas (con al menos dos categorías y menos de 1000), además de descartar las variables duplicadas payment_type y quantity_group, junto a recorded_by (con una única categoría), obteniendo la siguiente puntuación:
Train.accuracy Driven.Data
Num + Cat (> 1 & < 1000) sin duplicados 0.8168687 0.8128
  1. Por otro lado, los chunks con las ejecuciones de los modelos están sin ejecutar, para no ralentizar la elaboración del documento. No obstante, en cada apartado se incluyen tablas con las salidas obtenidas, así como el score del test resultante en DrivenData. Además, en el último apartado se dispone de un anexo con el código final ejecutable.

Durante toda la práctica los modelos se han realizado con la misma semilla: 1234.

1. Carga de datos y librerías

Para la realización de la práctica, los ficheros utilizados son los siguientes:

  1. train_values.csv: fichero con el conjunto de variables y observaciones con los que entrenar el modelo.

  2. train_labels.csv: fichero con las variables objetivo de cada una de las observaciones en train_values.csv:

  • functional

  • non functional

  • functional needs repair

  1. test_values.csv: fichero de prueba con el que maximizar las predicciones obtenidas con el correspondiente modelo

En relación con las librerías empleadas, debemos destacar las siguientes:

#--- Librerias empleadas
suppressPackageStartupMessages({
  library(dplyr)           # Manipulacion de datos 
  library(data.table)      # Lectura y escritura de ficheros
  library(ggplot2)         # Representacion grafica
  library(inspectdf)       # EDAs automaticos
  library(ranger)          # randomForest (+ rapido que caret)
  library(forcats)         # Tratamiento de variables categoricas
  library(tictoc)          # Calculo del tiempo de ejecucion
  library(missRanger)      # Imputacion de valores NA
  library(knitr)           # Generacion de documentos y formateo de tablas
  library(kableExtra)      # Formateo de tablas (colores)
  library(gmt)             # Calculo de la distancia geografica
  library(stringi)         # Tratamiento de strings
  library(FeatureHashing)  # Hashing Encoding
  library(dataPreparation) # Target encoding
  library(embed)           # Word-Embeddings
  library(recipes)         # Definicion de recetas (preprocesamiento de datos)
  library(lubridate)       # Tratamiento de fechas
  library(ggrepel)         # Añadir etiquetas (texto) a ggplot
  library(xgboost)         # Modelo XGboost
})

# Nota: junto al resto de paquetes debe añadirse h2o, aunque no se ha incluido inicialmente,
#       dado que sobreescribe algunas funciones como las disponibles en "lubridate"
#--- Carga de ficheros train y test
dattrainOr    <- fread(file = "./data/train_values.csv", data.table = FALSE )
dattrainLabOr <- fread(file = "./data/train_labels.csv", data.table = FALSE )
dattestOr     <- fread(file = "./data/test_values.csv",  data.table = FALSE )

2. Desarrollo durante el concurso

En esta sección se comentará las aportaciones que realicé en el grupo durante el concurso para mejorar la puntuación obtenida en clase, dado que cada miembro del grupo realizó sus pruebas de forma independiente, sin poner en común un mismo modelo, selección de variables o Feature Engineering.

2.1 Análisis de variables existentes

En primer lugar, en base a las funciones disponibles en la librería inspectdf, se extrajo un esquema inicial con las variables existentes en el conjunto de datos y su correspondiente tipo:

# Analizamos el tipo de dato en cada variable
show_plot(inspect_types(dattrainOr))

En conjunto, nos encontramos con:

  1. 27 variables de tipo carácter (categóricas)
  2. 10 variables numéricas (entre integer y numeric)
  3. 2 variables lógicas
  4. 1 variable tipo fecha

Por tanto, del mismo modo que lo realizado en el modelo 2:

  1. Conservamos únicamente las variables numéricas y categóricas con más de uno y menos de 1000 valores diferentes:
#-- Variables categoricas
datcat_df <- dattrainOr %>% select(where(is.character))
# Mediante un bucle for recuperamos el numero de categorias de cada variable caracter
numlev_df <- data.frame(
                          "vars" = names(datcat_df),
                           "levels" = apply(datcat_df, 2, 
                                  function(x) length(unique(x)))
                          
                        )
# Eliminamos los nombres de fila
rownames(numlev_df) <- NULL

#-- Ordenamos de menor a mayor numero de categorias
numlev_df_ordenado <- numlev_df %>% arrange(levels)

color_me <- which(numlev_df_ordenado$levels < 1000 & 
                    numlev_df_ordenado$levels > 1)
kable(numlev_df_ordenado) %>%
  kable_styling() %>%
  row_spec(color_me, bold = TRUE, background = "#4D934D")
vars levels
recorded_by 1
source_class 3
management_group 5
quantity 5
quantity_group 5
quality_group 6
waterpoint_type_group 6
extraction_type_class 7
payment 7
payment_type 7
source_type 7
waterpoint_type 7
water_quality 8
basin 9
source 10
management 12
scheme_management 13
extraction_type_group 13
extraction_type 18
region 21
lga 125
funder 1898
ward 2092
installer 2146
scheme_name 2697
subvillage 19288
wpt_name 37400

Por tanto, escogemos desde la variable source_class hasta lga, dado que recorded_by solo presenta una única categoría:

#-- Conservamos variables con categorias > 1 & < 1000
vars_gd <- numlev_df %>%
  filter(levels < 1000, levels > 1) %>% 
  select(vars)
datcat_gd <- datcat_df[ , vars_gd$vars]
#-- Variables numericas
datnum_df <- dattrainOr %>% select(where(is.numeric))
# Unificamos ambos tipos de variables...
datnumcat_df <- cbind(datnum_df, datcat_gd)
# ...Como tambien la variable objetivo
dattrainOrlab <- merge(
  datnumcat_df, dattrainLabOr,
  by.x = c('id'), by.y = c('id'),
  sort = FALSE
)

El resto de variables las almacenamos en un data.table auxiliar, con el objetivo de emplearlas en la sección 3:

#-- Almacenamos variables categoricas sin utilizar
datcat_mas_1000_mas_logicas_dr <- data.table(
  funder         = c(dattrainOr$funder,         dattestOr$funder         ),
  ward           = c(dattrainOr$ward,           dattestOr$ward           ),
  installer      = c(dattrainOr$installer,      dattestOr$installer      ),
  scheme_name    = c(dattrainOr$scheme_name,    dattestOr$scheme_name    ),
  subvillage     = c(dattrainOr$subvillage,     dattestOr$subvillage     ),
  wpt_name       = c(dattrainOr$wpt_name,       dattestOr$wpt_name       ),
  permit         = c(dattrainOr$permit,         dattestOr$permit         ),
  public_meeting = c(dattrainOr$public_meeting, dattestOr$public_meeting ),
  date_recorded  = c(dattrainOr$date_recored,   dattestOr$date_recorded  )
)
  1. Eliminación de variables con información duplicada, concretamente payment_type y quantity_group:
#-- payment_type: categorias muy similares con la variable "payment"
table(dattrainOrlab$payment_type)
## 
##   annually    monthly  never pay on failure      other per bucket    unknown 
##       3642       8300      25348       3914       1054       8985       8157
table(dattrainOrlab$payment)
## 
##             never pay                 other          pay annually 
##                 25348                  1054                  3642 
##           pay monthly        pay per bucket pay when scheme fails 
##                  8300                  8985                  3914 
##               unknown 
##                  8157
dattrainOrlab$payment_type   <- NULL
#-- quantity_group: mismas categorias que con la variable quantity
table(dattrainOrlab$quantity_group)
## 
##          dry       enough insufficient     seasonal      unknown 
##         6246        33186        15129         4050          789
table(dattrainOrlab$quantity)
## 
##          dry       enough insufficient     seasonal      unknown 
##         6246        33186        15129         4050          789
dattrainOrlab$quantity_group <- NULL

2.1 Tratamiento de valores anómalos

Una vez pre-procesado el conjunto de entrenamiento, una de las primeras ideas que se tuvo en cuenta fue el tratamiento de valores anómalos, concretamente a cero, de determinadas variables, anotados durante la primera fase de análisis exploratorio realizada en clase:

  1. construction_year

  2. gps_height: en el caso de la altura GPS, no se han considerado los valores negativos como valores anómalos, principalmente porque, trasladando el problema a un entorno real, es posible la existencia de alturas negativas (regiones situadas por debajo del nivel del mar), además de no existir una elevada acumulación (como si sucede con el cero).

  3. longitude

Sobre dichas variables podemos observar en los siguientes histogramas una concentración “anómala” de valores a cero:

# Graficos de distribucion
# Histograms for numeric columns
show_plot(
  inspect_num(dattrainOrlab[c("construction_year",
                              "gps_height",
                              "longitude")])
  
)

Por tanto, para eliminar y posteriormente imputar dichos ceros, haremos uso de una función denominada impute_zeros, la cual permite establecer inicialmente a NA los valores a cero, así como su posterior imputación haciendo uso de la librería missRanger.

Nota: a lo largo de la memoria se definen varias funciones que se emplean comunmente tanto en el concurso como en el desarrollo individual, evitando con ello la redundancia de código.

# Funcion para la imputacion de ceros con missRanger
impute_zeros <- function(data, column, k, num_trees, impute = TRUE, seed = 1234) {
  new_feature <- paste0("fe_", column)
  
  data[, new_feature] <- data[, column]
  data[, new_feature] <- ifelse( 
                            data[, new_feature] == 0, 
                            NA, 
                            data[, new_feature]
  )
  original_data       <- data[, new_feature]
  
  data[, column]      <- NULL
  
  if (impute == TRUE) {
    data_imp <- missRanger(
                            data,  
                            pmm.k = k,
                            num.trees = num_trees,
                            seed = seed,
                            verbose = 0
                           )
    # Finalmente, creamos una nueva columna binaria
    # 1 = el valor era NA; 0 = No
    new_feature_na      <- paste0("fe_", column, "_na")
    data_imp[, new_feature_na] <- ifelse(
                                    is.na(original_data),
                                    1,
                                    0
    )
  }
  return(data_imp)
}

Una vez definida la funcion, imputamos cada una de las columnas anteriores, no sin antes eliminar la columna status_group, con el objetivo de evitar que la imputación también dependa del valor de la variable objetivo, empleando k = 5 como número de vecinos:

# Almacenamos en un vector la variable objetivo
status_group_vector        <- dattrainOrlab$status_group
dattrainOrlab$status_group <- NULL
#--- construction_year
dattrainOrlab_imp <- impute_zeros(dattrainOrlab,
                                  column = "construction_year",
                                  k = 5,
                                  num_trees = 100)
#--- gps_heigth
dattrainOrlab_imp <- impute_zeros(dattrainOrlab_imp,
                                  column = "gps_height",
                                  k = 5,
                                  num_trees = 100)
#--- longitude
dattrainOrlab_imp <- impute_zeros(dattrainOrlab_imp,
                                  column = "longitude",
                                  k = 5,
                                  num_trees = 100)

Tras la imputación, añadimos nuevamente la variable objetivo al conjunto de entrenamiento:

dattrainOrlab_imp$status_group <- status_group_vector
dattrainOrlab_imp$status_group <- as.factor(dattrainOrlab_imp$status_group)

En relación con los datos test, dado que no debemos aplicar la imputación de los valores missing a partir de dicho conjunto, a cada valor NA (es decir, igual a cero), se tomó la decisión de imputarlo por un valor muy diferente al resto de la población, concretamente a 99.999, así como añadir una columna binaria adicional que indice si el valor imputado era o no un valor missing:

dattestOr_imp <- dattestOr

#-- Completar valores anomales a 99999 + columnas "_na_"
for(column in c("construction_year", "gps_height", "longitude")) {
  new_feature <- paste0("fe_", column)
  dattestOr_imp[, new_feature] <- ifelse(
                                          dattestOr_imp[, column] == 0,
                                          99999,
                                          dattestOr_imp[, column]
  )
  dattestOr_imp[, column] <- NULL
  
  new_feature_na <- paste0("fe_", column, "_na")
  dattestOr_imp[, new_feature_na] <- ifelse(
                                          dattestOr_imp[, new_feature] == 99999,
                                          1,
                                          0
  ) 
}

Una vez imputadas las variables, por medio de la función ranger creamos un primer modelo de Random Forest:

#-- Funcion que devuelve un modelo Random Forest entrenado
fit_random_forest <- function(formula, data, num_trees = 500, mtry = NULL, seed = 1234) {
  tic()
  my_model <- ranger( 
    formula, 
    importance = 'impurity',
    data       = data,
    num.trees = num_trees,
    mtry = mtry,
    verbose = FALSE,
    seed = seed
  )
  # "Estimacion" del error / acierto "esperado" (OOB accuracy)
  success <- 1 - my_model$prediction.error
  print(success)
  toc()
  
  return(my_model)
}

De forma adicional, creamos una función que permita:

  1. Crear un grafico con la importancia de las variables del modelo pasado como argumento

  2. Almacenar el gráfico en disco

#-- Funcion para pintar importancia de variables (ademas de guardar el grafico)
save_importance_ggplot <- function(model, path, title) {
  impor_df <- as.data.frame(model$variable.importance)
  names(impor_df)[1] <- c('Importance')
  impor_df$vars <- rownames(impor_df)
  rownames(impor_df) <- NULL
  
  print(ggplot(impor_df, aes(fct_reorder(vars, Importance), Importance)) +
      geom_col(group = 1, fill = "darkred") +
      coord_flip() + 
      labs(x = 'Variables', y = 'Importancia', title = title) +
      theme(axis.text.y = element_text(face = "bold", colour = "black"))
    )
  ggsave(path)
}
#-- Entrenamiento del primer modelo
formula    <- as.formula('status_group~.')

my_model_1 <- fit_random_forest(formula, dattrainOrlab_imp)
[1] 0.8101178
36.403 sec elapsed

Tras crear el modelo, guardamos el gráfico con la importancia de cada una de las variables:

save_importance_ggplot(my_model_1, "./charts/01_num_cat_menos1000_todo_imp.png",
                       title = "Importancia Variables Modelo 01")

Analizando el gráfico de importancia, observamos que las variables imputadas (año de construcción, altura GPS y longitud) presentan una elevada importancia. No obstante, las variables "_na" no parecen aportar una importancia relevante al modelo. Pese a ello, realizamos las predicciones y comparamos el score obtenido:

#-- Funcion para realizar prediccion sobre el modelo pasado como parametro
make_predictions <- function(model, test_data) {
  # Prediccion
  my_pred <- predict(model, test_data)
  
  # Submission
  my_sub <- data.table(
    id = test_data[, "id"],
    status_group = my_pred$predictions
  )
  
  return(my_sub)
}
my_sub_1 <- make_predictions(my_model_1, dattestOr_imp)
# guardo submission
fwrite(my_sub_1, file = "./submissions/01_num_cat_menos1000_todo_imp.csv")
Train.accuracy Driven.Data
Num + Cat (> 1 & < 1000) sin duplicados 0.8168687 0.8128
Num + Cat (> 1 & < 1000) sin duplicados longitude-construction_year-gps_height imp. 0.8101178 0.8096

Como podemos comprobar, el hecho de imputar las variables anteriores no ha mejorado el modelo. Sin embargo, ¿Y si eliminamos las variables “_na” anteriores?

dattrainOrlab_imp_sin_col_na <- dattrainOrlab_imp
# Eliminamos las columnas "_na"
dattrainOrlab_imp_sin_col_na$fe_construction_year_na <- NULL
dattrainOrlab_imp_sin_col_na$fe_gps_height_na        <- NULL
dattrainOrlab_imp_sin_col_na$fe_longitude_na         <- NULL

# Idem con los datos test
dattestOr_imp_sin_col_na    <- dattestOr_imp
# Eliminamos las columnas "_na"
dattestOr_imp_sin_col_na$fe_construction_year_na <- NULL
dattestOr_imp_sin_col_na$fe_gps_height_na        <- NULL
dattestOr_imp_sin_col_na$fe_longitude_na         <- NULL
my_model_1_sin_cols_na <- fit_random_forest(formula, dattrainOrlab_imp_sin_col_na)
[1] 0.8105724
39.914 sec elapsed

En relación con el conjunto train, parece mejorar ligeramente (de 0.8103535 a 0.811229). Analicemos el conjunto test:

my_sub_1_sin_cols_na <- make_predictions(my_model_1_sin_cols_na, dattestOr_imp_sin_col_na)
# guardo submission
fwrite(my_sub_1_sin_cols_na, file = "./submissions/01_num_cat_menos1000_todo_imp_sin_cols_na.csv")
Train.accuracy Driven.Data
Num + Cat (> 1 & < 1000) sin duplicados 0.8168687 0.8128
Num + Cat (> 1 & < 1000) sin duplicados longitude-construction_year-gps_height imp. 0.8101178 0.8096
Anterior, eliminando columnas ’_na’ 0.8105724 0.8081

Podemos comprobar que, pese a eliminar dichas variables, el score del modelo empeora.

2.2 Creación de nuevas variables

En vistas a los resultados obtenidos en el apartado anterior, se decidió no imputar dichas variables y buscar una alternativa a la mejora del modelo. Para ello, se recurrió a la creación de nuevas variables o feature engineering:

  1. construction_year: antigüedad del pozo, restando la fecha más reciente en el dataset (2014) a cada uno de los valores.

  2. longitude y latitude: distancia geográfica. En lugar de la raíz cuadrada de la suma de las coordenadas al cuadrado, R dispone de un paquete denominado gmt con el que calcular la distancia geográfica de cada bomba al punto de referencia (0,0), en lugar de emplear la distancia Euclídea:

\[ D = cos^{−1}[sin \theta_1 sin \theta_2 + cos\theta_1 cos \theta_2 cos(\phi_1−\phi_2)]) \]

#-- Ejemplo
geodist(55.75,37.63, 0, 0)  # Moscow - (0,0)
## [1] 7059.486
  1. amount_tsh y population: analizando el gráfico de importancia del primer modelo, se comprobó que variables como la población o la cantidad de agua disponible no son especialmente relevantes. Como consecuencia, se propuso la creación de una nueva variable, basada en la cantidad de agua que habría disponible por persona en cada pozo, dividiendo la cantidad total de agua en cada observación ( amount_tsh ) entre el total de población.
#---- Feature Engineering
#---  Creacion de nuevas variables
#--   Antiguedad del pozo (2014 - fecha)
# Train
dattrainOrlab$fe_cyear <- 2014 - dattrainOrlab$construction_year

# Test
dattestOr$fe_cyear     <- 2014 - dattestOr$construction_year
#--   Distancia geografica al punto (0,0)
# geodist (latitud_origen, longitud_origen, latitud_destino, longitud_destino)
# Por defecto, las unidades estan expresadas en kilometros
# Train
dattrainOrlab$fe_dist <- geodist(dattrainOrlab$latitude, dattrainOrlab$longitude, 0, 0)

# Test
dattestOr$fe_dist     <- geodist(dattestOr$latitude, dattestOr$longitude, 0, 0)
#--  Cantidad de agua disponible por persona
# Train
dattrainOrlab$fe_cant_agua <- ifelse(dattrainOrlab$population == 0,
                                     0,
                                    round(dattrainOrlab$amount_tsh /
                                          dattrainOrlab$population, 3)
                              )

# Test
dattestOr$fe_cant_agua    <- ifelse(dattestOr$population == 0,
                                    0,
                                    round(dattestOr$amount_tsh /
                                          dattestOr$population, 3)
                              )

Una vez creadas las variables, entrenamos nuevamente el modelo random forest:

# Incluimos la variable objetivo
dattrainOrlab$status_group <- status_group_vector
dattrainOrlab$status_group <- as.factor(dattrainOrlab$status_group)
my_model_2 <- fit_random_forest(formula, dattrainOrlab)

A simple vista, el valor obtenido en el conjunto de entrenamiento mejora de 0.810 en el caso anterior a 0.812. Veamos la importancia de las variables:

[1] 0.8122391
38.009 sec elapsed
save_importance_ggplot(my_model_2, "./charts/02_num_cat_menos1000_fe_cyear_dist_cant_agua.png",
                       title = "Importancia Variables Modelo 02")

De forma general, las variables creadas se sitúan entre las variables con mayor importancia (a excepción de fe_cant_agua). A continuación, si analizamos las predicciones obtenidas en el conjunto test:

my_sub_2 <- make_predictions(my_model_2, dattestOr)
# guardo submission
fwrite(my_sub_2, file = "./submissions/02_num_cat_menos1000_fe_cyear_dist_cant_agua.csv")
Train.accuracy Driven.Data
Num + Cat (> 1 & < 1000) sin duplicados 0.8168687 0.8128
Num + Cat (> 1 & < 1000) sin duplicados longitude-construction_year-gps_height imp. 0.8101178 0.8096
Anterior, eliminando columnas ’_na’ 0.8105724 0.8081
Num + Cat (> 1 & < 1000) fe cyear + dist + cant_agua 0.8122391 0.8174

2.3 Inclusión de date_recorded

Con la inclusión de las nuevas variables, el resultado mejoró con respecto a los dos submission anteriores, por lo que se decidió mantener cada una ellas. En relación con la siguiente prueba, se tomó la decisión de incluir una variable adicional que hasta el momento no se había tenido en cuenta: date_recorded. Dado que la variable se encuentra en formato fecha, en lugar de incluir directamente el campo se dividió en un total de tres nuevas variables, gracias al paquete lubridate:

  1. Año
  2. Mes
  3. Diferencia entre la fecha de observación y la fecha de construcción de la bomba, de forma que se pueda tener una columna que indique su antigüedad, en relación con su fecha de instalación:
#-- date_recorded --> año
#   Train
dattrainOrlab$fe_dr_year <- year(dattrainOr$date_recorded)

#   Test
dattestOr$fe_dr_year     <- year(dattestOr$date_recorded)

En relación a la diferencia entre date_recorded y construction_year, llamó la atención la diferencia negativa entre algunas observaciones, es decir, bombas de agua en las que la fecha de registro fue posterior a la fecha de su construcción:

#-- date_recorded --> año date_recorded - año construction_date
#   Train
dattrainOrlab$fe_dr_year_cyear_diff <- dattrainOrlab$fe_dr_year - dattrainOrlab$construction_year

#   Test
dattestOr$fe_dr_year_cyear_diff     <- dattestOr$fe_dr_year - dattestOr$construction_year
#-- date_recorded --> mes
#   Train
dattrainOrlab$fe_dr_month <- month(dattrainOr$date_recorded)

#   Test
dattestOr$fe_dr_month     <- month(dattestOr$date_recorded)
sum(dattrainOrlab$fe_dr_year_cyear_diff < 0) # Train
## [1] 9
sum(dattestOr$fe_dr_year_cyear_diff < 0)     # Test
## [1] 3

Nuevamente, tras crear las variables entrenamos el modelo random forest:

#-- Modelo 3
my_model_3 <- fit_random_forest(formula, dattrainOrlab)
[1] 0.8114983
46.337 sec elapsed

En relación al conjunto train anterior, ha mejorado considerablemente ¿Y en relación gráfico de importancia?

save_importance_ggplot(my_model_3, "./charts/03_num_cat_menos1000_fe_cyear_dist_cant_agua_dr_year_month_old.png",
                       title = "Importancia Variables Modelo 03")

Salvo la variable fe_dr_year (año de date_recorded), por lo general son variables relevantes. Además, si analizamos el conjunto test:

my_sub_3 <- make_predictions(my_model_3, dattestOr)
# guardo submission
fwrite(my_sub_3, 
    file = "./submissions/03_num_cat_menos1000_fe_cyear_dist_cant_agua_dr_year_month_old.csv")
Train.accuracy Driven.Data
Num + Cat (> 1 & < 1000) sin duplicados 0.8168687 0.8128
Num + Cat (> 1 & < 1000) sin duplicados longitude-construction_year-gps_height imp. 0.8101178 0.8096
Anterior, eliminando columnas ’_na’ 0.8105724 0.8081
Num + Cat (> 1 & < 1000) fe cyear + dist + cant_agua 0.8122391 0.8174
Num + Cat (> 1 & < 1000) fe cyear + dist + cant_agua + dr_year + dr_month + dr_year_cyear_diff 0.8114983 0.8176

Observamos que el modelo ha mejorado a 0.8176. Sin embargo, no todas las variables son de especial importancia, como es el caso de fe_dr_year (ver gráfica anterior) ¿Y si eliminamos dicha variable?

dattrainOrlab$fe_dr_year <- NULL
dattestOr$fe_dr_year     <- NULL
# Calculamos el nuevo modelo
my_model_3_sin_dr_year <- fit_random_forest(formula, dattrainOrlab)
[1] 0.8124579
41.523 sec elapsed
my_sub_3_sin_dr_year <- make_predictions(my_model_3_sin_dr_year, dattestOr)
# guardo submission
fwrite(my_sub_3_sin_dr_year, 
    file = "./submissions/03_num_cat_menos1000_fe_cyear_dist_cant_agua_month_old.csv")
Train.accuracy Driven.Data
Num + Cat (> 1 & < 1000) sin duplicados 0.8168687 0.8128
Num + Cat (> 1 & < 1000) sin duplicados longitude-construction_year-gps_height imp. 0.8101178 0.8096
Anterior, eliminando columnas ’_na’ 0.8105724 0.8081
Num + Cat (> 1 & < 1000) fe cyear + dist + cant_agua 0.8122391 0.8174
Num + Cat (> 1 & < 1000) fe cyear + dist + cant_agua + dr_year + dr_month + dr_year_cyear_diff 0.8114983 0.8176
Num + Cat (> 1 & < 1000) fe cyear + dist + cant_agua + dr_month + dr_year_cyear_diff 0.8124579 0.8176 
#-- Guardamos el mejor dataset para la mejora del modelo tras el concurso
fwrite(dattrainOrlab, "./data/train_values_concurso.csv")
fwrite(dattestOr    , "./data/test_values_concurso.csv" )

Curiosamente, eliminando el año de date_recorded no hace empeorar el score, por lo que se decidió no tener en cuenta dicha variable, al no aportar mejoría alguna.

2.4 Tuneo modelo final

Como última submission realizada en clase se realizó un tuneo de hiperparámetros del modelo random forest, empleando todos los parámetros utilizados hasta el momento.

#-- Tuneo del numero de arboles y mtry
my_ntree <- c(500, 600, 700, 800)
my_mtry  <- c(5, 6, 7)
my_pars  <- expand.grid(my_ntree, my_mtry)
names(my_pars) <- c('myntree', 'mymtry')
my_pars$accuracy <- 0

for (i in 1:nrow(my_pars)) {
   my_model_tuning <- fit_random_forest(formula, dattrainOrlab,
                                        num_trees = my_pars$myntree[[i]],
                                        mtry = my_pars$mymtry[i])
   my_pars$accuracy[i] <- (1 - my_model_tuning$prediction.error)
}
rm(my_model_tuning)

Tras realizar el tuneo de hiperparámetros, comprobamos qué modelo obtiene un mayor accuracy:

load("./rdata/00.RData")
my_pars[which(my_pars$accuracy == max(my_pars$accuracy)), ]
##    myntree mymtry  accuracy
## 12     800      7 0.8135017
#-- Entrenamos el modelo con 800 arboles y mtry = 7
my_model_4 <- fit_random_forest(formula, dattrainOrlab, num_trees = 800, mtry = 7)
[1] 0.8135017
88.765 sec elapsed

Finalmente, analizamos el score de las predicciones obtenidas:

my_sub_4 <- make_predictions(my_model_4, dattestOr)
# guardo submission
fwrite(my_sub_4, 
    file = "./submissions/04_num_cat_menos1000_fe_tunned.csv")
Train.accuracy Driven.Data
Num + Cat (> 1 & < 1000) sin duplicados 0.8168687 0.8128
Num + Cat (> 1 & < 1000) sin duplicados imp 0.8101178 0.8096
Num + Cat (> 1 & < 1000) sin duplicados imp sin cols na 0.8105724 0.8081
Num + Cat (> 1 & < 1000) fe cyear + dist + cant_agua 0.8122391 0.8174
Num + Cat (> 1 & < 1000) fe cyear + dist + cant_agua + dr_year + dr_month + dr_year_cyear_diff 0.8114983 0.8176
Num + Cat (> 1 & < 1000) fe cyear + dist + cant_agua + dr_month + dr_year_cyear_diff 0.8124579 0.8176
Num + Cat (> 1 & < 1000) fe + tunning random forest 0.8135017 0.8174

Realizando el tuneo del Random Forest, el score del modelo no consigue mejorar.

3. Mejora del modelo tras el concurso

3.1 Feature Engineering

Tras el concurso, por cuenta propia se ha decidido partir del mejor modelo obtenido hasta el momento, correspondiente con:

  1. Variables numéricas
  2. Variables categóricas (> 1 y < 1000 categorías)
  3. Feature Engineering: cyear, dist, cant_agua, dr_month y dr_year_cyear_diff
  4. Empleando un modelo Random Forest estándar: 500 árboles y mtry = 6 (raíz cuadrada del número de variables)

Sobre dicho modelo, se ha obtenido una mejor puntuación de 0.8124579 en el conjunto train y 0.8176 en el conjunto de prueba, por lo que el objetivo consistirá en superar dicho accuracy, así como la mejor puntuación obtenida por el grupo ganador (0.8248).

Importante: con el objetivo de evitar redundancia de código, el conjunto de datos obtenido en la fase anterior se ha almacenado en ficheros csv denominados train_values_concurso.csv y test_values_concurso.csv. Además, de ahora en adelante se trabajará con data.table en lugar de data.frames, con fines únicamente didácticos (aprender una nueva estructura de datos):

# Cargamos el conjunto de datos con el que se ha obtenido un mejor score en clase
dattrainOrlab    <- fread(file = "./data/train_values_concurso.csv", data.table = FALSE )
dattestOr        <- fread(file = "./data/test_values_concurso.csv",  data.table = FALSE )

Una vez cargados ambos datos (entrenamiento y prueba), se tomó la decisión de juntar ambos ficheros en una única variable (dat_completo), evitando de este modo la duplicidad en el código. Por tanto, y dado que el conjunto test no dispone de la columna status_group, almacenamos la variable objetivo del conjunto train en un vector:

vector_status_group        <- dattrainOrlab$status_group
dattrainOrlab$status_group <- NULL

# Creamos el dataset completo (train y test unidos)
columnas_test  <- names(dattestOr)[names(dattestOr) %in% names(dattrainOrlab)]
datcompleto    <- rbind(dattrainOrlab, dattestOr[, columnas_test])

#-- Nos traemos funder, ward (menos de 2100 categorias)
#   Para el siguiente apartado
datcompleto$funder <- datcat_mas_1000_mas_logicas_dr$funder
datcompleto$ward   <- datcat_mas_1000_mas_logicas_dr$ward

# Lo convertimos a data.table
datcompleto <- as.data.table(datcompleto)

Además, de cara a la posterior separación del conjunto train y test para el entrenamiento, almacenamos en una variable el número de fila a partir del cual comienza el conjunto test:

# El conjunto test empieza a partir de la fila n. 59401
# 50785: primer id del conjunto test
fila_test <- which(datcompleto$id == 50785)

3.1.1 Variables < 2100 valores únicos

Hasta ahora, hemos incluido variables con menos de 1000 categorias únicas, esto es, hasta la variable lga. Durante el concurso, no se tuvieron en cuenta, dado que el objetivo fue maximizar el score del modelo sin necesidad de recurrir a variables con un alto número de categorías. Por ello, de cara a la parte individual se tomó la primera decisión de trabajar sobre dichos campos.

Con el objetivo de seguir un orden y no incluir todas las variables “de forma abrupta”, comenzamos con aquellas con menos de 2100 categorías, esto es, funder y ward:

#-- Niveles de las categoricas.
datcat_df <- as.data.frame(datcompleto %>% select(where(is.character)))

numlev_df <- data.frame()
for (i in 1:ncol(datcat_df)) {
  col_tmp <- datcat_df[, i]
  num_lev <- length(unique(col_tmp))
  numlev_df[i, 1] <- names(datcat_df)[i]
  numlev_df[i, 2] <- num_lev
}
names(numlev_df) <- c('vars', 'levels')
numlev_df %>% arrange(levels)
##                     vars levels
## 1           source_class      3
## 2       management_group      5
## 3               quantity      5
## 4          quality_group      6
## 5  waterpoint_type_group      6
## 6  extraction_type_class      7
## 7                payment      7
## 8            source_type      7
## 9        waterpoint_type      7
## 10         water_quality      8
## 11                 basin      9
## 12                source     10
## 13            management     12
## 14     scheme_management     13
## 15 extraction_type_group     13
## 16       extraction_type     18
## 17                region     21
## 18                   lga    125
## 19                  ward   2098
## 20                funder   2141

Dado que existen demasiadas categorías, podemos preguntarnos si realmente existen 2098 y 2141 valores unicos en ward y funder, respectivamente. De hecho, podemos comprobar como algunas de las categorías se diferencian por tan solo una o dos palabras:

Lgcdg   - Lgcbg  # En el caso de funder
Magole  - Magoma # En el caso de ward

Es decir, se tratan de regiones muy similares desde un punto de vista léxico. Por tanto, ¿quién nos asegura que dos o más categorías se diferencien por tan solo un espacio en blanco, una palabra en mayúscula o por signos de puntuación? Como consecuencia, antes de procesar un modelo se procedió a realizar una limpieza sobre ambas columnas, con el propósito de comprobar si se reduce el número de valores únicos. Para ello, empleamos la librería stringi, creando una función que permita filtrar:

  1. Mayúsculas
  2. Espacios en blanco
  3. Signos de puntuación
#-- Funcion para limpieza de strings (vectorizado)
clean_text <- function(text) {
  stri_trans_tolower(
    stri_replace_all_regex(
      text, 
      pattern = "[ +\\p{Punct}]", 
      replacement = ""
      )
    )
}

Una vez creada la función, lo aplicamos sobre el data.table:

datcompleto[, fe_funder := clean_text(funder)][, fe_ward := clean_text(ward)]

Tras aplicar dicha función, comprobamos si el número de categorías se ha visto reducido:

datcat_df <- as.data.frame(datcompleto %>% select(where(is.character)))

numlev_df <- data.frame()
for (i in 1:ncol(datcat_df)) {
  col_tmp <- datcat_df[, i]
  num_lev <- length(unique(col_tmp))
  numlev_df[i, 1] <- names(datcat_df)[i]
  numlev_df[i, 2] <- num_lev
}
names(numlev_df) <- c('vars', 'levels')
numlev_df %>% arrange(levels)
##                     vars levels
## 1           source_class      3
## 2       management_group      5
## 3               quantity      5
## 4          quality_group      6
## 5  waterpoint_type_group      6
## 6  extraction_type_class      7
## 7                payment      7
## 8            source_type      7
## 9        waterpoint_type      7
## 10         water_quality      8
## 11                 basin      9
## 12                source     10
## 13            management     12
## 14     scheme_management     13
## 15 extraction_type_group     13
## 16       extraction_type     18
## 17                region     21
## 18                   lga    125
## 19               fe_ward   2096
## 20                  ward   2098
## 21             fe_funder   2110
## 22                funder   2141

Analizando la salida anterior, podemos comprobar como el número de categorías se ha visto reducido ligeramente en ambas columnas:

  1. ward : de 2098 a 2096
  2. funder: de 2141 a 2110

Una vez creadas las nuevas columnas, eliminamos las originales y lanzamos el modelo:

datcompleto[, c("ward", "funder") := NULL]
#-- Modelo
# Dividimos entre conjunto de entrenamiento y prueba
train <- datcompleto[c(1:fila_test-1),]
train$status_group <- vector_status_group
train$status_group <- as.factor(train$status_group)

test <- datcompleto[c(fila_test:nrow(datcompleto)),]

formula   <- as.formula("status_group~.")

# Empleamos la funcion fit_random_forest definida anteriormente
my_model_5 <- fit_random_forest(formula, train)
[1] 0.8149832
39.414 sec elapsed

Tras crear el modelo, echemos un vistazo a la importancia de las variables:

Observamos que ambas columnas se sitúan entre las variables con mayor relevancia (importancia en torno a 1000). Veamos si con el conjunto test obtenemos un mejor score:

my_sub_5 <- make_predictions(my_model_5, dattestOr)
# guardo submission
fwrite(my_sub_5, 
    file = "./submissions/05_04_mas_fe_inst_funder_scheme_ward.csv")
Train.accuracy Driven.Data
Mejor accuracy en el concurso
0.8176
Num + Cat (> 1 & < 2100) fe anteriores + fe_funder + fe_ward 0.8149832 0.8197

Efectivamente, hemos conseguido un mejor score que el obtenido en clase.

# Guardamos el dataset
fwrite(datcompleto, file = "./data/datcompleto_clean_text.csv")

3.1.2 Lumping sobre ward y funder (mediana)

Tras la mejora en el modelo anterior, todavía quedan demasiadas categorías en ward y funder, concretamente 2098 y 2110 categorías, respectivamente. Además, si echamos un vistazo a la proporción de aparición de las diferentes categorías de cada variable, mediante la función prop.table:

prop_categorias_ward   <- summary(c(prop.table(table(datcompleto[, fe_ward])))) # fe_ward
prop_categorias_ward
##      Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
## 1.347e-05 1.751e-04 3.367e-04 4.771e-04 6.330e-04 5.199e-03
prop_categorias_funder <- summary(c(prop.table(table(datcompleto[, fe_funder])))) # fe_funder
prop_categorias_funder
##      Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
## 1.347e-05 1.347e-05 1.347e-05 4.739e-04 9.428e-05 1.522e-01

Observamos que muchas de las categorías aperecen con muy poca frecuencia. Por ejemplo, si nos fijamos en la mediana vemos que la mitad de las categorías de fe_ward aparecen con una proporción igual o inferior a 3.37e-04, así como en fe_funder, con una proporción aún más pequeña (del orden de 1e-5).

Por tanto, una posibilidad sería aplicar lumping sobre aquellas categorías con una proporción de aparición muy pequeña, agrupándose en torno a una única categoría ( other ). En un primer comienzo, podemos aplicarlo sobre la mediana. Por ejemplo, en el caso de la variable ward, dado que la mediana se sitúa en torno a 3.367e-04, aplicamos lumping a partir de aquellas categorías con una proporción inferior a 4e-04, de forma aproximada, haciendo uso de la función fct_lump_prop del paquete forcats:

#-- fe_ward
datcompleto[, fe_ward := fct_lump_prop(datcompleto[,fe_ward], 4e-04, other_level = "other")]
datcompleto$fe_ward <- as.character(datcompleto$fe_ward)

# Pasamos de 2096 a 904 categorias (menos de la mitad de categorias)
sum(length(unique(datcompleto[, fe_ward])))
## [1] 904

En el caso de funder, aproximamos a 2e-05:

#-- fe_funder
datcompleto[, fe_funder := fct_lump_prop(datcompleto[,fe_funder], 2e-05, other_level = "other")]
datcompleto$fe_funder <- as.character(datcompleto$fe_funder)

# Pasamos de 2110 a 999 categorias (menos de la mitad de categorias)
sum(length(unique(datcompleto[, fe_funder])))
## [1] 999

Como podemos comprobar, conseguimos reducir significativamente el número de categorías en ambas variables. Una vez aplicado lumping sobre dichas columnas, lanzamos nuevamente un modelo Random Forest:

#-- Modelo
# Dividimos entre conjunto de entrenamiento y prueba
train <- datcompleto[c(1:fila_test-1),]
train$status_group <- vector_status_group
train$status_group <- as.factor(train$status_group)

test <- datcompleto[c(fila_test:nrow(datcompleto)),]

formula   <- as.formula("status_group~.")

# Empleamos la funcion fit_random_forest definida anteriormente
my_model_6 <- fit_random_forest(formula, train)
[1] 0.8159764
39.081 sec elapsed

En relación con el modelo anterior, el hecho de disminuir el número de categorías ha mejorado el resultado en el conjunto train (pasando de 0.8149 a 0.8159). A continuación, analicemos el conjunto test:

my_sub_6 <- make_predictions(my_model_6, dattestOr)
# guardo submission
fwrite(my_sub_6, 
    file = "./submissions/06_lumping_y_fe_sobre_funder_ward_mediana.csv")
Train.accuracy Driven.Data
Mejor accuracy en el concurso
0.8176
Num + Cat (> 1 & < 2100) fe anteriores + fe_funder + fe_ward 0.8149832 0.8197
Num + Cat (> 1 & < 2100) fe anteriores + lumping sobre funder + ward (mediana) 0.8159764 0.8212

Observamos que incluso en el conjunto test el score ha aumentado a 0.8212.

3.1.3 Lumping sobre ward y funder (primer y tercer cuartil)

En el apartado anterior, aplicamos lumping sobre aquellas categorías cuya proporción de aparición se situase por debajo de la mediana, aproximadamente. No obstante, ¿Y si aplicamos el punto de corte sobre el primer o tercer cuartil? Es decir, si recordamos la salida del summary anterior:

prop_categorias_ward   # fe_ward
##      Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
## 1.347e-05 1.751e-04 3.367e-04 4.771e-04 6.330e-04 5.199e-03
prop_categorias_funder # fe_funder
##      Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
## 1.347e-05 1.347e-05 1.347e-05 4.739e-04 9.428e-05 1.522e-01

Podríamos aplicar lumping sobre el primer cuartil, de forma aproximada, o bien sobre el tercer cuartil. De este modo, podemos comprobar si agrupando un mayor o menor número de categorías conseguimos mejorar la puntuación (“jugar con el punto de corte”). Para ello, recuperamos nuevamente el dataset obtenido tras aplicar la función clean_text:

datcompleto <- fread("./data/datcompleto_clean_text.csv")

3.1.3.1 Sobre el primer cuartil

#-- fe_ward (primer cuartil: 1.751e-04, aplicamos el punto de corte sobre 2e-04 aprox.)
datcompleto[, fe_ward := fct_lump_prop(datcompleto[,fe_ward], 2e-04, other_level = "other")]
datcompleto$fe_ward <- as.character(datcompleto$fe_ward)

# Pasamos de 2096 a 1483 categorias
sum(length(unique(datcompleto[, fe_ward])))
## [1] 1483
#-- fe_funder (primer cuartil: 1.347e-05, aplicamos el punto de corte sobre 2e-05 aprox.)
datcompleto[, fe_funder := fct_lump_prop(datcompleto[,fe_funder], 2e-05, other_level = "other")]
datcompleto$fe_funder <- as.character(datcompleto$fe_funder)

# Pasamos de 2110 a 999 categorias (primer cuartil = mediana)
sum(length(unique(datcompleto[, fe_funder])))
## [1] 999

Una vez aplicado lumping, recalculamos el modelo:

#-- Modelo
# Dividimos entre conjunto de entrenamiento y prueba
train <- datcompleto[c(1:fila_test-1),]
train$status_group <- vector_status_group
train$status_group <- as.factor(train$status_group)

test <- datcompleto[c(fila_test:nrow(datcompleto)),]

formula   <- as.formula("status_group~.")
# Resultado train: 0.8159259
my_model_7 <- fit_random_forest(formula,
                                train)

my_sub_7 <- make_predictions(my_model_7, test)

fwrite(my_sub_7, file = "./submissions/07_lumping_y_fe_sobre_funder_ward_primer_cuartil.csv")
# Score: 0.8196
Train.accuracy Driven.Data
lumping sobre mediana 0.8159764 0.8212
lumping sobre primer cuartil 0.8159259 0.8196

Observamos que manteniendo un mayor número de categorías, el score obtenido tanto en train como en test empeora.

3.1.3.2 Sobre el tercer cuartil

datcompleto <- fread("./data/datcompleto_clean_text.csv")
#-- fe_ward (primer cuartil: 6.330e-04, aplicamos el punto de corte sobre 7e-04 aprox.)
datcompleto[, fe_ward := fct_lump_prop(datcompleto[,fe_ward], 7e-04, other_level = "other")]
datcompleto$fe_ward <- as.character(datcompleto$fe_ward)

# Pasamos de 2096 a 452 categorias
sum(length(unique(datcompleto[, fe_ward])))
## [1] 452
#-- fe_funder (primer cuartil: 9.428-05, aplicamos el punto de corte sobre 1e-04 aprox.)
datcompleto[, fe_funder := fct_lump_prop(datcompleto[,fe_funder], 1e-04, other_level = "other")]
datcompleto$fe_funder <- as.character(datcompleto$fe_funder)

# Pasamos de 2110 a 502 categorias
sum(length(unique(datcompleto[, fe_funder])))
## [1] 502
#-- Modelo
# Dividimos entre conjunto de entrenamiento y prueba
train <- datcompleto[c(1:fila_test-1),]
train$status_group <- vector_status_group
train$status_group <- as.factor(train$status_group)

test <- datcompleto[c(fila_test:nrow(datcompleto)),]

formula   <- as.formula("status_group~.")
# Resultado train: 0.8146633
my_model_8 <- fit_random_forest(formula,
                                train)

my_sub_8 <- make_predictions(my_model_8, test)

fwrite(my_sub_8, file = "./submissions/08_lumping_y_fe_sobre_funder_ward_tercer_cuartil.csv")
# Score: 0.8203
Train.accuracy Driven.Data
lumping sobre mediana 0.8159764 0.8212
lumping sobre primer cuartil 0.8159259 0.8196
lumping sobre tercer cuartil 0.8146633 0.8203

Comprobamos que, pese a mejorar el score del modelo al aplicar el punto de corte sobre el tercer cuartil, el resultado obtenido sobre la mediana continua siendo mejor.

Por tanto, sobre las variables fe_ward y fe_funder aplicamos lumping sobre la mediana de proporción de aparición de cada una de sus categorías, reduciendo su número a 904 y 999, respectivamente.

Train.accuracy Driven.Data
Mejor accuracy en el concurso
0.8176
Num + Cat (> 1 & < 2100) fe anteriores + fe_funder + fe_ward 0.8149832 0.8197
Num + Cat (> 1 & < 2100) fe anteriores + lumping sobre funder + ward (mediana) 0.8159764 0.8212
Num + Cat (> 1 & < 2100) fe anteriores + lumping sobre funder + ward (primer cuartil) 0.8159259 0.8196
Num + Cat (> 1 & < 2100) fe anteriores + lumping sobre funder + ward (tercer cuartil) 0.8146633 0.8203

Pese a aplicar lumping sobre ambas variables, continuamos teniendo demasiadas categorías: 904 en el caso de fe_ward y 999 en el caso de fe_funder. Como consecuencia, se plantearon diversas alternativas con el fin de reducir nuevamente el número de valores únicos sobre dichas columnas, observando con ello si el score obtenido mejora. Concretamente, se probaron cuatro técnicas mencionadas a lo largo del módulo.

  1. Hashing encoding

  2. Frecuencias absolutas (realizada por el grupo ganador)

  3. Target encoding (realizada por el grupo ganador)

  4. Word Embedding

De este modo, en caso de que alguna de las técnicas anteriores, sobre las variables ward y funder, mejore el score (0.8212), se aplicará al resto de variables categóricas. Nuevamente, recuperamos el dataset sobre el que se ha obtenido el mejor score (clean text + lumping mediana):

datcompleto <- fread("./data/datcompleto_clean_text_lumping_mediana.csv")

3.1.4 Aplicando Hashing Encoding

También conocido como Feature Hashing, consiste en transformar un elevado conjunto de características en un vector, por medio de una función de transformación o hash function. No obstante, uno de los principales incovenientes de dicha técnica reside en el ratio de colisión o collision rate, es decir, cuando dos categorías presentan el mismo resultado en la función hash. Por tanto, cuanto mayor sea el tamaño de la función hash, menor será el ratio de colisión.

Para ello, haremos uso del paquete FeatureHashing disponible en CRAN. Dicho paquete dispone de la función hash.size, el cual indica el tamaño mínimo (teórico) que permite reducir al mínimo el ratio de colisión entre valores hash:

#-- Feature Hashing (tras aplicar lumping sobre funder y ward)
# Podemos elegir el tamaño minimo (teorico) que permite reducir el ratio de colision entre valores hash
tam_minimo <- hash.size(datcompleto[, c("fe_funder", "fe_ward")])
tam_minimo
## [1] 2048

La función nos indica que el tamaño mínimo es de 2048. A continuación, mediante la función hashed.model.matrix creamos la matriz hash correspondiente con la que indexar las categorías de las columnas ward y funder:

mat_hash   <- hashed.model.matrix(~., datcompleto[, c("fe_funder", "fe_ward")], tam_minimo, create.mapping = TRUE)

# Comprobamos el ratio de colisión
mean(duplicated(hash.mapping(mat_hash))) # 0.3520757
## [1] 0.3520757

Comprobamos que el ratio de colisión obtenido no es demasiado alto (alrededor de 0.35). Una vez elaborada la matriz, debemos sustituir los categorías originales por su valor hash correspondiente, haciendo uso de la función hash.mapping, encargada de devolver la equivalencia (mapeo) entre cada categoría y su valor hash:

# Sustituimos las columnas fe_funder y fe_ward por su valor hash correspondiente
vector_hash <- hash.mapping(mat_hash)
mat_hash_dt <- data.table("feature" = names(vector_hash), 
                          "values" = vector_hash)

head(mat_hash_dt)
##                    feature values
## 1:           fe_fundernddp   1122
## 2:           fe_wardmsindo   1392
## 3: fe_fundervillagecouncil    331
## 4:          fe_wardkighare    772
## 5:          fe_funderdonor   1114
## 6:       fe_wardurushimbwe   1273

Una vez obtenida la matriz de equivalencias, sustituimos en el data.table original:

# Por defecto, hashed.model.matrix añade el nombre de columna a la variable
# De forma que debemos incluirlo tambien tanto en fe_funder como en fe_ward
# para hacer coincidir las categorias
datcompleto[, fe_ward := paste0("fe_ward", fe_ward)]
datcompleto[, fe_funder := paste0("fe_funder", fe_funder)]

datcompleto[mat_hash_dt, fe_ward_hashed := i.values,  on = .(fe_ward = feature)]
datcompleto[mat_hash_dt, fe_funder_hashed := i.values,  on = .(fe_funder = feature)]

Además, observamos cómo el número de “categorías” (valores hash) se reduce a 727 y 803 valores únicos:

#- Numero de valores unicos en fe_ward_hashed
length(unique(datcompleto$fe_ward_hashed))   # 727 categoricas
## [1] 727
#- Numero de valores unicos en fe_funder_hashed
length(unique(datcompleto$fe_funder_hashed)) # 803 categorias
## [1] 803
# Eliminamos las columnas fe_ward y fe_funder originales
datcompleto[, `:=`(fe_funder = NULL, fe_ward = NULL)]

Una vez reemplazadas las categorías, entrenamos un nuevo modelo:

my_model_9 <- fit_random_forest(formula, train)
[1] 0.8158418
35.099 sec elapsed

En relación al conjunto de entrenamiento, el accuracy obtenido es similar al obtenido mediante lumping (del orden de 0.815). A continuación, guardamos la submission y comprobamos su score en DrivenData:

my_sub_9 <- make_predictions(my_model_9, test)
# guardo submission
fwrite(my_sub_9, file = "./submissions/09_lumping_fe_y_hashing_sobre_funder_ward.csv")
Train.accuracy Driven.Data
Mejor accuracy en el concurso
0.8176
Num + Cat (> 1 & < 2100) fe anteriores + fe_funder + fe_ward 0.8149832 0.8197
Num + Cat (> 1 & < 2100) fe anteriores + lumping sobre funder + ward (mediana) 0.8159764 0.8212
Num + Cat (> 1 & < 2100) fe anteriores + lumping sobre funder + ward (primer cuartil) 0.8159259 0.8196
Num + Cat (> 1 & < 2100) fe anteriores + lumping sobre funder + ward (tercer cuartil) 0.8146633 0.8203
Num + Cat (> 1 & < 2100) fe anteriores + lumping (mediana) + hashed sobre funder + ward 0.8158418 0.8214

Sorprendentemente, aplicando Feature Hashing conseguimos mejorar el score del modelo de 0.8212 a 0.8214.

Dado que 2048 es el tamaño mínimo que permite reducir el ratio de colisión entre valores hash, ¿Y si reducimos dicho tamaño a la mitad, 1024? ¿Mejorará el modelo al reducir el número de categogrías únicas?

mat_hash_2 <- hashed.model.matrix(~., datcompleto[, c("fe_funder", "fe_ward")], 2^10, create.mapping = TRUE, )
mean(duplicated(hash.mapping(mat_hash_2))) # 0.5491329
## [1] 0.5491329

Lógicamente, el ratio de colisión aumenta al disminuir el tamaño (de 0.35 a 0.55, aproximadamente). Si entrenamos nuevamente el modelo con la nueva matriz hash:

# Una vez realizados los mismos pasos que en el caso anterior...
my_model_9 <- fit_random_forest(formula, train)
[1] 0.8162121
35.918 sec elapsed

Aparentemente, el valor obtenido en el conjunto de entrenamiento mejora ligeramente (de 0.815 a 0.816). Sin embargo, al predecir el conjunto test, observamos que el score empeora:

my_sub_9 <- make_predictions(my_model_9, test)
# guardo submission
fwrite(my_sub_9, file = "./submissions/09_lumping_fe_y_hashing_sobre_funder_ward_2_up_10.csv")
Train.accuracy Driven.Data
Hashing (2048 hash size) 0.8158418 0.8214
Hashing (1024 hash size) 0.8162121 0.8197

3.1.5 Aplicando Frecuencias Absolutas

Otra alternativa a Feature Hashing consiste en sustituir cada categoría en fe_funder y fe_ward por su frecuencia de aparición correspondiente, un tipo de codificación con el que obtuvo buenos resultados el grupo ganador y, por ello, debemos comprobar si mejora el score obtenido en Feature Hashing (0.8214):

#- fe_ward
datcompleto[, fe_ward_freq := as.numeric(.N), by = fe_ward]
length(unique(datcompleto[, fe_ward_freq]))  # 129 categorias
## [1] 129
#-  fe_funder
datcompleto[, fe_funder_freq := as.numeric(.N), by = fe_funder]
length(unique(datcompleto[, fe_funder_freq])) # 169 categorias
## [1] 169
# Eliminamos las columnas originales
datcompleto[, `:=`(fe_funder = NULL, fe_ward = NULL)]

Observamos que el número de categorías en ambas variables se ve reducido de 904 y 999 a 129 y 169 valores únicos, respectivamente. A continuación, dividimos el conjunto de datos entre train y test, entrenando un nuevo modelo:

my_model_10 <- fit_random_forest(formula, train)
[1] 0.8154882
38.746 sec elapsed

En relación al accuracy obtenido en el entrenamiento, el porcentaje es ligeramente inferior con respecto al obtenido en el apartado anterior (0.8158418). No obstante, si realizamos la predicción del conjunto test:

Train.accuracy Driven.Data
Mejor accuracy en el concurso
0.8176
Num + Cat (> 1 & < 2100) fe anteriores + fe_funder + fe_ward 0.8149832 0.8197
Num + Cat (> 1 & < 2100) fe anteriores + lumping sobre funder + ward (mediana) 0.8159764 0.8212
Num + Cat (> 1 & < 2100) fe anteriores + lumping sobre funder + ward (tercer cuartil) 0.8146633 0.8203
Num + Cat (> 1 & < 2100) fe anteriores + lumping sobre funder + ward (primer cuartil) 0.8159259 0.8196
Num + Cat (> 1 & < 2100) fe anteriores + lumping (mediana) + hashed sobre funder + ward 0.8158418 0.8214
Num + Cat (> 1 & < 2100) fe anteriores + lumping (mediana) + freq. abs. sobre funder + ward 0.8154882 0.8216

Comprobamos que sustituyendo las categorías de fe_ward y fe_funder por sus frecuencias absolutas (en lugar de aplicar FeatureHashing), el score obtenido aumenta de 0.8214 a 0.8216.

3.1.6 Aplicando Target Encoding

3.1.6.1 Empleando la media

Otra técnica propuesta para la codificación de fe_ward y fe_funder consiste en reemplazar cada categoría por la agregación de la variable objetivo (reemplazando la variables objetivo por valores numéricos: 0, 1 y 2), sustituyendo cada categoría en fe_ward y fe_funder por el promedio de dicha variable. A modo de ejemplo:

   student grades grades_mean_by_student
1:   Marie      0                   0.75 = (0 + 0 + 1 + 2 / 4)
2:   Marie      0                   0.75
3:  Pierre      1                   1.00
4:   Louis      2                   1.75
5:   Louis      2                   1.75
6:   Marie      1                   0.75
7:   Marie      2                   0.75
8:   Louis      2                   1.75
9:   Louis      1                   1.75

Para aplicar target_encoding, recuperamos únicamente las variables categóricas fe_ward y fe_funder:

# Recuperamos las columnas categoricas
datcat_df <- as.data.frame(datcompleto %>% select(fe_funder, fe_ward))

En muchas ocasiones, se recomienda añadir un valor de “ruido aleatorio” a las columnas del conjunto de entrenamiento por lo que separaremos nuevamente datcompleto en entrenamiento y prueba:

# Dividimos datcompleto en entrenamiento y prueba
train <- datcompleto[c(1:fila_test-1),]
train$status_group <- vector_status_group
train$status_group <- as.factor(train$status_group)

test <- datcompleto[c(fila_test:nrow(datcompleto)),]

A continuación, de cara la codificación convertimos la variable objetivo en formato numérico:

# 0: functional ;  1: functional needs repair ; 2: non functional
train[, status_group := ifelse(status_group == "functional", 0, 
                               ifelse(status_group == "functional needs repair", 1
                                      , 2))]

Para realizar Target Encoding, empleamos la función build_target_encoding disponible en el paquete DataPreparation, sobre la que debemos indicar el conjunto de entrenamiento, el nombre las columnas a codificar, la variable objetivo, así como la función de agregación (en este caso la media):

#-- Target encoding (mean)
dat_encoded <- build_target_encoding(train, 
                                     cols_to_encode = names(datcat_df),
                                     target_col = "status_group",
                                     functions = "mean",
                                     verbose = TRUE)
## [1] "build_target_encoding: Start to compute encoding for target_encoding according to col: status_group."

Dicha función devuelve una lista con la correspondencia entre cada categoría y su valor codificado mediante Target Encoding. Veamos un ejemplo:

# Creamos un data.table para observar los resultados obtenidos tras aplicar Target Encoding
dat_encoded_dt          <- as.data.table(dat_encoded)
## Warning in as.data.table.list(dat_encoded): Item 2 has 904 rows but longest item
## has 989; recycled with remainder.
names(dat_encoded_dt)   <- stri_replace_all_regex(names(dat_encoded_dt), 
                                                  pattern = ".*\\.",
                                                  replacement = "")

dat_encoded_dt[1, ] # Veamos un ejemplo
##    fe_funder status_group_mean_by_fe_funder  fe_ward
## 1:     roman                      0.3236364 mundindi
##    status_group_mean_by_fe_ward
## 1:                         0.75

En la salida anterior, observamos que la categoría roman presenta un promedio de 0.32 con respecto a la variable objetivo, esto es, el estado de las bombas de agua creadas por dicha empresa, de media, es “no funcional” (0). Por el contrario, en la localidad de mundindi las bombas de agua presentan un promedio más cercano a 1 (0.75).

Tras obtener el Target Encoding de cada categoría, debemos sustituir el valor de cada categoría por su agregación correspondiente, tanto en el conjunto train como test. Para ello, empleamos la función target_encode:

train <- target_encode(train, target_encoding = dat_encoded)
## [1] "target_encode: Start to encode columns according to target."
test <- target_encode(test, target_encoding = dat_encoded)
## [1] "target_encode: Start to encode columns according to target."
sum(is.na(test)) # 21 valores missing en status_group_by_fe_funder
## [1] 21

En el caso del conjunto test, detectamos valores missing en la columna status_group_by_fe_funder. Como consecuencia, y dado que son pocos valores, se tomó la decisión de reeemplazarlos directamente por el valor medio:

# Reemplazamos los valores nulos por la media
media <- mean(test[, status_group_mean_by_fe_funder], na.rm = TRUE)
setnafill(test, cols="status_group_mean_by_fe_funder", fill = media)

A continuación, añadimos un valor de ruido al conjunto train, además de recodificar la variable objetivo y eliminar las columnas fe_ward y fe_funder originales:

# Añadimos ruido a los datos train
cat_cols <- paste0('status_group_mean_by_',names(datcat_df))

train[, cat_cols] <- train[, lapply(.SD, 
                                    function(x) x * rnorm(length(x), mean = 1, 
                                    sd = 0.05)), .SDcols = cat_cols]

# Recodificamos la variable objetivo
train[, status_group := ifelse(status_group == 0, "functional", 
                               ifelse(status_group == 1, "functional needs repair"
                                      , "non functional"))]

# Eliminamos las variables originales
train[, `:=`(fe_funder = NULL, fe_ward = NULL)]
test[, `:=`(fe_funder = NULL, fe_ward = NULL)]

Finalmente, entrenamos el modelo y realizamos las predicciones sobre el conjunto test:

#-- Modelo
formula   <- as.formula("status_group~.")
# 0.8154377
my_model_11 <- fit_random_forest(formula,
                                 train)
[1] 0.8154377
39.044 sec elapsed
my_sub_11 <- make_predictions(my_model_11, test)

# guardo submission
fwrite(my_sub_11, file = "./submissions/11_lumping_fe_sobre_funder_ward_target_encoding_media.csv")
Train.accuracy Driven.Data
Mejor accuracy en el concurso
0.8176
Num + Cat (> 1 & < 2100) fe anteriores + fe_funder + fe_ward 0.8149832 0.8197
Num + Cat (> 1 & < 2100) fe anteriores + lumping sobre funder + ward (mediana) 0.8159764 0.8212
Num + Cat (> 1 & < 2100) fe anteriores + lumping sobre funder + ward (tercer cuartil) 0.8146633 0.8203
Num + Cat (> 1 & < 2100) fe anteriores + lumping sobre funder + ward (primer cuartil) 0.8159259 0.8196
Num + Cat (> 1 & < 2100) fe anteriores + lumping (mediana) + hashed sobre funder + ward 0.8158418 0.8214
Num + Cat (> 1 & < 2100) fe anteriores + lumping (mediana) + freq. abs. sobre funder + ward 0.8154882 0.8216
Num + Cat (> 1 & < 2100) fe anteriores + lumping (mediana) + freq. abs. y target enc. (media) sobre funder + ward 0.8154377 0.8164

Como podemos comprobar en la tabla resultante, aplicando Target Encoding ha empeorado los resultados (de 0.8216 a 0.8164).

3.1.6.2 Empleando la suma

Revisando la documentación de la función build_target_encoding, me llamó la atención un ejemplo propuesto en el que se emplea la función de agregación sum en lugar de la media. A modo de ejemplo, supongamos el siguiente data.table con diferentes bombas de agua:

     fe_funder    fe_ward   status_group
1:       roman   mundindi              0
2:       other      natta              1
3:      unicef      other              1
4:      unicef      natta              2
5:       other   mundindi              0

0: non functional ;  1: functional needs repair ; 2: functional

¿Y si en lugar de aplicar la media realizamos la suma de los pesos de cada variable objetivo? Es decir, dado que cada categoría tiene un valor asignado (0, 1 o 2), una posibilidad sería calcular (sumar) el total de la variable objetivo por cada categoría, de este modo:

     fe_funder    fe_ward   status_group     fe_funder_sum     fe_ward_sum   
1:       roman   mundindi              0                 0       0 + 0 = 0   
2:       other      natta              1         1 + 0 = 1       1 + 2 = 3
3:      unicef      other              1         1 + 2 = 3               1
4:      unicef      natta              2         1 + 2 = 3       1 + 2 = 3
5:       other   mundindi              0         1 + 0 = 1       0 + 0 = 0

0: non functional ;  1: functional needs repair ; 2: functional

Podemos comprobar como empresas/entidades como unicef, así como regiones como es el caso de natta presentan un mayor número de bombas “defectuosas”, lo cual se refleja en la suma total. No obstante, un mayor valor de la suma no implica necesariamente un mayor riesgo en las bombas de agua. Un ejemplo sería el siguiente:

     fe_funder   status_group     fe_funder_sum   
1:       roman              2                 2 
2:      unicef              1         1 + 1 = 2
3:      unicef              1         1 + 1 = 2

0: non functional ;  1: functional needs repair ; 2: functional

En el caso anterior, la suma de las variables objetivo en ambos valores de fe_funder es el mismo (2). Sin embargo, mientras que en el caso de unicef las bombas continuan funcionando (aunque necesiten reparación), en el caso de roman, la única bomba que presenta no funciona. Por tanto, podemos asegurar que se trata de una representación con muy poca interpretabilidad.

No obstante, probaremos con un nuevo modelo, comprobando si los resultandos mejoran de cara al concurso:

#-- Target encoding (sum)
dat_encoded <- build_target_encoding(train, 
                                      cols_to_encode = names(datcat_df),
                                      target_col = "status_group",
                                      functions = "sum", # sustituimos "mean" por "sum"
                                      verbose = TRUE)
Train.accuracy Driven.Data
Mejor accuracy en el concurso
0.8176
Num + Cat (> 1 & < 2100) fe anteriores + fe_funder + fe_ward 0.8149832 0.8197
Num + Cat (> 1 & < 2100) fe anteriores + lumping sobre funder + ward (mediana) 0.8159764 0.8212
Num + Cat (> 1 & < 2100) fe anteriores + lumping sobre funder + ward (tercer cuartil) 0.8146633 0.8203
Num + Cat (> 1 & < 2100) fe anteriores + lumping sobre funder + ward (primer cuartil) 0.8159259 0.8196
Num + Cat (> 1 & < 2100) fe anteriores + lumping (mediana) + hashed sobre funder + ward 0.8158418 0.8214
Num + Cat (> 1 & < 2100) fe anteriores + lumping (mediana) + freq. abs. sobre funder + ward 0.8154882 0.8216
Num + Cat (> 1 & < 2100) fe anteriores + lumping (mediana) + freq. abs. y target enc. (media) sobre funder + ward 0.8154377 0.8164
Num + Cat (> 1 & < 2100) fe anteriores + lumping (mediana) + freq. abs. y target enc. (suma) sobre funder + ward 0.814899 0.8215

Curiosamente, el resultado obtenido mejora con respecto al cálculo mediante la media, aunque sigue siendo mejor (por poca diferencia), la codificación por frecuencias absolutas.

3.1.7 Aplicando Word Embedding

Sin duda, uno de los grandes avances en el campo del procesamiento del lenguaje han sido las técnicas de Word Embedding, que permiten la representación vectorial (n-dimensional) de una palabra o token. Por ello, de cara a la práctica se planteó la posibilidad de realizar una representación vectorial de las diferentes categorías en las columnas fe_ward y fe_funder, con el objetivo de observar si la precisión del modelo mejora o no.

Para ello, haremos uso de dos paquetes, pertenecientes al conjunto de librerías tidymodels:

  1. embed: permite la creación de modelos word-embedding, mediante el uso de la librería TensorFlow.

  2. recipes: se emplea fundamentalmente para el preprocesamiento de datos, mediante el uso de “recetas” definidas de forma secuencial por medios de pipelines.

En primer lugar, y una vez aplicado lumping sobre ward y funder, dividimos nuevamente los datos en entrenamiento y test:

# 1. Comprobamos que las variables categoricas estan codificadas como factor
emb_cols <- c("fe_ward", "fe_funder")
datcompleto[,(emb_cols):= lapply(.SD, as.factor), .SDcols = emb_cols]

# Dividimos en train y test
train <- datcompleto[c(1:fila_test-1),]
train$status_group <- vector_status_group
train$status_group <- as.factor(train$status_group)

test <- datcompleto[c(fila_test:nrow(datcompleto)),]

Una vez dividido el conjunto de datos, elaboramos el modelo embedded. Para ello, mediante la función recipe elaboramos una “receta” para la transformación de variables, concretamente sobre fe_ward, fe_funder y status_group. Sobre dicha receta o fórmula aplicamos la función step_embed del paquete embed, con la que realizar la transformación Word Embedding. Por defecto, el número de dimensiones es de 2, por lo que probaremos inicialmente con dicho valor:

# 2. Elaboramos el modelo embedded (empezando con dimensionalidad 2)
emb_cols_target <- c(emb_cols, "status_group")

base_recipe <- recipe(status_group ~ ., train[, ..emb_cols_target])
for(feat in emb_cols){
  base_recipe <- base_recipe %>% 
    step_embed({{feat}},
                num_terms = 2,
                outcome = vars(status_group),
                options = embed_control(epochs = 5, validation_split = 0.2)
    )
}

Finalmente, una vez elaborada la “receta” para la transformación de ambas columnas, lo aplicamos inicialmente sobre el conjunto train y posteriormente sobre el conjunto test, haciendo uso de las funciones prep y bake, respectivamente:

# 3. Creacion de las columnas Word Embeddings (train y test)
#    Sobre el conjunto train, elaboramos la "receta": word_embeddings
train_prepped <- prep(base_recipe, train[, ..emb_cols_target])
#    Sobre el conjunto test, lo aplicamos
test_prepped <-  bake(train_prepped, test[, ..emb_cols])

train_final  <-  cbind(as.data.table(juice(train_prepped)),
                      train[, setdiff(names(train), emb_cols_target), with = FALSE]
                 )

test_final   <-  cbind(test_prepped,
                     test[, setdiff(names(test), emb_cols), with = FALSE]
                 )
# Echemos un vistazo a las columnas transformadas
train_final[1, c("fe_ward_embed_1", "fe_ward_embed_2", "fe_funder_embed_1", "fe_funder_embed_2")]
##    fe_ward_embed_1 fe_ward_embed_2 fe_funder_embed_1 fe_funder_embed_2
## 1:     -0.02231506      -0.0347387        0.00927105       0.003461588

Como podemos comprobar, cada categoría se transforma en dos columnas con su representación bi-dimensional correspondiente ¿Mejorará el modelo si aplicamos Word Embedding sobre fe_ward y fe_funder?

# Train
my_model_12 <- fit_random_forest(formula, train_final)

# Predicciones
my_sub_12  <- make_predictions(my_model_12, test_final)
Train.accuracy Driven.Data
Mejor accuracy en el concurso
0.8176
Num + Cat (> 1 & < 2100) fe anteriores + fe_funder + fe_ward 0.8149832 0.8197
Num + Cat (> 1 & < 2100) fe anteriores + lumping sobre funder + ward (mediana) 0.8159764 0.8212
Num + Cat (> 1 & < 2100) fe anteriores + lumping sobre funder + ward (tercer cuartil) 0.8146633 0.8203
Num + Cat (> 1 & < 2100) fe anteriores + lumping sobre funder + ward (primer cuartil) 0.8159259 0.8196
Num + Cat (> 1 & < 2100) fe anteriores + lumping (mediana) + hashed sobre funder + ward 0.8158418 0.8214
Num + Cat (> 1 & < 2100) fe anteriores + lumping (mediana) + freq. abs. sobre funder + ward 0.8154882 0.8216
Num + Cat (> 1 & < 2100) fe anteriores + lumping (mediana) + freq. abs. y target enc. (media) sobre funder + ward 0.8154377 0.8164
Num + Cat (> 1 & < 2100) fe anteriores + lumping (mediana) + freq. abs. y target enc. (suma) sobre funder + ward 0.814899 0.8215
Num + Cat (> 1 & < 2100) fe anteriores + lumping (mediana) + word embed sobre funder y ward (dim = 2) 0.8155724 0.8204

Creando un modelo Word Embedding sobre ambas variables, el resultado no parece mejorar ¿Y si aumentamos el número de dimensiones en la representación vectorial? Normalmente, los modelos Word Embedding suelen emplear decenas e incluso cientos de dimensiones en su representación.

Train.accuracy Driven.Data
Word Embedding (Dim = 2) 0.8155724 0.8204
Word Embedding (Dim = 5) 0.8150337 0.8192
Word Embedding (Dim = 10) 0.8145286 0.8157

Pese a aumentar el número de dimensiones (y con ello la complejidad), la precisión del modelo disminuye.

Conclusión: dado el mejor score obtenido, elegimos como mejor alternativa la transformación por frecuencias absolutas, con respecto a las variables categóricas.

3.1.8 Transformar por frecuencias absolutas el resto de variables categóricas

A la vista de las transformaciones anteriores, la que mejor resultado ha obtenido ha sido mediante frecuencias absolutas. No obstante, y dado que solo las hemos aplicado sobre fe_ward y fe_funder ¿Y si la aplicamos sobre el resto de variables categóricas, concretamente de tipo carácter?

Para ello, y una vez aplicado lumping sobre ward y funder, realizamos la imputación de las variables categóricas (character) por sus correspondientes frecuencias absolutas. Inicialmente, almacenamos en una variable los nombres de las columnas a imputar:

#-- Imputacion de las variables categoricas por sus frecuencias absolutas
cat_cols <- names(datcompleto[, which(sapply(datcompleto, is.character)), with = FALSE])

Antes de realizar la imputación, guardamos en una variable el número de valores únicos de cada columna, con el fin de contrastar posteriormente la variación en el número de categorías:

#   Antes de imputar
freq_antes_fe <- apply(datcompleto[, ..cat_cols], 2, function(x) length(unique(x)))

A continuación, sustituimos el valor de cada columna por su correspondiente frecuencia absoluta, además de eliminar las columnas originales:

# Sustitucion de cada valor por su frecuencia absoluta
for (cat_col in cat_cols) {
  datcompleto[, paste0("fe_", cat_col) := as.numeric(.N), by = cat_col]
}
# fe_funder y fe_ward aparecen como fe_fe_ward y fe_fe_funder
# por lo que lo reemplazamos por fe_ward y fe_funder
names(datcompleto) <- stri_replace_all_fixed(names(datcompleto),
                                             "fe_fe_", "fe_")

# Eliminamos las columnas originales
for (cat_col in cat_cols) {
  datcompleto[, paste(cat_col) := NULL]
}

Una vez transformadas las columnas, creamos una nueva variable con el número de valores únicos de cada columna:

new_cat_cols <- paste0("fe_", stri_replace_all_fixed(cat_cols, "fe_", ""))

freq_despues_fe <- apply(datcompleto[, ..new_cat_cols], 2, function(x) length(unique(x)))

#-- Solo cambian funder, ward y lga en relacion al numero de categorias
data.table(rbind(freq_antes_fe, freq_despues_fe), keep.rownames = TRUE)
##                 rn basin region lga scheme_management extraction_type
## 1:   freq_antes_fe     9     21 125                13              18
## 2: freq_despues_fe     9     21 120                13              18
##    extraction_type_group extraction_type_class management management_group
## 1:                    13                     7         12                5
## 2:                    13                     7         12                5
##    payment water_quality quality_group quantity source source_type source_class
## 1:       7             8             6        5     10           7            3
## 2:       7             8             6        5     10           7            3
##    waterpoint_type waterpoint_type_group fe_funder fe_ward
## 1:               7                     6       999     904
## 2:               7                     6       169     129

Analizando la tabla resultante, podemos comprobar como tan solo funder, ward y lga se ven reducidos en cuanto al número de valores únicos se refiere. No obstante, veamos los resultados obtenidos tras entrenar el nuevo modelo:

my_model_13 <- fit_random_forest(formula, train)

my_sub_13   <- make_predictions(my_model_13, test)
Train.accuracy Driven.Data
Mejor accuracy en el concurso
0.8176
Num + Cat (> 1 & < 2100) fe anteriores + fe_funder + fe_ward 0.8149832 0.8197
Num + Cat (> 1 & < 2100) fe anteriores + lumping sobre funder + ward (mediana) 0.8159764 0.8212
Num + Cat (> 1 & < 2100) fe anteriores + lumping sobre funder + ward (tercer cuartil) 0.8146633 0.8203
Num + Cat (> 1 & < 2100) fe anteriores + lumping sobre funder + ward (primer cuartil) 0.8159259 0.8196
Num + Cat (> 1 & < 2100) fe anteriores + lumping (mediana) + hashed sobre funder + ward 0.8158418 0.8214
Num + Cat (> 1 & < 2100) fe anteriores + lumping (mediana) + freq. abs. sobre funder + ward 0.8154882 0.8216
Num + Cat (> 1 & < 2100) fe anteriores + lumping (mediana) + freq. abs. y target enc. (media) sobre funder + ward 0.8154377 0.8164
Num + Cat (> 1 & < 2100) fe anteriores + lumping (mediana) + freq. abs. y target enc. (suma) sobre funder + ward 0.814899 0.8215
Num + Cat (> 1 & < 2100) fe anteriores + lumping (mediana) + word embed sobre funder y ward (dim = 2) 0.8155724 0.8204
Num + Cat (> 1 & < 2100) fe anteriores + lumping (mediana) sobre funder y ward + freq. abs. cat. 0.8157071 0.8226

Codificando las variables categóricas por sus frecuencias absolutas mejoramos los resultados obtenidos en el conjunto test, pasando de 0.8216 a 0.8226.

3.1.9 Variables scheme_name e installer

Hasta ahora, hemos trabajado con variables categóricas que tuviesen menos de 2100 valores únicos. Sin embargo, existen variables con un mayor número de categorías como scheme_name, installer, subvillage y wpt_name. Con el propósito de seguir un orden, empezamos añadiendo a nuestro dataset aquellas variables con menos de 10000 categorías, esto es, scheme_name e installer.

Sobre dichos campos aplicaremos exactamente los mismos pasos que los realizados hasta ahora con las variables categóricas (y con los que se han obtenido mejores resultados):

  1. Limpieza de strings (conversión a minúsculas + eliminación de espacios en blanco + eliminación de signos de puntuación)

  2. Aplicar lumping sobre la mediana en la proporción de apariciones

  3. Reemplazar las categorías por sus frecuencias absolutas

# Incluimos las variables scheme_name e installer en datcompleto
datcompleto[, scheme_name := datcat_mas_1000_mas_logicas_dr[, scheme_name]]
datcompleto[, installer   := datcat_mas_1000_mas_logicas_dr[, installer]]
#-- fe_scheme_name
#   Numero de categorías antes de aplicar clean_text
length(unique(datcompleto[, scheme_name]))
## [1] 2869
#-- fe_installer
#   Numero de categorías antes de aplicar clean_text
length(unique(datcompleto[, installer]))
## [1] 2411
cols <- c('installer', 'scheme_name')
datcompleto[ , paste0('fe_',cols) := lapply(.SD, clean_text), .SDcols = cols]
rm(cols)

#   Numero de categorías tras aplicar clean_text
length(unique(datcompleto[, fe_scheme_name]))
## [1] 2615
length(unique(datcompleto[, fe_installer]))
## [1] 2069

Como podemos observar en las salidas anteriores, el número de categorías se ha visto reducido en ambas variables. A continuación, aplicamos lumping sobre dichas columnas:

#-- fe_scheme_name
summary(c(prop.table(table(datcompleto[, fe_scheme_name]))))
##      Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
## 0.0000135 0.0000269 0.0000808 0.0003824 0.0002020 0.4748552
#-- Aplicamos lumping sobre la mediana (50 % de categorias con una proporcion menor a 8.1e-05, aprox.)
datcompleto[, fe_scheme_name := fct_lump_prop(datcompleto[,fe_scheme_name], 8.1e-05, other_level = "other")]
datcompleto$fe_scheme_name <- as.character(datcompleto$fe_scheme_name)

datcompleto[, scheme_name := NULL]
# Pasamos de 2615 a 1205
length(unique(datcompleto[, fe_scheme_name]))
## [1] 1205
#-- fe_installer
summary(c(prop.table(table(datcompleto[, fe_installer]))))
##      Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
## 1.347e-05 1.347e-05 1.347e-05 4.833e-04 8.081e-05 2.934e-01
#-- Aplicamos lumping sobre la mediana (50 % de categorias con una proporcion menor a 2e-05, aprox.)
datcompleto[, fe_installer := fct_lump_prop(datcompleto[,fe_installer], 2e-05, other_level = "other")]
datcompleto$fe_installer <- as.character(datcompleto$fe_installer)

datcompleto[, installer := NULL]
# Pasamos de 2069 a 1029 categorias
length(unique(datcompleto[, fe_installer]))
## [1] 1029

Finalmente, una vez reducido el número de categorías al 50 %, aproximadamente, realizamos la transformación por frecuencias absolutas:

#-- Creamos una funcion para evitar redundancia de codigo
impute_freq_abs <- function(cat_cols) {

  #   Antes de imputar
  freq_antes_fe <- apply(datcompleto[, ..cat_cols], 2, function(x) length(unique(x)))
  
  for (cat_col in cat_cols) {
    datcompleto[, paste0("fe_", cat_col) := as.numeric(.N), by = cat_col]
  }
  names(datcompleto) <- stri_replace_all_fixed(names(datcompleto), "fe_fe_", "fe_")
  
  for (cat_col in cat_cols) {
    datcompleto[, paste(cat_col) := NULL]
  }
  new_cat_cols <- paste0("fe_", stri_replace_all_fixed(cat_cols, "fe_", ""))
  
  #  Despues de imputar
  freq_despues_fe <- apply(datcompleto[, ..new_cat_cols], 2, function(x) length(unique(x)))
  
  comparacion_freq_dt <- data.table(rbind(freq_antes_fe[cat_cols], freq_despues_fe[new_cat_cols]))
  
  print(comparacion_freq_dt)
  return(datcompleto)
}
# Imputacion de las variables categoricas por sus frecuencias absolutas
cat_cols <- names(datcompleto[, which(sapply(datcompleto, is.character)), with = FALSE])

datcompleto <- impute_freq_abs(cat_cols)
##    fe_installer fe_scheme_name
## 1:         1029           1205
## 2:          164            126

Comprobamos como el número de categorías en scheme_name e installer ha disminuido a 126 y 164 valores únicos, respectivamente.

Una vez depuradas ambas columnas, lanzamos el nuevo modelo:

my_model_14 <- fit_random_forest(formula, train)

my_sub_14 <- make_predictions(my_model_15, test)
Train.accuracy Driven.Data
Num + Cat (> 1 & < 10000) fe anteriores + lumping (mediana) + freq. abs. cat. 0.8158081 0.8231

Aparentemente, incluir ambas variables obtiene un mejor resultado (de 0.8226 a 0.8231). No obstante, ¿Merece la pena incluir ambas variables? Es decir, ¿Mejoraría el modelo si descartamos una de ellas? Veamos:

Train.accuracy Driven.Data
Empleando tanto scheme_name como installer 0.8158081 0.8231
Empleando unicamente installer 0.8160774 0.8211
Empleando unicamente scheme_name 0.8161616 0.8239

Curiosamente, si descartamos la variable installer, el score del modelo mejora de 0.8231 a 0.8239. Por tanto, en lugar en emplear ambas columnas añadiremos únicamente scheme_name:

Train.accuracy Driven.Data
Mejor accuracy en el concurso
0.8176
Num + Cat (> 1 & < 2100) fe anteriores + fe_funder + fe_ward 0.8149832 0.8197
Num + Cat (> 1 & < 2100) fe anteriores + lumping sobre funder + ward (mediana) 0.8159764 0.8212
Num + Cat (> 1 & < 2100) fe anteriores + lumping sobre funder + ward (tercer cuartil) 0.8146633 0.8203
Num + Cat (> 1 & < 2100) fe anteriores + lumping sobre funder + ward (primer cuartil) 0.8159259 0.8196
Num + Cat (> 1 & < 2100) fe anteriores + lumping (mediana) + hashed sobre funder + ward 0.8158418 0.8214
Num + Cat (> 1 & < 2100) fe anteriores + lumping (mediana) + freq. abs. sobre funder + ward 0.8154882 0.8216
Num + Cat (> 1 & < 2100) fe anteriores + lumping (mediana) + freq. abs. y target enc. (media) sobre funder + ward 0.8154377 0.8164
Num + Cat (> 1 & < 2100) fe anteriores + lumping (mediana) + freq. abs. y target enc. (suma) sobre funder + ward 0.814899 0.8215
Num + Cat (> 1 & < 2100) fe anteriores + lumping (mediana) + word embed sobre funder y ward (dim = 2) 0.8155724 0.8204
Num + Cat (> 1 & < 2100) fe anteriores + lumping (mediana) sobre funder y ward + freq. abs. cat. 0.8157071 0.8226
Num + Cat (> 1 & < 10000) fe anteriores + lumping (mediana) + freq. abs. cat. (sin installer) 0.8161616 0.8239 
# Guardamos datcompleto
datcompleto[, fe_installer := NULL]
fwrite(datcompleto, 
       file = "./data/datcompleto_clean_text_lumping_mediana_freq_abs_cat_menos_installer.csv")

3.1.10 Añadiendo subvillage

Hasta el momento, hemos empleado variables con menos de 10.000 valores únicos ¿Y si damos nuevamente un salto y añadimos la variable subvillage? Por supuesto, realizando los mismos pasos que en el resto de variables categóricas (limpieza de textos, lumping y recategorización a frecuencias absolutas):

# Incluimos la variable subvillage en dat_completo
datcompleto[, subvillage := datcat_mas_1000_mas_logicas_dr[, subvillage]]
#-- fe_subvillage
#   Numero de categorías antes de aplicar clean_text
length(unique(datcompleto[, subvillage]))
## [1] 21426
cols <- c('subvillage')
datcompleto[ , paste0('fe_',cols) := lapply(.SD, clean_text), .SDcols = cols]
rm(cols)

#   Numero de categorías tras aplicar clean_text
length(unique(datcompleto[, fe_subvillage]))
## [1] 21295

Nuevamente, el número de valores únicos se ve reducido de 21.426 a 21.295. A continuación, aplicamos lumping e imputación por frecuencias absolutas:

#-- fe_subvillage
summary(c(prop.table(table(datcompleto[, fe_subvillage]))))
##      Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
## 1.347e-05 1.347e-05 2.694e-05 4.696e-05 4.040e-05 8.714e-03
#-- Aplicamos lumping sobre la mediana (50 % de categorias con una proporcion menor a 3e-05, aprox.)
datcompleto[, fe_subvillage := fct_lump_prop(datcompleto[,fe_subvillage], 3e-05, other_level = "other")]
datcompleto$fe_subvillage <- as.character(datcompleto$fe_subvillage)

datcompleto[, subvillage := NULL]
# Pasamos de 21295 a 7511
length(unique(datcompleto[, fe_subvillage]))
## [1] 7511
#-- Finalmente, realizamos la imputacion por frecuencias absolutas
cat_cols <- names(datcompleto[, which(sapply(datcompleto, is.character)), with = FALSE])

datcompleto <- impute_freq_abs(cat_cols)
##    fe_subvillage
## 1:          7511
## 2:           109

Como primera impresión, el número de categorías en la variable subvillage se ha visto reducido drásticamente (hasta 109 valores únicos), de forma que podemos estar perdiendo información relevante que permita ganar en puntuación, aunque sea en el orden de milésimas. No obstante, probamos a entrenar el nuevo modelo y realizar la submission correspondiente:

# Train
my_model_14 <- fit_random_forest(formula, train)

# Submission
my_sub_14   <- make_predictions(my_model_14, test)
Train.accuracy Driven.Data
Mejor accuracy hasta el momento 0.8161616 0.8239
Num + Cat (> 1 & < 40000) fe anteriores + lumping (mediana) + freq. abs. cat. (sin installer) 0.8165993 0.8219

Como podemos comprobar en la tabla anterior, el hecho de añadir la variable subvillage ha provocado que el modelo empeore, pasando de 0.8239 a 0.8219. Esto puede ser debido al hecho de haber reducido de forma significativa el número de valores únicos en subvillage, por lo que veamos cómo quedaría el modelo sin aplicar ninguna técnica de depuración sobre dicha columna:

Train.accuracy Driven.Data
Mejor accuracy hasta el momento 0.8161616 0.8239
Num + Cat (> 1 & < 40000) fe anteriores + lumping (mediana) + freq. abs. cat. (sin installer) 0.8165993 0.8219
Sin depurar subvillage (sin installer) 0.8162963 0.8220

Incluso sin depurar la variable, el score obtenido continua siendo menor al mejor modelo obtenido, por lo que descartamos añadir dicha columna.

3.1.11 Añadiendo wpt_name

Junto con subvillage, nos queda por añadir una variable con un elevado número de categorías: wpt_name, con 45.684 valores únicos. Del mismo modo que en el apartado anterior, aplicamos el mismo proceso de depuración sobre dicha variable.

# Eliminamos fe_subvillage e incluimos la variable wpt_name en dat_completo
datcompleto[, subvillage := NULL]
## Warning in `[.data.table`(datcompleto, , `:=`(subvillage, NULL)): Column
## 'subvillage' does not exist to remove
datcompleto[, wpt_name   := datcat_mas_1000_mas_logicas_dr[, wpt_name]]
#-- fe_wpt_name
#   Numero de categorías antes de aplicar clean_text
length(unique(datcompleto[, wpt_name]))
## [1] 45684
cols <- c('wpt_name')
datcompleto[ , paste0('fe_',cols) := lapply(.SD, clean_text), .SDcols = cols]
rm(cols)

#   Numero de categorías tras aplicar clean_text
length(unique(datcompleto[, fe_wpt_name]))
## [1] 45042

Nuevamente, el número de categorías se reduce tras aplicar la función clean_text. Por último, si aplicamos lumping e imputamos por frecuencias absolutas:

#-- Aplicamos lumping sobre la mediana (50 % de categorias con una proporcion menor a 2e-05, aprox.)
datcompleto[, fe_wpt_name := fct_lump_prop(datcompleto[,fe_wpt_name], 2e-05, other_level = "other")]
datcompleto$fe_wpt_name <- as.character(datcompleto$fe_wpt_name)

datcompleto[, wpt_name := NULL]
# Pasamos de 45042 a 5924
length(unique(datcompleto[, fe_wpt_name]))
## [1] 5924
#-- Finalmente, realizamos la imputacion por frecuencias absolutas
cat_cols <- names(datcompleto[, which(sapply(datcompleto, is.character)), with = FALSE])

datcompleto <- impute_freq_abs(cat_cols)
##    fe_wpt_name
## 1:        5924
## 2:          87

Cabe destacar la considerable disminución en el número de categorías, por lo que es posible que estemos perdiendo información relevante de cara al modelo, del mismo modo que sucedía con subvillage. No obstante, si entrenamos el nuevo modelo:

# Train
my_model_15 <- fit_random_forest(formula, train)

# Submission
my_sub_15   <- make_predictions(my_model_18, test)
Train.accuracy Driven.Data
Mejor accuracy hasta el momento 0.8161616 0.8239
Num + Cat (> 1 & < 60000) fe anteriores + lumping (mediana) + freq. abs. cat. (sin installer y subvillage) 0.8165993 0.8236

Sorprendentemente, pese al haber reducido drásticamente el número de categorías, el score obtenido no se aleja tanto del mejor modelo hasta el momento (0.8236 frente a 0.8239) ¿Y si mantenemos la variable wpt_name original? Es decir, sin depurar:

Train.accuracy Driven.Data
Mejor accuracy hasta el momento 0.8161616 0.8239
Num + Cat (> 1 & < 60000) fe anteriores + lumping (mediana) + freq. abs. cat. (sin installer y subvillage) 0.8165993 0.8236
Sin depurar wpt_name 0.8155892 0.8220

Por el contrario, si mantenemos la variable original, el hecho de disponer de un elevado número de categorías provoca una disminución en el score del modelo.

Por tanto, de las variables con un elevado número de categorías nos quedaremos únicamente con ward, funder y scheme_name:

Train.accuracy Driven.Data
Num + Cat (> 1 & < 10000) fe anteriores + lumping (mediana) + freq. abs. cat. (sin installer) 0.8161616 0.8239

3.1.12 Añadiendo variables lógicas

Como últimas variables, nos quedan por añadir permit y public meeting, ambas variables lógicas. Sin embargo, nos encontramos un problema: ambas columnas presentan valores missing:

# Numero NAs permit y public_meeting
sum(is.na(datcompleto[, permit]))
## [1] 3793
sum(is.na(datcompleto[, public_meeting]))
## [1] 4155

Por tanto, dado que existen valores missing, debemos imputarlos por medio de la función missRanger (empleando los mismos parámetros que los utilizados durante el concurso, es decir, k = 5 y número de iteraciones = 100).

# Imputacion de permit y public_meeting por missRanger
datcompleto_imp <- missRanger(datcompleto,
                            pmm.k = 5,
                            seed = 1234,
                            maxiter = 100, 
                            verbose = 0)
## Growing trees.. Progress: 88%. Estimated remaining time: 4 seconds.
## Growing trees.. Progress: 85%. Estimated remaining time: 5 seconds.
# Comprobamos que no existen valores missing
sum(is.na(datcompleto_imp))
## [1] 0

A continuación, podría resultar de gran utilidad incluir dos variables adicionales ( is_na_public_meeting e is_na_permit ), que indiquen si el valor original, previamente de ser imputado, era un valor missing o no (1,0):

# Creamos dos columnas adicionales que indiquen si el valor original era o no NA
datcompleto_imp[, is_na_public_meeting := ifelse(is.na(datcompleto[, public_meeting]), 1, 0)]
datcompleto_imp[, is_na_permit := ifelse(is.na(datcompleto[, permit]), 1, 0)]

Finalmente, lanzamos el nuevo modelo:

# Train
my_model_15 <- fit_random_forest(formula, train)

# Submission
my_sub_15   <- make_predictions(my_model_15, test)
Train.accuracy Driven.Data
Mejor accuracy hasta el momento 0.8161616 0.8239
Num + Cat (> 1 & < 10000) fe anteriores + lumping (mediana) + freq. abs. cat. (sin installer) + variables logicas 0.8168855 0.8251

Añadiendo ambas variables lógicas, el modelo ha mejorado de 0.8239 a 0.8251, superando incluso a la mejor puntuación obtenida en clase (0.8248).

No obstante, si hacemos “zoom” sobre el gráfico de importancia:

Observamos que las variables binarias is_na apenas tienen importancia en el modelo ¿Y si las eliminamos?

Train.accuracy Driven.Data
Mejor accuracy hasta el momento 0.8161616 0.8239
Num + Cat (> 1 & < 10000) fe anteriores + lumping (mediana) + freq. abs. cat. (sin installer) + variables logicas 0.8168855 0.8251
Eliminando is_na_public_meeting e is_na_permit 0.8171886 0.8235

Pese a la poca importancia que presentaban, eliminando ambas columnas ha provocado que el modelo empeore en cuanto a score se refiere. Por otro lado, tanto public_meeting como permit son muy similares en cuanto a importancia ¿Mejoraría el modelo si descartamos una de ellas?

Train.accuracy Driven.Data
Mejor accuracy hasta el momento 0.8161616 0.8239
Num + Cat (> 1 & < 10000) fe anteriores + lumping (mediana) + freq. abs. cat. (sin installer) + variables logicas 0.8168855 0.8251
Eliminando is_na_public_meeting e is_na_permit 0.8171886 0.8235
Empleando solo permit (mas su columna is_na) 0.8168519 0.8234
Empleando solo public_meeting (mas su columna is_na) 0.8172559 0.8236

Nuevamente, en ninguno de los casos el modelo se ve mejorado. Por tanto, mantenemos ambas variables lógicas más sus correspondientes columnas is_na.

Train.accuracy Driven.Data
Mejor accuracy en el concurso
0.8176
Num + Cat (> 1 & < 2100) fe anteriores + fe_funder + fe_ward 0.8149832 0.8197
Num + Cat (> 1 & < 2100) fe anteriores + lumping sobre funder + ward (mediana) 0.8159764 0.8212
Num + Cat (> 1 & < 2100) fe anteriores + lumping sobre funder + ward (tercer cuartil) 0.8146633 0.8203
Num + Cat (> 1 & < 2100) fe anteriores + lumping sobre funder + ward (primer cuartil) 0.8159259 0.8196
Num + Cat (> 1 & < 2100) fe anteriores + lumping (mediana) + hashed sobre funder + ward 0.8158418 0.8214
Num + Cat (> 1 & < 2100) fe anteriores + lumping (mediana) + freq. abs. sobre funder + ward 0.8154882 0.8216
Num + Cat (> 1 & < 2100) fe anteriores + lumping (mediana) + freq. abs. y target enc. (media) sobre funder + ward 0.8154377 0.8164
Num + Cat (> 1 & < 2100) fe anteriores + lumping (mediana) + freq. abs. y target enc. (suma) sobre funder + ward 0.814899 0.8215
Num + Cat (> 1 & < 2100) fe anteriores + lumping (mediana) + word embed sobre funder y ward (dim = 2) 0.8155724 0.8204
Num + Cat (> 1 & < 2100) fe anteriores + lumping (mediana) sobre funder y ward + freq. abs. cat. 0.8157071 0.8226
Num + Cat (> 1 & < 10000) fe anteriores + lumping (mediana) + freq. abs. cat. (sin installer) 0.8161616 0.8239
Num + Cat (> 1 & < 10000) fe anteriores + lumping (mediana) + freq. abs. cat. (sin installer) + variables logicas 0.8168855 0.8251

3.1.13 Extracción de variables sobre date_recorded

Nuevas variables extraidas de date_recorded

Tras añadir las variables lógicas public_meeting y permit, hemos utilizado todas las variables disponibles hasta el momento. No obstante, una de los campos que, en mi opinión, no llegó a aprovecharse lo suficiente, fue date_recorded, del cual solo mantuvimos los campos:

  1. fe_month: mes de date_recorded
  2. fe_dr_year_cyear_diff: diferencia entre el año de date_recorded y el año de construcción.

No obstante, de dicho campo podemos aprovechar mucho más, extrayendo diferentes variables como las que se proponen a continuación:

  1. year: durante el concurso en clase, acabamos descartando dicha variable al no aportar mejoría alguna al modelo. Sin embargo, con el dataset depurado hasta ahora ¿Mejorará si lo incluimos?

  2. Podemos incluso añadir nuevas variables a partir de date_recorded como el día, el día de la semana, el día del cuatrimestre, la semana, el cuatrimestre o incluso si el día en el que se registró la bomba de agua era fin de semana o no.

Comenzando añadiendo nuevamente el campo year al dataset depurado hasta el momento, haciendo uso de la librería lubridate:

#-- Icluimos year
year <- c(year(datcat_mas_1000_mas_logicas_dr[, date_recorded]))
datcompleto_imp$fe_dr_year <- year
Train.accuracy Driven.Data
Mejor modelo hasta el momento 0.8168855 0.8251
Añadiendo year (date_recorded) 0.8173064 0.8236

Nuevamente, nos encontramos con que la variable year no aporta mejoría alguna al modelo.

En contraposición, tal y como se mencionó anteriormente, podemos añadir nuevas variables extraídas a partir de date_recorded, gracias al paquete lubridate de tidyverse:

#-- Incluimos otros campos, mediante el paquete lubridate
datcompleto_imp[, fe_dr_day :=  day(date_recorded)]        # Dia
datcompleto_imp[, fe_dr_wday := wday(date_recorded)]       # Dia de la semana
datcompleto_imp[, fe_dr_qday := qday(date_recorded)]       # Dia del cuatrimestre
datcompleto_imp[, fe_dr_week := week(date_recorded)]       # Semana
datcompleto_imp[, fe_dr_quarter := quarter(date_recorded)] # Cuatrimestre

#-- Incluimos si es o no fin de semana
datcompleto_imp[, fe_is_weekend := ifelse(fe_dr_wday %in% c(1,7), 1, 0)]
Train.accuracy Driven.Data
Mejor modelo hasta el momento 0.8168855 0.8251
Añadiendo year (date_recorded) 0.8173064 0.8236
Incluyendo otros campos 0.8165657 0.8224

Pese a incluir tantos campos, el modelo empeora. Analicemos el gráfico de importancia:

Nota: para una mejor visualización, se han eliminado el resto de variables.

Podemos observar que variables como el cuatrimestre o si es o no fin de semana (fe_is_weekend y fe_dr_quarter) no presentan una importancia demasiado relevante, lo que podría estar provocando una disminución del score final ¿Mejoraría el modelo si descartamos ambas columnas?

Train.accuracy Driven.Data
Mejor modelo hasta el momento 0.8168855 0.8251
Añadiendo year (date_recorded) 0.8173064 0.8236
Incluyendo otros campos 0.8165657 0.8224
Solo dr_qday, dr_day, dr_week y dr_wday 0.8171212 0.8226

Descartando ambas variables, el modelo mejora ligeramente, aunque continua lejos del mejor score (0.8251). No obstante, podemos dar un salto más y eliminar las siguientes variables con menor importancia como el día de la semana ( fe_dr_wday ) y la semana ( fe_dr_week ):

Train.accuracy Driven.Data
Mejor modelo hasta el momento 0.8168855 0.8251
Añadiendo year (date_recorded) 0.8173064 0.8236
Incluyendo otros campos 0.8165657 0.8224
Solo dr_qday, dr_day, dr_week y dr_wday 0.8171212 0.8226
Solo dr_qday y dr_day 0.8175084 0.8242

Aunque bien es cierto que el modelo mejora de 0.8226 a 0.8242, continua estando lejos de 0.8251, incluso si solo mantuviéramos fe_dr_qday:

Train.accuracy Driven.Data
Mejor modelo hasta el momento 0.8168855 0.8251
Añadiendo year (date_recorded) 0.8173064 0.8236
Incluyendo otros campos 0.8165657 0.8224
Solo dr_qday, dr_day, dr_week y dr_wday 0.8171212 0.8226
Solo dr_qday y dr_day 0.8175084 0.8242
Solo dr_qday 0.8170202 0.8243

El score apenas mejora. Por tanto, descartamos añadir cualquiera de las variables anteriores al modelo actual.

Day_count

Hemos comprobado como, de todas las posibles variables que podían extraerse de date_recorded, las únicas que hemos mantenido han sido fe_dr_year_cyear_diff y fe_dr_month. No obstante, y en relación a date_recorded, podríamos crear una nueva variable, similar a lo que realizamos en clase con construction_year (restando a la variable un año de referencia: 2014).

Para empezar, obtenemos una fecha de referencia, la fecha más reciente de la columna date_recorded:

#-- Obtenemos una fecha de referencia (mas reciente)
date_recorded    <- c(datcat_mas_1000_mas_logicas_dr$date_recorded)
fecha_referencia <- max(ymd(date_recorded)) # Tenemos 2013-12-03
## Warning: All formats failed to parse. No formats found.
fecha_referencia
## [1] NA

Sabemos que la máxima fecha registrada corresponde con el 3 de diciembre del año 2013 ¿Y si creamos una nueva variable que almacenen la diferencia (en días) entre dicha fecha y los valores de date_recorded?

datcompleto_imp[fe_dr_day_count := as.numeric(fecha_referencia - ymd(date_recorded))]
Train.accuracy Driven.Data
Mejor modelo hasta el momento 0.8168855 0.8251
Añadiendo day_count 0.8168855 0.8247

Pese al añadir dicha variable, el modelo continua sin mejorar.

Conclusión Feature Engineering

Tras realizar todos los pasos anteriores, dado el mejor score obtenido, nos decantamos por:

  1. Selección de variables numéricas

  2. Selección de variables categóricas (con menos de 10000 categorías), a excepción de installer.

  3. Aplicar la función clean_text sobre las variables categoricas (conversión a minúsculas, eliminación de espacios en blanco y signos de puntuación).

  4. Aplicar lumping sobre las variables categóricas, concretamente sobre la mediana en proporción de aparición.

  5. Añadir las variables lógicas permit y public_meeting, imputando NAs por missRanger.

Mejor score obtenido en DrivenData: 0.8251.

3.2 Tuneo de modelos

A continuación, y en base al dataset modelado por los pasos anteriores, se propone el tuneo de tres modelos Machine Learning:

  1. Random Forest

  2. XGboost

  3. autoML (h2o)

3.2.1 Random Forest

Durante la etapa de Feature Engineering, empleamos la siguiente configuración:

  1. Número de árboles: 500

  2. mtry: por defecto (raíz cuadrada del número de variables, alrededor de 6)

El tuneo mediante Random Forest se divide en dos etapas:

3.2.1.1. Selección del número de árboles, empleando el valor mtry por defecto:

#-- Tuneo numero de arboles
n_trees <- c(400, 500, 600, 700, 800, 900)

prediction_errors <- c()
for (n_tree in n_trees) {
  rf  <- fit_random_forest(formula, train, num_trees = n_tree)
  prediction_error <- rf$prediction.error
  prediction_errors <- c(prediction_errors, 1 - prediction_error)
  
  submission <- make_predictions(rf, test)
  fwrite(submission, paste0("./submissions/tunning_models/random_forest/ntrees/random_forest_with_",n_tree,"_ntrees.csv"))
}
rm(n_tree); rm(prediction_error); rm(rf); rm(submission)

Aparentemente, podemos comprobar que 400-500 arboles es una buena opcion. Sin embargo, por cuestiones asociadas a la reproducibilidad de la semilla, el score del modelo aumenta de 0.8251 a 0.8253.

3.2.1.2. Selección del mejor valor mtry:

Dado que con 400-500 árboles se obtiene un score superior al obtenido en la fase de Feature Engineering, realizamos el tuneo del valor mtry sobre ambos:

Empleando 400 árboles:

Empleando 500 árboles:

En ambos casos, pese a aumentar o disminuir el valor mtry, el modelo no mejora en cuanto a score se refiere.

3.2.1.3. Tuneado de la profundidad del arbol (max_depth)

Por último, procedemos a tunear la profundidad máxima del árbol, cuyo valor por defecto es NULL, es decir, profundidad ilimitada:

#-- Veamos si tuneando la profundidad maxima del arbol mejora
max_depth_vector           <- c(200, 150, 80, 50)
prediction_errors <- c()
for (max_depth in max_depth_vector) {
  # Empleamos 500 arboles y mtry = 6
  rf  <- fit_random_forest(formula, train, num_trees = 500, mtry = 6, max_depth = max_depth)
  prediction_error <- rf$prediction.error
  prediction_errors <- c(prediction_errors, 1 - prediction_error)
  submission <- make_predictions(rf, test)
  fwrite(submission, paste0("./submissions/tunning_models/random_forest/max_depth/random_forest_with_",max_depth,"_max_depth.csv"))
}
rm(n_tree); rm(prediction_error); rm(rf); rm(submission)

Curiosamente, empleando una profundidad máxima de 50-80, el score del modelo se estabiliza en torno a 0.8253, mismo score que el obtenido con profudidad ilimitada. Por tanto, con Random Forest mantendremos la configuración original:

Driven.Data
Mejor score Feature Egineering 0.8251
Mejor score Random Forest 0.8253

3.2.2 XGboost

Continuando con XGboost, nos encontramos con los siguientes parámetros a tunear:

  1. nrounds o número de iteraciones.
  2. max_depth o profundidad máxima del árbol.
  3. eta o parámetro de regularización.
  4. colsample_bytree o número de variables a sortear por cada iteración.
  5. colsample_bylevel o número de variables a sortear por cada nuevo nivel.
  6. colsample_bynode o número de variables a sortear cada vez que se tenga que realizar una división en el árbol.
  7. subsample o porcentaje de “submuestras” a escoger durante el entrenamiento.

Nuevamente, para el tuneo del XGboost comenzamos tuneando el número de iteraciones o nrounds. No obstante, y dado que con un mtry bajo obtuvimos buenos resultados en Random Forest, de cara a XGboost realizamos un sorteo de columnas a través de la variable colsample_bytree, concretamente un porcentaje bajo (0.3), de forma que se sortee un pequeño subconjunto de variables (30 %) en lugar de todas las columnas (por defecto). Además, de cara a obtener un accuracy “aproximado” del modelo, empleamos mlogloss como métrica de evaluación (por defecto):

params = list(
  objective = "multi:softmax",
  num_class = 3,
  colsample_bytree = 0.3,  # Sorteamos el 30 % de las variables
  eval_metric = "mlogloss" # multiclass logloss
)

3.2.2.1. Tuneo del número de iteraciones

Comenzamos tuneando el número de iteraciones o nrounds:

#-- Modelo 1: parametrizando el numero de iteraciones
nrounds     <- c(400, 500, 600, 700)
lista_accuracy <- c()
for(nround in nrounds) {
  
  my_model <- fit_xgboost_model(params, xgb.train, xgb.train, nrounds = nround)
  xgb_pred <- make_predictions_xgboost(my_model, test)
  fwrite(xgb_pred, 
      file = paste0("./submissions/tunning_models/xgboost/rounds/xgboost_with_",nround,"_iters.csv")
  )

  # "accuracy" = 1 - mlogloss
  lista_accuracy <- c(lista_accuracy, 1 - tail(my_model$evaluation_log$val1_mlogloss, 1))
}

Como podemos comprobar en la salida anterior, tuneando el número de iteraciones obtenemos una puntuación máxima de 0.816 en el conjunto test, un score bajo en comparación con los obtenidos tanto en Feature Engineering como en Random Forest.

3.2.2.2. Tuneo de max_depth y eta:

A continuación, y utilizando el mejor número de iteraciones (500), procedemos a tunear otros parámetros fundamentales en el modelo, concretamente:

  1. Profundidad máxima del árbol o max_depth.

  2. Parámetro de regularización o eta.

search_grid <- expand.grid(colsample_bytree = c(0.3),
                           max_depth = c(20, 15, 10, 8, 6, 3),
                           eta = c(0.3, 0.4, 0.5)
)

En primera instancia, empleando una profundidad máxima de 15 en el árbol, así como un parametro de regularización en torno a 0.3, conseguimos mejorar el resultado, aunque no lo suficiente en comparación con el mejor score obtenido hasta el momento: 0.8253.

Sin embargo, de la gráfica anterior debemos detallar un aspecto fundamental: a medida que se reduce el parametro de regularización eta, el score del test aumenta. Por tanto, ¿Y si reducimos aún más dicho parámetro? Veamos los resultados, empleando tanto 500 iteraciones como una profundidad máxima de 15 (mejores valores hasta el momento):

Como podemos observar, reduciendo parámetro eta hasta 0.02 ha permitido mejorar el score del conjunto test a 0.8257, superior al obtenido por Random Forest. Sin embargo, debemos recordar un detalle relevante: a medida que disminuye el parámetro de regularización, el tiempo de cómputo aumenta, por lo que puede ser necesario añadir más iteraciones. Por tanto, probamos a continuación modelos XGboost variando no solo el parámetro eta (entre 0.01 y 0.04, donde se han obtenido los resultados más altos), sino además variando el número de iteraciones (entre 500 y 700):

Comprobamos como con 600 iteraciones y un parametro de regularización de 0.02 conseguimos mejorar aun más el modelo, de 0.8257 a 0.8260.

3.2.2.3 Tuneo de colsample_bytree

A continuación, y empleando 600 iteraciones, max_depth = 15 y eta = 0.02, procedemos a tunear el parámetro colsample_bytree, el cual hemos mantenido por defecto a 0.3. De este modo, comprobamos si sorteando un mayor o menor número de variables mejora el score del modelo:

Pese a aumentar o disminuir el sorteo de variables, el score del modelo no mejora, por lo que mantenemos la variable colsample_bytree a 0.3.

3.2.2.4 Tuneo de subsample

Por último, tunearemos un últimos parámetro: subsample, cuyo valor por defecto es 1, es decir, se emplean todas las observaciones durante en el entrenamiento. Veamos si mejora el modelo al variar el tamaño de la muestra de entrenamiento:

Nuevamente, ajustando el tamaño de la muestra de entrenamiento, el modelo continua sin mejorar.

3.2.2.5 Tuneo de colsample_bylevel y colsample_bynode

En último lugar, nos queda por tunear los parámetros colsample_bylevel y colsamle_bynode, los cuales permiten sortear las variables a utilizar por cada nivel y nodo del árbol, respectivamente. Dichos parámetros son acumulativos. Pongamos un ejemplo:

  1. colsample_bytree: 0.3. Sorteo un 30 % de las variables en la construcción de cada árbol.

  2. colsample_bylevel: 0.7. Del 30 % anterior, se sortean un 70 % de las variables por cada nuevo nivel de profundidad, es decir, total_variables * 0.3 * 0.7.

  3. colsample_bynode: 0.6. Del total_variables * 0.3 * 0.7 anterior, se sortea un 60 % de las variables cada vez que se tenga que realizar una división en el árbol.

Por ello, del 30 % de variables sorteadas en cada arbol, probamos a sortear un porcentaje de variables en cada nivel, tuneando colsample_bylevel a 0.9, 0.8 o 0.7:

Driven.Data
colsample_bytree = 0.3 0.8260
colsample_bytree = 0.3 & colsample_bylevel = 0.9 0.8265
colsample_bytree = 0.3 & colsample_bylevel = 0.8 0.8257
colsample_bytree = 0.3 & colsample_bylevel = 0.7 0.8260

Curiosamente, sorteando además el número de observaciones en cada nivel, el modelo mejora a 0.8265 ¿Y si damos un más y sorteamos también el número de observaciones en cada división del árbol?

Driven.Data
colsample_bytree = 0.3 & colsample_bylevel = 0.9 0.8265
colsample_bytree = 0.3 & colsample_bylevel = 0.9 & colsample_bynode = 0.9 0.8256
colsample_bytree = 0.3 & colsample_bylevel = 0.9 & colsample_bynode = 0.8 0.8244
colsample_bytree = 0.3 & colsample_bylevel = 0.9 & colsample_bynode = 0.7 0.8272

Incluso sorteando las variables a tan bajo nivel, ¡Sigue mejorando el modelo! ¿Y si reducimos aún más colsample_bynode?

Driven.Data
colsample_bytree = 0.3 & colsample_bylevel = 0.9 0.8265
colsample_bytree = 0.3 & colsample_bylevel = 0.9 & colsample_bynode = 0.9 0.8256
colsample_bytree = 0.3 & colsample_bylevel = 0.9 & colsample_bynode = 0.8 0.8244
colsample_bytree = 0.3 & colsample_bylevel = 0.9 & colsample_bynode = 0.7 0.8272
colsample_bytree = 0.3 & colsample_bylevel = 0.9 & colsample_bynode = 0.6 0.8261
colsample_bytree = 0.3 & colsample_bylevel = 0.9 & colsample_bynode = 0.5 0.8266

Aparentemente, 0.7 es una buena opción. Como conclusión, con XGboost empleamos la siguiente configuración:

  1. colsample_bytree: 0.3
  2. max_depth: 15
  3. eta: 0.02
  4. nrounds: 600
  5. subsample: 1
  6. colsample_bylevel: 0.9
  7. colsample_bynode: 0.7
Driven.Data
Mejor score Feature Egineering 0.8251
Mejor score Random Forest 0.8253
Mejor score XGboost 0.8272

3.2.1 XGboost con H2o

Tras tunear el modelo con el paquete XGboost original, ¿Que ocurriría si trasladamos los parámetros del mejor modelo a otra librería? Concretamente, H2o dispone de una función denominada h2o.xgboost, el cual permite tunear prácticamente los mismos parámetros que en el paquete anterior:

#-- Conversion train y test a objeto h2o 
train_h <- as.h2o(train)
test_h  <- as.h2o(test)

xgboost_model <- h2o.xgboost(x = pred, y = y,
                             training_frame = train_h, ntrees = 600, 
                             seed = 1234, distribution = "multinomial",
                             eta = 0.02, max_depth = 15, colsample_bytree = 0.3, 
                             colsample_bylevel = 0.9, colsample_bynode = 0.7)
Driven.Data
Mejor score XGboost (paquete xgboost) 0.8272
Mejor score XGboost (paquete H2o) 0.8255

Dado que se tratan se diferentes paquetes, los resultados del modelo son ligeramente diferentes. En este caso, el resultado con h2o es menor, por lo que nos decantamos por el paquete XGboost anterior.

3.2.2 Autoxgboost

Antes de finalizar con XGboost realizamos una última prueba, concretamente con un paquete denominado autoxgboost, desarrollado por Janek Thomas. Dicho paquete permite realizar, de un modo similar a h2o, el tuneo automático del mejor modelo a través de una optimización bayesiana:

#-- autoxgboost
bomb_task <- makeClassifTask(data = as.data.frame(train), target = "status_group")

#  max.rounds:  numero maximo de iteraciones por modelo (indicamos un maximo de 600, el mejor
#               parametro obtenido hasta el momento)
#  time.budget: tiempo maximo de tuneo (lo establecemos a media hora: 1800 segundos)
my_model  <- autoxgboost(bomb_task, max.nrounds = 600, time.budget = 1800)

submission <- make_predictions(model = my_model, test_data = as.data.frame(test))
Driven.Data
Mejor score Feature Egineering 0.8251
Mejor score Random Forest 0.8253
Mejor score XGboost 0.8272
autoxgboost 0.8206

Pese al tuneo realizado, continua siendo mejor el modelo obtenido en el apartado anterior.

3.3 AutoML con H2o

Finalmente, con los datos depurados en el apartado Feature Engineering elaboramos modelos automl empleando el paquete h2o:

library(h2o)  # AutoML
#-- Iniciamos el cluster h2o
h2o.init()

#-- Conversion train y test a objeto h2o 
train_h <- as.h2o(train)
test_h  <- as.h2o(test)

#-- Lanzamos el modelo
#   1. Probamos con 10 modelos
#   2. Excluiremos los algoritmos mas costosos
#      2.1 Deep Learning
#      2.2 Ensamblado
#   3. Extraemos los 10 mejores modelos
#   4. Indicamos un balanceo de clases en la validacion cruzada
aml <- h2o.automl(x = pred, y = y,
                 training_frame = train_h,
                 max_models = 5,
                 seed = 1234,
                 exclude_algos = c("DeepLearning", "StackedEnsemble"),
                 balance_classes = TRUE
)

#-- Realizamos las predicciones y guardamos submissions
model_id <- as.vector(aml@leaderboard$model_id)
for(i in 1:length(model_id)) {
  aml_aux <- h2o.getModel(aml@leaderboard[i, 1])
  prediccion    <- h2o.predict(aml_aux, newdata = test_h)
  prediccion_df <- as.data.frame(prediccion)
  
  prediccion_df    <- prediccion_df$predict
  submission       <- data.table(id = test$id, status_group = prediccion_df)
  path             <- paste0("./submissions/tunning_models/h2o/",model_id[i],".csv")
  fwrite(submission, path)
}

Tras elaborar los modelos automl, analicemos las submissions:

Driven.Data
Mejor score Feature Egineering 0.8251
Mejor score Random Forest 0.8253
Mejor score XGboost 0.8272
autoxgboost 0.8206
XGBoost_2_AutoML 0.8101
XGBoost_1_AutoML 0.8118
XGBoost_3_AutoML 0.8020
GBM_5_AutoML 0.8079
GBM_4_AutoML 0.8137

Analizando las salidas anterioes, ninguno de los modelos automl obtenidos mejora el score del apartado anterior.

Conclusiones

Finalmente, y en base a los resultados obtenidos a lo largo de la práctica, el mejor score obtenido ha sido de 0.8272:

4. Anexo: código final

A continuación, se incluye el código final con el mejor score:

suppressPackageStartupMessages({
  library(dplyr)                # Manipulacion de datos 
  library(data.table)           # Lectura y escritura de ficheros
  library(ranger)               # randomForest (+ rapido que caret)
  library(forcats)              # Tratamiento de variables categoricas
  library(tictoc)               # Calculo de tiempo de ejecucion
  library(missRanger)           # Imputacion de valores NA
  library(knitr)                # Generacion de informes (formateo de tablas)
  library(gmt)                  # Calculo de la distancia geografica
  library(stringi)              # Tratamiento de strings
  library(missRanger)           # Tratamiento de valores missing (mediante random forest)
  library(xgboost)              # XGboost
})

# -----------------------------
# ----- funciones propias -----
# -----------------------------

#-- Funcion para limpieza de textos
#   1. Conversion a minusculas
#   2. Eliminacion de espacios en blanco
#   3. Eliminacion de signos de puntuacion
clean_text <- function(text) {
  stri_trans_tolower(
    stri_replace_all_regex(
      text, 
      pattern = "[ +\\p{Punct}]", 
      replacement = ""
    )
  )
}

#-- Funcion para entrenar un modelo XGboost, en base a los parametros proporcionados
fit_xgboost_model <- function(params, train, val, nrounds, early_stopping_rounds = 20, show_log_error = TRUE, seed = 1234) {
  set.seed(seed)
  my_model <- xgb.train(
    data   = train,
    params = params,
    watchlist=list(val1=val),
    verbose = 1,
    nrounds= nrounds,
    early_stopping_rounds = early_stopping_rounds,
    nthread=4
  )
  return(my_model)
}

#-- Funcion para realizar la prediccion de un modelo XGboost
make_predictions_xgboost <- function(my_model, test) {
  xgb_pred <- predict(my_model,as.matrix(test),reshape=T)
  xgb_pred <- ifelse(xgb_pred == 0, "functional", 
                     ifelse(xgb_pred == 1, 
                            "functional needs repair", 
                            "non functional"
                            )
                     )
  
  xgb_pred <- data.table(id = test$id, status_group = xgb_pred)
  return(xgb_pred)
}

# -----------------------------
# -------- script final -------
# -----------------------------

#---------------------- Carga de ficheros train y test ---------------------
dattrainOr    <- fread(file = "./data/train_values.csv", data.table = FALSE )
dattrainLabOr <- fread(file = "./data/train_labels.csv", data.table = FALSE )
dattestOr     <- fread(file = "./data/test_values.csv", data.table  = FALSE )

#-------------------------- Variables categoricas ---------------------------
datcat_df <- dattrainOr %>% select(where(is.character))

# Mediante un bucle for... 
numlev_df <- data.frame(
  "vars" = names(datcat_df),
  "levels" = apply(datcat_df, 2, 
                   function(x) length(unique(x)))
  
)

# Eliminamos los nombres de fila
rownames(numlev_df) <- NULL

numlev_df %>% arrange(levels)

# Unimos dattrainOr con la variable objetivo
dattrainOrlab <- merge(
  dattrainOr, dattrainLabOr,
  by.x = c('id'), by.y = c('id'),
  sort = FALSE
)

#-- Eliminamos las variables no empleadas
dattrainOrlab$recorded_by    <- NULL; dattestOr$recorded_by    <- NULL
dattrainOrlab$payment_type   <- NULL; dattestOr$payment_type   <- NULL
dattrainOrlab$quantity_group <- NULL; dattestOr$quantity_group <- NULL
dattrainOrlab$installer      <- NULL; dattestOr$installer      <- NULL
dattrainOrlab$wpt_name       <- NULL; dattestOr$wpt_name       <- NULL
dattrainOrlab$subvillage     <- NULL; dattestOr$subvillage     <- NULL

vector_status_group <- dattrainOrlab$status_group

dattrainOrlab$status_group <- NULL
datcompleto <- as.data.table(rbind(dattrainOrlab, dattestOr))
#-- Guardamos el indice de fila donde comienza el conjunto "test",
#   concretamente la posicion 59401
fila_test <- which(datcompleto$id == 50785)

#----------------------------- Feature Engineering ------------------------------
#-- fe_cyear: 2014 - construction_year
datcompleto$fe_cyear     <- 2014 - datcompleto$construction_year

#-- fe_dist: geodist(latitude, longitude) al (0,0)
datcompleto$fe_dist <- geodist(datcompleto$latitude, datcompleto$longitude, 0, 0)

#-- fe_cant_agua: cantidad de agua / hab.
datcompleto$fe_cant_agua <- ifelse(datcompleto$population == 0,
                                         0,
                                         round(datcompleto$amount_tsh /
                                                 datcompleto$population, 3)
)

#-- month: mes date_recorded
datcompleto$fe_dr_month           <- month(datcompleto$date_recorded)

#-- fe_dr_year_cyear_diff: año date_recorded - construction_year
datcompleto$fe_dr_year_cyear_diff <- year(datcompleto$date_recorded) - datcompleto$construction_year

#-- Eliminamos date_recorded
datcompleto$date_recorded <- NULL

#-- Limpieza de variables categoricas mediante clean_text
cols <- c('funder', 'ward', 'scheme_name')
datcompleto[ , paste0('fe_',cols) := lapply(.SD, clean_text), .SDcols = cols]
rm(cols)

#-- lumping sobre la mediana de la proporcion de aparicion
#-  fe_funder
summary(c(prop.table(table(datcompleto[, fe_funder]))))
datcompleto[, fe_funder := fct_lump_prop(datcompleto[,fe_funder], 2e-05, other_level = "other")]
datcompleto$fe_funder <- as.character(datcompleto$fe_funder)

datcompleto[, funder := NULL]

#- fe_ward
summary(c(prop.table(table(datcompleto[, fe_ward]))))
datcompleto[, fe_ward := fct_lump_prop(datcompleto[,fe_ward], 4e-04, other_level = "other")]
datcompleto$fe_ward <- as.character(datcompleto$fe_ward)

datcompleto[, ward := NULL]

#- fe_scheme_name
summary(c(prop.table(table(datcompleto[, fe_scheme_name]))))
datcompleto[, fe_scheme_name := fct_lump_prop(datcompleto[,fe_scheme_name], 8.1e-05, other_level = "other")]
datcompleto$fe_scheme_name <- as.character(datcompleto$fe_scheme_name)

datcompleto[, scheme_name := NULL]

#-- Imputacion de las variables categoricas por sus frecuencias absolutas
cat_cols <- names(datcompleto[, which(sapply(datcompleto, is.character)), with = FALSE])

#-  Antes de imputar
freq_antes_fe <- apply(datcompleto[, ..cat_cols], 2, function(x) length(unique(x)))

for (cat_col in cat_cols) {
  datcompleto[, paste0("fe_", cat_col) := as.numeric(.N), by = cat_col]
}
names(datcompleto) <- stri_replace_all_fixed(names(datcompleto),
                                             "fe_fe_", "fe_")
#-- Eliminamos las variables originales
for (cat_col in cat_cols) {
  datcompleto[, paste(cat_col) := NULL]
}

new_cat_cols <- paste0("fe_", stri_replace_all_fixed(cat_cols, "fe_", ""))

#- Despues de imputar
freq_despues_fe <- apply(datcompleto[, ..new_cat_cols], 2, function(x) length(unique(x)))

cat("ANTES DE IMPUTAR")
freq_antes_fe
cat("------------------------------")
cat("DESPUES DE IMPUTAR")
freq_despues_fe

#-- Variables logicas (debemos imputarlas, dado que en el conjunto test existen missings)
sum(is.na(dattestOr$permit))
sum(is.na(dattestOr$public_meeting))

# Ordenamos las columnas para obtener la misma imputacion de las variables logicas con missRanger
# dado que al variar el orden de las columnas, la imputacion tambien varia
cols_orden <- c(names(datcompleto)[c(1:9)],  "construction_year", "fe_cyear", 
                "fe_dist", "fe_cant_agua", "fe_dr_year_cyear_diff", "fe_dr_month", 
                "public_meeting", "permit", names(datcompleto)[c(18:38)])

setcolorder(datcompleto, cols_orden)

#-- Imputamos por missRanger
datcompleto_imp <- missRanger(datcompleto,
                              pmm.k = 5,
                              seed = 1234,
                              maxiter = 100)

# Comprobamos que las variables logicas estan imputadas
sum(is.na(datcompleto_imp))

# Las convertimos a variables numericas
datcompleto_imp[, public_meeting := as.numeric(public_meeting)]
datcompleto_imp[, permit := as.numeric(permit)]

#-- Creamos dos columnas adicionales que indiquen si la variable logica era o no NA:
#   is_na_public_meeting
#   is_na_permit
datcompleto_imp[, is_na_public_meeting := ifelse(is.na(datcompleto[, public_meeting]), 1, 0)]
datcompleto_imp[, is_na_permit := ifelse(is.na(datcompleto[, permit]), 1, 0)]
#---------------------------- Fin Feature Engineering -----------------------------

#------------------------------------ Modelo --------------------------------------
formula <- as.formula("status_group~.")

train <- datcompleto_imp[c(1:fila_test-1),]

test  <- datcompleto_imp[c(fila_test:nrow(datcompleto_imp)),]

# Transformamos la variable objetivo en numerica
vector_status_group <- ifelse(vector_status_group == "functional", 0
                              , ifelse(vector_status_group == "functional needs repair", 1
                              , 2))
xgb.train <- xgb.DMatrix(data=as.matrix(train), label=vector_status_group)


params = list(
  objective = "multi:softmax",
  num_class = 3,
  colsample_bytree  = 0.3,
  colsample_bylevel = 0.9,
  colsample_bynode  = 0.7,
  max_depth         = 15,
  eta               = 0.02
)

tic()
my_model <- fit_xgboost_model(params, train = xgb.train, val = xgb.train, nrounds = 600)
toc()

# accuracy = 1 - mlogloss
accuracy <- 1 - tail(my_model$evaluation_log$val1_mlogloss, 1)
accuracy

xgb_pred <- make_predictions_xgboost(my_model, test)

# Guardamos la submission
fwrite(xgb_pred, file = "final_submission.csv")