Práctica guiada - Clasificación con Regresión logística

Introducción

Vamos a realizar una regresión logística. La regresión logística se utiliza para clasificar, más allá que el nombre de “regresión” pueda parecer engañoso, lo usamos en ML para hacer clasificaciones ya sean binarias (2 clases) o multiclase (más de una clase). El ejemplo que vamos a utilizar está disponible en este link, PERO se optimiza otro algoritmo, Random Forest.

https://juliasilge.com/blog/water-sources/

Esto nos va a ser útil para poder comparar y ver qué modelo tiene menor error.

Variables del dataset

Los datos están disponibles en este link: https://github.com/rfordatascience/tidytuesday/blob/master/data/2021/2021-05-04/readme.md Se ofrece una descripción del dataset además de detallarse los pasos de limpieza que se realizan a la base de datos original. Para una descripción del problema visitar https://www.waterpointdata.org/

La variable a predecir es status_id y nos dice si hay o no disponibles agua, según diferentes variables.

Estas variables son:

row_id: id
lat_deg: latitud
lon_deg: longitud
report_date: fecha en que se reporto esa fuente de agua.
status_id: identifica si está disponible esa fuente de agua, el dia de la visita (si/no).
water_source_clean: describe la fuente de agua (variable categórica)
water_tech: describe si el sistema utilizado para transportar el agua desde la fuente al punto de recolección (incluye la posibilidad de q sea un distribuidor).
facility_type:
country_name: pais
install_year: año de la instalación
installer: instalador
pago: como se realiza el pago.
status: estado de la condición física / mecánica del punto de agua.

Ingresamos las librerias

library(tidymodels)
library(tidyverse)
library(themis) #esta libreria me permite trabajar con datos desbalanceados

Ingreso el dataset

water_raw <- read_csv("https://raw.githubusercontent.com/rfordatascience/tidytuesday/master/data/2021/2021-05-04/water.csv")

## Rows: 473293 Columns: 13

## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (9): report_date, status_id, water_source, water_tech, facility_type, co...
## dbl (4): row_id, lat_deg, lon_deg, install_year

## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

water_raw %>%
  glimpse()

## Rows: 473,293
## Columns: 13
## $ row_id        <dbl> 3957, 33512, 35125, 37760, 38118, 38501, 46357, 46535, 4…
## $ lat_deg       <dbl> 8.0731360, 7.3737842, 0.7734576, 0.7805757, 0.7792664, 0…
## $ lon_deg       <dbl> 38.61704, 40.50382, 34.92951, 34.96364, 34.97112, 34.122…
## $ report_date   <chr> "04/06/2017", "08/04/2020", "03/18/2015", "03/18/2015", …
## $ status_id     <chr> "y", "y", "y", "y", "y", "y", "y", "y", "y", "y", "y", "…
## $ water_source  <chr> NA, "Protected Spring", "Protected Shallow Well", "Boreh…
## $ water_tech    <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, "Hand Pump -…
## $ facility_type <chr> NA, "Improved", "Improved", "Improved", "Improved", "Imp…
## $ country_name  <chr> "Ethiopia", "Ethiopia", "Kenya", "Kenya", "Kenya", "Keny…
## $ install_year  <dbl> NA, 2019, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 2011, …
## $ installer     <chr> "Private-CRS", "WaterAid", NA, NA, NA, NA, NA, NA, NA, N…
## $ pay           <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
## $ status        <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, "Yes pump_pr…

Graficamos

water_raw %>%
  filter(
    country_name == "Sierra Leone",
    lat_deg > 0, lat_deg < 15, lon_deg < 0,
    status_id %in% c("y", "n")
  ) %>%
  ggplot(aes(lon_deg, lat_deg, color = status_id)) +
  geom_point(alpha = 0.1) +
  coord_fixed() +
  guides(color = guide_legend(override.aes = list(alpha = 1)))

Preprocesamiento de datos (limpieza)

Este paso se realiza sobre el dataset entero, ya que recodifico variables, elijo filas, etc. Es una etapa de limpieza previo al modelado, por eso lo hago sobre el dataset entero.

Este preprocesamiento es específico para este set de datos en particular.

water <- water_raw %>%
  filter(
    country_name == "Sierra Leone",
    lat_deg > 0, lat_deg < 15, lon_deg < 0,
    status_id %in% c("y", "n")
  ) %>%
  mutate(pay = case_when(
    str_detect(pay, "^No") ~ "no",
    str_detect(pay, "^Yes") ~ "yes",
    is.na(pay) ~ pay,
    TRUE ~ "it's complicated"
  )) %>%
  select(-country_name, -status, -report_date) %>%
  mutate_if(is.character, as.factor)

Luego de realizar la limpieza de datos, nos queda este dataset. La variable a predecir es status_id

Finalmente nos quedamos con 10 variables para modelar.

water %>%
  glimpse()

## Rows: 54,793
## Columns: 10
## $ row_id        <dbl> 252544, 36202, 34613, 36006, 36063, 36291, 34364, 34977,…
## $ lat_deg       <dbl> 9.445763, 8.369090, 8.483420, 8.378530, 8.483860, 8.3853…
## $ lon_deg       <dbl> -12.48694, -13.23060, -13.23500, -13.24050, -13.24310, -…
## $ status_id     <fct> n, y, y, y, n, n, y, y, y, y, y, y, y, n, y, n, y, y, y,…
## $ water_source  <fct> NA, NA, Surface Water (River/Stream/Lake/Pond/Dam), NA, …
## $ water_tech    <fct> Hand Pump, Rope and Bucket, Hydram, Hydram, Hand Pump - …
## $ facility_type <fct> NA, NA, No facilities, NA, NA, NA, Improved, Unimproved,…
## $ install_year  <dbl> 1997, 2017, 2011, 2014, 1994, 2017, NA, 2003, 2017, 2016…
## $ installer     <fct> NA, "Community", "Guma", "no", "private owner", "Communi…
## $ pay           <fct> NA, yes, yes, yes, yes, yes, NA, yes, yes, yes, yes, yes…

Veamos cuantos casos tenemos en SI y en NO.

water %>%
  count(status_id)

## # A tibble: 2 × 2
##   status_id     n
##   <fct>     <int>
## 1 n         13773
## 2 y         41020

Modelado

Divido los datos

#seteo la semilla
set.seed(123)
#divido los datos
water_split <- initial_split(water, strata = status_id) #estratificamos segun las clases 
water_train <- training(water_split)
water_test <- testing(water_split)

#folds de validación cruzada
set.seed(234)
water_folds <- vfold_cv(water_train, strata = status_id)
water_folds

## #  10-fold cross-validation using stratification 
## # A tibble: 10 × 2
##    splits               id    
##    <list>               <chr> 
##  1 <split [36984/4110]> Fold01
##  2 <split [36984/4110]> Fold02
##  3 <split [36984/4110]> Fold03
##  4 <split [36984/4110]> Fold04
##  5 <split [36984/4110]> Fold05
##  6 <split [36985/4109]> Fold06
##  7 <split [36985/4109]> Fold07
##  8 <split [36985/4109]> Fold08
##  9 <split [36985/4109]> Fold09
## 10 <split [36986/4108]> Fold10

water_train

## # A tibble: 41,094 × 10
##    row_id lat_deg lon_deg status_id water_source     water_tech    facility_type
##     <dbl>   <dbl>   <dbl> <fct>     <fct>            <fct>         <fct>        
##  1 252544    9.45   -12.5 n         <NA>             Hand Pump     <NA>         
##  2  36063    8.48   -13.2 n         <NA>             Hand Pump - … <NA>         
##  3  36291    8.39   -13.3 n         <NA>             <NA>          <NA>         
##  4  34838    8.28   -13.2 n         <NA>             Rope and Buc… <NA>         
##  5  36468    8.29   -13.2 n         <NA>             Rope and Buc… <NA>         
##  6  35939    8.42   -13.3 n         <NA>             Hydram        <NA>         
##  7  35149    8.38   -13.3 n         <NA>             <NA>          <NA>         
##  8  36779    8.34   -13.2 n         Borehole         Mechanized P… Improved     
##  9  35375    8.47   -13.2 n         Protected Shall… Hand Pump     Improved     
## 10  36789    8.34   -13.2 n         <NA>             <NA>          <NA>         
## # … with 41,084 more rows, and 3 more variables: install_year <dbl>,
## #   installer <fct>, pay <fct>

water_test

## # A tibble: 13,699 × 10
##    row_id lat_deg lon_deg status_id water_source     water_tech    facility_type
##     <dbl>   <dbl>   <dbl> <fct>     <fct>            <fct>         <fct>        
##  1  36202    8.37   -13.2 y         <NA>             Rope and Buc… <NA>         
##  2  36006    8.38   -13.2 y         <NA>             Hydram        <NA>         
##  3  34977    8.47   -13.2 y         Unprotected Spr… <NA>          Unimproved   
##  4  36340    8.31   -13.2 y         <NA>             Rope and Buc… <NA>         
##  5  34645    8.47   -13.2 y         Undefined Shall… <NA>          Unknown      
##  6  35418    8.41   -13.3 y         Borehole         <NA>          Improved     
##  7  35960    8.40   -13.3 y         Borehole         Mechanized P… Improved     
##  8  35913    8.47   -13.2 n         <NA>             Hand Pump     <NA>         
##  9  35430    8.47   -13.2 y         <NA>             Hand Pump - … <NA>         
## 10  36125    8.46   -13.2 y         Undefined Shall… <NA>          Unknown      
## # … with 13,689 more rows, and 3 more variables: install_year <dbl>,
## #   installer <fct>, pay <fct>

water_train %>%
  count(status_id)

## # A tibble: 2 × 2
##   status_id     n
##   <fct>     <int>
## 1 n         10329
## 2 y         30765

water_test %>%
  count(status_id)

## # A tibble: 2 × 2
##   status_id     n
##   <fct>     <int>
## 1 n          3444
## 2 y         10255

Como podemos ver las proporciones en las clases se conservan.

Receta

En esta etapa voy a PREPROCESAR los datos para el modelo de machine learning.

La variable a predecir es status_id

Funciones:

recipe(): aca ponemos la formula que vamos a modelar, es decir, el status_id según todas las variables presentes en el dataset, los datos que uso son los de train.
update_role(): asignar un nuevo rol a una nueva variable.
step_unknown(): asignar “unknown” a valores perdidos.
step_other(): para trabajar con las categorias que tienen poca ocurrencia.
step_impute_linear(): imputa un valor numérico segun un modelo lineal.
step_downsample(): para balancear las clases (disminuimos la clase mayoritaria).

library(themis) #esta libreria me permite trabajar con datos desbalanceados
logistic_recipe <-
  recipe(formula = status_id ~ ., data = water_train) %>%
  update_role(row_id, new_role = "id") %>%
  step_unknown(all_nominal_predictors()) %>%
  step_other(all_nominal_predictors(), threshold = 0.03) %>%
  step_impute_linear(install_year) %>%
  step_downsample(status_id)

library(themis) #esta libreria me permite trabajar con datos desbalanceados
logistic_recipe <-
  recipe(formula = status_id ~ lat_deg + lon_deg + row_id, data = water_train) %>%
  update_role(row_id, new_role = "id") %>%
  step_downsample(status_id)

Modelo propiamente dicho

Especificamos el modelo, en este caso, una regresión logística. El set_mode() puede omitirse ya que con regresion logistica siempre hacemos clasificación.

logistic_spec <- logistic_reg() %>%
  set_mode("classification") # esta linea puede omitirse

Workflow

Procedemos igual que las veces anteriores, agregando la receta y el modelo.

log_workflow <- workflow() %>%
  add_recipe(logistic_recipe) %>%
  add_model(logistic_spec)

Entrenamiento del algoritmo

Con la función fit_resamples() hacemos el ajuste de la función.

doParallel::registerDoParallel()
set.seed(7443)
logistic_rs <-
  fit_resamples(log_workflow,
    resamples = water_folds,
    control = control_resamples(save_pred = TRUE)
  )

Performance

Mediante la función collect_metrics() vamos a ver la performance del modelo, en TRAIN.

logistic_rs %>%
  collect_metrics()

## # A tibble: 2 × 6
##   .metric  .estimator  mean     n std_err .config             
##   <chr>    <chr>      <dbl> <int>   <dbl> <fct>               
## 1 accuracy binary     0.536    10 0.00165 Preprocessor1_Model1
## 2 roc_auc  binary     0.523    10 0.00272 Preprocessor1_Model1

collect_predictions(logistic_rs) %>%
  group_by(id) %>%
  roc_curve(status_id, .pred_n) %>%
  autoplot()

Predicción

En esta etapa con la función last_fit elijo el mejor modelo y lo utilizo para predecir en el set de datos de TEST. * last_fit(): toma como argumento el workflow q cree anteriormente y el split q tiene los datos de TRAIN y de TEST.

final_fitted <- last_fit(log_workflow, water_split)
final_fitted #imprimo

## # Resampling results
## # Manual resampling 
## # A tibble: 1 × 6
##   splits                id               .metrics  .notes .predictions .workflow
##   <list>                <chr>            <list>    <list> <list>       <list>   
## 1 <split [41094/13699]> train/test split <tibble … <tibb… <tibble [13… <workflo…

Métricas

Voy a ver las métricas del modelo con collect_metrics().

final_fitted %>%
  collect_metrics()

## # A tibble: 2 × 4
##   .metric  .estimator .estimate .config             
##   <chr>    <chr>          <dbl> <fct>               
## 1 accuracy binary         0.541 Preprocessor1_Model1
## 2 roc_auc  binary         0.525 Preprocessor1_Model1

Tengo un accuracy del modelo, igual a 87.4 %, para una clasificación binaria, es un valor aceptable. El área bajo la curva es de 0.90, y es un valor muy bueno también.

Predicciones del modelo

Las dos columnas que interesa mirar son:

.pred_class: es la clase predicha por el modelo
status_id: es la clase verdadera.

final_fitted %>%
  collect_predictions()

## # A tibble: 13,699 × 7
##    id               .pred_n .pred_y  .row .pred_class status_id .config         
##    <chr>              <dbl>   <dbl> <int> <fct>       <fct>     <fct>           
##  1 train/test split   0.470   0.530     2 y           y         Preprocessor1_M…
##  2 train/test split   0.470   0.530     4 y           y         Preprocessor1_M…
##  3 train/test split   0.473   0.527     8 y           y         Preprocessor1_M…
##  4 train/test split   0.470   0.530    10 y           y         Preprocessor1_M…
##  5 train/test split   0.473   0.527    11 y           y         Preprocessor1_M…
##  6 train/test split   0.471   0.529    12 y           y         Preprocessor1_M…
##  7 train/test split   0.470   0.530    13 y           y         Preprocessor1_M…
##  8 train/test split   0.473   0.527    14 y           n         Preprocessor1_M…
##  9 train/test split   0.473   0.527    17 y           y         Preprocessor1_M…
## 10 train/test split   0.473   0.527    28 y           y         Preprocessor1_M…
## # … with 13,689 more rows

Matriz de confusión

Para ello, vamos a utilizar la función conf_mat().

La matriz de confusión nos permite ver de manera resumida en cuantos casos el modelo “acertó” en la clase (realizando una predicción correcta), y que en que casos, no (predicción incorrecta).

final_fitted %>%
  collect_predictions() %>%
  conf_mat(status_id, .pred_class)

##           Truth
## Prediction    n    y
##          n 1824 4671
##          y 1620 5584

Como podemos ver los casos clasificados correctamente estan en la diagonal principal de la matriz de confusión, quedando por fuera de la misma, los FP y los FN.

Mas información sobre la matriz de confusión en las slides teóricas de la clase 1.

Actividad Práctica

El modelo anterior de regresión logística utiliza todas las variables del dataset para realizar una clasificación. Esto lo vemos en este chunk de código. En R, eso lo indicamos escribiendo

status_id ~ .

Esto nos indica que estamos utilizando todas las variables para modelar la regresión logística.

Consigna

En vez de utilizar todas las variables, elija 2 variables (lat y long) y realice los mismos cálculos de regresión logística. La pregunta a contestar es ¿obtengo los mismos resultados utilizando solo los datos espaciales?

Para ello reemplazar con el siguiente chunk de R que corresponde a la forma en que se preprocesan los datos.

library(themis) 
logistic_recipe <-
  recipe(formula = status_id ~ lat_deg + lon_deg + row_id, data = water_train) %>%
  update_role(row_id, new_role = "id") %>%
  step_downsample(status_id)