Pump it Up - Driven Data - Machine Learning - R

Librerías

require(dplyr)
require(lubridate)
require(fastDummies)
require(ranger)
require(data.table)
require(tidyverse)
require(tidyr)
require(VIM)
require(DMwR)

Ingesta de datos iniciales

Se descarga el conjunto de datos y se guarda en las variables myTrain, myTest y myLabel.

#Por reproducibilidad, se descargan los archivos directamente de la página web de la competición. 
#Si se rompen los links, pueden ser extraídos de: https://www.drivendata.org/competitions/7/pump-it-up-data-mining-the-water-table/data/

myTrain <- fread("https://drivendata-prod.s3.amazonaws.com/data/7/public/4910797b-ee55-40a7-8668-10efd5c1b960.csv?X-Amz-Algorithm=AWS4-HMAC-SHA256&X-Amz-Credential=AKIARVBOBDCYVI2LMPSY%2F20210222%2Fus-east-1%2Fs3%2Faws4_request&X-Amz-Date=20210222T094545Z&X-Amz-Expires=86400&X-Amz-SignedHeaders=host&X-Amz-Signature=ea5642952302945bedaa9575f80accd501ddf3b21ec782d60b1f84947415418a")
myTest <- fread("https://drivendata-prod.s3.amazonaws.com/data/7/public/702ddfc5-68cd-4d1d-a0de-f5f566f76d91.csv?X-Amz-Algorithm=AWS4-HMAC-SHA256&X-Amz-Credential=AKIARVBOBDCYVI2LMPSY%2F20210222%2Fus-east-1%2Fs3%2Faws4_request&X-Amz-Date=20210222T094545Z&X-Amz-Expires=86400&X-Amz-SignedHeaders=host&X-Amz-Signature=0942443c72b4c39d886209712aa314f592b8a984fd786d0d6bd74cefeebc2ae2")
myLabel <- fread("https://drivendata-prod.s3.amazonaws.com/data/7/public/0bf8bc6e-30d0-4c50-956a-603fc693d966.csv?X-Amz-Algorithm=AWS4-HMAC-SHA256&X-Amz-Credential=AKIARVBOBDCYVI2LMPSY%2F20210222%2Fus-east-1%2Fs3%2Faws4_request&X-Amz-Date=20210222T094545Z&X-Amz-Expires=86400&X-Amz-SignedHeaders=host&X-Amz-Signature=4ae10b39b2f99fe2c5edc6650d65025cd22cfefc90392e31bb9f702006a54b4b")

Preprocesado 1. Tratamiento de valores no válidos y nulos.

# Anexamos la variable objetivo al conjunto de train mediante su "id".
trainIn <- left_join(myTrain, myLabel, by = "id")

# Etiquetamos pertenencia a grupo train o test
trainIn$Train <- 1
myTest$status_group <- "Conjunto myTest"
myTest$Train <- 0

# Sacamos las filas que tienen gps_height a 0 y longitude distinto de 0, puesto que si no tenemos las coordenadas correctas, calcularemos una elevación errónea.
my_coord <- trainIn %>% 
        filter(gps_height == 0 & longitude != 0) %>%
                select(id, gps_height, longitude, latitude) %>%
                        arrange(id)

# Dividimos el dataset en grupos de 2000 observaciones. Esto es debido a que la API que calcula elevaciones a través de longitudes y latitudes permite un número limitado de consultas a la hora y al día.
n <- 2000
nr <- nrow(my_coord)
splited_df <- split(my_coord, rep(1:ceiling(nr/n), each=n, length.out=nr))
list2env(setNames(splited_df, paste0("coord_", 1:10)), environment())

# Calculamos la elevación para cada uno de los 10 bloques de coordenadas (se necesita un username para ejecutarlo).
coord_1 <- coord_1 %>%
        mutate(gps_height = elevation(latitude = latitude, 
                                      longitude = longitude, 
                                      username = "username"))

coord_2 <- coord_2 %>%
        mutate(gps_height = elevation(latitude = latitude, 
                                      longitude = longitude, 
                                      username = "username"))

coord_3 <- coord_3 %>%
        mutate(gps_height = elevation(latitude = latitude, 
                                      longitude = longitude, 
                                      username = "username"))

coord_4 <- coord_4 %>%
        mutate(gps_height = elevation(latitude = latitude, 
                                      longitude = longitude, 
                                      username = "username"))

coord_5 <- coord_5 %>%
        mutate(gps_height = elevation(latitude = latitude, 
                                      longitude = longitude, 
                                      username = "username"))

coord_6 <- coord_6 %>%
        mutate(gps_height = elevation(latitude = latitude, 
                                      longitude = longitude, 
                                      username = "username"))

coord_7 <- coord_7 %>%
        mutate(gps_height = elevation(latitude = latitude, 
                                      longitude = longitude, 
                                      username = "username"))

coord_8 <- coord_8 %>%
        mutate(gps_height = elevation(latitude = latitude, 
                                      longitude = longitude, 
                                      username = "username"))

coord_9 <- coord_9 %>%
        mutate(gps_height = elevation(latitude = latitude, 
                                      longitude = longitude, 
                                      username = "username"))

coord_10 <- coord_10 %>%
        mutate(gps_height = elevation(latitude = latitude, 
                                      longitude = longitude, 
                                      username = "username"))

# Juntamos los splits de coordenadas en un unico dataframe.
coord_height <- do.call("rbind", list(coord_1, 
                                      coord_2, 
                                      coord_3, 
                                      coord_4, 
                                      coord_5, 
                                      coord_6, 
                                      coord_7, 
                                      coord_8, 
                                      coord_9, 
                                      coord_10))

# La función `elevation` devuelve un df de tres columnas, con lo que vamos a hacer un "subset" para extraer la columna que nos interesa y tirar el resto.
coord_height$gps_height <- coord_height$gps_height$elevation_geonames

my_coord_test <- myTest %>% 
                    filter(gps_height == 0 & longitude != 0) %>%
                                select(id, gps_height, longitude, latitude) %>%
                                        arrange(id)

# Dividimos el dataset en grupos de 2000 observaciones por lo que hemos comentado previamente.
n <- 2000
nr <- nrow(my_coord_test)
splited_df_test <- split(my_coord_test, rep(1:ceiling(nr/n), each=n, length.out=nr))
list2env(setNames(splited_df_test, paste0("coord_", 1:3)), environment())

# Calculamos la elevación para cada uno de los 3 bloques de coordenadas.
coord_1 <- coord_1 %>%
        mutate(gps_height = elevation(latitude = latitude, 
                                      longitude = longitude, 
                                      username = "username"))

coord_2 <- coord_2 %>%
        mutate(gps_height = elevation(latitude = latitude, 
                                      longitude = longitude, 
                                      username = "username"))

coord_3 <- coord_3 %>%
        mutate(gps_height = elevation(latitude = latitude, 
                                      longitude = longitude, 
                                      username = "username"))

# Juntamos los paquetitos de coordenadas en un gran dataframe
coord_height_test <- do.call("rbind", list(coord_1, 
                                           coord_2, 
                                           coord_3))

# La función `elevation` devuelve un df de tres columnas, con lo que vamos a hacer un subset para extraer la columna que nos interesa y tirar el resto.
coord_height_test$gps_height <- coord_height_test$gps_height$elevation_geonames

#Unimos los dos datasets (los tenemos bien identificados con el id).
coord_height_full <- rbind(coord_height_train, coord_height_test)

# Le ponemos una etiqueta a las elevaciones sacadas con la API.
coord_height_full$geoname <- TRUE

# Quitamos las columnas que no necesitamos de coord_height
coord_height_full <- coord_height_full[, c(1, 2, 5)]

# Juntamos train + test en un df
tt_df <- rbind(trainIn, myTest)

# Se añaden las elevaciones calculadas con la API mediante un left_join.
tt_df <- tt_df %>% left_join(coord_height_full, by = "id", 
                             suffix = c("", "_geoname"))

tt_df$gps_height <- ifelse((tt_df$gps_height == 0 & tt_df$longitude != 0), 
                                tt_df$gps_height_geoname,
                                tt_df$gps_height)

# Preparamos los archivos para seguir con su procesado.
myTrain <- tt_df %>%
                filter(Train == 1) %>%
                        select(-c(status_group, Train, geoname, gps_height_geoname))

myTest <- tt_df %>%
                filter(Train == 0) %>%
                        select(-c(status_group, Train, geoname, gps_height_geoname))

Preprocesado 2. Tratamiento de NA y unión de conjuntos.

Se combinan los conjuntos de myTrain y myLabel mediante la variable id en myTrain y luego se crea una nueva variable myTrainTest combinando myTrain con myTest, exluyendo la variable objetivo del primer conjunto.

# fread lee los "" como NA, así que hay que reemplazarlos manualmente.
myTrain$permit[is.na(myTrain$permit)] <- ""
myTrain$public_meeting[is.na(myTrain$public_meeting)] <- ""
myTest$permit[is.na(myTest$permit)] <- ""
myTest$public_meeting[is.na(myTest$public_meeting)] <- ""

# Anexamos la variable objetivo al conjunto de entrenamiento.
myTrain     <- left_join(myTrain, myLabel, by = "id")

# Guardamos la variable objetivo para más tarde.
myVarObj <- myTrain$status_group

# Unimos los conjuntos de entrenamiento y test para asegurarnos que aplicamos las mismas transformaciones a ambos.
myTrainTest <- myTrain %>%
                  select(-status_group) %>%
                     bind_rows(myTest)

Feature Engineering

Se crean nuevas variables a partir de otras, así como se modifican algunas existentes para intentar sacarles un mayor rendimiento.

# Extraemos valor de la variable `date_recorded` puesto que la variable con formato fecha no tiene potencial.
ref_date <- max(ymd(myTrainTest$date_recorded))
myTrainTest$day_count <- as.numeric(ref_date - ymd(myTrainTest$date_recorded))

# Modificamos las variables `region_code` y `district_code` a categóricas para que el modelo no le de un sentido que probablemente no tiene.
myTrainTest$region_code <- paste("RegCode", myTrainTest$region_code, sep = "_")
myTrainTest$district_code <- paste("DistCode", myTrainTest$district_code, sep = "_")

# Creamos las variables mes y año por si la época del año o el mismo año en que se registró tiene alguna relación con el estado de las bombas.
myTrainTest$month <- as.numeric(month(ymd(myTrainTest$date_recorded)))
myTrainTest$year  <- as.numeric(year(ymd(myTrainTest$date_recorded)))

# Convertimos la variable de water_quality a dummy.
wq_dum <- dummy_cols(myTrainTest, select_columns = 'water_quality')
wq_dum <- wq_dum[,44:51]
colnames(wq_dum) <- c("wq1", "wq2", "wq3", "wq4", "wq5", "wq6", "wq7", "wq8")

# Convertimos la variable de quantity a dummy.
quan_dum <- dummy_cols(myTrainTest, select_columns = 'quantity')
quan_dum <- quan_dum[,44:48]
colnames(quan_dum) <- c("q1", "q2", "q3", "q4", "q5")

# Anexamos las columnas recientemente creadas.
myTrainTest <- cbind(myTrainTest, wq_dum)
myTrainTest <- cbind(myTrainTest, quan_dum)
myTrainTest <- as.data.table(myTrainTest)

# Convertimos todas las columnas de tipo character a tipo factor.
myTrainTest <- myTrainTest %>% 
                  mutate_if(is.character, as.factor)

# Seleccionamos los nombres de las columnas que son factor.
cfactor <- myTrainTest %>%
              select_if(is.factor) %>%
                  names()

# Creamos un bucle que utiliza los nombres de las variables de tipo factor para anexar columnas de frecuencia para estas variables.
for (i in 1:length(cfactor)) {
        myTrainTest[, paste(cfactor[i], 'count', sep = '_') := as.numeric(.N), by = eval(cfactor[i])]
}

# Guardamos los nombres de las columnas que acabamos de crear en una variable.
cCount <- myTrainTest %>%
                select(contains("count")) %>%
                        names()

# Seteamos a -10000 aquellas observaciones cuya frecuencia es inferior a 5 para que el modelo no las considere relevantes.
for (i in 1:length(cCount)) {
        myTrainTest[get(cCount[i]) < 5 , c(cCount[i]) := -10000]
}

# Quitamos las columnas de tipo factor puesto que ya no nos sirven.
myTrainTest <- myTrainTest %>%
                        select_if(negate(is.factor))

Particionado

Particionado de la variable myTrainTest en train_num y test_num.

# Particionamos en train_num y asignamos le asignamos la variable objetivo.
train_num   <- myTrainTest[1:nrow(myTrain),]
train_num$myVarObj <- as.factor(myVarObj)

# Parcicionamos asimismo el test_num.
test_num    <- myTrainTest[(nrow(myTrain) + 1):nrow(myTrainTest),]

Modelado y predicción

Modelado de random forest con Ranger.

# Semilla para reproducibilidad.
set.seed(1)

# Modelo Random Forest con Ranger.
my_model <- ranger(
        myVarObj ~.,
        importance = "impurity",
        num.trees = 1000,
        mtry = 7,
        data = train_num
)

# Hacemos la predicción sobre test_num con el modelo.
my_predict <- predict(my_model, data = test_num)
pred_status_group <- my_predict$predictions

# Creamos el data.frame para la submission.
my_sub <- data.frame(id = test_num$id, status_group = pred_status_group)

# Escribimos el data.frame para la submission en un .csv.
fwrite(my_sub, "./submissions/Modelo_RFR_S20022021_F.csv")