Machine learning con R - Preparación de datos

Author

Mateo Vega

Capitulo 2 - Proceso de modelamiento

Paquetes que se van a usar

# datasets principales

ames = AmesHousing::make_ames()
attrition =  modeldata::attrition
mnist = dslabs::read_mnist()
url <- "https://koalaverse.github.io/homlr/data/my_basket.csv"
basket = readr::read_csv(url)

library(dplyr)
library(ggplot2)

# para modelamiento
library(rsample)
library(caret)
library(h2o)

# establecer h2o
h2o.no_progress()
h2o.init()

Mucha parte de machine learning está hecha con el paquete H2o, por lo tanto es mejor convertir los data sets a objetos H2O, H2O no maneja dactores ordenados, por lo tanto hay que arreglarlos antes de convertirlos.

ames.h2o = as.h2o(ames)

churn = attrition %>%
  mutate_if(is.ordered, .funs = factor, ordered = F)
churn.h2o = as.h2o(churn)

Separar datos

Necesitamos un algoritmo \(f(x)\) que prediga con precición valores futuros \(\hat{Y}\) basado en algunos features \(X\), queremos un algoritmo que funcione bien en nuestros datos pasados y que prediga el futuro con exactitud, esto se llama la generalizabilidad de nuestro algoritmo, para esto separamos nuestros datos en entrenamiento y prueba:

Entrenamiento: se usan para generar el conjunto de categorias, entrenar el algoritmo, sintonizar hiperparámetros, comparar modelos y el resto de actividades para escoger un modelo final.
Prueba: Si ya se tiene el modelo final, se usan para generar una evaluación del rendimiento del modelo.

Se recomienda separar en 70-30 o 80-20, mucho entrenamiento no nos daría buenas evaluaciones y mucha prueba no nos dejaría encontrar los mejores parámetros del modelo, las dos formas más comunes de separar los datos son el muestreo aleatorio simple y el muestreo estratificados.

Muestreo aleatorio simple

Este método no tiene en cuenta ningún atributo de los datos,

# R base
set.seed(1234)
index_1 = sample(1:nrow(ames), round(nrow(ames) * 0.7))
train_1 = ames[index_1, ]
test_1 = ames[-index_1, ]

# Usando Caret
index_2 = createDataPartition(ames$Sale_Price, p = 0.7,
                              list = F)
train_2 = ames[index_2, ]
test_2 = ames[-index_2, ]

# usando rsample
split_1 = initial_split(ames, prop = 0.7)
train_3 = training(split_1)
test_3 = testing(split_1)

# Usando h2o
split_2 = h2o.splitFrame(ames.h2o, ratios = 0.7,
                         seed = 123)
train_4 = split_2[[1]]
test_4 = split_2[[2]]

Con el suficiente tamaño de muestra las distribuciones de los métodos va a ser casi la misma.

Muestreo estratificado

Se usa si queremos controlar que nuestros datos de entrenamiento y prueba tengan distribuciones similares de \(Y\), común en problemas de clasificación donde la variable respuesta está muy desbalanceada (90% si 10% no), igualmente para regresión si en una muestra pequeña la variable respuesta se aleja de la normalidad.

# distribución de la variable attrition

table(churn$Attrition) %>% prop.table()


       No       Yes 
0.8387755 0.1612245

# muestreo estratificado

set.seed(123)
split_strat = initial_split(churn, prop = 0.7,
                            strat = 'Attrition')
train_strat = training(split_strat)
test_strat = testing(split_strat)

# distribución
table(train_strat$Attrition) %>% prop.table()


       No       Yes 
0.8394942 0.1605058

table(test_strat$Attrition) %>% prop.table()


       No       Yes 
0.8371041 0.1628959

Clases desbalanceadas

Se tiene cuando una clase tiene una proporción muy baja de observaciones (5% default vs 95% no default), esto se puede solucionar con muchos métodos, aquí se presentan up-sampling y down-sampling.

Down-sampling balancea el conjunto de datos reduciendo el tamaño de la clase abundante, este método se usa cuando la cantidad de datos son suficientes.

Up-sampling se usa cuando la cantidad de datos es insuficiente, balancea el conjunto de datos incrementando el tamaño de las muestras raras, se generan nuevas muestras raras usando repetición o bootstraping.

También se puede juntar los dos, Synthetic Minority Over-Sampling (SMOTE), ver ?caret::trainControl(), h2o::weights_column, h2o::balance_classes.

Creando modelos con R

Formas de crear una formula

Suponiendo que model_fn() es una función que recibe una formula, hay varias formas de crear un modelo:

# precio en funcion del barrio y del año
model_fn(Sale_Price ~ Neighborhood + Year_Sold,
         data = ames)

# variables + interacción
model_fn(Sale_Price ~ Neighborhood + Year_Sold +
           Neighborhood:Year_Sold, data = ames)

# todos los predictores
model_fn(Sale_Price ~ ., data = ames)

# Funciones en linea
model_fn(log10(Sale_Price) ~ ns(Longitude, df = 3) +
           ns(Latitude, df = 3), data = ames)

Métodos de re muestreo

Para medir el desempeño de nuestro modelo en la fase de entrenamiento usamos un método de validación, este se basa en separar los datos de entrenamiento para crear dos partes: entrenamiento y validación. Entrenamos el modelo en los datos de entrenamiento y estimamos el rendimiento en los datos de validación, esto puede ser malo si tenemos una muestra pequeña. Para esto usamos metodos de re muestreo, los métodos mas comúnes son validación cruzada k-fold y bootstraping.

Validación cruzada k-fold

Este método divide al azar el conjunto de entrenamiento en k grupos (folds) de casi el mismo tamaño, el modelo se ajusta en k - 1 grupos y el restante se usa para medir el rendimiento, este proceso se repite k veces, cada vez el grupo de validación cambia, esto resulta en k estimaciones de la generalización del error, se estima el promedio de estos y obtenemos una aproximación al error, se ha encontrado que k = 10 es de los más óptimos, este método tiene más variabilidad que bootstraping.

h2o.cv = h2o.glm(
  x = x,
  y = y,
  training_frame = ames.h2o,
  nfolds = 10
)

O se puede hacer a parte,

vfold_cv(ames, v = 10)

#  10-fold cross-validation 
# A tibble: 10 × 2
   splits             id    
   <list>             <chr> 
 1 <split [2637/293]> Fold01
 2 <split [2637/293]> Fold02
 3 <split [2637/293]> Fold03
 4 <split [2637/293]> Fold04
 5 <split [2637/293]> Fold05
 6 <split [2637/293]> Fold06
 7 <split [2637/293]> Fold07
 8 <split [2637/293]> Fold08
 9 <split [2637/293]> Fold09
10 <split [2637/293]> Fold10

Bootstraping

Es un muestreo aleatorio de los datos con reemplazo, el tamaño de muestra del bootstrap es el mismo que el conjunto de donde se sacó, los datos que no queden en la muestra se consideran out-of-bag (OOB), estos se pueden usar como datos de validación.

En muestras pequeñas debido a la repetición se puede generar sesgo, para \(n \geq 1000\) el sesgo es mucho menor.

bootstraps(ames, times = 10)

# Bootstrap sampling 
# A tibble: 10 × 2
   splits              id         
   <list>              <chr>      
 1 <split [2930/1062]> Bootstrap01
 2 <split [2930/1076]> Bootstrap02
 3 <split [2930/1066]> Bootstrap03
 4 <split [2930/1045]> Bootstrap04
 5 <split [2930/1087]> Bootstrap05
 6 <split [2930/1108]> Bootstrap06
 7 <split [2930/1075]> Bootstrap07
 8 <split [2930/1078]> Bootstrap08
 9 <split [2930/1053]> Bootstrap09
10 <split [2930/1067]> Bootstrap10

Sesgo contra varianza

Sesgo es la diferencia entre la predicción esperada de nuestro modelo y el valor correcto que estámos tratando de predecir, mide que tan lejos estamos de los valores correctos.

Varianza es la variabilidad de la predicción de un modelo para unos datos dados, modelos con mucha varianza son más probables de sufrir sobrejauste.

Sintonizar hiperparámetros

Los hiperparámetros son los ajustes para controlar la complejidad de los algoritmos de machine learning, no todos los algoritmos tienen hiperparámetros pero la mayoria si.

La mejor forma de sintonizar estos hiperparámetros es con una grid search, una busqueda automática de muchas combinaciones de hiperparámetros.

Evaluación de modelos

La manera más efectiva de evaluar un modelo es a través de la función de perdida, las funciones de perdida son metricas que comparan el valor predicho con el valor real.

Métricas

MSE: \(MSE = \frac{1}{n}\sum_{i=1}^n(y_i - \hat{y}_i)^2\), es el promedio de los errores al cuadrado. Objetivo: minimizar.
RMSE = \(RMSE = \sqrt{\frac{1}{n}\sum_{i=1}^n(y_i - \hat{y}_i)^2}\), toma la raiz del anterior para que la métrica tenga la misma unidad que la variable respuesta. Objetivo: minimizar.
Devianza: Presenta el grado al cual un modelo explica la variación de un conjunto de datos usando estimación de máxima verosimilitúd, compara el modelo saturado con el modelo no saturado, Objetivo: minimizar.
MAE: Error medio absoluto, \(MAE = \frac{1}{n}\sum_{i=1}^n(|y_i - \hat{y}_i|)\), menos enfasis en errores muy grandes, Objetivo: minimizar.
RMSLE: Raiz del cuadrado medio del error logarítmico, \(RMSLE = \sqrt{\frac{1}{n}\sum_{i=1}^n(log(y_i + 1) - log(\hat{y}_i + 1))^2}\), se usa cuando la variable respuesta tiene un rango muy grande, Objetivo: minimizar.
\(R^2\): proporción de la varianza en la variable dependiente explicada por el modelo, no hay que poner mucho cuidado en esta metrica, Objetivo: maximizar.

Para modelos de clasificación

mala clasificación: El error generalel porcentaje de observaciones mal clasificadas, por ejemplo se tienen 3 clases con 25, 30, 35 observaciones, si se clasifican mal 3, 6, 4, respectivamente la métrica es \(\frac{13}{90}\), 14%, Objetivo: minimizar.
Error medio por clase: El error promedio para cada clase, por ejemplo, la media de \(\frac{3}{25}, \frac{6}{30}, \frac{4}{35}\), 14.5%, Objetivo: minimizar.
MSE: Error cuadrado medio, computa la distancia de 1 a la probabilidad sugerida, por ejemplo, tenemos tres clases y el modelo predice probabilidades de 0.91, 0.07, 0.02, si la respuesta correcta es la primer clase entonces \(MSE = 0.09^2\), si la respuesta correcta es la segunda clase, entonces \(MSE = 0.93^2\), Objetivo: minimizar.
Cross-entropy: Similar al MSE pero incopora logaritmo a la probabilidad predicha multiplicada por la clase verdadera, Objetivo: minimizar.
Indice de Gini: Usada principalmente en métodos de árboles, se refiere a la pureza donde un valor pequeño indica que un nodo contiene observaciones predominantes de una sola clase, Objetivo: minimizar.

Ahora, dada la matriz de confusión se pueden obtener las siguientes métricas:

Accuracy (Exactitud): Que tan seguido es el modelo correcto? \(\frac{TP+TN}{total}\), Objetivo: maximizar.
Precision: Que tan exacto el modelo predice eventos? Se preocupa por maximizar la proporción de positivos verdaderos (TP) frente a falsos positivos (FP), de las predicciones que hicimos, cuantas fueron correctas? \(\frac{TP}{TP + FP}\), Objetivo: maximizar.
Specificity (recall): Que tan exacto el modelo clasifica eventos reales? se busca maximizar la proporción de verdaderos positivos frente a falsos negativos, de los eventos que ocurrieron, cuantos predijimos? \(\frac{TP}{TP+FN}\), Objetivo: maximizar.
Specificity: Que tan exacto el modelo clasifica no eventos reales, \(\frac{TN}{TN+FP}\), Objetivo: maximizar.
AUC: Area debajo de la curva, un buen clasificador tienen alta precisión y sensitividad, el modelo clasifica bien cuando predice un evento y cuando no va a suceder, minimiza falsos positivos y falsos negativos, para capturar este balance usamos una curva ROC que grafica el porcentaje de falsos positivos en el eje x y los verdaderos positivos en el eje y, la linea diagonal dice que nuestro modelo no es mejor que adivinar, entre más cerca la linea a la esquina izquierda de arriba, mejor.

Ejemplo del proceso

Primero hacemos un muestreo estratificado con el paquete rsample

set.seed(123)
split = initial_split(ames, prop = 0.7, 
                      strata = 'Sale_Price')
ames_train = training(split)
ames_test = testing(split)

Se hace una regresión KNN a nuestros datos usando caret usando los siguientes pasos:

Metodo de re muestreo: usamos un CV - 10 fold repetido 5 veces.
Grid search: especificamos los hiperparámetros a buscar \((k = 2, 3, 4, \dots, 25)\).
Entrenamiento de modelo y validación: entrenamos el modelo knn y usamos la metrica RMSE.

# especificar estrategia de re muestreo
cv = trainControl(
  method = 'repeatedcv',
  number = 10,
  repeats = 5
)

# crear una cuadricula de hiperparámetros
hyper_grid = expand.grid(k = seq(2, 25, by = 1))

# sintonizar un modelo knn usando la cuadricula
knn_fit = train(
  Sale_Price ~ ., 
  data = ames_train,
  method = 'knn',
  trControl = cv,
  tuneGrid = hyper_grid,
  metric = 'RMSE'
)

knn_fit

k-Nearest Neighbors 

2049 samples
  80 predictor

No pre-processing
Resampling: Cross-Validated (10 fold, repeated 5 times) 
Summary of sample sizes: 1844, 1844, 1843, 1844, 1844, 1845, ... 
Resampling results across tuning parameters:

  k   RMSE      Rsquared   MAE     
   2  47206.74  0.6596344  31133.25
   3  45684.91  0.6773635  30007.71
   4  44771.97  0.6901038  29276.09
   5  44064.54  0.7005790  28996.93
   6  43846.05  0.7045166  28895.48
   7  43858.13  0.7059401  28883.74
   8  44181.13  0.7033657  29055.05
   9  44352.03  0.7028597  29109.00
  10  44332.46  0.7053884  29129.62
  11  44282.81  0.7083442  29081.39
  12  44486.34  0.7075253  29155.41
  13  44647.15  0.7076206  29256.68
  14  44790.79  0.7077073  29307.11
  15  45041.02  0.7063767  29423.86
  16  45119.37  0.7073844  29484.00
  17  45264.01  0.7070891  29586.22
  18  45366.02  0.7072968  29641.92
  19  45537.84  0.7066304  29766.52
  20  45746.63  0.7052851  29907.89
  21  45983.35  0.7031524  30058.55
  22  46187.92  0.7017539  30192.39
  23  46406.95  0.7001361  30329.98
  24  46605.01  0.6986611  30483.94
  25  46824.70  0.6971044  30617.33

RMSE was used to select the optimal model using the smallest value.
The final value used for the model was k = 6.

Encontramos un k óptimo igual a 6, con un \(RMSE = 43846.05\).

ggplot(knn_fit)

Capitulo 3 - Feature engineering

En esta parte se pre procesan los datos, diferentes modelos tienen diferentes necesidades para funcionar mejor. En este capitulo se aplican estas técnicas.

# Helper packages
library(dplyr)    # for data manipulation
library(ggplot2)  # for awesome graphics
library(visdat)   # for additional visualizations

# Feature engineering packages
library(caret)    # for various ML tasks
library(recipes)  # for feature engineering tasks

# adicionales
library(purrr)
library(tidyr)
library(readr)
library(kableExtra)

Transformaciones

Por ejemplo, la regresión lineal ordinal asume que los errores y la variable respuesta está distribuida normal, a veces una simple transformación logaritmica puede transformar la variable a una distribución normal,

modelos = c('Residuales del modelo sin transformar', 
            'Residuales del modelo con transformación logaritmica')
list(
  m1 = lm(Sale_Price ~ Year_Built, data = ames_train),
  m2 = lm(log(Sale_Price ) ~ Year_Built, data = ames_train)
) %>%
  map2_dfr(modelos, ~ broom::augment(.x) %>% mutate(model = .y)) %>%
  ggplot(aes(.resid)) +
    geom_histogram(bins = 75) +
    facet_wrap(~ model, scales = 'free_x') +
    ylab(NULL) +
    xlab('Residuales')

Vemos que después de la transformación logaritmica los residuos parecen distribuirse de forma más normal, para datos con sesgo positivo hay dos opciones:

Opcion 1: Normalizar con una transformación log,

variable_transformada = log(ames_train$Sale_Price)

Sin embargo se quiere representar el pre procesamiento como un plano de lo que hacemos para que se puede re utilizar, para esto usamos la libreria recipe,

ames_recipe = recipe(Sale_Price ~ ., data = ames_train) %>%
  step_log(all_outcomes())
ames_recipe

── Recipe ──────────────────────────────────────────────────────────────────────

── Inputs

Number of variables by role

outcome:    1
predictor: 80

── Operations

• Log transformation on: all_outcomes()

Esta transformación logartimica solo se puede hacer para valores positivos, si los valores toman valores entre (-0.99 y 0), podemos sumar 1 y hacer la transformación log1p(), si los valores son menores a -1, se puede usar la transformación Yeo-Johnson.

log(-0.5)

Warning in .Primitive("log")(x, ...): Se han producido NaNs

[1] NaN

log1p(-0.5)

[1] -0.6931472

Opción 2: Usar la transformación Box Cox, esta trata de aproximar los datos a una distribución normal

# Log transformation
train_log_y <- log(ames_train$Sale_Price)
test_log_y  <- log(ames_train$Sale_Price)

# Box Cox transformation
lambda  <- forecast::BoxCox.lambda(ames_train$Sale_Price)

Registered S3 method overwritten by 'quantmod':
  method            from
  as.zoo.data.frame zoo

train_bc_y <- forecast::BoxCox(ames_train$Sale_Price, lambda)
test_bc_y  <- forecast::BoxCox(ames_test$Sale_Price, lambda)

# Plot differences
levs <- c("Normal", "Log_Transform", "BoxCox_Transform")
data.frame(
  Normal = ames_train$Sale_Price,
  Log_Transform = train_log_y,
  BoxCox_Transform = train_bc_y
) %>%
  gather(Transform, Value) %>%
  mutate(Transform = factor(Transform, levels = levs)) %>% 
  ggplot(aes(Value, fill = Transform)) +
    geom_histogram(show.legend = FALSE, bins = 40) +
    facet_wrap(~ Transform, scales = "free_x")

Warning: attributes are not identical across measure variables; they will be
dropped

Valores perdidos

Hay dos grandes razones para los datos perdidos:

Valores perdidos informativos: Aca el valor perdido tiene sentido, se podría pensar en crear una categoria para estos datos.
Valores perdidos aleatorios: Estos datos son independientes de la recolección, son valores perdidos de verdad y se deberían imputar.

Visualizar valores perdidos

Es importante entender la distribución de los datos perdidos

sum(is.na(AmesHousing::ames_raw))

[1] 13997

Ahora para ver la dsitribución,

AmesHousing::ames_raw %>%
  is.na() %>%
  reshape2::melt() %>%
  ggplot(aes(Var2, Var1, fill = value)) +
    geom_raster() +
    coord_flip() + 
    scale_y_continuous(NULL, expand = c(0, 0)) +
    scale_fill_grey(name = '',
                    labels = c('Presente',
                               'Perdido')) +
    xlab('Observaciones') +
    theme(axis.text.y = element_text(size = 4))

En este caso la variable Garage_xx tiene muchos NA, tal vez no sabían como clasificar las casas sin garage, en este caso sería mejor reemplazar los NA con su propia categoria ‘None’.

AmesHousing::ames_raw %>%
  filter(is.na(`Garage Type`)) %>%
  select(`Garage Type`, `Garage Cars`, `Garage Area`)

# A tibble: 157 × 3
   `Garage Type` `Garage Cars` `Garage Area`
   <chr>                 <int>         <int>
 1 <NA>                      0             0
 2 <NA>                      0             0
 3 <NA>                      0             0
 4 <NA>                      0             0
 5 <NA>                      0             0
 6 <NA>                      0             0
 7 <NA>                      0             0
 8 <NA>                      0             0
 9 <NA>                      0             0
10 <NA>                      0             0
# ℹ 147 more rows

Otra forma de visualizarlo,

vis_miss(AmesHousing::ames_raw, cluster = T)

Imputación

Es el proceso de reemplazar valores perdidos con un sustituto, hay varios métodos.

Con un valor fijo

La forma más sencilla es con una valor fijo como la media, la mediana o la moda, aunque fácil este método no tiene enc cuenta los atributos de los datos, no es aconsejable usarlo, importante recordar que la imputación deber ser hecha dentro del proceso de re muestreo.

ames_recipe %>%
  step_impute_median(Gr_Liv_Area)

── Recipe ──────────────────────────────────────────────────────────────────────

── Inputs

Number of variables by role

outcome:    1
predictor: 80

── Operations

• Log transformation on: all_outcomes()

• Median imputation for: Gr_Liv_Area

KNN

Identifica los valores más parecidos al valor perdido con respecto a otras categorias, este método es bastante pesado para conjuntos de datos muy grandes.

ames_recipe %>%
  step_impute_knn(all_predictors(), neighbors = 5)

── Recipe ──────────────────────────────────────────────────────────────────────

── Inputs

Number of variables by role

outcome:    1
predictor: 80

── Operations

• Log transformation on: all_outcomes()

• K-nearest neighbor imputation for: all_predictors()

Basado en árboles

ames_recipe %>% 
  step_impute_bag(all_predictors())

── Recipe ──────────────────────────────────────────────────────────────────────

── Inputs

Number of variables by role

outcome:    1
predictor: 80

── Operations

• Log transformation on: all_outcomes()

• Bagged tree imputation for: all_predictors()

Filtrado de features

Es muy importante filtrar las variables que realmente sirven para el modelo, aunque algunos modelos no sufren tanto por variables no informativas la computación puede volverse pesada.

Las primeras variables que debemos eliminar son las que tengan varianza cero o casi cero, no aportan información importante al modelo, para detectar variables con varianza casi cero se puede tener en cuenta:

La fracción de valores únicos es baja (menos de 10%).
El radio de la frecuencia del valor más prevalente al segundo más prevalente es grande (mas de 20%).

Si ambos criterios son verdad es mejor quitar las variables del modelo, por ejemplo,

caret::nearZeroVar(ames_train, saveMetrics = T) %>%
  tibble::rownames_to_column() %>%
  filter(nzv)

              rowname  freqRatio percentUnique zeroVar  nzv
1              Street  226.66667    0.09760859   FALSE TRUE
2               Alley   24.25316    0.14641288   FALSE TRUE
3        Land_Contour   19.50000    0.19521718   FALSE TRUE
4           Utilities 1023.00000    0.14641288   FALSE TRUE
5          Land_Slope   22.15909    0.14641288   FALSE TRUE
6         Condition_2  202.60000    0.34163006   FALSE TRUE
7           Roof_Matl  144.35714    0.39043436   FALSE TRUE
8           Bsmt_Cond   20.24444    0.29282577   FALSE TRUE
9      BsmtFin_Type_2   25.85294    0.34163006   FALSE TRUE
10       BsmtFin_SF_2  453.25000    9.37042460   FALSE TRUE
11            Heating  106.00000    0.29282577   FALSE TRUE
12    Low_Qual_Fin_SF 1010.50000    1.31771596   FALSE TRUE
13      Kitchen_AbvGr   21.23913    0.19521718   FALSE TRUE
14         Functional   38.89796    0.39043436   FALSE TRUE
15     Enclosed_Porch  102.05882    7.41825281   FALSE TRUE
16 Three_season_porch  673.66667    1.12249878   FALSE TRUE
17       Screen_Porch  169.90909    4.63640800   FALSE TRUE
18          Pool_Area 2039.00000    0.53684724   FALSE TRUE
19            Pool_QC  509.75000    0.24402147   FALSE TRUE
20       Misc_Feature   34.18966    0.24402147   FALSE TRUE
21           Misc_Val  180.54545    1.56173743   FALSE TRUE

Podemos usar step_nzv o step_zv para quitar estas variables.

Variables numericas

Las variables numericas pueden crear problemas cuando tienen sesgo, contienen outliers o tienen muchos rangos, esto se puede arreglar en parte normalizando y estandarizando variables con mucho sesgo/

Sesgo

Los modelos parámetricos necesitan el supuesto de normalidad, cuando se normalizan variables es mejor usa Box Cox o Yeo Johnson.

recipe(Sale_Price ~ ., data = ames_train) %>%
  step_YeoJohnson(all_numeric())

── Recipe ──────────────────────────────────────────────────────────────────────

── Inputs

Number of variables by role

outcome:    1
predictor: 80

── Operations

• Yeo-Johnson transformation on: all_numeric()

Estandarización

Para algunos modelos es mejor estandarizar las variables para que queden en una unidad comparable con media cero y varianza unitarias, se debe estandarizar en el plano para que los datos de entrenamiento y prueba tengan la misma media y varianza y se minimize la perdida de datos.

ames_recipe %>% 
  step_center(all_numeric(), -all_outcomes()) %>%
  step_scale(all_numeric(), -all_outcomes())

── Recipe ──────────────────────────────────────────────────────────────────────

── Inputs

Number of variables by role

outcome:    1
predictor: 80

── Operations

• Log transformation on: all_outcomes()

• Centering for: all_numeric(), -all_outcomes()

• Scaling for: all_numeric(), -all_outcomes()

Variables categóricas

Agrupar

Cuando se tienen variables con niveles con pocas observaciones, por ejemplo en ames hay 28 barrios únicos,

count(ames_train, Neighborhood) %>% arrange(n)

# A tibble: 28 × 2
   Neighborhood            n
   <fct>               <int>
 1 Landmark                1
 2 Green_Hills             2
 3 Greens                  3
 4 Blueste                 8
 5 Veenker                15
 6 Northpark_Villa        17
 7 Bloomington_Heights    18
 8 Meadow_Village         22
 9 Briardale              23
10 Clear_Creek            26
# ℹ 18 more rows

count(ames_train, Screen_Porch) %>% arrange(n)

# A tibble: 95 × 2
   Screen_Porch     n
          <int> <int>
 1           40     1
 2           53     1
 3           60     1
 4           64     1
 5           80     1
 6           84     1
 7           88     1
 8           94     1
 9           99     1
10          104     1
# ℹ 85 more rows

Es mejor agrupar estas categorias con pocos valores en una nueva categoria, por ejemplo los niveles con menos del 10% de la muestra.

ames_train$Screen_Porch = as.factor(ames_train$Screen_Porch)
lumping <- recipe(Sale_Price ~ ., data = ames_train) %>%
    step_other(Neighborhood, threshold = 0.01, 
               other = "other")

# aplicarlo el plano 

apply_2_training <- prep(lumping, training = ames_train) %>%
    bake(ames_train)

count(apply_2_training, Neighborhood) %>% arrange(n)

# A tibble: 22 × 2
   Neighborhood                                n
   <fct>                                   <int>
 1 Meadow_Village                             22
 2 Briardale                                  23
 3 Clear_Creek                                26
 4 South_and_West_of_Iowa_State_University    32
 5 Stone_Brook                                40
 6 Northridge                                 51
 7 Timberland                                 52
 8 other                                      64
 9 Brookside                                  74
10 Crawford                                   75
# ℹ 12 more rows

One hot y codificación dummy

Muchos modelos necesitan que todas las variables sean numericas, por lo tanto se transforman las variables categóricas a numericas, el one hot transforma la variable categorica para que cada nivel de la variable sea representada por un booleano, esto causa colinealidad, por lo tanto se quita un nivel, esto es la codificacion dummy.

recipe(Sale_Price ~., data = ames_train) %>%
  step_dummy(all_nominal(), one_hot = T)

── Recipe ──────────────────────────────────────────────────────────────────────

── Inputs

Number of variables by role

outcome:    1
predictor: 80

── Operations

• Dummy variables from: all_nominal()

Si se tienen variables con muchas categorias la dimensión del conjunto de datos se aumenta mucho, ahí sería mejor una codificación ordianl u otras alternativa.

Codificación de etiquetas

Es una organización numerica de una variable categorica, puede ser en orden u orden alfabetico,

count(ames_train, MS_SubClass)

# A tibble: 15 × 2
   MS_SubClass                                   n
   <fct>                                     <int>
 1 One_Story_1946_and_Newer_All_Styles         756
 2 One_Story_1945_and_Older                     94
 3 One_Story_with_Finished_Attic_All_Ages        4
 4 One_and_Half_Story_Unfinished_All_Ages       13
 5 One_and_Half_Story_Finished_All_Ages        203
 6 Two_Story_1946_and_Newer                    404
 7 Two_Story_1945_and_Older                     90
 8 Two_and_Half_Story_All_Ages                  14
 9 Split_or_Multilevel                          85
10 Split_Foyer                                  37
11 Duplex_All_Styles_and_Ages                   76
12 One_Story_PUD_1946_and_Newer                132
13 Two_Story_PUD_1946_and_Newer                 86
14 PUD_Multilevel_Split_Level_Foyer             12
15 Two_Family_conversion_All_Styles_and_Ages    43

# codificado

recipe(Sale_Price ~ ., data = ames_train) %>%
  step_integer(MS_SubClass) %>%
  prep(ames_train) %>%
  bake(ames_train) %>%
  count(MS_SubClass)

# A tibble: 15 × 2
   MS_SubClass     n
         <int> <int>
 1           1   756
 2           2    94
 3           3     4
 4           4    13
 5           5   203
 6           6   404
 7           7    90
 8           8    14
 9           9    85
10          10    37
11          11    76
12          12   132
13          14    86
14          15    12
15          16    43

Hay que tenre cuidado ya que se transforma a una variable ordinal, si la categorica no tiene orden esto no tendría sentido, si la variable categórica tiene un orden entonces es directa la transformación.

ames_train %>% select(contains('Qual'))

# A tibble: 2,049 × 6
   Overall_Qual  Exter_Qual Bsmt_Qual   Low_Qual_Fin_SF Kitchen_Qual Garage_Qual
   <fct>         <fct>      <fct>                 <int> <fct>        <fct>      
 1 Average       Typical    Typical                   0 Typical      Typical    
 2 Above_Average Typical    Typical                   0 Typical      Typical    
 3 Good          Typical    Good                      0 Typical      Typical    
 4 Above_Average Typical    Good                      0 Typical      Typical    
 5 Average       Typical    Typical                   0 Typical      Typical    
 6 Above_Average Typical    No_Basement               0 Typical      Typical    
 7 Average       Typical    Typical                   0 Typical      Typical    
 8 Average       Typical    Typical                   0 Typical      Typical    
 9 Average       Typical    Typical                   0 Typical      No_Garage  
10 Average       Typical    Typical                   0 Typical      Typical    
# ℹ 2,039 more rows

Estas variables pueden transformarse a una variable ordinal.

count(ames_train, Overall_Qual)

# A tibble: 10 × 2
   Overall_Qual       n
   <fct>          <int>
 1 Very_Poor          3
 2 Poor               9
 3 Fair              33
 4 Below_Average    147
 5 Average          580
 6 Above_Average    521
 7 Good             408
 8 Very_Good        242
 9 Excellent         82
10 Very_Excellent    24

recipe(Sale_Price ~ ., data = ames_train) %>%
  step_integer(Overall_Qual) %>%
  prep(ames_train) %>%
  bake(ames_train) %>%
  count(Overall_Qual)

# A tibble: 10 × 2
   Overall_Qual     n
          <int> <int>
 1            1     3
 2            2     9
 3            3    33
 4            4   147
 5            5   580
 6            6   521
 7            7   408
 8            8   242
 9            9    82
10           10    24

Reducción de dimensión

# pca
recipe(Sale_Price ~ ., data = ames_train) %>%
  step_center(all_numeric()) %>%
  step_scale(all_numeric()) %>%
  step_pca(all_numeric(), threshold = 0.95)

── Recipe ──────────────────────────────────────────────────────────────────────

── Inputs

Number of variables by role

outcome:    1
predictor: 80

── Operations

• Centering for: all_numeric()

• Scaling for: all_numeric()

• PCA extraction with: all_numeric()

Implementación

Si se va a usar Box Cox que sea antes de centrar los datos ya que eso volveria los datos negativos, o se puede usar Yeo Johnson.
One hot o dummy resulta en datos sparse, si se estandariza despues se queda con datos densos, y se pierde la eficiencia computacional, es mejor estandarizar primero y despues hacer la codificacion.
Si se agrupan categoriasinfrecuentes juntas, hacerlo antes de la one hot / dummy.
Primero se hace la reducción de dimensión en variables numéricas.

Un paso a paso puede ser:

Filtrar variables con varianza cero o casi cero.
Hacer imputación.
Normalizar para resolver sesgo numerico.
Estandarizar variables numericas.
Hace reducción de dimensión en variables numericas.
One hot o Dummy codificación en variables categoricas.

Fuga de datos

La fuga de datos pasa en el pre procesamiento de los datos, para minimizar esto se debe hacer este feature engineering isolado de cada iteración de remuestreo, por ejemplo, cada datos de entrenamiento re muestreado debe usar su propia media y varianza y estos valores aplicarlos a los datos de prueba.

Ejemplo del proceso

Primero una introducción del paquete recipe, hay tres pasos:

recipe: defines tus pasos de feature engineering para crear el plano.
prep: estimar parámetros de feature engineering basado en los datos de entrenamiento.
bake: aplicar el plano a nuevos datos.

Por ejemplo:

Eliminar variables con varianza cero.
Codificar nuestras variables ordinales (las que tienen sentido).
Centrar y escalar (estandarizar) todas las variables numericas.
Reducir dimensión aplicando PCA a las variables numericas.

plano = recipe(Sale_Price ~ ., data = ames_train) %>%
  step_nzv(all_nominal()) %>%
  step_integer(matches('Qual|Cond|QC|Qu')) %>%
  step_center(all_numeric(), - all_outcomes()) %>%
  step_scale(all_numeric(), -all_outcomes()) %>%
  step_pca(all_numeric(), -all_outcomes())

plano

── Recipe ──────────────────────────────────────────────────────────────────────

── Inputs

Number of variables by role

outcome:    1
predictor: 80

── Operations

• Sparse, unbalanced variable filter on: all_nominal()

• Integer encoding for: matches("Qual|Cond|QC|Qu")

• Centering for: all_numeric(), -all_outcomes()

• Scaling for: all_numeric(), -all_outcomes()

• PCA extraction with: all_numeric(), -all_outcomes()

Luego se entrena el plano en los datos de entrenamiento

prepare = prep(plano, training = ames_train) 
prepare

── Recipe ──────────────────────────────────────────────────────────────────────

── Inputs

Number of variables by role

outcome:    1
predictor: 80

── Training information

Training data contained 2049 data points and no incomplete rows.

── Operations

• Sparse, unbalanced variable filter removed: Street, Alley, ... | Trained

• Integer encoding for: Condition_1, Overall_Qual, Overall_Cond, ... | Trained

• Centering for: Lot_Frontage, Lot_Area, Condition_1, ... | Trained

• Scaling for: Lot_Frontage, Lot_Area, Condition_1, Overall_Qual, ... | Trained

• PCA extraction with: Lot_Frontage, Lot_Area, Condition_1, ... | Trained

Por ultimos aplicamos el plano a datos nuevos

baked_train = bake(prepare, new_data = ames_train)
baked_test = bake(prepare, new_data = ames_test)

Warning:  There was 1 column that was a factor when the recipe was prepped:
 'Screen_Porch'.
 This may cause errors when processing new data.

baked_train

# A tibble: 2,049 × 27
   MS_SubClass MS_Zoning Lot_Shape Lot_Config Neighborhood Bldg_Type House_Style
   <fct>       <fct>     <fct>     <fct>      <fct>        <fct>     <fct>      
 1 Two_Story_… Resident… Regular   Inside     Briardale    Twnhs     Two_Story  
 2 Two_Story_… Resident… Regular   Inside     Briardale    Twnhs     Two_Story  
 3 One_Story_… Resident… Regular   FR2        Northpark_V… Twnhs     One_Story  
 4 One_Story_… Resident… Regular   Inside     Sawyer_West  TwnhsE    One_Story  
 5 One_Story_… Resident… Regular   Corner     Sawyer_West  OneFam    One_Story  
 6 Duplex_All… Resident… Regular   Inside     Sawyer       Duplex    One_and_Ha…
 7 One_Story_… Resident… Slightly… Inside     Sawyer       OneFam    One_Story  
 8 One_Story_… Resident… Slightly… Corner     Sawyer       OneFam    One_Story  
 9 Duplex_All… Resident… Slightly… Inside     North_Ames   Duplex    One_Story  
10 One_Story_… Resident… Regular   Inside     North_Ames   OneFam    One_Story  
# ℹ 2,039 more rows
# ℹ 20 more variables: Roof_Style <fct>, Exterior_1st <fct>,
#   Exterior_2nd <fct>, Mas_Vnr_Type <fct>, Foundation <fct>,
#   Bsmt_Exposure <fct>, BsmtFin_Type_1 <fct>, Central_Air <fct>,
#   Electrical <fct>, Garage_Type <fct>, Garage_Finish <fct>,
#   Paved_Drive <fct>, Fence <fct>, Sale_Type <fct>, Sale_Price <int>,
#   PC1 <dbl>, PC2 <dbl>, PC3 <dbl>, PC4 <dbl>, PC5 <dbl>

Ahora, la meta es desarrollar el plano y con cada iteracion del remuestreo aplicar prep y bake a nuestros datos de entrenamiento remuestreados, el paquete caret hace esto por nosotros.

plano = recipe(Sale_Price ~ ., data = ames_train) %>%
  step_nzv(all_nominal()) %>%
  step_integer(matches('Qual|Cond|QC|Qu')) %>%
  step_center(all_numeric(), - all_outcomes()) %>%
  step_scale(all_numeric(), -all_outcomes()) %>%
  step_dummy(all_nominal(), -all_outcomes(), one_hot = T)

Luego aplicamos el metodo de remuestreo y la busqueda de hiperparámetros,

# especificar el plan de remuestreo
cv = trainControl(
  method = 'repeatedcv',
  number = 10,
  repeats = 5
)

# Construir la cuadricula de los valores de hiperparámetros
hyper_grid = expand.grid(k = seq(2, 25, by = 1))

# sintonizar un modelo knn usando la busqueda de cuadricula
knn_fit2 = train(
  plano, 
  data = ames_train,
  method = 'knn',
  trControl = cv,
  tuneGrid = hyper_grid,
  metric = 'RMSE'
)

Mostramos ahora los resultados,

knn_fit2

k-Nearest Neighbors 

2049 samples
  80 predictor

Recipe steps: nzv, integer, center, scale, dummy 
Resampling: Cross-Validated (10 fold, repeated 5 times) 
Summary of sample sizes: 1844, 1844, 1845, 1843, 1844, 1844, ... 
Resampling results across tuning parameters:

  k   RMSE      Rsquared   MAE     
   2  35771.61  0.8080363  22998.57
   3  34678.41  0.8222120  21991.80
   4  34410.74  0.8275060  21614.02
   5  34391.61  0.8292146  21379.88
   6  34142.04  0.8333363  21177.03
   7  33939.18  0.8364766  20935.74
   8  33620.35  0.8406797  20738.86
   9  33482.66  0.8427131  20702.59
  10  33494.32  0.8430202  20747.83
  11  33492.38  0.8436009  20834.21
  12  33495.78  0.8442750  20851.08
  13  33551.10  0.8444398  20872.45
  14  33752.07  0.8432682  20965.94
  15  33814.27  0.8432751  20978.73
  16  33927.15  0.8429191  21040.96
  17  34011.68  0.8428927  21113.15
  18  34141.46  0.8420149  21161.50
  19  34269.05  0.8410037  21234.79
  20  34445.50  0.8394513  21334.53
  21  34619.40  0.8383019  21454.17
  22  34741.25  0.8375805  21514.74
  23  34810.59  0.8375506  21574.76
  24  34950.92  0.8366352  21662.61
  25  34986.06  0.8365953  21684.43

RMSE was used to select the optimal model using the smallest value.
The final value used for the model was k = 9.

ggplot(knn_fit2)

Se redujo el error por más de 10000 dólares.