Número de observaciones y valores ausentes

Junto con el estudio del tipo de variables, es básico conocer el número de observaciones disponibles y si todas ellas están completas. Los valores ausentes son muy importantes a la hora de crear modelos, algunos algoritmos no aceptan observaciones incompletas o bien se ven muy influenciados por ellas. Aunque la imputación de valores ausentes es parte del preprocesado y, por lo tanto, debe de aprenderse únicamente con los datos de entrenamiento, su identificación se tiene que realizar antes de separar los datos para asegurar que se establecen todas las estrategias de imputación necesarias

# Número de datos ausentes por variable
datos %>% map_dbl(.f = function(x){sum(is.na(x))})

##             precio     metros_totales         antiguedad     precio_terreno 
##                  0                  0                  0                  0 
##  metros_habitables     universitarios        dormitorios           chimenea 
##                  0                  0                  0                  0 
##             banyos       habitaciones        calefaccion consumo_calefacion 
##                  0                  0                  0                  0 
##            desague        vistas_lago nueva_construccion aire_acondicionado 
##                  0                  0                  0                  0

plot_missing(
  data    = datos, 
  title   = "Porcentaje de valores ausentes",
  ggtheme = theme_bw(),
  theme_config = list(legend.position = "none")
)

Variable respuesta

Cuando se crea un modelo, es muy importante estudiar la distribución de la variable respuesta, ya que, a fin de cuentas, es lo que interesa predecir. La variable precio tiene una distribución asimétrica con una cola positiva debido a que, unas pocas viviendas, tienen un precio muy superior a la media. Este tipo de distribución suele visualizarse mejor tras aplicar el logarítmica o la raíz cuadrada.

p1 <- ggplot(data = datos, aes(x = precio)) +
      geom_density(fill = "steelblue", alpha = 0.8) +
      geom_rug(alpha = 0.1) +
      scale_x_continuous(labels = scales::comma) +
      labs(title = "Distribución original") +
      theme_bw() 

p2 <- ggplot(data = datos, aes(x = sqrt(precio))) +
      geom_density(fill = "steelblue", alpha = 0.8) +
      geom_rug(alpha = 0.1) +
      scale_x_continuous(labels = scales::comma) +
      labs(title = "Transformación raíz cuadrada") +
      theme_bw() 

p3 <- ggplot(data = datos, aes(x = log(precio))) +
      geom_density(fill = "steelblue", alpha = 0.8) +
      geom_rug(alpha = 0.1) +
      scale_x_continuous(labels = scales::comma) +
      labs(title = "Transformación logarítmica") +
      theme_bw() 

ggarrange(p1, p2, p3, ncol = 1, align = "v")

# Tabla de estadísticos principales 
summary(datos$precio)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    5000  145000  189900  211967  259000  775000

Algunos modelos de machine learning y aprendizaje estadístico requieren que la variable respuesta se distribuya de una forma determinada. Por ejemplo, para los modelos de regresión lineal (LM), la distribución tiene que ser de tipo normal. Para los modelos lineales generalizados (GLM), la distribución tiene que ser de la familia exponencial.

Existen varios paquetes en R que permiten identificar a qué distribución se ajustan mejor los datos, uno de ellos es univariateML. Para conocer más sobre cómo identificar distribuciones consultar Ajuste de distribuciones con R.

# Se comparan únicamente las distribuciones con un dominio [0, +inf)
# Cuanto menor el valor AIC mejor el ajuste
comparacion_aic <- AIC(
                    mlbetapr(datos$precio),
                    mlexp(datos$precio),
                    mlinvgamma(datos$precio),
                    mlgamma(datos$precio),
                    mllnorm(datos$precio),
                    mlrayleigh(datos$precio),
                    mlinvgauss(datos$precio),
                    mlweibull(datos$precio),
                    mlinvweibull(datos$precio),
                    mllgamma(datos$precio)
                   )
comparacion_aic %>% rownames_to_column(var = "distribucion") %>% arrange(AIC)

Las dos distribuciones con mejor ajuste acorde al valor AIC son la normal y la gamma.

Variables continuas

plot_density(
  data    = datos %>% select(-precio),
  ncol    = 3,
  title   = "Distribución variables continuas",
  ggtheme = theme_bw(),
  theme_config = list(
                  plot.title = element_text(size = 16, face = "bold"),
                  strip.text = element_text(colour = "black", size = 12, face = 2)
                 )
  )

La variable chimenea, aunque es de tipo numérico, apenas toma unos pocos valores y la gran mayoría de observaciones pertenecen a solo dos de ellos. En casos como este, suele ser conveniente tratar la variable como cualitativa.

# Valores observados de chimenea.
table(datos$chimenea)

## 
##   0   1   2   3   4 
## 740 942  42   2   2

# Se convierte la variable chimenea a factor.
datos <- datos %>%
         mutate(chimenea = as.factor(chimenea))

Como el objetivo del estudio es predecir el precio de las viviendas, el análisis de cada variable se hace también en relación a la variable respuesta precio. Analizando los datos de esta forma, se pueden empezar a extraer ideas sobre qué variables están más relacionadas con el precio y de qué forma. Las variables antiguedad, metros_totales y precio_terreno se convierten a escala logarítmica para mejorar su representación gráfica.

datos <- datos %>%
         mutate(
           # Se le añade +0.1 a la antigüedad para que cuando la antigüedad es
           # 0 no de -Inf
           log_antiguedad     = log10(antiguedad + 0.1),
           log_metros_totales = log10(metros_totales),
           log_precio_terreno = log10(precio_terreno)
         )

custom_corr_plot <- function(variable1, variable2, df, alpha=0.3){
  p <- df %>%
       mutate(
         # Truco para que se ponga el título estilo facet
        title = paste(toupper(variable2), "vs", toupper(variable1))
       ) %>%
       ggplot(aes(x = !!sym(variable1), y = !!sym(variable2))) + 
       geom_point(alpha = alpha) +
       # Tendencia no lineal
       geom_smooth(se = FALSE, method = "gam", formula =  y ~ splines::bs(x, 3)) +
       # Tendencia lineal
       geom_smooth(se = FALSE, method = "lm", color = "firebrick") +
       facet_grid(. ~ title) +
       theme_bw() +
       theme(strip.text = element_text(colour = "black", size = 8, face = 2),
             axis.title = element_blank())
  return(p)
}

variables_continuas <- c("metros_habitables", "universitarios", "dormitorios",
                         "banyos", "habitaciones", "log_antiguedad", "log_metros_totales",
                         "log_precio_terreno")

plots <- map(
            .x = variables_continuas,
            .f = custom_corr_plot,
            variable2 = "precio",
            df = datos
         )

ggarrange(plotlist = plots, ncol = 3, nrow = 3) %>%
  annotate_figure(
    top = text_grob("Correlación con precio", face = "bold", size = 16,
                    x = 0.20)
  )

datos <- datos %>%
         select(-c(log_antiguedad, log_metros_totales, log_precio_terreno))

Correlación variables continuas

Algunos modelos (LM, GLM, …) se ven perjudicados si incorporan predictores altamente correlacionados. Por esta razón, es conveniente estudiar el grado de correlación entre las variables disponibles.

plot_correlation(
  data = datos,
  type = "continuous",
  title = "Matriz de correlación variables continuas",
  theme_config = list(legend.position = "none",
                      plot.title = element_text(size = 12, face = "bold"),
                      axis.title = element_blank(),
                      axis.text.x = element_text(angle = -45, hjust = +0.1)
                     )
)

GGally::ggscatmat(
  data = datos %>% select_if(is.numeric),
  alpha = 0.1) +
theme_bw() +
labs(title = "Correlación por pares") +
theme(
  plot.title = element_text(size = 16, face = "bold"),
  axis.text = element_blank(),
  strip.text = element_text(colour = "black", size = 5, face = 1)
)

Variables cualitativas

plot_bar(
  datos,
  ncol    = 3,
  title   = "Número de observaciones por grupo",
  ggtheme = theme_bw(),
  theme_config = list(
                   plot.title = element_text(size = 16, face = "bold"),
                   strip.text = element_text(colour = "black", size = 10, face = 2),
                   legend.position = "none"
                  )
)

Si alguno de los niveles de una variable cualitativa tiene muy pocas observaciones en comparación a los otros niveles, puede ocurrir que, durante la validación cruzada o bootstrapping, algunas particiones no contengan ninguna observación de dicha clase (varianza cero), lo que puede dar lugar a errores. En estos casos, suele ser conveniente:

• Eliminar las observaciones del grupo minoritario si es una variable multiclase.
• Eliminar la variable si solo tiene dos niveles.
• Agrupar los niveles minoritarios en un único grupo.
• Asegurar que, en la creación de las particiones, todos los grupos estén representados en cada una de ellas.

Para este caso, hay que tener precaución con la variable chimenea. Se unifican los niveles de 2, 3 y 4 en un nuevo nivel llamado “2_mas”.

table(datos$chimenea)

## 
##   0   1   2   3   4 
## 740 942  42   2   2

datos <- datos %>%
         mutate(
           chimenea = recode_factor(
                        chimenea,
                        `2` = "2_mas",
                        `3` = "2_mas",
                        `4` = "2_mas"
                      )
         )

table(datos$chimenea)

## 
## 2_mas     0     1 
##    46   740   942

custom_box_plot <- function(variable1, variable2, df, alpha=0.3){
  p <- df %>%
       mutate(
         # Truco para que se ponga el título estilo facet
        title = paste(toupper(variable2), "vs", toupper(variable1))
       ) %>%
       ggplot(aes(x = !!sym(variable1), y = !!sym(variable2))) + 
       geom_violin(alpha = alpha) +
       geom_boxplot(width = 0.1, outlier.shape = NA) +
       facet_grid(. ~ title) +
       theme_bw() +
       theme(strip.text = element_text(colour = "black", size = 8, face = 2),
             axis.text.x = element_text(size = 7),
             axis.title = element_blank())
  return(p)
}

variables_cualitativas <- c("chimenea", "calefaccion", "consumo_calefacion", "desague",
                            "vistas_lago", "nueva_construccion", "aire_acondicionado")

plots <- map(
            .x = variables_cualitativas,
            .f = custom_box_plot,
            variable2 = "precio",
            df = datos
         )

ggarrange(plotlist = plots, ncol = 3, nrow = 3) %>%
  annotate_figure(
    top = text_grob("Correlación con precio", face = "bold", size = 16,
                    x = 0.16)
  )

Task

El objeto task encapsula los datos, información sobre el tipo de variables, número de observaciones, el tipo de modelo (regresión, clasificación…) e identifica el rol de cada variable (predictores, respuesta, estratificación, la clase positiva para calcular métricas de clasificación binaria, pesos…).

task_datos <- TaskRegr$new(
                id      = "task_datos",
                backend = datos,
                target  = "precio"
              )
task_datos

## <TaskRegr:task_datos> (1728 x 16)
## * Target: precio
## * Properties: -
## * Features (15):
##   - fct (7): aire_acondicionado, calefaccion, chimenea,
##     consumo_calefacion, desague, nueva_construccion, vistas_lago
##   - int (6): antiguedad, dormitorios, habitaciones, metros_habitables,
##     precio_terreno, universitarios
##   - dbl (2): banyos, metros_totales

mlr3 incorpora toda una serie de atributos y métodos con los que extraer o modificar información de un objeto task: $col_info, $col_roles, $missings(), $levels(), $select(), filter$(). Es importante tener en cuenta que, $select() y filter$(), modifican el objeto task.

# Información de todas las variables.
task_datos$col_info

# Role de todas las variables.
task_datos$col_roles

## $feature
##  [1] "aire_acondicionado" "antiguedad"         "banyos"            
##  [4] "calefaccion"        "chimenea"           "consumo_calefacion"
##  [7] "desague"            "dormitorios"        "habitaciones"      
## [10] "metros_habitables"  "metros_totales"     "nueva_construccion"
## [13] "precio_terreno"     "universitarios"     "vistas_lago"       
## 
## $target
## [1] "precio"
## 
## $name
## character(0)
## 
## $order
## character(0)
## 
## $stratum
## character(0)
## 
## $group
## character(0)
## 
## $weight
## character(0)

# Número de valores NA por variable.
task_datos$missings()

##             precio aire_acondicionado         antiguedad             banyos 
##                  0                  0                  0                  0 
##        calefaccion           chimenea consumo_calefacion            desague 
##                  0                  0                  0                  0 
##        dormitorios       habitaciones  metros_habitables     metros_totales 
##                  0                  0                  0                  0 
## nueva_construccion     precio_terreno     universitarios        vistas_lago 
##                  0                  0                  0                  0

# Niveles de cada variable de tipo factor.
task_datos$levels()

## $aire_acondicionado
## [1] "Yes" "No" 
## 
## $calefaccion
## [1] "hot air"         "hot water/steam" "electric"       
## 
## $chimenea
## [1] "2_mas" "0"     "1"    
## 
## $consumo_calefacion
## [1] "gas"      "electric" "oil"     
## 
## $desague
## [1] "septic"            "public/commercial" "none"             
## 
## $nueva_construccion
## [1] "Yes" "No" 
## 
## $vistas_lago
## [1] "Yes" "No"

# Mostrar datos por número de fila
task_datos$data(rows = 1:3)

División train-test

Evaluar la capacidad predictiva de un modelo consiste en comprobar cómo de próximas son sus predicciones a los verdaderos valores de la variable respuesta. Para poder cuantificarlo de forma correcta, se necesita disponer de un conjunto de observaciones, de las que se conozca la variable respuesta, pero que el modelo no haya “visto”, es decir, que no hayan participado en su ajuste. Con esta finalidad, se dividen los datos disponibles en un conjunto de entrenamiento y un conjunto de test. El tamaño adecuado de las particiones depende en gran medida de la cantidad de datos disponibles y la seguridad que se necesite en la estimación del error, 80%-20% suele dar buenos resultados. El reparto debe hacerse de forma aleatoria o aleatoria-estratificada. Evaluación de modelos.

El paquete mlr3 dispone del objeto resampling = rsmp("holdout") para repartir en dos partes (train-test) los datos contenidos en un task. Con $train_set() y $test_set() se accede a los índices de las observaciones que caen en cada parte.

set.seed(1234)
rsmp_holdout <- rsmp("holdout", ratio = 0.8)
rsmp_holdout

## <ResamplingHoldout> with 1 iterations
## * Instantiated: FALSE
## * Parameters: ratio=0.8

# Para que el reparto sea estratificado, es necesario indicar en el task qué
# columna o columnas se emplea como stratum. Solo para variables cualitativas.
# task_datos$col_roles$stratum <- "precio"

rsmp_holdout$instantiate(task = task_datos)

# Índices train
id_train <- rsmp_holdout$train_set(i = 1)
# Índices test
id_test <- rsmp_holdout$test_set(i = 1)

# Se crean dos nuevas task, una con los datos de train y otra con los de test.
# Dado que se va a aplicar un filtrado, y para no alterar task_datos, se emplea
# antes del filtro el método $clone() para hacer una copia.
task_train <- task_datos$clone()$filter(id_train)
task_test  <- task_datos$clone()$filter(id_test)

Es importante verificar que la distribución de la variable respuesta es similar en el conjunto de entrenamiento y en el de test.

summary(task_train$data()[["precio"]])

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   10300  145000  189000  211594  256838  775000

summary(task_test$data()[["precio"]])

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    5000  142250  194000  213455  260000  760000

Cuando se emplea un resampling, si existe una o varias columnas con el role stratum, se realiza un reparto estratificado, garantizando así una distribución aproximada (por el momento, esto solo es válido para variables cualitativas). El reparto estratificado es muy importante cuando hay grupos minoritarios en alguna variable. Por ejemplo, en un set de datos con 100 observaciones, un predictor binario que tenga 90 observaciones de un grupo y solo 10 de otro, tiene un alto riesgo de que, en alguna de las particiones, el grupo minoritario no tenga representantes. Si esto ocurre en el conjunto de entrenamiento, algunos algoritmos darán error al aplicarlos al conjunto de test, ya que no entenderán el valor que se les está pasando. Este problema puede evitarse eliminando variables con varianza próxima a cero (ver más adelante).

Preprocesado-mlr3pipelines

El preprocesado engloba todas aquellas transformaciones realizadas sobre los datos con el objetivo que puedan ser interpretados por el algoritmo de machine learning lo más eficientemente posible. Todo preprocesado de datos debe aprenderse con las observaciones de entrenamiento y luego aplicarse al conjunto de entrenamiento y al de test. Esto es muy importante para no violar la condición de que ninguna información procedente de las observaciones de test participe o influya en el ajuste del modelo. Este principio debe aplicarse también si se emplea validación cruzada (ver más adelante). En tal caso, el preporcesado debe realizarse dentro de cada iteración de validación, para que las estimaciones que se hacen con cada partición de validación no contengan información del resto de particiones. A continuación, se resumen algunos de los pasos de preprocesado que más se suelen necesitar.

Imputación de valores ausentes

La gran mayoría de algoritmos no aceptan observaciones incompletas, por lo que, cuando el set de datos contiene valores ausentes, se puede:

• Eliminar aquellas observaciones que estén incompletas.

• Eliminar aquellas variables que contengan valores ausentes.

• Tratar de estimar los valores ausentes empleando el resto de información disponible (imputación).

Las primeras dos opciones, aunque sencillas, suponen perder información. La eliminación de observaciones solo puede aplicarse cuando se dispone de muchas y el porcentaje de registros incompletos es muy bajo. En el caso de eliminar variables, el impacto dependerá de cuánta información aporten dichas variables al modelo. Cuando se emplea imputación, es muy importante tener en cuenta el riesgo que se corre al introducir valores en predictores que tengan mucha influencia en el modelo. Supóngase un estudio médico en el que, cuando uno de los predictores es positivo, el modelo predice casi siempre que el paciente está sano. Para un paciente cuyo valor de este predictor se desconoce, el riesgo de que la imputación sea errónea es muy alto, por lo que es preferible obtener una predicción basada únicamente en la información disponible. Esta es otra muestra de la importancia que tiene que el analista conozca el problema al que se enfrenta y pueda así tomar la mejor decisión.

Exclusión de variables con varianza próxima a cero

No se deben incluir en el modelo predictores que contengan un único valor (cero-varianza) ya que no aportan información. Tampoco es conveniente incluir predictores que tengan una varianza próxima a cero, es decir, predictores que toman solo unos pocos valores, de los cuales, algunos aparecen con muy poca frecuencia. El problema con estos últimos es que pueden convertirse en predictores con varianza cero cuando se dividen las observaciones por validación cruzada o bootstrap.

Si bien la eliminación de predictores no informativos podría considerarse un paso propio del proceso de selección de predictores, dado que consiste en un filtrado por varianza, tiene que realizarse antes de estandarizar los datos, ya que después, todos los predictores tienen varianza 1.

Estandarización y escalado de variables numéricas

Cuando los predictores son numéricos, la escala en la que se miden, así como la magnitud de su varianza pueden influir en gran medida en el modelo. Muchos algoritmos de machine learning (SVM, redes neuronales, lasso…) son sensibles a esto, de forma que, si no se igualan de alguna forma los predictores, aquellos que se midan en una escala mayor o que tengan más varianza dominarán el modelo aunque no sean los que más relación tienen con la variable respuesta. Existen principalmente 2 estrategias para evitarlo:

• Centrado: consiste en restarle a cada valor la media del predictor al que pertenece. Si los datos están almacenados en un dataframe, el centrado se consigue restándole a cada valor la media de la columna en la que se encuentra. Como resultado de esta transformación, todos los predictores pasan a tener una media de cero, es decir, los valores se centran en torno al origen.

• Normalización (estandarización): consiste en transformar los datos de forma que todos los predictores estén aproximadamente en la misma escala. Hay dos formas de lograrlo:

  • Normalización Z-score: dividir cada predictor entre su desviación típica después de haber sido centrado, de esta forma, los datos pasan a tener una distribución normal. z=x−μσ
  • Estandarización max-min: transformar los datos de forma que estén dentro del rango [0, 1]. Xnorm=X−XminXmax−Xmin

Nunca se debe estandarizar las variables después de ser binarizadas (ver a continuación).

Binarización de las variables cualitativas

La binarización consiste en crear nuevas variables dummy con cada uno de los niveles de las variables cualitativas. A este proceso también se le conoce como one hot encoding. Por ejemplo, una variable llamada color que contenga los niveles rojo, verde y azul, se convertirá en tres nuevas variables (color_rojo, color_verde, color_azul), todas con el valor 0 excepto la que coincide con la observación, que toma el valor 1. Muchos modelos de machine learning y aprendizaje estadístico no pueden interpretar variables de tipo character o factor, por lo que es necesario hacer esta transformación.

La forma de preprocesar los datos dentro del ecosistema mlr3 es mediante los paquetes mlr3pipelines y mlrfilters. Estos paquetes permite encadenar operadores PipeOp y Filters, cada uno aplica una transformación concreta, y ejecutarlos todos de forma secuencial. Una vez que han sido definidos, pueden ser entrenados con $train() para aprender las tranformaciones y luego aplicados a nuevos conjuntos de datos con $pred(). Son muchos los operadores disponibles y se continúan añadiendo más.

# Operadores PipeOp disponibles
mlr_pipeops$keys()

##  [1] "boxcox"                "branch"                "chunk"                
##  [4] "classbalancing"        "classifavg"            "classweights"         
##  [7] "colapply"              "collapsefactors"       "colroles"             
## [10] "copy"                  "datefeatures"          "encode"               
## [13] "encodeimpact"          "encodelmer"            "featureunion"         
## [16] "filter"                "fixfactors"            "histbin"              
## [19] "ica"                   "imputeconstant"        "imputehist"           
## [22] "imputelearner"         "imputemean"            "imputemedian"         
## [25] "imputemode"            "imputeoor"             "imputesample"         
## [28] "kernelpca"             "learner"               "learner_cv"           
## [31] "missind"               "modelmatrix"           "multiplicityexply"    
## [34] "multiplicityimply"     "mutate"                "nmf"                  
## [37] "nop"                   "ovrsplit"              "ovrunite"             
## [40] "pca"                   "proxy"                 "quantilebin"          
## [43] "randomprojection"      "randomresponse"        "regravg"              
## [46] "removeconstants"       "renamecolumns"         "replicate"            
## [49] "scale"                 "scalemaxabs"           "scalerange"           
## [52] "select"                "smote"                 "spatialsign"          
## [55] "subsample"             "targetinvert"          "targetmutate"         
## [58] "targettrafoscalerange" "textvectorizer"        "threshold"            
## [61] "tunethreshold"         "unbranch"              "vtreat"               
## [64] "yeojohnson"

# Filtros disponibles
mlr_filters$keys()

##  [1] "anova"            "auc"              "carscore"         "cmim"            
##  [5] "correlation"      "disr"             "find_correlation" "importance"      
##  [9] "information_gain" "jmi"              "jmim"             "kruskal_test"    
## [13] "mim"              "mrmr"             "njmim"            "performance"     
## [17] "permutation"      "variance"

Para este ejemplo se normalizan los datos y se hace una codificación one-hot de las variables cualitativas.

preproces_pipeline <- po("scale", param_vals = list(center = TRUE, scale = TRUE)) %>>%
                      po("encode",param_vals = list(method = "one-hot"))

preproces_pipeline$train(task_train)

## $encode.output
## <TaskRegr:task_datos> (1382 x 27)
## * Target: precio
## * Properties: -
## * Features (26):
##   - dbl (26): aire_acondicionado.No, aire_acondicionado.Yes,
##     antiguedad, banyos, calefaccion.electric, calefaccion.hot.air,
##     calefaccion.hot.water.steam, chimenea.0, chimenea.1,
##     chimenea.2_mas, consumo_calefacion.electric,
##     consumo_calefacion.gas, consumo_calefacion.oil, desague.none,
##     desague.public.commercial, desague.septic, dormitorios,
##     habitaciones, metros_habitables, metros_totales,
##     nueva_construccion.No, nueva_construccion.Yes, precio_terreno,
##     universitarios, vistas_lago.No, vistas_lago.Yes

Los resultados (output) de un operador mlr3pipelines siempre son una lista, en la que el primer elemento ([[1]]) contiene los datos transformados.

task_train_prerp <- preproces_pipeline$predict(task_train)[[1]]$clone()
task_test_prerp  <- preproces_pipeline$predict(task_test)[[1]]$clone()
task_train_prerp$data() %>% head()

Modelo-Learner

El siguiente paso tras definir los datos de entrenamiento, es seleccionar el algoritmo que se va a emplear. En mlr3, esto se hace mediante la creación de un objeto learner. En concreto, este objeto almacena el nombre del algoritmo, sus parámetros e hiperparámetros y el tipo de resultado que se muestra al predecir nuevas observaciones. El paquete mlr3learners contiene solo los principales algoritmos empleados para clasificación y regresión (listado).

# Listado de learners disponibles en mlr3learners
mlr_learners$keys()

##  [1] "classif.cv_glmnet"   "classif.debug"       "classif.featureless"
##  [4] "classif.glmnet"      "classif.kknn"        "classif.lda"        
##  [7] "classif.log_reg"     "classif.multinom"    "classif.naive_bayes"
## [10] "classif.qda"         "classif.ranger"      "classif.rpart"      
## [13] "classif.svm"         "classif.xgboost"     "regr.cv_glmnet"     
## [16] "regr.featureless"    "regr.glmnet"         "regr.kknn"          
## [19] "regr.km"             "regr.lm"             "regr.ranger"        
## [22] "regr.rpart"          "regr.svm"            "regr.xgboost"       
## [25] "surv.cv_glmnet"      "surv.glmnet"         "surv.ranger"        
## [28] "surv.xgboost"

Pero existen toda una serie de extensiones con otros algoritmos (listado). Por ejemplo, los modelos SVM del paquete kernlab.

#install.packages("mlr3learners.kernlab")
# Si no funciona probar:
#install.packages("remotes")
#remotes::install_github("mlr3learners/mlr3learners.kernlab")
library(mlr3learners.kernlab)
mlr_learners$keys()

##  [1] "classif.cv_glmnet"   "classif.debug"       "classif.featureless"
##  [4] "classif.glmnet"      "classif.kknn"        "classif.ksvm"       
##  [7] "classif.lda"         "classif.log_reg"     "classif.multinom"   
## [10] "classif.naive_bayes" "classif.qda"         "classif.ranger"     
## [13] "classif.rpart"       "classif.svm"         "classif.xgboost"    
## [16] "regr.cv_glmnet"      "regr.featureless"    "regr.glmnet"        
## [19] "regr.kknn"           "regr.km"             "regr.ksvm"          
## [22] "regr.lm"             "regr.ranger"         "regr.rpart"         
## [25] "regr.svm"            "regr.xgboost"        "surv.cv_glmnet"     
## [28] "surv.glmnet"         "surv.ranger"         "surv.xgboost"

Para crear un nuevo learner se puede emplear a la función lrn() o al método mlr_learners$get().

# Se crea un learner con el algoritmo de SVM
learner_svm <- lrn("regr.svm")
# learner_svm <- mlr_learners$get(key = "regr.svm")

Imprimiendo el learner se muestra información general sobre el algoritmo:

• feature_types: tipo de variables que el algoritmo es capaz de aceptar (importante para identificar aquellos que no pueden tratar con factores.
• packages: el paquete al que pertenece el algoritmo.
• properties: propiedades adicionales del algoritmo. Por ejemplo, “missings” significa que puede tratar con valores ausentes, y “importance” significa que el algoritmo es capaz de calcular la importancia de los predictores.
• predict_types: tipo de predicción.
  • Clasificación: “response” devuelve únicamente la clase con mayor probabilidad y “prob” devuelve además la probabilidad de cada clase.
  • Regresión: “response” devuelve el valor predicho y “se” devuelve además el error estándar de la predicción.
  • Supervivencia (Pendiente de implementar): “response” devuelve una medida de “risk” y “prob” devuelve la probabilidad en función del tiempo.
  • Clustering (Pendiente de implementar): “response” devuelve el id del cluster al que se ha asignado cada observación y “prob” devuelve la probabilidad de asignación a cada cluster (solo para fuzzy clustering).

# Información general del learner
print(learner_svm)

## <LearnerRegrSVM:regr.svm>
## * Model: -
## * Parameters: list()
## * Packages: e1071
## * Predict Type: response
## * Feature types: logical, integer, numeric
## * Properties: -

Para conocer en detalle sus parámetros (hiperparámetros), el rango de valores válido y los valores por defecto, se emplea el atributo $param_set.

# Tabla con toda la información
knitr::kable(as.data.table(learner_svm$param_set), format = "markdown")

# Valores de los parámetros cualitativos
learner_svm$param_set$levels %>% purrr::discard(is.null)

## $type
## [1] "eps-regression" "nu-regression" 
## 
## $kernel
## [1] "linear"     "polynomial" "radial"     "sigmoid"   
## 
## $shrinking
## [1]  TRUE FALSE
## 
## $fitted
## [1]  TRUE FALSE

# Valores por defecto
learner_svm$param_set$default %>% as.data.frame()

Los parámetros (hiperparámetros) del learner pueden configurarse en el momento de su creación o modificarlos posteriormente.

learner_svm <- lrn("regr.svm", type = "eps-regression", kernel = "linear", cost = 0.5)
learner_svm

## <LearnerRegrSVM:regr.svm>
## * Model: -
## * Parameters: type=eps-regression, kernel=linear, cost=0.5
## * Packages: e1071
## * Predict Type: response
## * Feature types: logical, integer, numeric
## * Properties: -

learner_svm$param_set$values <- list(type = "eps-regression", kernel = "linear", cost = 0.5)

Además de los parámetros, otra característica importante que se puede modificar es predict_type, que determina tipo de resultado devuelto en las predicciones.

learner_svm$predict_type <- "response"

Entrenamiento

Una vez creado el task y el learner, se puede entrenar el modelo. Se procede a ajustar un modelo máquinas vector soporte lineal (SVM) que prediga el precio de la vivienda en función de todos los predictores. A excepción del tipo de SVM y de kernel, todos los demás parámetros se dejan por defecto.

# Definición del learner
learner_svm <- lrn("regr.svm", type = "eps-regression", kernel = "linear")

El entrenamiento se realiza únicamente con las observaciones que forman parte de la partición de train, que están almacenadas en task_train.

# Entrenamiento del learner
learner_svm$train(task = task_train_prerp)

El modelo generado se almacena en dentro del learner en el elemento $model. El modelo sigue manteniendo la misma clase que en el paquete original, por lo que puede manipularse acorde a todos sus métodos y funciones originales. Esto es importante a la hora de decidir qué hacer con el modelo final. Por ejemplo, si este modelo se quiere poner en producción, no sería necesario instalar todos los paquetes de mlr3, únicamente el paquete e1071,

learner_svm$model

## 
## Call:
## svm.default(x = data, y = task$truth(), type = "eps-regression", 
##     kernel = "linear")
## 
## 
## Parameters:
##    SVM-Type:  eps-regression 
##  SVM-Kernel:  linear 
##        cost:  1 
##       gamma:  0.03846154 
##     epsilon:  0.1 
## 
## 
## Number of Support Vectors:  1161

Validación

La finalidad última de un modelo es predecir la variable respuesta en observaciones futuras o en observaciones que el modelo no ha “visto” antes. El error mostrado por defecto tras entrenar un modelo suele ser el error de entrenamiento, el error que comete el modelo al predecir las observaciones que ya ha “visto”. Si estos errores son útiles para entender cómo está aprendiendo el modelo (estudio de residuos), no es una estimación realista de cómo se comporta el modelo ante nuevas observaciones (el error de entrenamiento suele ser demasiado optimista). Para conseguir una estimación más certera, y antes de recurrir al conjunto de test, se pueden emplear estrategias de validación basadas en resampling. El paquete mlr3 incorpora las siguientes estrategias: Bootstrap (mlr_resamplings_bootstrap), V-Fold Cross-Validation (mlr_resamplings_cv), Repeated V-Fold Cross-Validation (mlr_resamplings_repeated_cv) y Monte Carlo Cross-Validation (mlr_resamplings_subsampling) Anexo 1. Cada uno funciona internamente de forma distinta, pero todos ellos se basan en la idea: ajustar y evaluar el modelo de forma repetida, empleando cada vez distintos subconjuntos creados a partir de los datos de entrenamiento y obteniendo en cada repetición una estimación del error. El promedio de todas las estimaciones tiende a converger en el valor real del error de test.

# Métodos de resampling disponibles
as.data.table(mlr_resamplings)

Se ajusta de nuevo el modelo, esta vez con validación cruzada repetida para estimar su error.

En primer lugar, se crea un objeto Resampling en el que se define el tipo de validación y el número de repeticiones.

# Validación cruzada con 10 particiones
set.seed(123)
resampling_cv <- rsmp("cv", folds = 10)
resampling_cv

## <ResamplingCV> with 10 iterations
## * Instantiated: FALSE
## * Parameters: folds=10

Con la función resample() se combinan un task, un learner y la estrategia de Resampling. De esta forma el modelo se ajusta y evalúa con cada una de las particiones de forma automática, calculando y almacenando en cada iteración la métrica de interés.

# Métricas disponibles
mlr_measures

## <DictionaryMeasure> with 54 stored values
## Keys: classif.acc, classif.auc, classif.bacc, classif.bbrier,
##   classif.ce, classif.costs, classif.dor, classif.fbeta, classif.fdr,
##   classif.fn, classif.fnr, classif.fomr, classif.fp, classif.fpr,
##   classif.logloss, classif.mbrier, classif.mcc, classif.npv,
##   classif.ppv, classif.prauc, classif.precision, classif.recall,
##   classif.sensitivity, classif.specificity, classif.tn, classif.tnr,
##   classif.tp, classif.tpr, debug, oob_error, regr.bias, regr.ktau,
##   regr.mae, regr.mape, regr.maxae, regr.medae, regr.medse, regr.mse,
##   regr.msle, regr.pbias, regr.rae, regr.rmse, regr.rmsle, regr.rrse,
##   regr.rse, regr.rsq, regr.sae, regr.smape, regr.srho, regr.sse,
##   selected_features, time_both, time_predict, time_train

# Se define la métrica de evaluación
metrica <- msr("regr.rmse")
  
# Se aplica el resampling
resultados_resampling <- resample(
                           task         = task_train_prerp,
                           learner      = learner_svm,
                           resampling   = resampling_cv,
                           store_models = FALSE
                         )

## INFO  [21:20:15.900] Applying learner 'regr.svm' on task 'task_datos' (iter 6/10) 
## INFO  [21:20:16.357] Applying learner 'regr.svm' on task 'task_datos' (iter 4/10) 
## INFO  [21:20:16.816] Applying learner 'regr.svm' on task 'task_datos' (iter 5/10) 
## INFO  [21:20:17.253] Applying learner 'regr.svm' on task 'task_datos' (iter 8/10) 
## INFO  [21:20:17.724] Applying learner 'regr.svm' on task 'task_datos' (iter 7/10) 
## INFO  [21:20:18.168] Applying learner 'regr.svm' on task 'task_datos' (iter 10/10) 
## INFO  [21:20:18.632] Applying learner 'regr.svm' on task 'task_datos' (iter 2/10) 
## INFO  [21:20:19.082] Applying learner 'regr.svm' on task 'task_datos' (iter 9/10) 
## INFO  [21:20:19.518] Applying learner 'regr.svm' on task 'task_datos' (iter 3/10) 
## INFO  [21:20:19.981] Applying learner 'regr.svm' on task 'task_datos' (iter 1/10)

El objeto ResampleResult creado almacena toda la información sobre cada iteración.

resultados_resampling

## <ResampleResult> of 10 iterations
## * Task: task_datos
## * Learner: regr.svm
## * Warnings: 0 in 0 iterations
## * Errors: 0 in 0 iterations

# Métricas promedio de todas las particiones
resultados_resampling$aggregate(measures = msr("regr.rmse"))

## regr.rmse 
##  58186.11

# Métricas individuales de cada una de las particiones
resultados_resampling$score(msr("regr.rmse")) %>%
  as.data.table() 

# Distribución del error de validación 
metricas <- resultados_resampling$score(msr("regr.rmse")) %>%
            as.data.table() %>% 
            select(iteration, regr.rmse)

p1 <- ggplot(data = metricas, aes(x = regr.rmse)) +
      geom_density(alpha = 0.5, fill = "gray50") +
      geom_rug() + 
      geom_vline(xintercept = mean(metricas$regr.rmse),
                 linetype = "dashed") +
      theme_bw() 

p2 <- ggplot(data = metricas, aes(x = 1, y = regr.rmse)) +
      geom_boxplot(outlier.shape = NA, alpha = 0.5, fill = "gray50") +
      geom_jitter(width = 0.05) +
      coord_flip() +
      labs(x = "") +
      theme_bw() +
      theme(axis.text.x = element_blank(), axis.ticks.x = element_blank())
  
ggarrange(p1, p2, nrow = 2, align = "v") %>%
annotate_figure(
    top = text_grob("Accuracy obtenido en la validación cruzada",
                    size = 15)
)

#autoplot(resultados_resampling)

Almacenar las predicciones de cada partición es útil para poder evaluar los residuos del modelo y diagnosticar su comportamiento. Si se emplea validación cruzada repetida o bootstrapping, una misma observación puede formar parte de la partición de validación varias veces.

# Predicciones individuales de cada observación. Si 
as.data.table(resultados_resampling$prediction()) %>% head(3)

predicciones_validacion <- resultados_resampling$prediction() %>%
                           as.data.table() %>%
                           mutate(
                             residuo = response - truth
                           )
  
p1 <- ggplot(data = predicciones_validacion, aes(x = truth, y = response)) +
      geom_point(alpha = 0.3) +
      geom_abline(slope = 1, intercept = 0, color = "firebrick") +
      labs(title = "Valor predicho vs valor real") +
      theme_bw()


p2 <- ggplot(data = predicciones_validacion, aes(x = row_id, y = residuo)) +
      geom_point(alpha = 0.3) +
      geom_hline(yintercept =  0, color = "firebrick") +
      labs(title = "Residuos del modelo") +
      theme_bw()

p3 <- ggplot(data = predicciones_validacion, aes(x = residuo)) +
      geom_density() + 
      labs(title = "Distribución residuos del modelo") +
      theme_bw()

p4 <- ggplot(data = predicciones_validacion, aes(sample = residuo)) +
      geom_qq() +
      geom_qq_line(color = "firebrick") +
      labs(title = "Q-Q residuos del modelo") +
      theme_bw()

ggarrange(plotlist = list(p1, p2, p3, p4)) %>%
annotate_figure(
  top = text_grob("Distribución residuos", size = 15, face = "bold")
)

Como se comentó anteriormente, cuando se realizan validaciones por resampling, el preprocesado debe ocurrir dentro de cada iteración. Para seguir este principio, se emplea el objeto GraphLearner, que permite combinar un mlr3pipeline de preprocesado con un learner. Una vez creado el objeto GraphLearner, este puede emplearse para como si de un learner se tratase (entrenamiento, predicción, validación…).

# Pipeline de preprocesado
preproces_pipeline <- po("scale", param_vals = list(center = TRUE, scale = TRUE)) %>>%
                      po("encode", param_vals = list(method = "one-hot"))

# Learner
# learner_svm <- lrn("regr.svm", type = "eps-regression", kernel = "linear", cost = 0.5)

# Unión del pipeline con el learner
graph <- preproces_pipeline %>>% learner_svm
graph$plot()

graph_learner <- GraphLearner$new(graph)

Es importante tener en cuenta que, el objeto GraphLearner, incorpora las acciones de preprocesado. A la hora de entrenarlo o estimar su error mediante resample(), se tienen que utilizar los datos sin haber sido preprocesado.

# Entrenamiento del GraphLearner
graph_learner$train(
    task = task_train
)

# Validación del GraphLearner
set.seed(123)
resampling_cv <- rsmp("cv", folds = 10)
metrica <- msr("regr.rmse")
lgr::get_logger("mlr3")$set_threshold("warn")
resultados_resampling <- resample(
                           task         = task_train,
                           learner      = graph_learner,
                           resampling   = resampling_cv,
                           store_models = FALSE
                         )

A efectos prácticos, cuando se aplican métodos de resampling hay que tener en cuenta dos cosas: el coste computacional que implica ajustar múltiples veces un modelo, cada vez con un subconjunto de datos distinto, y la reproducibilidad en la creación de las particiones. resample() permite paralelizar el proceso si se registra un backed paralelo.

future::plan(future::multiprocess())
lgr::get_logger("mlr3")$set_threshold("warn")
resultados_resampling <- resample(
                           task         = task_train,
                           learner      = graph_learner,
                           resampling   = resampling_cv,
                           store_models = FALSE
                         )

resultados_resampling$aggregate(measures = msr("regr.rmse"))

## regr.rmse 
##  58918.14

El rmse promedio estimado mediante validación cruzada es de 5.891814210^{4}. Este valor será contrastado más adelante cuando se calcule el rmse del modelo con el conjunto de test.

Predicción

Una vez que el modelo ha sido entrenado, bien empleando un learner o un GraphLearner, con el método $predict() se pueden predecir nuevas observaciones. El objeto devuelto almacena toda la información relativa a la predicción junto con el valor real (en caso de estar disponible).

# Si se emplea un learner, los datos de test tienen que haber sido preprocesados
# con antelación
predicciones <- learner_svm$predict(
                  task = task_test_prerp
                )
as.data.table(predicciones) %>% head(3)

# Si se emplea un GraphLearner, los datos de test no tienen que haber sido
# preprocesados con antelación
predicciones <- graph_learner$predict(
                  task = task_test
                )
as.data.table(predicciones) %>% head(3)

Error de test

Aunque mediante los métodos de validación (CV, bootstrapping…) se consiguen buenas estimaciones del error que tiene un modelo al predecir nuevas observaciones, la mejor forma de evaluar un modelo final es prediciendo un conjunto test, es decir, un conjunto de observaciones que se ha mantenido al margen del proceso de entrenamiento y optimización. Dependiendo del problema en cuestión, pueden ser interesantes unas métricasAnexo 2 u otras. El paquete mlr3 incorpora una variedad considerable de métricas para evaluar la calidad de las predicciones y las hace accesibles a través del objeto Measure combinado con un objeto Prediction.

# Métricas disponibles
mlr_measures$keys()

##  [1] "classif.acc"         "classif.auc"         "classif.bacc"       
##  [4] "classif.bbrier"      "classif.ce"          "classif.costs"      
##  [7] "classif.dor"         "classif.fbeta"       "classif.fdr"        
## [10] "classif.fn"          "classif.fnr"         "classif.fomr"       
## [13] "classif.fp"          "classif.fpr"         "classif.logloss"    
## [16] "classif.mbrier"      "classif.mcc"         "classif.npv"        
## [19] "classif.ppv"         "classif.prauc"       "classif.precision"  
## [22] "classif.recall"      "classif.sensitivity" "classif.specificity"
## [25] "classif.tn"          "classif.tnr"         "classif.tp"         
## [28] "classif.tpr"         "debug"               "oob_error"          
## [31] "regr.bias"           "regr.ktau"           "regr.mae"           
## [34] "regr.mape"           "regr.maxae"          "regr.medae"         
## [37] "regr.medse"          "regr.mse"            "regr.msle"          
## [40] "regr.pbias"          "regr.rae"            "regr.rmse"          
## [43] "regr.rmsle"          "regr.rrse"           "regr.rse"           
## [46] "regr.rsq"            "regr.sae"            "regr.smape"         
## [49] "regr.srho"           "regr.sse"            "selected_features"  
## [52] "time_both"           "time_predict"        "time_train"

Se procede a calcular el rmse para el conjunto de test.

# rmse de test
predicciones$score(measures = msr("regr.rmse"))

## regr.rmse 
##  58863.05

En el apartado de validación, se estimó, mediante validación cruzada repetida, que el rmse del modelo era de 58918.14, un valor muy próximo al obtenido con el conjunto de test (5.886305410^{4}).

Hiperparámetros (tuning)

Muchos modelos, entre ellos las máquinas vector soporte (SVM), contienen parámetros que no pueden aprenderse a partir de los datos de entrenamiento y, por lo tanto, deben de ser establecidos por el analista. A estos se les conoce como hiperparámetros. Los resultados de un modelo pueden depender en gran medida del valor que tomen sus hiperparámetros, sin embargo, no se puede conocer de antemano cuál es el adecuado. Aunque con la práctica, los especialistas en machine learning ganan intuición sobre qué valores pueden funcionar mejor en cada problema, no hay reglas fijas. La forma más común de encontrar los valores óptimos es probando diferentes posibilidades.

1. Escoger un conjunto de valores para el o los hiperparámetros.

2. Para cada valor (combinación de valores si hay más de un hiperparámetro), entrenar el modelo y estimar su error mediante un método de validación Anexo 2.

3. Finalmente, ajustar de nuevo el modelo, esta vez con todos los datos de entrenamiento y con los mejores hiperparámetros encontrados.

El modelo learner_svm empleado hasta ahora tienen un hiperparámetro llamado coste. Este hiperparámetro controla la penalización que se aplica a las clasificaciones erróneas cuando se entrena el modelo. Si su valor es alto, el modelo resultante es más flexible y se ajusta mejor a las observaciones de entrenamiento pero tiene más riesgo de sufrir overfitting. mlr3 permite explorar diferentes valores de hiperparámetros mediante las extensiones mlr3tuning y paradox.

Se vuelve a ajustar un modelo learner_svm con diferentes valores de cost, y empleando validación cruzada repetida para identificar con cuáles de ellos se obtienen mejores resultados.

# LEARNER
# =============================================================================
learner_svm <- lrn("regr.svm", type = "eps-regression", kernel = "linear")
 
# PIPELINE PREPROCESADO
# =============================================================================
preproces_pipeline <- po("scale",
                         param_vals = list(center = TRUE, scale = TRUE)) %>>%
                      po("encode",
                         param_vals = list(method = "one-hot"))

# GRAPHLEARNER: PIPELINE + LEARNER
# =============================================================================
# Unión del pipeline con el learner
graph <- preproces_pipeline %>>%
         learner_svm
graph_learner <- GraphLearner$new(graph)

# ESTRATEGIA RESAMPLING
# =============================================================================
resampling_cv <- rsmp("cv", folds = 5)

# MÉTRICA VALIDACIÓN
# =============================================================================
metrica <- msr("regr.rmse")

# HIPERPARÁMETROS
# =============================================================================
# Espacio de búsqueda
hiperparam <- ParamSet$new(
                  list(
                    ParamDbl$new("regr.svm.cost", lower = 1e-5, upper = 1)
                  )
                )

# CONDICIÓN DE STOP
# =============================================================================
# En este caso, al ser un "random_search" se evalúan 15 valores aleatorios
# dentro del espacio de búsqueda establecido
terminator <- trm("evals", n_evals = 20)

# INSTANCIA TuningInstanceSingleCrit
# =============================================================================
grid <- TuningInstanceSingleCrit$new(
          task         = task_train,
          learner      = graph_learner,
          resampling   = resampling_cv,
          measure      = metrica,
          search_space = hiperparam,
          terminator   = terminator
        )

# ESTRATEGIA DE BÚSQUEDA
# =============================================================================
# Estrategia para extraer valores: design_points, gensa, grid_search, random_search
set.seed(123)
tuner_grid <- tnr("random_search")

# OPTIMIZACIÓN
# =============================================================================
future::plan(future::multiprocess())
lgr::get_logger("mlr3")$set_threshold(level = "warn")
lgr::get_logger("bbotk")$set_threshold("warn")
silent_results <- tuner_grid$optimize(grid)

Dentro del objeto grid$archive$data() quedan almacenados los resultados obtenidos con cada valor de hiperparámetros.

grid$archive$data(unnest = "x_domain") %>%
  select(regr.rmse,dplyr::starts_with("x_domain")) %>%
  arrange(regr.rmse)

ggplot(
  data = grid$archive$data(unnest = "x_domain"),
  aes(x = x_domain_regr.svm.cost, y = regr.rmse)
) +
geom_point() +
geom_line() +
labs(title = "Evolución del error de validación") +
theme_bw()

También se almacenan el valor de los hiperparámetros con los que se ha conseguido el mejor resultado.

# Mejores hiperparámetros encontrados
grid$result_learner_param_vals %>% as.data.frame()

Una vez identificados los mejores hiperparámetros, se reentrena el modelo indicando los valores óptimos en sus argumentos.

# Learner final
graph_learner$param_set$values = grid$result_learner_param_vals
graph_learner$train(task_train)

predicciones_test <- graph_learner$predict(task_test)
predicciones_test$score(metrica)

## regr.rmse 
##  58850.12

GLM

Los modelos GLM (John Nelder y Robert Wedderburn) son una generalización de los modelos de regresión lineal que permiten superar limitación de que la variable respuesta se distribuya de forma normal. En los modelos GLM la variable respuesta puede seguir cualquier distribución de la familia exponencial (Gaussian, Poisson, binomial…). Para conseguirlo, la media de la variable respuesta ya no está relacionada directamente con los predictores, sino que se hace a través de una función link g(.). Por ejemplo, el modelo lineal generalizado es equivale al modelo lineal por mínimos cuadrados cuando la función link es la identidad, o a la regresión logística cuando es binomial (logit).

# LEARNER
# =============================================================================
learner_glm <- lrn("regr.glmnet")
knitr::kable(as.data.table(learner_glm$param_set), format = "markdown")

# PIPELINE PREPROCESADO
# =============================================================================
preproces_pipeline <- po("scale",
                         param_vals = list(center = TRUE, scale = TRUE)) %>>%
                      po("encode",
                         param_vals = list(method = "one-hot"))

# GRAPHLEARNER: PIPELINE + LEARNER
# =============================================================================
# Unión del pipeline con el learner
graph <- preproces_pipeline %>>%
         learner_glm
graph_learner <- GraphLearner$new(graph)

# ESTRATEGIA RESAMPLING
# =============================================================================
resampling_cv <- rsmp("cv", folds = 5)

# MÉTRICA VALIDACIÓN
# =============================================================================
metrica <- msr("regr.rmse")

# HIPERPARÁMETROS
# =============================================================================
# Espacio de búsqueda
hiperparam <- ParamSet$new(
                list(
                  ParamDbl$new("regr.glmnet.alpha", lower = 0, upper = 1)
                )
              )

# CONDICIÓN DE STOP
# =============================================================================
# En este caso, al ser un "grid_search" se evalúan todos los valores.
terminator <- trm("none")

# INSTANCIA TuningInstanceSingleCrit
# =============================================================================
grid <- TuningInstanceSingleCrit$new(
          task         = task_train,
          learner      = graph_learner,
          resampling   = resampling_cv,
          measure      = metrica,
          search_space = hiperparam,
          terminator   = terminator
        )

# ESTRATEGIA DE BÚSQUEDA
# =============================================================================
# Estrategia para extraer valores: design_points, gensa, grid_search, random_search
set.seed(123)
tuner_grid <- tnr("grid_search", resolution = 20)

# OPTIMIZACIÓN
# =============================================================================
future::plan(future::multiprocess())
lgr::get_logger("mlr3")$set_threshold(level = "warn")
lgr::get_logger("bbotk")$set_threshold("warn")
silent_results <- tuner_grid$optimize(grid)

# RESULTADOS
# =============================================================================
grid$archive$data(unnest = "x_domain") %>%
  select(regr.rmse,dplyr::starts_with("x_domain")) %>%
  arrange(regr.rmse)

# Mejores hiperparámetros encontrados
grid$result_learner_param_vals

## $scale.center
## [1] TRUE
## 
## $scale.scale
## [1] TRUE
## 
## $scale.robust
## [1] FALSE
## 
## $encode.method
## [1] "one-hot"
## 
## $regr.glmnet.family
## [1] "gaussian"
## 
## $regr.glmnet.alpha
## [1] 0.3157895

ggplot(
  data = grid$archive$data(unnest = "x_domain"),
  aes(x = x_domain_regr.glmnet.alpha, y = regr.rmse)
) +
geom_point() +
geom_line() +
labs(title = "Evolución del error de validación") +
theme_bw()

# MODELO FINAL
# =============================================================================
graph_learner$param_set$values = grid$result_learner_param_vals
graph_learner$train(task_train)

# PREDICCIONES TEST
# =============================================================================
predicciones_test <- graph_learner$predict(task_test)
predicciones_test$score(metrica)

## regr.rmse 
##  58488.96

XGBOOST

Un modelo Gradient Boosting está formado por un conjunto de árboles de decisión individuales, entrenados de forma secuencial, de forma que cada nuevo árbol trata de mejorar los errores de los árboles anteriores. La predicción de una nueva observación se obtiene agregando las predicciones de todos los árboles individuales que forman el modelo.

# LEARNER
# =============================================================================
learner_xgboost <- lrn("regr.xgboost")
knitr::kable(as.data.table(learner_xgboost$param_set), format = "markdown")

# PIPELINE PREPROCESADO
# =============================================================================
preproces_pipeline <- po("scale",
                         param_vals = list(center = TRUE, scale = TRUE)) %>>%
                      po("encode",
                         param_vals = list(method = "one-hot"))

# GRAPHLEARNER: PIPELINE + LEARNER
# =============================================================================
# Unión del pipeline con el learner
graph <- preproces_pipeline %>>% learner_xgboost
graph_learner <- GraphLearner$new(graph)

# ESTRATEGIA RESAMPLING
# =============================================================================
resampling_cv <- rsmp("cv", folds = 3)

# MÉTRICA VALIDACIÓN
# =============================================================================
metrica <- msr("regr.rmse")

# HIPERPARÁMETROS
# =============================================================================
# Espacio de búsqueda
hiperparam <- ParamSet$new(
                list(
                  ParamDbl$new(id = "regr.xgboost.eta", lower = 0.2, upper = .4),
                  ParamDbl$new(id = "regr.xgboost.min_child_weight", lower = 1, upper = 20),
                  ParamDbl$new(id = "regr.xgboost.subsample", lower = .7, upper = .8),
                  ParamDbl$new(id = "regr.xgboost.colsample_bytree",  lower = .9, upper = 1),
                  ParamDbl$new(id = "regr.xgboost.colsample_bylevel", lower = .5, upper = .7),
                  ParamInt$new(id = "regr.xgboost.nrounds", lower = 1L, upper = 25)
                )
              )

# CONDICIÓN DE STOP
# =============================================================================
# En este caso, al ser un "random_search" se evalúan 15 valores aleatorios
# dentro del espacio de búsqueda establecido
terminator <- trm("evals", n_evals = 20)

# INSTANCIA TuningInstanceSingleCrit
# =============================================================================
grid <- TuningInstanceSingleCrit$new(
          task         = task_train,
          learner      = graph_learner,
          resampling   = resampling_cv,
          measure      = metrica,
          search_space = hiperparam,
          terminator   = terminator
        )

# ESTRATEGIA DE BÚSQUEDA
# =============================================================================
# Estrategia para extraer valores: design_points, gensa, grid_search, random_search
set.seed(123)
tuner_points <- tnr("random_search")

# OPTIMIZACIÓN
# =============================================================================
future::plan(future::multiprocess())
lgr::get_logger("mlr3")$set_threshold(level = "warn")
lgr::get_logger("bbotk")$set_threshold("warn")
silent_results <- tuner_grid$optimize(grid)

# RESULTADOS
# =============================================================================
grid$archive$data(unnest = "x_domain") %>%
  select(regr.rmse,dplyr::starts_with("x_domain")) %>%
  arrange(regr.rmse)

# Mejores hiperparámetros encontrados
grid$result_learner_param_vals

## $scale.center
## [1] TRUE
## 
## $scale.scale
## [1] TRUE
## 
## $scale.robust
## [1] FALSE
## 
## $encode.method
## [1] "one-hot"
## 
## $regr.xgboost.nrounds
## [1] 22
## 
## $regr.xgboost.verbose
## [1] 0
## 
## $regr.xgboost.eta
## [1] 0.2
## 
## $regr.xgboost.min_child_weight
## [1] 4
## 
## $regr.xgboost.subsample
## [1] 0.7263158
## 
## $regr.xgboost.colsample_bytree
## [1] 0.9526316
## 
## $regr.xgboost.colsample_bylevel
## [1] 0.6473684

# MODELO FINAL
# =============================================================================
graph_learner$param_set$values = grid$result_learner_param_vals
graph_learner$train(task_train)

# PREDICCIONES TEST
# =============================================================================
predicciones_test <- graph_learner$predict(task_test)
predicciones_test$score(metrica)

## regr.rmse 
##  54611.76

KNN

K-Nearest Neighbor es uno de los algoritmos de machine learning más simples. Se fundamenta en la idea de identificar observaciones en el conjunto de entrenamiento que se asemejen a la observación de test (observaciones vecinas) y asignarle como valor predicho la clase predominante entre dichas observaciones. A pesar de su sencillez, en muchos escenarios consigue resultados aceptables.

# LEARNER
# =============================================================================
learner_kknn <- lrn("regr.kknn")
knitr::kable(as.data.table(learner_kknn$param_set), format = "markdown")

# PIPELINE PREPROCESADO
# =============================================================================
preproces_pipeline <- po("scale",
                         param_vals = list(center = TRUE, scale = TRUE)) %>>%
                      po("encode",
                         param_vals = list(method = "one-hot"))

# GRAPHLEARNER: PIPELINE + LEARNER
# =============================================================================
# Unión del pipeline con el learner
graph <- preproces_pipeline %>>% learner_kknn
graph_learner <- GraphLearner$new(graph)

# ESTRATEGIA RESAMPLING
# =============================================================================
resampling_cv <- rsmp("cv", folds = 5)

# MÉTRICA VALIDACIÓN
# =============================================================================
metrica <- msr("regr.rmse")

# HIPERPARÁMETROS
# =============================================================================
# Espacio de búsqueda
hiperparam <- ParamSet$new(
                list(
                  ParamInt$new("regr.kknn.k", lower = 3, upper = 50),
                  ParamDbl$new("regr.kknn.distance", lower = 1, upper = 3),
                  ParamFct$new("regr.kknn.kernel",
                               c("rectangular", "gaussian", "rank", "optimal")),
                  ParamLgl$new("regr.kknn.scale")
                )
              )

# CONDICIÓN DE STOP
# =============================================================================
# En este caso, al ser un "random_search" se evalúan 15 valores aleatorios
# dentro del espacio de búsqueda establecido
terminator <- trm("evals", n_evals = 20)

# INSTANCIA TuningInstanceSingleCrit
# =============================================================================
grid <- TuningInstanceSingleCrit$new(
          task         = task_train,
          learner      = graph_learner,
          resampling   = resampling_cv,
          measure      = metrica,
          search_space = hiperparam,
          terminator   = terminator
        )

# ESTRATEGIA DE BÚSQUEDA
# =============================================================================
# Estrategia para extraer valores: design_points, gensa, grid_search, random_search
set.seed(123)
tuner_points <- tnr("random_search")

# OPTIMIZACIÓN
# =============================================================================
future::plan(future::multiprocess())
lgr::get_logger("mlr3")$set_threshold(level = "warn")
lgr::get_logger("bbotk")$set_threshold("warn")
silent_results <- tuner_grid$optimize(grid)

# RESULTADOS
# =============================================================================
grid$archive$data(unnest = "x_domain") %>%
  select(regr.rmse,dplyr::starts_with("x_domain")) %>%
  arrange(regr.rmse)

# Mejores hiperparámetros encontrados
grid$result_learner_param_vals

## $scale.center
## [1] TRUE
## 
## $scale.scale
## [1] TRUE
## 
## $scale.robust
## [1] FALSE
## 
## $encode.method
## [1] "one-hot"
## 
## $regr.kknn.k
## [1] 13
## 
## $regr.kknn.distance
## [1] 2.473684
## 
## $regr.kknn.kernel
## [1] "rank"
## 
## $regr.kknn.scale
## [1] FALSE

# MODELO FINAL
# =============================================================================
graph_learner$param_set$values = grid$result_learner_param_vals
graph_learner$train(task_train)

# PREDICCIONES TEST
# =============================================================================
predicciones_test <- graph_learner$predict(task_test)
predicciones_test$score(metrica)

## regr.rmse 
##  61809.16

Anexo1: Métodos de validación

Los métodos de validación, también conocidos como resampling, son estrategias que permiten estimar la capacidad predictiva de los modelos cuando se aplican a nuevas observaciones, haciendo uso únicamente de los datos de entrenamiento. La idea en la que se basan todos ellos es la siguiente: el modelo se ajusta empleando un subconjunto de observaciones del conjunto de entrenamiento y se evalúa (calcular una métrica que mida cómo de bueno es el modelo, por ejemplo, accuracy) con las observaciones restantes. Este proceso se repite múltiples veces y los resultados se agregan y promedian. Gracias a las repeticiones, se compensan las posibles desviaciones que puedan surgir por el reparto aleatorio de las observaciones. La diferencia entre métodos suele ser la forma en la que se generan los subconjuntos de entrenamiento/validación.

k-Fold-Cross-Validation (CV)

Las observaciones de entrenamiento se reparten en k folds (conjuntos) del mismo tamaño. El modelo se ajusta con todas las observaciones excepto las del primer fold y se evalúa prediciendo las observaciones del fold que ha quedado excluido, obteniendo así la primera métrica. El proceso se repite k veces, excluyendo un fold distinto en cada iteración. Al final, se generan k valores de la métrica, que se agregan (normalmente con la media y la desviación típica) generando la estimación final de validación.

Leave-One-Out Cross-Validation (LOOCV)

LOOCV es un caso especial de k-Fold-Cross-Validation en el que el número k de folds es igual al número de observaciones disponibles en el conjunto de entrenamiento. El modelo se ajusta cada vez con todas las observaciones excepto una, que se emplea para evaluar el modelo. Este método supone un coste computacional muy elevado, el modelo se ajusta tantas veces como observaciones de entrenamiento, por lo que, en la práctica, no suele compensar emplearlo.

Repeated k-Fold-Cross-Validation (repeated CV)

Es exactamente igual al método k-Fold-Cross-Validation pero repitiendo el proceso completo n veces. Por ejemplo, 10-Fold-Cross-Validation con 5 repeticiones implica a un total de 50 iteraciones ajuste-validación, pero no equivale a un 50-Fold-Cross-Validation.

Leave-Group-Out Cross-Validation (LGOCV)

LGOCV, también conocido como repeated train/test splits o Monte Carlo Cross-Validation, consiste simplemente en generar múltiples divisiones aleatorias entrenamiento-test (solo dos conjuntos por repetición). La proporción de observaciones que va a cada conjunto se determina de antemano, 80%-20% suele dar buenos resultados. Este método, aunque más simple de implementar que CV, requiere de muchas repeticiones (>50) para generar estimaciones estables.

Bootstrapping

Una muestra bootstrap es una muestra obtenida a partir de la muestra original por muestreo aleatorio con reposición, y del mismo tamaño que la muestra original. Muestreo aleatorio con reposición (resampling with replacement) significa que, después de que una observación sea extraída, se vuelve a poner a disposición para las siguientes extracciones. Como resultado de este tipo de muestreo, algunas observaciones aparecerán múltiples veces en la muestra bootstrap y otras ninguna. Las observaciones no seleccionadas reciben el nombre de out-of-bag (OOB). Por cada iteración de bootstrapping se genera una nueva muestra bootstrap, se ajusta el modelo con ella y se evalúa con las observaciones out-of-bag.

1. Obtener una nueva muestra del mismo tamaño que la muestra original mediante muestreo aleatorio con reposición.

2. Ajustar el modelo empleando la nueva muestra generada en el paso 1.

3. Calcular el error del modelo empleando aquellas observaciones de la muestra original que no se han incluido en la nueva muestra. A este error se le conoce como error de validación.

4. Repetir el proceso n veces y calcular la media de los n errores de validación.

5. Finalmente, y tras las n repeticiones, se ajusta el modelo final empleando todas las observaciones de entrenamiento originales.

La naturaleza del proceso de bootstrapping genera cierto bias en las estimaciones que puede ser problemático cuando el conjunto de entrenamiento es pequeño. Existen ciertas modificaciones del algoritmo original para corregir este problema, algunos de ellos son: 632 method y 632+ method.

No existe un método de validación que supere al resto en todos los escenarios, la elección debe basarse en varios factores.

• Si el tamaño de la muestra es pequeño, se recomienda emplear repeated k-Fold-Cross-Validation, ya que consigue un buen equilibrio bias-varianza y, dado que no son muchas observaciones, el coste computacional no es excesivo.

• Si el objetivo principal es comparar modelos mas que obtener una estimación precisa de las métricas, se recomienda bootstrapping ya que tiene menos varianza.

• Si el tamaño muestral es muy grande, la diferencia entre métodos se reduce y toma más importancia la eficiencia computacional. En estos casos, 10-Fold-Cross-Validation simple es suficiente.

Puede encontrarse un estudio comparativo de los diferentes métodos en Comparing Different Species of Cross-Validation.

Anexo2: Métricas

Existe una gran variedad de métricas que permiten evaluar cómo de bueno es un algoritmo realizando predicciones. La idoneidad de cada una depende completamente del problema en cuestión, y su correcta elección dependerá de cómo de bien entienda el analista el problema al que se enfrenta. A continuación, se describen algunas de las más utilizadas.

Accuracy y Kappa

Estas dos métricas son las más empleadas en problemas de clasificación binaria y multiclase. Accuracy es el porcentaje de observaciones correctamente clasificadas respecto al total de predicciones. Kappa o Cohen’s Kappa es el valor de accuracy normalizado respecto del porcentaje de acierto esperado por azar. A diferencia de accuracy, cuyo rango de valores puede ser [0, 1], el de kappa es [-1, 1]. En problemas con clases desbalanceadas, donde el grupo mayoritario supera por mucho a los otros, Kappa es más útil porque evita caer en la ilusión de creer que un modelo es bueno cuando realmente solo supera por poco lo esperado por azar.

MSE, RMSE, MAE

Estas son las métricas más empleadas en problemas de regresión.

MSE (Mean Squared Error) es la media de los errores elevados al cuadrado. Suele ser muy buen indicativo de cómo funciona el modelo en general, pero tiene la desventaja de estar en unidades cuadradas. Para mejorar la interpretación, suele emplearse RMSE (Root Mean Squared Error), que es la raíz cuadrada del MSE y por lo tanto sus unidades son las mismas que la variable respuesta.

MAE (Mean Absolute Error) es la media de los errores en valor absoluto. La diferencia respecto a MSE es que, este último, eleva al cuadrado los errores, lo significa que penaliza mucho más las desviaciones grandes. A modo general, MSE favorece modelos que se comportan aproximadamente igual de bien en todas las observaciones, mientras que MAE favorece modelos que predicen muy bien la gran mayoría de observaciones aunque en unas pocas se equivoque por mucho.