Un primer modelo de Machine Learning

Dataset

El dataset corresponde a Boston, que se encuentra disponible en Kaggle; una plataforma de competencia de machine learning. Más info en: https://www.cs.toronto.edu/~delve/data/boston/bostonDetail.html

Cargo las librerias

#si no tengo instalado la libreria, la instalo
#install.packages(tidymodels)
library(tidymodels)
library(tidyverse)
library(magrittr)
library(corrr)
library(MASS) #el dataset se encuentra en esta librería 
library(corrplot)

Ingreso los datos

#ingreso el dataset que esta en la libreria 
data(Boston)

Con la función head() de rbase, nos permite ver las primeras filas del dataset

head(Boston) #veamos el encabezado del dataset
##      crim zn indus chas   nox    rm  age    dis rad tax ptratio  black lstat
## 1 0.00632 18  2.31    0 0.538 6.575 65.2 4.0900   1 296    15.3 396.90  4.98
## 2 0.02731  0  7.07    0 0.469 6.421 78.9 4.9671   2 242    17.8 396.90  9.14
## 3 0.02729  0  7.07    0 0.469 7.185 61.1 4.9671   2 242    17.8 392.83  4.03
## 4 0.03237  0  2.18    0 0.458 6.998 45.8 6.0622   3 222    18.7 394.63  2.94
## 5 0.06905  0  2.18    0 0.458 7.147 54.2 6.0622   3 222    18.7 396.90  5.33
## 6 0.02985  0  2.18    0 0.458 6.430 58.7 6.0622   3 222    18.7 394.12  5.21
##   medv
## 1 24.0
## 2 21.6
## 3 34.7
## 4 33.4
## 5 36.2
## 6 28.7

Con la función glimpse() de dplyr podemos ver una breve descripción de las variables del dataset, ver que tipo de dato forman parte del dataset, etc.

glimpse(Boston)
## Rows: 506
## Columns: 14
## $ crim    <dbl> 0.00632, 0.02731, 0.02729, 0.03237, 0.06905, 0.02985, 0.08829,…
## $ zn      <dbl> 18.0, 0.0, 0.0, 0.0, 0.0, 0.0, 12.5, 12.5, 12.5, 12.5, 12.5, 1…
## $ indus   <dbl> 2.31, 7.07, 7.07, 2.18, 2.18, 2.18, 7.87, 7.87, 7.87, 7.87, 7.…
## $ chas    <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,…
## $ nox     <dbl> 0.538, 0.469, 0.469, 0.458, 0.458, 0.458, 0.524, 0.524, 0.524,…
## $ rm      <dbl> 6.575, 6.421, 7.185, 6.998, 7.147, 6.430, 6.012, 6.172, 5.631,…
## $ age     <dbl> 65.2, 78.9, 61.1, 45.8, 54.2, 58.7, 66.6, 96.1, 100.0, 85.9, 9…
## $ dis     <dbl> 4.0900, 4.9671, 4.9671, 6.0622, 6.0622, 6.0622, 5.5605, 5.9505…
## $ rad     <int> 1, 2, 2, 3, 3, 3, 5, 5, 5, 5, 5, 5, 5, 4, 4, 4, 4, 4, 4, 4, 4,…
## $ tax     <dbl> 296, 242, 242, 222, 222, 222, 311, 311, 311, 311, 311, 311, 31…
## $ ptratio <dbl> 15.3, 17.8, 17.8, 18.7, 18.7, 18.7, 15.2, 15.2, 15.2, 15.2, 15…
## $ black   <dbl> 396.90, 396.90, 392.83, 394.63, 396.90, 394.12, 395.60, 396.90…
## $ lstat   <dbl> 4.98, 9.14, 4.03, 2.94, 5.33, 5.21, 12.43, 19.15, 29.93, 17.10…
## $ medv    <dbl> 24.0, 21.6, 34.7, 33.4, 36.2, 28.7, 22.9, 27.1, 16.5, 18.9, 15…

Como podemos ver, todas las variables son cuantitativas continuas, lo que nos facilita no tener que codificarlas (en el caso que sean categóricas). Las variables (columnas) en total son 14; de las cuales una es la variable que nos interesa predecir, en este caso, medv. Las otras variables van a ser las variables para nuestro modelo.

Breve descripción de las variables

(https://stat.ethz.ch/R-manual/R-devel/library/MASS/html/Boston.html)

crim: tasa de criminalidad per cápita por ciudad.

zn: proporción de terreno residencial dividido en zonas para lotes de más de 25,000 pies cuadrados.

indus: proporción de acres comerciales no minoristas por ciudad.

chas: Variable ficticia de Charles River (= 1 si el tramo limita con el río; 0 en caso contrario).

nox: concentración de óxidos de nitrógeno (partes por 10 millones).

rm: número medio de habitaciones por vivienda.

age: Proporción de unidades ocupadas por sus propietarios construidas antes de 1940.

dis: media ponderada de las distancias a cinco centros de empleo de Boston.

rad:índice de accesibilidad a carreteras radiales.

tax: Tasa de impuesto a la propiedad de valor total por  $ 10,000.

ptratio: Proporción alumno-profesor por ciudad.

black: 1000 (Bk - 0.63) ^ 2 donde Bk es la proporción de negros por ciudad.

lstat: estatus más bajo de la población (porcentaje).

medv: valor medio de las viviendas ocupadas por sus propietarios en  $ 1000

Pequeño análisis exploratorio de toda la base de datos

Esta etapa nos ayuda a visualizar distribuciones de probabilidad de las variables, si hay algun outlier y demás.

cor(Boston)
##                crim          zn       indus         chas         nox
## crim     1.00000000 -0.20046922  0.40658341 -0.055891582  0.42097171
## zn      -0.20046922  1.00000000 -0.53382819 -0.042696719 -0.51660371
## indus    0.40658341 -0.53382819  1.00000000  0.062938027  0.76365145
## chas    -0.05589158 -0.04269672  0.06293803  1.000000000  0.09120281
## nox      0.42097171 -0.51660371  0.76365145  0.091202807  1.00000000
## rm      -0.21924670  0.31199059 -0.39167585  0.091251225 -0.30218819
## age      0.35273425 -0.56953734  0.64477851  0.086517774  0.73147010
## dis     -0.37967009  0.66440822 -0.70802699 -0.099175780 -0.76923011
## rad      0.62550515 -0.31194783  0.59512927 -0.007368241  0.61144056
## tax      0.58276431 -0.31456332  0.72076018 -0.035586518  0.66802320
## ptratio  0.28994558 -0.39167855  0.38324756 -0.121515174  0.18893268
## black   -0.38506394  0.17552032 -0.35697654  0.048788485 -0.38005064
## lstat    0.45562148 -0.41299457  0.60379972 -0.053929298  0.59087892
## medv    -0.38830461  0.36044534 -0.48372516  0.175260177 -0.42732077
##                  rm         age         dis          rad         tax    ptratio
## crim    -0.21924670  0.35273425 -0.37967009  0.625505145  0.58276431  0.2899456
## zn       0.31199059 -0.56953734  0.66440822 -0.311947826 -0.31456332 -0.3916785
## indus   -0.39167585  0.64477851 -0.70802699  0.595129275  0.72076018  0.3832476
## chas     0.09125123  0.08651777 -0.09917578 -0.007368241 -0.03558652 -0.1215152
## nox     -0.30218819  0.73147010 -0.76923011  0.611440563  0.66802320  0.1889327
## rm       1.00000000 -0.24026493  0.20524621 -0.209846668 -0.29204783 -0.3555015
## age     -0.24026493  1.00000000 -0.74788054  0.456022452  0.50645559  0.2615150
## dis      0.20524621 -0.74788054  1.00000000 -0.494587930 -0.53443158 -0.2324705
## rad     -0.20984667  0.45602245 -0.49458793  1.000000000  0.91022819  0.4647412
## tax     -0.29204783  0.50645559 -0.53443158  0.910228189  1.00000000  0.4608530
## ptratio -0.35550149  0.26151501 -0.23247054  0.464741179  0.46085304  1.0000000
## black    0.12806864 -0.27353398  0.29151167 -0.444412816 -0.44180801 -0.1773833
## lstat   -0.61380827  0.60233853 -0.49699583  0.488676335  0.54399341  0.3740443
## medv     0.69535995 -0.37695457  0.24992873 -0.381626231 -0.46853593 -0.5077867
##               black      lstat       medv
## crim    -0.38506394  0.4556215 -0.3883046
## zn       0.17552032 -0.4129946  0.3604453
## indus   -0.35697654  0.6037997 -0.4837252
## chas     0.04878848 -0.0539293  0.1752602
## nox     -0.38005064  0.5908789 -0.4273208
## rm       0.12806864 -0.6138083  0.6953599
## age     -0.27353398  0.6023385 -0.3769546
## dis      0.29151167 -0.4969958  0.2499287
## rad     -0.44441282  0.4886763 -0.3816262
## tax     -0.44180801  0.5439934 -0.4685359
## ptratio -0.17738330  0.3740443 -0.5077867
## black    1.00000000 -0.3660869  0.3334608
## lstat   -0.36608690  1.0000000 -0.7376627
## medv     0.33346082 -0.7376627  1.0000000

Vamos a dividir los datos

En esta instancia se deben dividir los datos en train y test.

set.seed(123)
p_split <- Boston %>%
  initial_split(prop=0.75) # divido en 75%
p_train <- training(p_split)
p_test <- testing(p_split)
p_split
## <Analysis/Assess/Total>
## <379/127/506>

Esta división es arbitraria, podemos hacer divisiones de 75/25, 80/20, 90/10, etc.

#Ejecutamos este comando para corrobar que la data de TRAIN haya sido randomizada al hacer el sampling. 
p_train %>%
  head()
##         crim zn indus chas   nox    rm   age    dis rad tax ptratio  black
## 415 45.74610  0 18.10    0 0.693 4.519 100.0 1.6582  24 666    20.2  88.27
## 463  6.65492  0 18.10    0 0.713 6.317  83.0 2.7344  24 666    20.2 396.90
## 179  0.06642  0  4.05    0 0.510 6.860  74.4 2.9153   5 296    16.6 391.27
## 14   0.62976  0  8.14    0 0.538 5.949  61.8 4.7075   4 307    21.0 396.90
## 195  0.01439 60  2.93    0 0.401 6.604  18.8 6.2196   1 265    15.6 376.70
## 426 15.86030  0 18.10    0 0.679 5.896  95.4 1.9096  24 666    20.2   7.68
##     lstat medv
## 415 36.98  7.0
## 463 13.99 19.5
## 179  6.92 29.9
## 14   8.26 20.4
## 195  4.38 29.1
## 426 24.39  8.3

Lo mismo para TEST, vemos si han sido tomados al azar los datos que forman parte de TEST.

p_test %>%
  head()
##       crim   zn indus chas   nox    rm   age    dis rad tax ptratio  black
## 1  0.00632 18.0  2.31    0 0.538 6.575  65.2 4.0900   1 296    15.3 396.90
## 6  0.02985  0.0  2.18    0 0.458 6.430  58.7 6.0622   3 222    18.7 394.12
## 8  0.14455 12.5  7.87    0 0.524 6.172  96.1 5.9505   5 311    15.2 396.90
## 9  0.21124 12.5  7.87    0 0.524 5.631 100.0 6.0821   5 311    15.2 386.63
## 15 0.63796  0.0  8.14    0 0.538 6.096  84.5 4.4619   4 307    21.0 380.02
## 17 1.05393  0.0  8.14    0 0.538 5.935  29.3 4.4986   4 307    21.0 386.85
##    lstat medv
## 1   4.98 24.0
## 6   5.21 28.7
## 8  19.15 27.1
## 9  29.93 16.5
## 15 10.26 18.2
## 17  6.58 23.1

Receta (Preprocesamiento de los datos)

Con este nombre agrupamos todos los pre-procesamientos que vamos a realizar al dataset de train, con lo cual haremos el modelado. Estas mismas operaciones se realizan al final en el dataset de test.

En este paso lo que vamos a hacer es eliminar las correlaciones, centrar, escalar los datos.

recipe_dt <- p_train %>%
  recipe(medv~.) %>%
  step_corr(all_predictors()) %>% #elimino las correlaciones
  step_center(all_predictors(), -all_outcomes()) %>% #centrado
  step_scale(all_predictors(), -all_outcomes()) %>% #escalado
  prep()

Vemos la receta que hemos creado

recipe_dt
## Recipe
## 
## Inputs:
## 
##       role #variables
##    outcome          1
##  predictor         13
## 
## Training data contained 379 data points and no missing data.
## 
## Operations:
## 
## Correlation filter removed tax [trained]
## Centering for crim, zn, indus, chas, nox, rm, age, dis, rad, ... [trained]
## Scaling for crim, zn, indus, chas, nox, rm, age, dis, rad, ... [trained]

Modelo

En esta instancia vamos a elegir el modelo para el entrenamiento, tambien seleccionamos de que libreria vamos a tomar las funciones de R y también el tipo de tarea. En este caso vamos a hacer, un árbol de decisión, de la libreria rpart y una tarea de regresión.

set.seed(123)
tree_spec <- decision_tree() %>% #arboles de decisión
  set_engine("rpart") %>% #librería rpart
  set_mode("regression") #regresión
tree_spec
## Decision Tree Model Specification (regression)
## 
## Computational engine: rpart

Inicio el workflow

En workflow nos permite trabajar de manera ordenada en MACHINE LEARNING, y es una buena práctica utilizarlo siempre que tengamos más de un paso de pre-procesamiento.

Luego imprimo el workflow

tree_wf <- workflow() %>%
  add_recipe(recipe_dt) %>% #agrego la receta
  add_model(tree_spec) #agrego el modelo
# Imprimo el workflow
tree_wf
## ══ Workflow ════════════════════════════════════════════════════════════════════
## Preprocessor: Recipe
## Model: decision_tree()
## 
## ── Preprocessor ────────────────────────────────────────────────────────────────
## 3 Recipe Steps
## 
## • step_corr()
## • step_center()
## • step_scale()
## 
## ── Model ───────────────────────────────────────────────────────────────────────
## Decision Tree Model Specification (regression)
## 
## Computational engine: rpart

Error de estimacion del modelo en TRAIN

Esta es la primer etapa de nuestro modelo, vamos entrenarlo en esta etapa, antes SOLO LO DEFINIMOS NO HEMOS CALCULADO NADA. Con la función fit(), se realiza el entrenamiento en este caso.

Depende de como decidamos modelar los datos, esta función fit(), puede cambiar por otra, por ejemplo, fit_resamples(), si hacemos validación cruzada.

set.seed(123) 

tree_fit <- tree_spec %>%
  fit(medv ~ ., data = p_train) 

#imprimo el modelo
tree_fit
## parsnip model object
## 
## Fit time:  13ms 
## n= 379 
## 
## node), split, n, deviance, yval
##       * denotes terminal node
## 
##  1) root 379 32157.2400 22.51214  
##    2) rm< 6.941 326 13439.7400 20.02577  
##      4) lstat>=14.805 120  2071.2040 14.27167  
##        8) crim>=7.46495 54   589.6593 11.40370 *
##        9) crim< 7.46495 66   673.9782 16.61818 *
##      5) lstat< 14.805 206  5080.9170 23.37767  
##       10) rm< 6.542 159  3543.5980 22.17233  
##         20) dis>=1.6009 152  1086.2130 21.64342 *
##         21) dis< 1.6009 7  1491.5570 33.65714 *
##       11) rm>=6.542 47   524.8362 27.45532 *
##    3) rm>=6.941 53  4305.8880 37.80566  
##      6) rm< 7.437 33  1523.8810 32.52424  
##       12) lstat>=6.97 7   577.0171 25.95714 *
##       13) lstat< 6.97 26   563.6985 34.29231 *
##      7) rm>=7.437 20   342.7320 46.52000 *

Error en TRAIN

En este paso lo que vamos a hacer es agregar la columna de la predicción a los datos de TRAIN, mediante la función bind_cols. Esto nos va a permitir ver la diferencia entre la predicción y el valor verdadero.

results <- p_train %>%
    bind_cols(predict(tree_fit, p_train) %>%
                  rename(.pred_tree = .pred))

Ahora vamos a ver las primeras filas de este nuevo dataset que hemos armado.

#veamos nuestra nueva tabla
head(results)
##       crim zn indus chas   nox    rm   age    dis rad tax ptratio  black lstat
## 1 45.74610  0 18.10    0 0.693 4.519 100.0 1.6582  24 666    20.2  88.27 36.98
## 2  6.65492  0 18.10    0 0.713 6.317  83.0 2.7344  24 666    20.2 396.90 13.99
## 3  0.06642  0  4.05    0 0.510 6.860  74.4 2.9153   5 296    16.6 391.27  6.92
## 4  0.62976  0  8.14    0 0.538 5.949  61.8 4.7075   4 307    21.0 396.90  8.26
## 5  0.01439 60  2.93    0 0.401 6.604  18.8 6.2196   1 265    15.6 376.70  4.38
## 6 15.86030  0 18.10    0 0.679 5.896  95.4 1.9096  24 666    20.2   7.68 24.39
##   medv .pred_tree
## 1  7.0   11.40370
## 2 19.5   21.64342
## 3 29.9   27.45532
## 4 20.4   21.64342
## 5 29.1   27.45532
## 6  8.3   11.40370

Vamos a calcular las metricas correspondientes del modelo. Para ello utilizamos la función metrics() disponible tambien en tidymodels. Agregamos los el dataset, y le decimos cual es el valor verdadero y cual es la predicción.

metrics(results, truth = medv, estimate = .pred_tree)
## # A tibble: 3 × 3
##   .metric .estimator .estimate
##   <chr>   <chr>          <dbl>
## 1 rmse    standard       3.93 
## 2 rsq     standard       0.818
## 3 mae     standard       2.83

Error de estimacion del modelo en TEST

Esta etapa es la predicción del modelo de ML, y es el resultado que generalmente se reporta. IMPORTANTE: los resultados se deben reportar SI O SI, puede o no ir acompañado de los resultados de TRAIN.

results_test <- p_test %>%
    bind_cols(predict(tree_fit, p_test) %>%
                  rename(.pred_tree = .pred))

Vemos la primera fila del dataset de test con las predicciones al final.

head(results_test)
##      crim   zn indus chas   nox    rm   age    dis rad tax ptratio  black lstat
## 1 0.00632 18.0  2.31    0 0.538 6.575  65.2 4.0900   1 296    15.3 396.90  4.98
## 2 0.02985  0.0  2.18    0 0.458 6.430  58.7 6.0622   3 222    18.7 394.12  5.21
## 3 0.14455 12.5  7.87    0 0.524 6.172  96.1 5.9505   5 311    15.2 396.90 19.15
## 4 0.21124 12.5  7.87    0 0.524 5.631 100.0 6.0821   5 311    15.2 386.63 29.93
## 5 0.63796  0.0  8.14    0 0.538 6.096  84.5 4.4619   4 307    21.0 380.02 10.26
## 6 1.05393  0.0  8.14    0 0.538 5.935  29.3 4.4986   4 307    21.0 386.85  6.58
##   medv .pred_tree
## 1 24.0   27.45532
## 2 28.7   21.64342
## 3 27.1   16.61818
## 4 16.5   16.61818
## 5 18.2   21.64342
## 6 23.1   21.64342

Ahora vemos las métricas.

metrics(results_test, truth = medv, estimate = .pred_tree)
## # A tibble: 3 × 3
##   .metric .estimator .estimate
##   <chr>   <chr>          <dbl>
## 1 rmse    standard       5.22 
## 2 rsq     standard       0.689
## 3 mae     standard       3.67

Este es el modelo más simple que podemos obtener con tidymodels. La próxima clase este modelo va a complejizarse mucho mas. :)

Vamos a plotear el árbol de decisión

Ingreso las librerias para hacer el plot del árbol.

library(rpart)
library(rpart.plot)

Plot del árbol, mediante dos funciones: extract_fit() y rpart.plot().

tree_fit %>%
  extract_fit_engine() %>%
  rpart.plot(roundint = FALSE, extra=1)

Ploteo las predicciones y el valor verdadero de los resultados en TEST.

results_test %>%
    ggplot(aes(medv, .pred_tree)) +
    geom_abline(lty = 2, color = "gray50") +
    geom_point(size = 1.5, alpha = 0.3, show.legend = FALSE)