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

Introducción

En este notebook vamos a implementar el Spatial resampling o remuestreo espacial y ver como varia con el anterior notebook presentado. En muchas partes es muy similar al notebook anterior, ya que tenemos el mismo algoritmo (Regresión logística) y el mismo objetivo (clasificar).

Lo que varía es la forma en que vamos a hacer el resampleo, es decir, en vez de hacer validación cruzada, los folds que vamos a crear van a tener en cuenta la variable espacial.

Recordemos que: 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
library(spatialsample)

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 con 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

Folds con spatial resampling

set.seed(123)
water_sp_folds <- spatial_clustering_cv(water_train, coords = c("lat_deg", "lon_deg"), v = 10)
water_sp_folds

## #  10-fold spatial cross-validation 
## # A tibble: 10 × 2
##    splits               id    
##    <list>               <chr> 
##  1 <split [38042/3052]> Fold01
##  2 <split [33486/7608]> Fold02
##  3 <split [34936/6158]> Fold03
##  4 <split [37804/3290]> Fold04
##  5 <split [35832/5262]> Fold05
##  6 <split [34298/6796]> Fold06
##  7 <split [38499/2595]> Fold07
##  8 <split [38431/2663]> Fold08
##  9 <split [38723/2371]> Fold09
## 10 <split [39795/1299]> Fold10

A modo de inspección visual

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).

Resampleo común

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)

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.

Notar aca que en vez de darle los water_folds (folds con validación cruzada) utilizo water_sp_folds (los folds que tienen en cuenta el resampleo espacial)

doParallel::registerDoParallel()
set.seed(7443)
logistic_rs <-
  fit_resamples(log_workflow,
    resamples = water_sp_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.856    10  0.0133 Preprocessor1_Model1
## 2 roc_auc  binary     0.891    10  0.0102 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.878 Preprocessor1_Model1
## 2 roc_auc  binary         0.906 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.0172    0.983     2 y           y         Preprocessor1_M…
##  2 train/test split 0.127     0.873     4 y           y         Preprocessor1_M…
##  3 train/test split 0.0112    0.989     8 y           y         Preprocessor1_M…
##  4 train/test split 0.0264    0.974    10 y           y         Preprocessor1_M…
##  5 train/test split 0.00414   0.996    11 y           y         Preprocessor1_M…
##  6 train/test split 0.109     0.891    12 y           y         Preprocessor1_M…
##  7 train/test split 0.160     0.840    13 y           y         Preprocessor1_M…
##  8 train/test split 0.106     0.894    14 y           n         Preprocessor1_M…
##  9 train/test split 0.124     0.876    17 y           y         Preprocessor1_M…
## 10 train/test split 0.00476   0.995    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 2557  785
##          y  887 9470

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.

Si observamos atentamente la matriz de confusión cambia muy poco a la matriz de confusión obtenida anteriormente. La decisión de incluir o no el resampleo espacial en un análisis de Machine learning, dependerá de los datos, de las métricas a optimizar, entre otras cosas