Etapas de un problema de machine learning

El siguiente es un listado de las etapas que suelen formar parte de la mayoría de problemas de machine learning.

• Definir el problema: ¿Qué se pretende predecir? ¿De qué datos se dispone? o ¿Qué datos es necesario conseguir?

• Explorar y entender los datos que se van a emplear para crear el modelo.

• Métrica de existo: definir una forma apropiada de cuantificar cómo de buenos son los resultados obtenidos.

• Preparar la estrategia para evaluar el modelo: separar las observaciones en un conjunto de entrenamiento, un conjunto de validación (este último suele ser un subconjunto del de entrenamiento) y un conjunto de test. Ninguna información del conjunto de test debe participar en el proceso de entrenamiento del modelo.

• Preprocesar los datos: aplicar las transformaciones necesarias para que los datos puedan ser interpretados por el algoritmo de machine learning seleccionado.

• Ajustar un primer modelo capaz de superar unos resultados mínimos. Por ejemplo, en problemas de clasificación, el mínimo a superar es el porcentaje de la clase mayoritaria (la moda).

• Gradualmente, mejorar el modelo optimizando sus hiperparámetros.

• Evaluar la capacidad del modelo final con el conjunto de test para tener una estimación de la capacidad que tiene el modelo cuando predice nuevas observaciones.
  
  

Paquete caret

R es uno de los lenguajes de programación que domina dentro del ámbito de la estadística, data mining y machine learning. Al tratarse de un software libre, innumerables usuarios han podido implementar sus algoritmos, dando lugar a un número muy elevado de paquetes/librerías donde encontrar prácticamente todas las técnicas de machine learning existentes. Sin embargo, esto tiene un lado negativo, cada paquete tiene una sintaxis, estructura e implementación propia, lo que dificulta su aprendizaje. El paquete caret, desarrollado por Max Kuhn, es una interfaz que unifica bajo un único marco cientos de funciones de distintos paquetes, facilitando en gran medida todas las etapas de de preprocesado, entrenamiento, optimización y validación de modelos predictivos.

El paquete caret ofrece tal cantidad de posibilidades que, difícilmente, pueden ser mostradas con un único ejemplo. En este documento, se emplean solo algunas de sus funcionalidades. Si en algún caso se requiere una explicación detallada, para que no interfiera con la narrativa del análisis, se añadirá un anexo. Aun así, para conocer bien todas sus funcionalidades se recomienda leer su documentación.

Otros proyectos similares a caret pero más recientes son: mlr3, tidymodels y H2O

# Instalación de los paquetes que unifica caret. Esta instalación puede tardar.
# Solo es necesario ejecutarla si no se dispone de los paquetes.
install.packages("caret", dependencies = c("Depends", "Suggests"))

library(caret)

Dataset Titanic

A poco que el lector haya buscado ejemplos sobre machine learning o análisis predictivo, seguro que ha oído hablar del dataset Titanic disponible en la plataforma Kaggle. Este set de datos contiene información sobre los pasajeros del RMS Titanic, el transatlántico británico que se hundió en abril de 1912 durante su viaje inaugural desde Southampton a Nueva York. Entre la información almacenada se encuentran la edad, género, características socio-económicas de los pasajeros y si sobrevivieron o no al naufragio. Aunque pueda resultar un clásico poco original, estos datos tienen una serie de características que los hacen idóneos para ser utilizados como ejemplo introductorio al machine learning:

• Contiene suficientes observaciones para entrenar modelos que consigan un poder predictivo alto.

• Incluye tanto variables continuas como cualitativas, lo que permite mostrar diferentes análisis exploratorios.

• La variable respuesta es binaria. Aunque la mayoría de los algoritmos de clasificación mostrados en este capítulo se generalizan para múltiples clases, su interpretación suele ser más sencilla cuando solo hay dos.

• Contiene valores ausentes. La forma en que se manejan estos registros (eliminación o imputación) influye en gran medida en el modelo final.

• Aunque a primera vista no lo parezca, las variables contienen información adicional que puede ser extraída mediante feature engineering.

• Requiere ciertos pasos de limpieza y conversión de datos.

y lo más importante de todo, se trata de un problema y de unos datos cuyas variables pueden entenderse de forma sencilla. Es intuitivo comprender el impacto que puede tener la edad, el sexo, la localización del camarote… en la supervivencia de los pasajeros. Aunque no lo parezca, comprender a fondo el problema que se pretende modelar es lo más importante para lograr buenos resultados.

En el paquete titanic se pueden encontrar todos los datos proporcionados por kaggle. Están divididos en dos partes, un conjunto de entrenamiento y un conjunto de test. En los datos de test, se desconoce si el pasajero sobrevivió o no, esto es así porque la plataforma Kaggle evalúa con ellos la capacidad predictiva de los modelos presentados, es decir, el competidor crea un modelo con los datos de entrenamiento, predice la variable respuesta (superviviente o no superviviente) en el conjunto de test y envía sus predicciones a la plataforma. La plataforma contrasta las predicciones con los verdaderos resultados y devuelve una puntuación. En este ejemplo, para evitar tener que estar enviado los resultados de cada modelo, se emplea únicamente el conjunto de entrenamiento. Además, esto refleja mejor lo que ocurre en la práctica cuando un analista se enfrenta a un nuevo set de datos.

library(tidyverse)
library(titanic)
datos <- titanic_train

Tipo de variables

Una de las primeras comprobaciones que hay que hacer tras cargar los datos, es verificar que cada variable se ha almacenado con el tipo de valor que le corresponde, es decir, que las variables numéricas sean números y las cualitativas factor, character o booleanas. En el lenguaje de programación R, cuando la variable es cualitativa, conviene almacenarla con el tipo factor.

# Resumen del set de datos
glimpse(datos)

## Rows: 891
## Columns: 12
## $ PassengerId <int> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17…
## $ Survived    <int> 0, 1, 1, 1, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 1, 0, 1, 0, …
## $ Pclass      <int> 3, 1, 3, 1, 3, 3, 1, 3, 3, 2, 3, 1, 3, 3, 3, 2, 3, 2, 3, …
## $ Name        <chr> "Braund, Mr. Owen Harris", "Cumings, Mrs. John Bradley (F…
## $ Sex         <chr> "male", "female", "female", "female", "male", "male", "ma…
## $ Age         <dbl> 22, 38, 26, 35, 35, NA, 54, 2, 27, 14, 4, 58, 20, 39, 14,…
## $ SibSp       <int> 1, 1, 0, 1, 0, 0, 0, 3, 0, 1, 1, 0, 0, 1, 0, 0, 4, 0, 1, …
## $ Parch       <int> 0, 0, 0, 0, 0, 0, 0, 1, 2, 0, 1, 0, 0, 5, 0, 0, 1, 0, 0, …
## $ Ticket      <chr> "A/5 21171", "PC 17599", "STON/O2. 3101282", "113803", "3…
## $ Fare        <dbl> 7.2500, 71.2833, 7.9250, 53.1000, 8.0500, 8.4583, 51.8625…
## $ Cabin       <chr> "", "C85", "", "C123", "", "", "E46", "", "", "", "G6", "…
## $ Embarked    <chr> "S", "C", "S", "S", "S", "Q", "S", "S", "S", "C", "S", "S…

En la tabla anterior puede verse el tipo de cada variable, así como un ejemplo de los valores que toma. A priori, el único caso en el que el tipo de valor no se corresponde con la naturaleza de la variable es Survived. Aunque esta variable está codificada como 1 si el pasajero sobrevivió y 0 si murió, no conviene almacenarla en formato numérico, ya que esto puede llevar a errores como el de tratar de calcular su media. Para evitar este tipo de problemas, se recodifica la variable para que sus dos posibles niveles sean “Si”-“No” y se convierte a factor.

datos$Survived <- if_else(datos$Survived == 1, "Si", "No")
datos$Survived <- as.factor(datos$Survived)

La variable Pclass es cualitativa ordinal, es decir, toma distintos valores cualitativos ordenados siguiendo una escala establecida, aunque no es necesario que el intervalo entre mediciones sea uniforme. Por ejemplo, se asume que la diferencia entre primera y segunda clase es menor que la diferencia entre primera y tercera, sin embargo, las diferencias entre primera-segunda y segunda-tercera no tiene por qué ser iguales. Dado que las propiedades matemáticas no se cumplen \((2-1 \neq 3-2)\) es preferible no almacenarlas como números.

datos$Pclass <- as.factor(datos$Pclass)

Las variables SibSp y Parch son cuantitativas discretas, pueden tomar únicamente determinados valores numéricos. En este caso, al tratarse de número de personas (familiares e hijos), solo pueden ser números enteros. No existe una norma clara sobre como almacenar estas variables. Para este estudio exploratorio, dado que solo toman unos pocos valores, se decide almacenarlas como factor.

datos$SibSp <- as.factor(datos$SibSp)
datos$Parch <- as.factor(datos$Parch)

Las variables Sex y Embarked también se convierten a tipo factor.

datos$Sex      <- as.factor(datos$Sex)
datos$Embarked <- as.factor(datos$Embarked)

El análisis exploratorio de variables suele caracterizarse por el cálculo de sumatorios, reestructuración de los datos y representaciones gráficas. Los paquetes dplyr, tidyr, ggplot… englobados dentro de la filosofía tydiverse facilitan en gran medida todos estos pasos. Muchos de los cálculos y representaciones que se realizan a lo largo de este capítulo se consiguen de forma más rápida si los datos están almacenados siguiendo la estructura: observación, variable, valor. Como el set de datos no es lo suficientemente grande como para dar problemas de memoria, se crea un segundo dataframe con esta estructura.

datos_long <- datos %>% gather(key = "variable", value = "valor", -PassengerId)
head(datos_long)

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, es decir, verificar si para cada observación se ha registrado el valor de cada una de las variables.

# Número de observaciones del set de datos
nrow(datos)

## [1] 891

# Detección si hay alguna fila incompleta
any(!complete.cases(datos))

## [1] TRUE

Una vez detectado que existen valores ausentes, se estudia la distribución de los mismos.

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

## PassengerId    Survived      Pclass        Name         Sex         Age 
##           0           0           0           0           0         177 
##       SibSp       Parch      Ticket        Fare       Cabin    Embarked 
##           0           0           0           0           0           0

Viendo esta tabla, aparece una contradicción con el glimse() mostrado en el apartado anterior. El primer valor de la variable Cabin es "“, es decir, no se dispone del valor, sin embargo, el contaje de los valores ausentes no lo muestra. Esto ocurre porque R interpreta el valor”" como un character, no como valor ausente NA.

Se procede a identificar qué variables contienen valores "".

datos %>% map_lgl(.f = function(x){any(!is.na(x) & x == "")})

## PassengerId    Survived      Pclass        Name         Sex         Age 
##       FALSE       FALSE       FALSE       FALSE       FALSE       FALSE 
##       SibSp       Parch      Ticket        Fare       Cabin    Embarked 
##       FALSE       FALSE       FALSE       FALSE        TRUE        TRUE

Las variables Cabin y Embarked contienen al menos un valor "", se sustituyen por NA.

# La variable Cabin está almacenada como character
datos$Cabin[datos$Cabin == ""] <- NA

La sustitución de valores en factors debe de hacerse con más cautela que cuando se hace en vectores de otro tipo. Esto se debe a que los factors almacenan los niveles originales y no se eliminan, aunque el factor ya no contenga ese valor. Una forma sencilla de evitar este problema es convertir el factor a character, hacer la modificación y volver a convertirlo a factor.

levels(datos$Embarked)

## [1] ""  "C" "Q" "S"

datos$Embarked <- as.character(datos$Embarked)
datos$Embarked[datos$Embarked == ""] <- NA
datos$Embarked <- as.factor(datos$Embarked)
levels(datos$Embarked)

## [1] "C" "Q" "S"

# Este cambio también se aplica al dataframe datos_long
datos_long$valor[datos_long$valor == ""] <- NA

Una vez sustituidos se repite el contaje de valores ausentes.

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

## PassengerId    Survived      Pclass        Name         Sex         Age 
##           0           0           0           0           0         177 
##       SibSp       Parch      Ticket        Fare       Cabin    Embarked 
##           0           0           0           0         687           2

# Representación gráfica de los datos ausentes
datos_long <- datos_long %>%  mutate(ausente = is.na(valor))
ggplot(data = datos_long, aes(x = variable, y = PassengerId, fill = ausente)) +
  geom_raster() +
  scale_fill_manual(values = c("gray60", "orangered2")) +
  theme_bw() +
  labs(title = "Valores ausentes por variable") +
  theme(legend.position = "bottom")

datos_long %>%
  group_by(variable) %>% 
  summarize(porcentaje_NA = 100 * sum(is.na(valor)) / length(valor)) %>%
  ggplot(aes(x = reorder(variable, desc(porcentaje_NA)), y = porcentaje_NA)) +
    geom_col() +
    labs(title = "Porcentaje valores ausentes por variable",
         x = "Variable", y = "Porcentaje NAs") +
    theme_bw()

La distribución de los valores ausentes muestra que, para poco más del 20% de los pasajeros, se dispone de información sobre la cabina. La variable Age también está ausente en un número considerable de pasajeros. Solo para dos pasajeros se desconoce el puerto desde el que embarcaron.

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. Por ejemplo, solo dos observaciones tienen ausente la variable Embarked, si esas dos observaciones caen en el conjunto de test, al estudiar el conjunto de entrenamiento, no se identificaría la necesidad de imputar esta variable.

Distribución de variables 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 nos interesa predecir.

ggplot(data = datos, aes(x = Survived, y = ..count.., fill = Survived)) +
  geom_bar() +
  scale_fill_manual(values = c("gray50", "orangered2")) +
  labs(title = "Supervivencia") +
  theme_bw() +
  theme(legend.position = "bottom")

# Tabla de frecuencias 
table(datos$Survived)

## 
##  No  Si 
## 549 342

prop.table(table(datos$Survived)) %>% round(digits = 2)

## 
##   No   Si 
## 0.62 0.38

Para que un modelo predictivo sea útil, debe de tener un porcentaje de acierto superior a lo esperado por azar o a un determinado nivel basal. En problemas de clasificación, el nivel basal es el que se obtiene si se asignan todas las observaciones a la clase mayoritaria (la moda). En el naufragio del Titanic, dado que el 62% de los pasajeros fallecieron, si siempre se predice Survived = No, el porcentaje de aciertos será aproximadamente del 62%. Este es el porcentaje mínimo que hay que intentar superar con los modelos predictivos. (Siendo estrictos, este porcentaje tendrá que ser recalculado únicamente con el conjunto de entrenamiento).

# Porcentaje de aciertos si se predice para todas las observaciones que no sobrevivieron.
n_observaciones <- nrow(datos)
predicciones <- rep(x = "No",  n_observaciones)
mean(predicciones == datos$Survived) * 100

## [1] 61.61616

Distribución de variables continuas

Como el objetivo del estudio es predecir qué pasajeros sobrevivieron y cuáles no, el análisis de cada variable se hace en relación a la variable respuesta Survived. Analizando los datos de esta forma, se pueden empezar a extraer ideas sobre qué variables están más relacionadas con la supervivencia.

library(ggpubr)
p1 <- ggplot(data = datos, aes(x = Age, fill = Survived)) +
      geom_density(alpha = 0.5) +
      scale_fill_manual(values = c("gray50", "orangered2")) +
      geom_rug(aes(color = Survived), alpha = 0.5) +
      scale_color_manual(values = c("gray50", "orangered2")) +
      theme_bw()
p2 <- ggplot(data = datos, aes(x = Survived, y = Age, color = Survived)) +
      geom_boxplot(outlier.shape = NA) +
      geom_jitter(alpha = 0.3, width = 0.15) +
      scale_color_manual(values = c("gray50", "orangered2")) +
      theme_bw()
final_plot <- ggarrange(p1, p2, legend = "top")
final_plot <- annotate_figure(final_plot, top = text_grob("Age", size = 15))
final_plot

# Estadísticos de la edad de los supervivientes y fallecidos
datos %>% filter(!is.na(Age)) %>% group_by(Survived) %>%
          summarise(media = mean(Age),
                    mediana = median(Age),
                    min = min(Age),
                    max = max(Age))

La distribución de la edad de los pasajeros parece ser muy similar entre el grupo de supervivientes y fallecidos, con dos excepciones: en el rango de edad aproximado de 0 a 10 años, el porcentaje de supervivencia es mucho mayor, mientras que, en el extremo opuesto, a partir de los 60 años, la tendencia se invierte. Dos hipótesis que podrían explicar estos patrones son: que, según los registros, en el protocolo de evacuación del Titanic se priorizó que mujeres y niños subiesen a los botes salvavidas, y que los ancianos tuviesen menor movilidad para alcanzar las zonas de evacuación.

Cuando la información de una variable continua reside en si se superan o no determinados límites, los modelos predictivos suelen conseguir mejores resultados si la variable se discretiza en intervalos. En este caso, una posible aproximación es crear una nueva variable que clasifique a los pasajeros en niño, adulto o anciano.

datos <- datos %>%
         mutate(Age_grupo = case_when(Age <= 10  ~ "niño",
                                      Age > 10 & Age <= 60  ~ "adulto",
                                      Age > 60 ~ "anciano"))
datos$Age_grupo <- as.factor(datos$Age_grupo)

Esta nueva variable se analizará junto con el resto de variables cualitativas.

p1 <- ggplot(data = datos, aes(x = Fare, fill = Survived)) +
      geom_density(alpha = 0.5) +
      scale_fill_manual(values = c("gray50", "orangered2")) +
      geom_rug(aes(color = Survived), alpha = 0.5) +
      scale_color_manual(values = c("gray50", "orangered2")) +
      theme_bw()
p2 <- ggplot(data = datos, aes(x = Survived, y = Fare, color = Survived)) +
      geom_boxplot(outlier.shape = NA) +
      geom_jitter(alpha = 0.3, width = 0.15) +
      scale_color_manual(values = c("gray50", "orangered2")) +
      theme_bw()
final_plot <- ggarrange(p1, p2, legend = "top")
final_plot <- annotate_figure(final_plot, top = text_grob("Fare", size = 15))
final_plot

# Estadísticos del precio del billete de los supervivientes y fallecidos
datos %>% filter(!is.na(Fare)) %>% group_by(Survived) %>%
          summarise(media = mean(Fare),
                    mediana = median(Fare),
                    min = min(Fare),
                    max = max(Fare))

La variable Fare tiene una distribución asimétrica, muchos billetes tenían un coste bajo y unos pocos un coste alto. Este tipo de distribución suele visualizarse mejor tras una trasformación logarítmica.

p1 <- ggplot(data = datos, aes(x = log(Fare), fill = Survived)) +
      geom_density(alpha = 0.5) +
      scale_fill_manual(values = c("gray50", "orangered2")) +
      geom_rug(aes(color = Survived), alpha = 0.5) +
      scale_color_manual(values = c("gray50", "orangered2")) +
      theme_bw()
p2 <- ggplot(data = datos, aes(x = Survived, y = log(Fare), color = Survived)) +
      geom_boxplot(outlier.shape = NA) +
      geom_jitter(alpha = 0.3, width = 0.15) +
      scale_color_manual(values = c("gray50", "orangered2")) +
      theme_bw()
final_plot <- ggarrange(p1, p2, legend = "top")
final_plot <- annotate_figure(final_plot, top = text_grob("Log(Fare)", size =15))
final_plot

Los datos indican que el precio medio de los billetes de las personas que sobrevivieron era superior al de los que fallecieron. ¿Significa esto que, haber pagado más por el billete, garantizaba una mayor probabilidad de supervivencia? No parece muy verosímil que esto fuese así. Este es un caso típico en el que cabe recordar que correlación/asociación no son sinónimos de causalidad. Es probable que exista otra variable subyacente que sí tenga un vínculo de causalidad con la probabilidad de supervivencia y que, a su vez, esté correlacionada con el precio del billete. Una hipótesis podría ser: los pasajeros de primera clase tuvieron ciertas facilidades para acceder a los botes de evacuación por lo que consiguieron salvarse con más éxito que los de segunda y tercera clase. Los billetes de primera clase costaban mucho más dinero que los de las clases inferiores, de ahí que el precio medio de billetes entre los supervivientes fuese superior.

Distribución de variables cualitativas

ggplot(data = datos, aes(x = Pclass, y = ..count.., fill = Survived)) +
  geom_bar() +
  scale_fill_manual(values = c("gray50", "orangered2")) +
  labs(title = "Pclass") +
  theme_bw() +
  theme(legend.position = "bottom")

# Tabla de frecuencias relativas de supervivientes por clase
prop.table(table(datos$Pclass, datos$Survived), margin = 1) %>% round(digits = 2)

##    
##       No   Si
##   1 0.37 0.63
##   2 0.53 0.47
##   3 0.76 0.24

ggplot(data = datos, aes(x = Sex, y = ..count.., fill = Survived)) +
  geom_bar() +
  labs(title = "Sex") +
  scale_fill_manual(values = c("gray50", "orangered2")) +
  theme_bw() +
  theme(legend.position = "bottom")

# Tabla de frecuencias relativas de supervivientes por sexo
prop.table(table(datos$Sex, datos$Survived), margin = 1) %>% round(digits = 2)

##         
##            No   Si
##   female 0.26 0.74
##   male   0.81 0.19

ggplot(data = datos, aes(x = SibSp, y = ..count.., fill = Survived)) +
  geom_bar() +
  scale_fill_manual(values = c("gray50", "orangered2")) +
  labs(title = "SibSp") +
  theme_bw() +
  theme(legend.position = "bottom")

# Tabla de frecuencias relativas de supervivientes por número de familiares
prop.table(table(datos$SibSp, datos$Survived), margin = 1) %>% round(digits = 2)

##    
##       No   Si
##   0 0.65 0.35
##   1 0.46 0.54
##   2 0.54 0.46
##   3 0.75 0.25
##   4 0.83 0.17
##   5 1.00 0.00
##   8 1.00 0.00

ggplot(data = datos, aes(x = Parch, y = ..count.., fill = Survived)) +
  geom_bar() +
  scale_fill_manual(values = c("gray50", "orangered2")) +
  labs(title = "Parch") +
  theme_bw() +
  theme(legend.position = "bottom")

# Tabla de frecuencias relativas de supervivientes por Parch
prop.table(table(datos$Parch, datos$Survived), margin = 1) %>% round(digits = 2)

##    
##       No   Si
##   0 0.66 0.34
##   1 0.45 0.55
##   2 0.50 0.50
##   3 0.40 0.60
##   4 1.00 0.00
##   5 0.80 0.20
##   6 1.00 0.00

Las variables SibSp y Parch, ambas relacionadas con el número de familiares a bordo, se han tratado como una variable cualitativa para el análisis exploratorio. Al agrupar las observaciones en los diferentes niveles, algunos de ellos apenas contienen unas pocas (Parch = 4, 5, 6). Que niveles de una variable cualitativa tengan muy poca representación (distribución desbalanceada) puede suponer un problema a la hora de ajustar los modelos u optimizarlos mediante validación cruzada. Por esta razón, a la hora de ajustar el modelo, conviene seguir una de estas dos estrategias:

• Agrupar todos los grupos minoritarios en uno solo, por ejemplo, combinar los grupos Parch = 4, 5, 6 en uno nuevo que sea Parch > 3.

• Tratar el predictor como una variable continua.

En este caso, se opta por la segunda opción.

# Para pasar de factor a numeric primero se convierte a character
datos$SibSp <- as.character(datos$SibSp)
datos$SibSp <- as.numeric(datos$SibSp)
datos$Parch <- as.character(datos$Parch)
datos$Parch <- as.numeric(datos$Parch)

ggplot(data = datos, aes(x = Embarked, y = ..count.., fill = Survived)) +
  geom_bar() +
  scale_fill_manual(values = c("gray50", "orangered2")) +
  labs(title = "Embarked") +
  theme_bw() +
  theme(legend.position = "bottom")

# Tabla de frecuencias relativas de supervivientes por puerto de embarque
prop.table(table(datos$Embarked, datos$Survived), margin = 1) %>% round(digits = 2)

##    
##       No   Si
##   C 0.45 0.55
##   Q 0.61 0.39
##   S 0.66 0.34

ggplot(data = datos, aes(x = Age_grupo, y = ..count.., fill = Survived)) +
  geom_bar() +
  scale_fill_manual(values = c("gray50", "orangered2")) +
  labs(title = "Age_grupo") +
  theme_bw() +
  theme(legend.position = "bottom")

# Tabla de frecuencias relativas de supervivientes por grupo de edad
prop.table(table(datos$Age_grupo, datos$Survived), margin = 1) %>% round(digits = 2)

##          
##             No   Si
##   adulto  0.61 0.39
##   anciano 0.77 0.23
##   niño    0.41 0.59

Las variables PassengerId, Name y Cabin no se han incluido en la exploración porque, a priori, no aportan información relevante sobre la supervivencia. En la sección feature engineering se analizan con más detenimiento.

Importancia de las variables

La representación gráfica de la distribución de las variables en función de si los pasajeros sobrevivieron o no, ayuda a tener una idea de qué variables pueden ser buenos predictores para el modelo y cuales no aportan información o la que aportan es redundante. Aunque la creación de un buen modelo debe entenderse como un proceso iterativo, en el que se van ajustando y probando distintos modelos, existen ciertas pistas que pueden ayudar a realizar una selección inicial adecuada.

• Si dos variables numéricas están muy correlacionadas, añaden información redundante al modelo, por lo tanto, no conviene incorporar ambas. Si esto ocurre, se puede: excluir aquella que, acorde al criterio del analista, no está realmente asociada con la variable respuesta; o combinarlas para recoger toda su información en una única nueva variable, por ejemplo, con un PCA.

• Si una variable tiene varianza igual o próxima a cero (su valor es el mismo o casi el mismo para todas las observaciones) añade al modelo más ruido que información, por lo que suele ser conveniente excluirla.

• 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) o asegurar que, en la creación de las particiones, se garantice que todos los grupos estén representados en cada una de ellas.
  
  

cor.test(x = datos$Age, y = datos$Fare, method = "pearson")

## 
##  Pearson's product-moment correlation
## 
## data:  datos$Age and datos$Fare
## t = 2.5753, df = 712, p-value = 0.01022
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.02285549 0.16825304
## sample estimates:
##        cor 
## 0.09606669

ggplot(data = datos, aes(x = Age, y = log(Fare))) +
  geom_point(color = "gray30") +
  geom_smooth(color = "firebrick") +
  theme_bw()

ggplot(data = datos, aes(x = Age_grupo, y = ..count.., fill = Survived)) +
  geom_bar() +
  scale_fill_manual(values = c("gray50", "orangered2")) +
  labs(title = "Age_grupo") +
  theme_bw() +
  theme(legend.position = "bottom")

# Se excluyen las variables continuas y las cualitativas que no agrupan a los
# pasajeros. También la variable Cabin por su alto % de valores NA.
datos_cualitativos <- datos %>%
                      select(-Age, -Fare, -Name, -Ticket, -Cabin, -PassengerId)

datos_cualitativos_tidy <- datos_cualitativos %>%
                           gather(key = "variable", value = "grupo",-Survived)

# Se eliminan los valores NA para que no se interpreten como un grupo
datos_cualitativos_tidy <- datos_cualitativos_tidy %>% filter(!is.na(grupo))
  
# Se añade un identificador formado por el nombre de la variable y el grupo 
datos_cualitativos_tidy <- datos_cualitativos_tidy %>%
                           mutate(variable_grupo = paste(variable, grupo, sep = "_"))

# Función que calcula el test de proporciones para la columna "Survived" de un df
test_proporcion <- function(df){
  n_supervivientes <- sum(df$Survived == "Si") 
  n_fallecidos     <- sum(df$Survived == "No")
  n_total <- n_supervivientes + n_fallecidos
  test <- prop.test(x = n_supervivientes, n = n_total, p = 0.3838)
  prop_supervivientes <- n_supervivientes / n_total
  return(data.frame(p_value = test$p.value, prop_supervivientes))
}

# Se agrupan los datos por "variable_grupo" y se aplica a cada grupo la función
# test_proporcion()
analisis_prop <- datos_cualitativos_tidy %>%
                 group_by(variable_grupo) %>%
                 nest() %>%
                 arrange(variable_grupo) %>%
                 mutate(prop_test = map(.x = data, .f = test_proporcion)) %>%
                 unnest(prop_test) %>%
                 arrange(p_value) %>% 
                 select(variable_grupo,p_value, prop_supervivientes)
analisis_prop 

# Representación gráfica de la distribución de los 6 grupos con menor p-value
top6_grupos <- analisis_prop %>% pull(variable_grupo) %>% head(6)

# Se crea una función que, dados un dataframe y el nombre de un grupo, genere la
# representación gráfica de supervivientes y no supervivientes.
plot_grupo <- function(grupo, df, threshold_line = 0.3838){

  p <- ggplot(data = df, aes(x = 1, y = ..count.., fill = Survived)) +
            geom_bar() +
            scale_fill_manual(values = c("gray50", "orangered2")) +
            # Se añade una línea horizontal en el nivel basal
            geom_hline(yintercept = nrow(df) * threshold_line,
                       linetype = "dashed") +
            labs(title = grupo) +
            theme_bw() +
            theme(legend.position = "bottom",
                  axis.text.x = element_blank(),
                  axis.title.x = element_blank(),
                  axis.ticks.x = element_blank())
  return(p)
}

datos_graficos <- datos_cualitativos_tidy %>%
                  filter(variable_grupo %in% top6_grupos) %>%
                  group_by(variable_grupo) %>%
                  nest() %>%
                  arrange(variable_grupo)

plots <- map2(datos_graficos$variable_grupo, .y = datos_graficos$data,
              .f = plot_grupo)

ggarrange(plotlist = plots, common.legend = TRUE)

library(randomForest)
datos_rf <- datos %>%
            select(-PassengerId, -Name, -Ticket, -Cabin, -Age, -Age_grupo) %>%
            na.omit()
datos_rf <- map_if(.x = datos_rf, .p = is.character, .f = as.factor) %>%
            as.data.frame()
modelo_randforest <- randomForest(formula = Survived ~ . ,
                                  data = na.omit(datos_rf),
                                  mtry = 5,
                                  importance = TRUE, 
                                  ntree = 1000) 
importancia <- as.data.frame(modelo_randforest$importance)
importancia <- rownames_to_column(importancia,var = "variable")

p1 <- ggplot(data = importancia, aes(x = reorder(variable, MeanDecreaseAccuracy),
                               y = MeanDecreaseAccuracy,
                               fill = MeanDecreaseAccuracy)) +
      labs(x = "variable", title = "Reducción de Accuracy") +
      geom_col() +
      coord_flip() +
      theme_bw() +
      theme(legend.position = "bottom")

p2 <- ggplot(data = importancia, aes(x = reorder(variable, MeanDecreaseGini),
                               y = MeanDecreaseGini,
                               fill = MeanDecreaseGini)) +
      labs(x = "variable", title = "Reducción de pureza (Gini)") +
      geom_col() +
      coord_flip() +
      theme_bw() +
      theme(legend.position = "bottom")
ggarrange(p1, p2)

Conclusión análisis exploratorio

La exploración de los datos, el estudio de su distribución, y su posible relación con la variable respuesta parecen indicar que los factores que más influyeron en la supervivencia de los pasajeros fueron: el sexo, la clase a la que pertenecían y si tenían o no al menos un hijo a bordo. También se ha detectado que las variables continuas no están correlacionadas y que las variables Age, Cabin y Embarked tienen valores ausentes.

Imputación de valores ausentes

Tal y como se ha identificado en la exploración de datos, las variables Cabin, Age y Embarked contienen 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 cuanta 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.

En el set de datos Titanic, si se eliminan las observaciones incompletas, se pasa de 891 observaciones a solo 183, por lo que esta no es una opción.

nrow(datos)

## [1] 891

nrow(na.omit(datos))

## [1] 183

La variable Cabin está ausente para casi un 80% de las observaciones, con un porcentaje tan alto de valores ausentes, no es conveniente imputarla, se excluye directamente del modelo. La variable Age_grupo se ha obtenido discretizando en intervalos la variable Age, por lo que contiene el mismo número de valores ausentes que esta última. Esto deja al modelo con dos variables que requieren de imputación: Age_grupo y Embarked.

La imputación es un proceso complejo que debe de realizarse con detenimiento, identificando cuidadosamente qué variables son las adecuadas para cada imputación. Para más detalles sobre por qué se realizan de la siguiente forma, ver \(^{\text{Anexo 7}}\).

El paquete recipes permite 4 métodos de imputación distintos:

• step_bagimpute(): imputación vía Bagged Trees.

• step_knnimpute(): imputación vía K-Nearest Neighbors.

• step_meanimpute(): imputación vía media del predictor (predictores continuos).

• step_modeimpute(): imputación vía moda del predictor (predictores cualitativos).

Cabe destacar también los paquetes Hmisc missForest y MICE que permiten aplicar otros métodos.

Se imputa la variable Embarked con el valor C. Como no existe una función step() que haga sustituciones por valores concretos, se realiza una sustitución de forma externa al recipe.

datos_train <- datos_train %>%
               mutate(Embarked = replace(Embarked, is.na(Embarked), "C"))
datos_test <-  datos_test %>%
               mutate(Embarked = replace(Embarked, is.na(Embarked), "C"))

La variable Age_grupo se imputa con el método bagging empleando todos los otros predictores.

library(recipes)
# Se crea un objeto recipe() con la variable respuesta y los predictores. 
# Las variables *PassengerId*, *Name*, *Ticket* no parecen aportar información
# relevante sobre la supervivencia de los pasajeros. Excluyendo todas estas
# variables, se propone como modelo inicial el formado por los predictores:
#  Pclass + Sex + SibSp + Parch + Fare + Embarked + Age_grupo.

objeto_recipe <- recipe(formula = Survived ~ Pclass + Sex + SibSp + Parch +
                                  Fare + Embarked + Age_grupo,
                        data =  datos_train)
objeto_recipe

## Data Recipe
## 
## Inputs:
## 
##       role #variables
##    outcome          1
##  predictor          7

objeto_recipe <- objeto_recipe %>% step_bagimpute(Age_grupo)
objeto_recipe

## Data Recipe
## 
## Inputs:
## 
##       role #variables
##    outcome          1
##  predictor          7
## 
## Operations:
## 
## Bagged tree imputation for Age_grupo

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.

La función nearZeroVar() del paquete caret y step_nzv() del paquete recipe identifican como predictores potencialmente problemáticos aquellos que tienen un único valor (cero varianza) o que cumplen las dos siguientes condiciones:

• Ratio de frecuencias: ratio entre la frecuencia del valor más común y la frecuencia del segundo valor más común. Este ratio tiende a 1 si las frecuencias están equidistribuidas y a valores grandes cuando la frecuencia del valor mayoritario supera por mucho al resto (el denominador es un número decimal pequeño). Valor por defecto freqCut = 95/5.

• Porcentaje de valores únicos: número de valores únicos dividido entre el total de muestras (multiplicado por 100). Este porcentaje se aproxima a cero cuanto mayor es la variedad de valores. Valor por defecto uniqueCut = 10.

datos %>% select(Pclass, Sex, SibSp, Parch, Fare, Embarked, Age_grupo) %>%
          nearZeroVar(saveMetrics = TRUE)

Entre los predictores incluidos en el modelo, no se detecta ninguno con varianza cero o próxima a cero.

objeto_recipe <- objeto_recipe %>% step_nzv(all_predictors())

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

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 = \frac{x - \mu}{\sigma}\]
  • Estandarización max-min: transformar los datos de forma que estén dentro del rango [0, 1]. \[X_{norm} = \frac{X - X_{min}}{X_{max}-X_{min}}\]

Para este análisis se normalizan todas las variables numéricas.

objeto_recipe <- objeto_recipe %>% step_center(all_numeric())
objeto_recipe <- objeto_recipe %>% step_scale(all_numeric())

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

Binarización de 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.

Por defecto, la función step_dummy(all_nominal()) binariza todas las variables almacenadas como tipo factor o character. Además, elimina uno de los niveles para evitar redundancias. Volviendo al ejemplo anterior, no es necesario almacenar las tres variables, ya que, si color_rojo y color_verde toman el valor 0, la variable color_azul toma necesariamente el valor 1. Si color_rojo o color_verde toman el valor 1, entonces color_azul es necesariamente 0.

objeto_recipe <- objeto_recipe %>% step_dummy(all_nominal(), -all_outcomes())

Una vez que se ha creado el objeto recipe con todas las transformaciones de preprocesado, se aprenden con los datos de entrenamiento y se aplican a los dos conjuntos.

# Se entrena el objeto recipe
trained_recipe <- prep(objeto_recipe, training = datos_train)
trained_recipe

## Data Recipe
## 
## Inputs:
## 
##       role #variables
##    outcome          1
##  predictor          7
## 
## Training data contained 714 data points and 145 incomplete rows. 
## 
## Operations:
## 
## Bagged tree imputation for Age_grupo [trained]
## Sparse, unbalanced variable filter removed no terms [trained]
## Centering for SibSp, Parch, Fare [trained]
## Scaling for SibSp, Parch, Fare [trained]
## Dummy variables from Pclass, Sex, Embarked, Age_grupo [trained]

# Se aplican las transformaciones al conjunto de entrenamiento y de test
datos_train_prep <- bake(trained_recipe, new_data = datos_train)
datos_test_prep  <- bake(trained_recipe, new_data = datos_test)

glimpse(datos_train_prep)

## Rows: 714
## Columns: 11
## $ SibSp             <dbl> -0.4679339, -0.4679339, -0.4679339, -0.4679339, 2.4…
## $ Parch             <dbl> -0.4618939, -0.4618939, -0.4618939, -0.4618939, 0.7…
## $ Fare              <dbl> -0.521565287, -0.518757526, -0.509586256, 0.4653626…
## $ Survived          <fct> Si, No, No, No, No, Si, Si, Si, No, No, No, Si, Si,…
## $ Pclass_X2         <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 1, 1, …
## $ Pclass_X3         <dbl> 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 0, 0, 1, 1, 0, 0, …
## $ Sex_male          <dbl> 0, 1, 1, 1, 1, 0, 0, 0, 1, 1, 0, 0, 1, 0, 0, 1, 1, …
## $ Embarked_Q        <dbl> 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ Embarked_S        <dbl> 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, …
## $ Age_grupo_anciano <dbl> 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ Age_grupo_niño    <dbl> 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …

Tras el preprocesado de los datos, se han generado un total de 11 variables (10 predictores y la variable respuesta).

Métodos wrapper

Los métodos wrapper evalúan múltiples modelos, generados mediante la incorporación o eliminación de predictores, con la finalidad de identificar la combinación óptima que consigue maximizar la capacidad del modelo. Pueden entenderse como algoritmos de búsqueda que tratan a los predictores disponibles como valores de entrada y utilizan una métrica del modelo, por ejemplo, su error de predicción, como objetivo de la optimización. El paquete caret incorpora métodos wrapper basados en eliminación recursiva, algoritmos genéticos y simulated annealing.

# ELIMINACIÓN RECURSIVA MEDIANTE RANDOM FOREST Y BOOTSTRAPPING
# =============================================================================

# Se paraleliza el proceso para que sea más rápido. El número de cores debe 
# seleccionarse en función del ordenador que se está empleando.
library(doMC)
registerDoMC(cores = 4)

# Tamaño de los conjuntos de predictores analizados
subsets <- c(3:11)

# Número de resamples para el proceso de bootstrapping
repeticiones <- 30

# Se crea una semilla para cada repetición de validación. Esto solo es necesario si
# se quiere asegurar la reproducibilidad de los resultados, ya que la validación
# cruzada y el bootstrapping implican selección aleatoria.

# El número de semillas necesarias depende del número total de repeticiones: 
# Se necesitan B+1 elementos donde B es el número total de particiones (CV) o
# resampling (bootstrapping). Los primeros B elementos deben ser vectores formados
# por M números enteros, donde M es el número de modelos ajustados, que en este caso
# se corresponde con el número de tamaños. El último elemento solo necesita un único
# número para ajustar el modelo final.
set.seed(123)
seeds <- vector(mode = "list", length = repeticiones + 1)
for (i in 1:repeticiones) {
  seeds[[i]] <- sample.int(1000, length(subsets))
} 
seeds[[repeticiones + 1]] <- sample.int(1000, 1)

# Se crea un control de entrenamiento donde se define el tipo de modelo empleado
# para la selección de variables, en este caso random forest, la estrategia de
# resampling, en este caso bootstrapping con 30 repeticiones, y las semillas para
# cada repetición. Con el argumento returnResamp = "all" se especifica que se
# almacene la información de todos los modelos generados en todas las repeticiones.
ctrl_rfe <- rfeControl(functions = rfFuncs, method = "boot", number = repeticiones,
                       returnResamp = "all", allowParallel = TRUE, verbose = FALSE,
                       seeds = seeds)

# Se ejecuta la eliminación recursiva de predictores
set.seed(342)
rf_rfe <- rfe(Survived ~ ., data = datos_train_prep,
              sizes = subsets,
              metric = "Accuracy",
              # El accuracy es la proporción de clasificaciones correctas
              rfeControl = ctrl_rfe,
              ntree = 500)
# Dentro de rfe() se pueden especificar argumentos para el modelo empleado, por
# ejemplo, el hiperparámetro ntree=500.

# Se muestra una tabla resumen con los resultados
rf_rfe

## 
## Recursive feature selection
## 
## Outer resampling method: Bootstrapped (30 reps) 
## 
## Resampling performance over subset size:
## 
##  Variables Accuracy  Kappa AccuracySD KappaSD Selected
##          3   0.7973 0.5577    0.02728 0.05625         
##          4   0.8199 0.6041    0.02259 0.04918         
##          5   0.8189 0.6031    0.01894 0.04365         
##          6   0.8196 0.6048    0.01942 0.04339         
##          7   0.8223 0.6109    0.01874 0.04255         
##          8   0.8211 0.6068    0.01868 0.04255         
##          9   0.8240 0.6169    0.02015 0.04524         
##         10   0.8290 0.6262    0.02123 0.04838        *
## 
## The top 5 variables (out of 10):
##    Sex_male, Fare, Age_grupo_niño, Pclass_X3, SibSp

# El objeto rf_rfe almacena en optVariables las variables del mejor modelo.
rf_rfe$optVariables

##  [1] "Sex_male"          "Fare"              "Age_grupo_niño"   
##  [4] "Pclass_X3"         "SibSp"             "Embarked_S"       
##  [7] "Parch"             "Pclass_X2"         "Age_grupo_anciano"
## [10] "Embarked_Q"

# Valores de accuracy y kappa para cada tamaño de modelo en cada resample.
rf_rfe$resample %>% select(1, 2, 3, 8) %>% head(8)

# Métricas promedio de cada tamaño
rf_rfe$resample %>% group_by(Variables) %>%
                    summarise(media_accuracy = mean(Accuracy),
                              media_kappa = mean(Kappa)) %>%
                    arrange(desc(media_accuracy))

ggplot(data = rf_rfe$results, aes(x = Variables, y = Accuracy)) +
  geom_line() +
  scale_x_continuous(breaks  = unique(rf_rfe$results$Variables)) +
  geom_point() +
  geom_errorbar(aes(ymin = Accuracy - AccuracySD, ymax = Accuracy + AccuracySD),
                width = 0.2) +
  geom_point(data = rf_rfe$results %>% slice(which.max(Accuracy)),
             color = "red") +
  theme_bw()

head(rf_rfe$variables, 10)

rf_rfe$variables %>% filter(Variables == 10) %>% group_by(var) %>%
                     summarise(media_influencia = mean(Overall),
                               sd_influencia = sd(Overall)) %>%
                     arrange(desc(media_influencia))

# Paralelización
# =============================================================================
library(doMC)
registerDoMC(cores = 4)

# Control de entrenamiento
# =============================================================================
ga_ctrl <- gafsControl(functions = rfGA,
                       method = "cv",
                       number = 5,
                       allowParallel = TRUE,
                       genParallel = TRUE, 
                       verbose = FALSE)

# Selección de predictores
# =============================================================================
set.seed(10)
rf_ga <- gafs(x = datos_train_prep %>% select(-Survived),
              y = datos_train_prep$Survived,
              iters = 10, 
              popSize = 10,
              gafsControl = ga_ctrl,
              ntree = 100)

rf_ga

## 
## Genetic Algorithm Feature Selection
## 
## 714 samples
## 10 predictors
## 2 classes: 'No', 'Si' 
## 
## Maximum generations: 10 
## Population per generation: 10 
## Crossover probability: 0.8 
## Mutation probability: 0.1 
## Elitism: 0 
## 
## Internal performance values: Accuracy, Kappa
## Subset selection driven to maximize internal Accuracy 
## 
## External performance values: Accuracy, Kappa
## Best iteration chose by maximizing external Accuracy 
## External resampling method: Cross-Validated (5 fold) 
## 
## During resampling:
##   * the top 5 selected variables (out of a possible 10):
##     Age_grupo_niño (100%), Fare (100%), Pclass_X3 (100%), Sex_male (100%), Embarked_S (80%)
##   * on average, 8.2 variables were selected (min = 6, max = 10)
## 
## In the final search using the entire training set:
##    * 10 features selected at iteration 2 including:
##      SibSp, Parch, Fare, Pclass_X2, Pclass_X3 ... 
##    * external performance at this iteration is
## 
##    Accuracy       Kappa 
##      0.8221      0.6096

rf_ga$optVariables

##  [1] "SibSp"             "Parch"             "Fare"             
##  [4] "Pclass_X2"         "Pclass_X3"         "Sex_male"         
##  [7] "Embarked_Q"        "Embarked_S"        "Age_grupo_anciano"
## [10] "Age_grupo_niño"

# Accuracy media en cada generación
rf_ga$external %>% group_by(Iter) %>% summarize(accuracy_media = mean(Accuracy))

crear_poblacion <- function(popSize, n_variables, n_max){
  # Esta función crea una matriz binaria en la que el número de 1s por
  # fila no supera un valor máximo.
  # Argumentos:
  #   popSize:     número total de individuos (número de filas de la matriz).
  #   n_variables: longitud de los individuos (número de columnas de la matriz).
  #   n_max:       número máximo de 1 que puede contener un individuo.
  
  # Matriz donde almacenar los individuos generados.
  poblacion <- matrix(data = NA, nrow = popSize, ncol = n_variables)
  
  # Bucle para crear cada individuo.   
  for(i in 1:popSize){
    # Se selecciona con (igual probabilidad ) el número de valores = 1 que puede
    # tener el individuo.
    numero_1s <- sample(x = 1:n_max, size = 1)
    
    # Se crea un vector con todo ceros que representa el individuo.
    individuo <- rep(0, times = n_variables)
    
    # Se sustituyen (numero_1s) posiciones aleatorias por unos.
    individuo[sample(x = 1:n_variables, size = numero_1s)] <- 1
    
    # Se añade el nuevo individuo a la población.
    poblacion[i,] <- individuo
  }
    return(poblacion)
}

poblacion_inicial <- crear_poblacion(popSize = 10,
                                     n_variables = ncol(datos_train_prep) - 1 ,
                                     n_max = 4)
poblacion_inicial

##       [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
##  [1,]    1    0    0    0    0    0    0    0    0     1
##  [2,]    0    0    0    0    0    0    1    1    0     0
##  [3,]    0    0    0    1    0    0    0    0    0     1
##  [4,]    1    0    0    1    0    0    0    0    0     0
##  [5,]    0    0    0    1    0    0    0    0    0     0
##  [6,]    0    0    0    1    0    0    0    0    0     0
##  [7,]    0    1    0    0    0    0    0    1    0     0
##  [8,]    1    0    0    1    0    1    0    1    0     0
##  [9,]    0    0    0    1    0    0    0    0    1     0
## [10,]    0    0    0    1    1    0    0    1    0     0

set.seed(10)
rf_ga <- gafs(x = datos_train_prep %>% select(-Survived),
              y = datos_train_prep$Survived,
              iters = 10, 
              popSize = 10,
              suggestions = poblacion_inicial,
              gafsControl = ga_ctrl,
              ntree = 100)
rf_ga

## 
## Genetic Algorithm Feature Selection
## 
## 714 samples
## 10 predictors
## 2 classes: 'No', 'Si' 
## 
## Maximum generations: 10 
## Population per generation: 10 
## Crossover probability: 0.8 
## Mutation probability: 0.1 
## Elitism: 0 
## 
## Internal performance values: Accuracy, Kappa
## Subset selection driven to maximize internal Accuracy 
## 
## External performance values: Accuracy, Kappa
## Best iteration chose by maximizing external Accuracy 
## External resampling method: Cross-Validated (5 fold) 
## 
## During resampling:
##   * the top 5 selected variables (out of a possible 10):
##     Pclass_X2 (100%), Sex_male (100%), Age_grupo_niño (80%), Embarked_S (80%), SibSp (60%)
##   * on average, 4.8 variables were selected (min = 4, max = 6)
## 
## In the final search using the entire training set:
##    * 5 features selected at iteration 8 including:
##      SibSp, Parch, Pclass_X2, Pclass_X3, Age_grupo_niño  
##    * external performance at this iteration is
## 
##    Accuracy       Kappa 
##      0.8067      0.5836

Métodos de filtrado

Los métodos basados en filtrado evalúan la relevancia de los predictores fuera del modelo para, posteriormente, incluir únicamente aquellos que pasan un determinado criterio. Se trata por lo tanto de analizar la relación que tiene cada predictor con la variable respuesta. Por ejemplo, en problemas de clasificación con predictores continuos, se puede aplicar un ANOVA a cada predictor para identificar aquellos que varían dependiendo de la variable respuesta. Finalmente, se incorporan al modelo aquellos predictores con un p-value inferior a un determinado límite o los n mejores. Al igual que con los métodos wrapper, para evitar que la selección esté excesivamente influenciada por los datos de entrenamiento (overfitting), es conveniente repetir el proceso varias veces mediante validación cruzada o bootstrapping.

Al igual que con los métodos de eliminación recursiva, caret incorpora modelos predefinidos para la evaluación de los predictores filtrados (lmSBF, rfSBF, treebagSBF, ldaSBF y nbSBF), además, ofrece la de la posibilidad de emplear cualquier otro modelo (caretFuns). \(^{\text{Anexo 2}}\)

Los test estadísticos que sbf() emplea por defecto para cuantificar la relación entre los predictores y la variable respuesta son: ANOVA (cuando la variable respuesta es cualitativa) y GAMS (cuando la variable respuesta es continua). En ambos casos, el límite para la selección es p-value = 0.05. \(^{\text{Anexo 3}}\)

# FILTRADO DE PREDICTORES MEDIANTE ANOVA, RANDOM FOREST Y CV-REPETIDA
# =============================================================================

# Se paraleliza para que sea más rápido
library(doMC)
registerDoMC(cores = 4)

# Se crea una semilla para cada partición y cada repetición: el vector debe
# tener B+1 semillas donde B = particiones * repeticiones.
particiones = 10
repeticiones = 5
set.seed(123)
seeds <- sample.int(1000, particiones * repeticiones + 1)

# Control del filtrado
ctrl_filtrado <- sbfControl(functions = rfSBF, method = "repeatedcv",
                            number = particiones, repeats = repeticiones,
                            seeds = seeds, verbose = FALSE, 
                            saveDetails = TRUE, allowParallel = TRUE)
set.seed(234)
rf_sbf <- sbf(Survived ~ ., data = datos_train_prep,
              sbfControl = ctrl_filtrado,
              # argumentos para el modelo de evaluación
              ntree = 500)

rf_sbf

## 
## Selection By Filter
## 
## Outer resampling method: Cross-Validated (10 fold, repeated 5 times) 
## 
## Resampling performance:
## 
##  Accuracy  Kappa AccuracySD KappaSD
##    0.8199 0.5986    0.04298 0.09951
## 
## Using the training set, 8 variables were selected:
##    Parch, Fare, Pclass_X2, Pclass_X3, Sex_male...
## 
## During resampling, the top 5 selected variables (out of a possible 8):
##    Age_grupo_anciano (100%), Age_grupo_niño (100%), Embarked_S (100%), Fare (100%), Pclass_X3 (100%)
## 
## On average, 7.6 variables were selected (min = 6, max = 8)

rf_sbf$optVariables

## [1] "Parch"             "Fare"              "Pclass_X2"        
## [4] "Pclass_X3"         "Sex_male"          "Embarked_S"       
## [7] "Age_grupo_anciano" "Age_grupo_niño"

Empleando la selección por filtrado, se identifican como óptimos 8 predictores: Parch, Fare, Pclass_X2, Pclass_X3, Sex_male, Embarked_S, Age_grupo_anciano y Age_grupo_niño. El porcentaje que acompaña a las top 5 selected variables se corresponde con el porcentaje de repeticiones (resamples) en las que la variable ha sido seleccionada.

Se almacenan en una variable los predictores filtrados para emplearlos posteriormente en el ajuste de los modelos.

predictores_filtrados <- rf_sbf$optVariables

Comparación de métodos

Ambas estrategias, wrapper y filtrado, tienen ventajas y desventajas. Los métodos de filtrado son computacionalmente más rápidos por lo que suelen ser la opción factible cuando hay cientos o miles de predictores, sin embargo, el criterio de selección no está directamente relacionado con la efectividad del modelo. Además, en la mayoría de casos, los métodos de filtrado evalúan cada predictor de forma individual, por lo que no contemplan interacciones y pueden incorporar predictores redundantes (correlacionados). Los métodos wrapper, además de ser computacionalmente más costosos, para evitar overfitting, necesitan recurrir a validación cruzada o bootstrapping, por lo que requieren un número alto de observaciones. A pesar de ello, si se cumplen las condiciones, suelen conseguir una mejor selección.

Entrenamiento del modelo

Todos los modelos incorporados en el paquete caret se entrenan con la función train(). Entre los argumentos de esta función destacan:

• formula: la fórmula del modelo que se quiere crear.

• x, y: en lugar de una fórmula, se pueden pasar por separado los valores de los predictores y de la variable respuesta.

• method: el nombre del algoritmo que se desea emplear (listado).

• metric: las métricas empleadas para evaluar la capacidad predictiva del modelo. Por defecto, se emplea accuracy para problemas de clasificación y RMSE para regresión \(^{\text{Anexo 5}}\).

• trControl: especificaciones adicionales sobre la forma de llevar a cabo el entrenamiento del modelo.

• ...: argumentos propios del algoritmo empleado.
  
  

En primer lugar, se ajusta un modelo basado en una máquina vector soporte lineal que prediga la supervivencia del pasajero en función de todos los predictores disponibles tras el preprocesado. A excepción del nombre del algoritmo, todos los demás argumentos de la función train() y los propios de svmLinear se dejan por defecto.

# Predictores: Survived, SibSp, Parch, Fare, Pclass_X2, Pclass_X3, Sex_male, Embarked_Q,
#              Embarked_S, Age_grupo_anciano, Age_grupo_niño
modelo_svmlineal <- train(Survived ~ ., method = "svmLinear", data = datos_train_prep)

El objeto devuelto por train() contiene múltiples elementos y mucha información. Por el momento, se analiza únicamente $finalModel, que es el modelo final obtenido tras el ajuste.

modelo_svmlineal$finalModel

## Support Vector Machine object of class "ksvm" 
## 
## SV type: C-svc  (classification) 
##  parameter : cost C = 1 
## 
## Linear (vanilla) kernel function. 
## 
## Number of Support Vectors : 315 
## 
## Objective Function Value : -296.459 
## Training error : 0.207283

Al imprimir el objeto $finalModel se muestra el tipo de modelo creado, el valor de los hiperparámetros (más detalles sobre esto en la siguiente sección) e información adicional, entre ella, el error de entrenamiento (training error). Este error se corresponde con el error que comete el modelo al predecir las mismas observaciones con las que se ha entrenado, un 20% en este caso.

Evaluación del modelo mediante resampling

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 en $finalModel se corresponde con el error que comete el modelo al predecir observaciones que sí ha “visto”, por lo tanto, 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. caret incorpora los métodos boot, boot632, optimism_boot, boot_all, cv, repeatedcv, LOOCV, LGOCV \(^{\text{Anexo 4}}\). Cada uno funciona internamente de forma distinta, pero todos ellos se basan en la idea de ajustar y evaluar el modelo múltiples veces con distintos subconjuntos creados a partir de los datos de entrenamiento, 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.

Para especificar el tipo de validación, así como el número de repeticiones, se crea un control de entrenamiento mediante la función trainControl(), que se pasa al argumento trControl de la función train(). 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. caret permite paralelizar el proceso para que sea más rápido y establecer semillas para asegurar que cada resampling o partición pueda crearse de nuevo con exactamente las mismas observaciones.

Se ajusta de nuevo una máquina vector soporte lineal, esta vez con validación cruzada repetida para estimar su error.

# PARALELIZACIÓN DE PROCESO
#===============================================================================
# Se paraleliza para que sea más rápido. El número de cores activados depende del
# las características del ordenador donde se ejecute el código.
library(doMC)
registerDoMC(cores = 4)

# NÚMERO DE REPETICIONES Y SEMILLAS PARA CADA REPETICIÓN
#===============================================================================
# Es este caso se recurre a validación cruzada repetida como método de validación.
# Número de particiones y repeticiones
particiones  <- 10
repeticiones <- 5

# Las semillas solo son necesarias si se quiere asegurar la reproducibilidad de
# los resultados, ya que la validación cruzada y el bootstrapping implican selección
# aleatoria. Las semillas se almacenan en una lista con B+1 elementos donde B depende
# del método de validación empleado:
# 
#   "cv": B es el número total de particiones
#   "repeatedcv": B es el número de particiones por el número de repeticiones.
#   "boot": B es el número de resamples.
#   "LOGOCV":   B es el número de repeticiones.
#   
# Los primeros B elementos deben ser vectores formados por M números enteros,
# donde M es el número de modelos ajustados en cada partición o repetición,
# es decir, el total de hiperparámetros comparados. El último elemento (B+1) solo
# necesita un único número para ajustar el modelo final.

set.seed(123)
seeds <- vector(mode = "list", length = (particiones * repeticiones) + 1)
for (i in 1:(particiones * repeticiones)) {
  # Cada elemento de la lista, excepto el último tiene que tener tantas semillas
  # como hiperparámetros analizados. En este caso, se emplea el valor por 
  # defecto C = 1, por lo que cada elemento de seeds está formada por un
  # único valor.
  seeds[[i]] <- sample.int(1000, 1) 
}
# La última semilla se emplea para ajustar el modelo final con todas las observaciones.
seeds[[(particiones * repeticiones) + 1]] <- sample.int(1000, 1)

# DEFINICIÓN DEL ENTRENAMIENTO
#===============================================================================
control_train <- trainControl(method = "repeatedcv", number = particiones,
                              repeats = repeticiones, seeds = seeds,
                              returnResamp = "all", verboseIter = FALSE,
                              allowParallel = TRUE)

# AJUSTE DEL MODELO
# ==============================================================================
set.seed(342)
modelo_svmlineal <- train(Survived ~ ., data = datos_train_prep,
                          method = "svmLinear",
                          metric = "Accuracy",
                          trControl = control_train)
modelo_svmlineal

## Support Vector Machines with Linear Kernel 
## 
## 714 samples
##  10 predictor
##   2 classes: 'No', 'Si' 
## 
## No pre-processing
## Resampling: Cross-Validated (10 fold, repeated 5 times) 
## Summary of sample sizes: 642, 643, 643, 643, 643, 642, ... 
## Resampling results:
## 
##   Accuracy  Kappa    
##   0.792766  0.5554307
## 
## Tuning parameter 'C' was held constant at a value of 1

Una validación cruzada con 10 particiones y 5 repeticiones, implica ajustar y evaluar el modelo 10 x 5 = 50 veces, cada vez con una partición distinta, más un último ajuste con todos los datos de entrenamiento para crear el modelo final. Al imprimir el objeto devuelto por train() se muestra información sobre el algoritmo empleado, el método de resampling y el valor promedio de las métricas de validación, en este caso, al ser un problema de clasificación, se emplea por defecto accuracy. Además, si se indica en el control de entrenamiento returnResamp = "all", se almacenan los resultados obtenidos en cada una de las iteraciones, lo que permite un análisis más detallado. Véase la distribución de los valores de accuracy.

# Valores de validación (Accuracy y Kappa) obtenidos en cada partición y repetición.
modelo_svmlineal$resample %>% head(10)

summary(modelo_svmlineal$resample$Accuracy)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.6806  0.7639  0.8028  0.7928  0.8188  0.8889

p1 <- ggplot(data = modelo_svmlineal$resample, aes(x = Accuracy)) +
      geom_density(alpha = 0.5, fill = "gray50") +
      geom_vline(xintercept = mean(modelo_svmlineal$resample$Accuracy),
                 linetype = "dashed") +
      theme_bw() 
p2 <- ggplot(data = modelo_svmlineal$resample, aes(x = 1, y = Accuracy)) +
      geom_boxplot(outlier.shape = NA, alpha = 0.5, fill = "gray50") +
      geom_jitter(width = 0.05) +
      labs(x = "") +
      theme_bw() +
      theme(axis.text.x = element_blank(), axis.ticks.x = element_blank())
  
final_plot <- ggarrange(p1, p2)
final_plot <- annotate_figure(
                final_plot,
                top = text_grob("Accuracy obtenido en la validación", size = 15))
final_plot

El accuracy promedio estimado mediante validación cruzada repetida es de 0.79, el modelo predice correctamente la supervivencia de los pasajeros un 79% de las veces. Este valor será contrastado más adelante cuando se calcule el accuracy del modelo con el conjunto de test.

Optimización de hiperparámetros (parameter tuning)

Muchos modelos, entre ellos 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 \(^{\text{Anexo 4}}\).

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

El modelo svmLinear empleado hasta ahora tienen un hiperparámetro llamado coste (C). 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 con el riesgo de overfitting. caret permite explorar diferentes valores de hiperparámetros en el momento de entrenar el modelo. Los valores de los hiperparámetros se pasan en formato de dataframe al argumento tuneGrid de la función train().

Se vuelve a ajustar un modelo svmLinear con diferentes valores de C y se emplea validación cruzada repetida para identificar con cuál de ellos se obtienen mejores resultados.

# PARALELIZACIÓN DE PROCESO
#===============================================================================
library(doMC)
registerDoMC(cores = 4)

# HIPERPARÁMETROS, NÚMERO DE REPETICIONES Y SEMILLAS PARA CADA REPETICIÓN
#===============================================================================
particiones  <- 10
repeticiones <- 5

hiperparametros <- data.frame(C = c(0.001, 0.01, 0.1, 0.5, 1, 10))

set.seed(123)
seeds <- vector(mode = "list", length = (particiones * repeticiones) + 1)
for (i in 1:(particiones * repeticiones)) {
  seeds[[i]] <- sample.int(1000, nrow(hiperparametros)) 
}
seeds[[(particiones * repeticiones) + 1]] <- sample.int(1000, 1)

# DEFINICIÓN DEL ENTRENAMIENTO
#===============================================================================
control_train <- trainControl(method = "repeatedcv", number = particiones,
                              repeats = repeticiones, seeds = seeds,
                              returnResamp = "all", verboseIter = FALSE,
                              allowParallel = TRUE)

# AJUSTE DEL MODELO
# ==============================================================================
set.seed(342)
modelo_svmlineal <- train(Survived ~ ., data = datos_train_prep,
                          method = "svmLinear",
                          tuneGrid = hiperparametros,
                          metric = "Accuracy",
                          trControl = control_train)
modelo_svmlineal

## Support Vector Machines with Linear Kernel 
## 
## 714 samples
##  10 predictor
##   2 classes: 'No', 'Si' 
## 
## No pre-processing
## Resampling: Cross-Validated (10 fold, repeated 5 times) 
## Summary of sample sizes: 642, 643, 643, 643, 643, 642, ... 
## Resampling results across tuning parameters:
## 
##   C      Accuracy   Kappa     
##   1e-03  0.6308451  0.04826884
##   1e-02  0.7927660  0.55543074
##   1e-01  0.7927660  0.55543074
##   5e-01  0.7927660  0.55543074
##   1e+00  0.7927660  0.55543074
##   1e+01  0.7927660  0.55543074
## 
## Accuracy was used to select the optimal model using the largest value.
## The final value used for the model was C = 0.01.

Como se ha empleado validación cruzada con 10 particiones y 5 repeticiones, caret ajusta el modelo 10 x 5 = 50 veces por cada valor de C. Promediando los 50 valores de accuracy obtenidos para cada valor del hiperparámetro, se identifica cual es el mejor (valores mostrados en la tabla). Finalmente, se reajusta el modelo empleando todas las observaciones de entrenamiento y el mejor valor del hiperparámetro, este modelo se almacena en $finalModel. En este caso, de entre los valores de C probados, C = 0.01 consigue los mejores resultados con un accuaracy de 0.7927. Si varios valores consiguen el mismo resultado, se selecciona el que da lugar al modelo más sencillo.

Si se especifica en el control de entrenamiento returnResamp = "all", se almacenan los resultados obtenidos en cada iteración para cada valor del hiperparámetro, lo que permite un análisis más detallado.

modelo_svmlineal$resample %>% head()

ggplot(data = modelo_svmlineal$resample,
       aes(x = as.factor(C), y = Accuracy, color = as.factor(C))) +
  geom_boxplot(outlier.shape = NA, alpha = 0.6) +
  geom_jitter(width = 0.2, alpha = 0.6) +
  # Línea horizontal en el accuracy basal
  geom_hline(yintercept = 0.62, linetype = "dashed") +
  labs(x = "C") +
  theme_bw() + theme(legend.position = "none")

Además, caret incorpora un método basado en ggplot2 para visualizar la evolución de los modelos según el valor de los hiperparámetros.

ggplot(modelo_svmlineal, highlight = TRUE) +
  labs(title = "Evolución del accuracy del modelo en función de C") +
  theme_bw()

# PARALELIZACIÓN DE PROCESO
#===============================================================================
library(doMC)
registerDoMC(cores = 4)

# HIPERPARÁMETROS, NÚMERO DE REPETICIONES Y SEMILLAS PARA CADA REPETICIÓN
#===============================================================================
particiones  <- 10
repeticiones <- 5

# Hiperparámetros aleatorios
set.seed(123)
hiperparametros <- data.frame(C = runif(n = 5, min = 0.001, max = 20))

set.seed(123)
seeds <- vector(mode = "list", length = (particiones * repeticiones) + 1)
for (i in 1:(particiones * repeticiones)) {
  seeds[[i]] <- sample.int(1000, nrow(hiperparametros)) 
}
seeds[[(particiones * repeticiones) + 1]] <- sample.int(1000, 1)

# DEFINICIÓN DEL ENTRENAMIENTO
#===============================================================================
control_train <- trainControl(method = "repeatedcv", number = particiones,
                              repeats = repeticiones, seeds = seeds,
                              returnResamp = "all", verboseIter = FALSE,
                              allowParallel = TRUE)

# AJUSTE DEL MODELO
# ==============================================================================
set.seed(342)
modelo_svmlineal_random <- train(Survived ~ ., data = datos_train_prep,
                               method = "svmLinear",
                               tuneGrid = hiperparametros,
                               metric = "Accuracy",
                               trControl = control_train)

modelo_svmlineal_random

## Support Vector Machines with Linear Kernel 
## 
## 714 samples
##  10 predictor
##   2 classes: 'No', 'Si' 
## 
## No pre-processing
## Resampling: Cross-Validated (10 fold, repeated 5 times) 
## Summary of sample sizes: 642, 643, 643, 643, 643, 642, ... 
## Resampling results across tuning parameters:
## 
##   C          Accuracy  Kappa    
##    5.752263  0.792766  0.5554307
##    8.180129  0.792766  0.5554307
##   15.766314  0.792766  0.5554307
##   17.660465  0.792766  0.5554307
##   18.809405  0.792766  0.5554307
## 
## Accuracy was used to select the optimal model using the largest value.
## The final value used for the model was C = 5.752263.

# PARALELIZACIÓN DE PROCESO
#===============================================================================
library(doMC)
registerDoMC(cores = 4)

# HIPERPARÁMETROS, NÚMERO DE REPETICIONES Y SEMILLAS PARA CADA REPETICIÓN
#===============================================================================
particiones  <- 10
repeticiones <- 5

# Número de hiperparámetros aleatorios
n_hiperparametros <- 10

set.seed(123)
seeds <- vector(mode = "list", length = (particiones * repeticiones) + 1)
for (i in 1:(particiones * repeticiones)) {
  seeds[[i]] <- sample.int(1000, n_hiperparametros)
}
seeds[[(particiones * repeticiones) + 1]] <- sample.int(1000, 1)

# DEFINICIÓN DEL ENTRENAMIENTO
#===============================================================================
control_train <- trainControl(method = "adaptive_cv", number = particiones,
                              repeats = repeticiones,
                              adaptive = list(min = 5, alpha = 0.05, 
                                             method = "gls", complete = TRUE),
                              returnResamp = "all", verboseIter = FALSE,
                              allowParallel = TRUE, search = "random")

# AJUSTE DEL MODELO
# ==============================================================================
set.seed(342)
modelo_svmlineal_adaptative <- train(Survived ~ ., data = datos_train_prep,
                                 method = "svmLinear", trControl = control_train,
                                 tuneLength = 5)

modelo_svmlineal_adaptative

## Support Vector Machines with Linear Kernel 
## 
## 714 samples
##  10 predictor
##   2 classes: 'No', 'Si' 
## 
## No pre-processing
## Resampling: Adaptively Cross-Validated (10 fold, repeated 5 times) 
## Summary of sample sizes: 642, 643, 643, 643, 643, 642, ... 
## Resampling results across tuning parameters:
## 
##   C            Accuracy   Kappa      Resamples
##     0.2060751  0.7942816  0.5592503  49       
##     2.7584259  0.7942816  0.5592503  49       
##    22.9018044  0.7942816  0.5592503  49       
##   337.3203713  0.7949210  0.5555385  11       
##   357.7548716  0.7935994  0.5567193  50       
## 
## Accuracy was used to select the optimal model using the largest value.
## The final value used for the model was C = 357.7549.

Predicción

Una vez que el modelo ha sido ajustado, con la función predict() del paquete caret, se pueden predecir nuevos datos utilizando el objeto devuelto por train(). En concreto, la función predict() utiliza el modelo almacenado en $finalModel, que es el mejor modelo de entre los ajustados en el proceso de optimización. Los argumentos de la función son:

• models: un modelo o lista de modelos creados mediante train().

• newdata: un dateframe con nuevas observaciones.

• type: tipo de predicción. En modelos de clasificación puede ser “raw” para obtener directamente la clase predicha o “prob” para obtener la probabilidad de cada clase. Algunos algoritmos no incorporan de forma nativa un cálculo de probabilidades, para asegurar que caret las calcula, se debe indicar en el control de entrenamiento classProbs = TRUE.

predicciones_raw <- predict(modelo_svmlineal, newdata = datos_test_prep,
                            type = "raw")
predicciones_raw

##   [1] No Si Si Si No No Si No Si Si Si No Si No Si No No Si Si Si No Si No No Si
##  [26] No No Si No No Si Si Si No No No No No Si Si No No No Si Si No Si No No No
##  [51] No No Si No No No Si No No No No Si No Si No No Si Si No No No Si Si Si No
##  [76] Si No No No Si No No No No Si No Si No No No No No No No No No Si Si No No
## [101] Si No No Si No No No No No No Si No No No Si No No No No No Si No No No Si
## [126] No Si No No No No No No No No No Si No Si No No Si Si No No No No No No No
## [151] Si No Si No Si No Si No No No Si Si No No No No Si No No Si Si Si Si No No
## [176] No No
## Levels: No Si

# El algoritmo svmLinear no calcula de forma nativa probabilidades, para obtenerlas
# se reajusta el modelo indicando `classProbs = TRUE`.
particiones  <- 10
repeticiones <- 5
hiperparametros <- expand.grid(C = c(1))

set.seed(123)
seeds <- vector(mode = "list", length = (particiones * repeticiones) + 1)
for (i in 1:(particiones * repeticiones)) {
  seeds[[i]] <- sample.int(1000, nrow(hiperparametros)) 
}
seeds[[(particiones * repeticiones) + 1]] <- sample.int(1000, 1)

control_train <- trainControl(method = "repeatedcv", number = particiones,
                              repeats = repeticiones, seeds = seeds,
                              returnResamp = "all", verboseIter = FALSE,
                              classProbs = TRUE, allowParallel = TRUE)

set.seed(342)
modelo_svmlineal <- train(Survived ~ ., data = datos_train_prep,
                          method = "svmLinear",
                          tuneGrid = hiperparametros,
                          metric = "Accuracy",
                          trControl = control_train)

predicciones_prob <- predict(modelo_svmlineal, newdata = datos_test_prep,
                             type = "prob")
predicciones_prob %>% head()

Obtener las probabilidades para cada clase es más informativo, ya que, además de conocer la clase predicha (la que tiene mayor probabilidad), se puede cuantificar la confianza de dicha predicción. Se tiene mucha más seguridad de que la predicción es correcta si las probabilidades son 0.9-0.1, que si son 0.49-0.51. Para problemas de clasificación binaria, se suele emplear un cutoff de probabilidad del 0.5, aquella clase que supera este límite es a la que se asigna la clasificación. Sin embargo, en determinados escenarios en los que se quiere minimizar el riesgo de falsos positivos o falsos negativos, este cutoff puede ser modificado. Esta es otra ventaja notable de predecir probabilidades en lugar de clases.

Tal y como se describe en secciones posteriores, caret permite obtener las predicciones de un listado de modelos a la vez. En estos casos, es preferible emplear la función extractPrediction(), ya que devuelve los resultados en formato de dataframe en lugar de list.

predicciones <- extractPrediction(
                  models = list(svm = modelo_svmlineal),
                  testX = datos_test_prep %>% select(-Survived),
                  testY = datos_test_prep$Survived
                  )
predicciones %>% head()

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étricas u otras, la función confusionMatrix() calcula las principales. También hay funciones específicas para cada métrica: sensitivity(), specificity(), posPredValue(), negPredValue(), precision(), recall(), y F_meas().

confusionMatrix(data = predicciones_raw, reference = datos_test_prep$Survived,
                positive = "Si")

## Confusion Matrix and Statistics
## 
##           Reference
## Prediction No Si
##         No 92 25
##         Si 17 43
##                                           
##                Accuracy : 0.7627          
##                  95% CI : (0.6931, 0.8233)
##     No Information Rate : 0.6158          
##     P-Value [Acc > NIR] : 2.393e-05       
##                                           
##                   Kappa : 0.4872          
##                                           
##  Mcnemar's Test P-Value : 0.2801          
##                                           
##             Sensitivity : 0.6324          
##             Specificity : 0.8440          
##          Pos Pred Value : 0.7167          
##          Neg Pred Value : 0.7863          
##              Prevalence : 0.3842          
##          Detection Rate : 0.2429          
##    Detection Prevalence : 0.3390          
##       Balanced Accuracy : 0.7382          
##                                           
##        'Positive' Class : Si              
## 

En el apartado de optimización, se estimó, mediante validación cruzada repetida, que el accuracy de modelo svmLinear con un C=0.01 era de 0.7927, es decir, el modelo predice correctamente el 79.27% de las observaciones, un valor casi muy próximo al obtenido con el conjunto de test (0.7627).

En problemas de clasificación, se emplea con frecuencia el porcentaje de predicciones erróneas (1-accuracy), también llamado error de clasificación. Esta medida no se muestra en el resultado de confusionMatrix(), pero que es muy sencilla de calcular.

# Error de test
error_test <- mean(predicciones_raw != datos_test_prep$Survived)
paste("El error de test del modelo:", round(error_test*100, 2), "%")

## [1] "El error de test del modelo: 23.73 %"

K-Nearest Neighbor (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.

Ajuste, optimización y validación del modelo

El método knn de caret emplea la función knn3() con un único hiperparámetro:

• k: número de observaciones vecinas empleadas.

# PARALELIZACIÓN DE PROCESO
#===============================================================================
library(doMC)
registerDoMC(cores = 4)

# HIPERPARÁMETROS, NÚMERO DE REPETICIONES Y SEMILLAS PARA CADA REPETICIÓN
#===============================================================================
particiones  <- 10
repeticiones <- 5

# Hiperparámetros
hiperparametros <- data.frame(k = c(1, 2, 5, 10, 15, 20, 30, 50))

set.seed(123)
seeds <- vector(mode = "list", length = (particiones * repeticiones) + 1)
for (i in 1:(particiones * repeticiones)) {
  seeds[[i]] <- sample.int(1000, nrow(hiperparametros)) 
}
seeds[[(particiones * repeticiones) + 1]] <- sample.int(1000, 1)

# DEFINICIÓN DEL ENTRENAMIENTO
#===============================================================================
control_train <- trainControl(method = "repeatedcv", number = particiones,
                              repeats = repeticiones, seeds = seeds,
                              returnResamp = "final", verboseIter = FALSE,
                              allowParallel = TRUE)

# AJUSTE DEL MODELO
# ==============================================================================
set.seed(342)
modelo_knn <- train(Survived ~ ., data = datos_train_prep,
                    method = "knn",
                    tuneGrid = hiperparametros,
                    metric = "Accuracy",
                    trControl = control_train)
modelo_knn

## k-Nearest Neighbors 
## 
## 714 samples
##  10 predictor
##   2 classes: 'No', 'Si' 
## 
## No pre-processing
## Resampling: Cross-Validated (10 fold, repeated 5 times) 
## Summary of sample sizes: 642, 643, 643, 643, 643, 642, ... 
## Resampling results across tuning parameters:
## 
##   k   Accuracy   Kappa    
##    1  0.7834977  0.5392593
##    2  0.7916784  0.5601159
##    5  0.7863576  0.5457335
##   10  0.7888419  0.5505144
##   15  0.7857864  0.5433513
##   20  0.7700900  0.5020490
##   30  0.7667371  0.4857700
##   50  0.7440258  0.4255430
## 
## Accuracy was used to select the optimal model using the largest value.
## The final value used for the model was k = 2.

# REPRESENTACIÓN GRÁFICA
# ==============================================================================
ggplot(modelo_knn, highlight = TRUE) +
  scale_x_continuous(breaks = hiperparametros$k) +
  labs(title = "Evolución del accuracy del modelo KNN", x = "K") +
  theme_bw()

Con un modelo kNN con k=2 se consigue un accuracy de validación promedio del 79%.

Naive Bayes

Empleando el teorema de Bayes, el algoritmo de Naive Bayes calcula las probabilidades condicionales de que una observación pertenezca a cada una de las clases dadas unas evidencias (valores de los predictores). El término naive (ingenuo en inglés) se debe a que el algoritmo asume que las variables son independientes, es decir, que el valor que toman no está influenciado por las demás. Aunque en la práctica, esta asunción raramente es cierta, permite calcular la probabilidad cuando hay múltiples predictores simplemente multiplicando las probabilidades individuales de cada uno (regla de eventos independientes). A pesar de que es solo una aproximación, puede resultar muy efectiva. Más información sobre Naive Bayes en este link.

Ajuste, optimización y validación del modelo

El método nb de caret emplea la función NaiveBayes() del paquete klaR con tres hiperparámetros:

• usekernel: TRUE para emplear un kernel que estime la densidad o FALSE para asumir una distribución de densidad gaussiana.

• fL: factor de corrección de Laplace, 0 para no aplicar ninguna corrección.

• adjust: parámetro pasado a la función density si usekernel = TRUE.

# PARALELIZACIÓN DE PROCESO
#===============================================================================
library(doMC)
registerDoMC(cores = 4)

# HIPERPARÁMETROS, NÚMERO DE REPETICIONES Y SEMILLAS PARA CADA REPETICIÓN
#===============================================================================
particiones  <- 10
repeticiones <- 5

# Hiperparámetros
hiperparametros <- data.frame(usekernel = FALSE, fL = 0 , adjust = 0)

set.seed(123)
seeds <- vector(mode = "list", length = (particiones * repeticiones) + 1)
for (i in 1:(particiones * repeticiones)) {
  seeds[[i]] <- sample.int(1000, nrow(hiperparametros))
}
seeds[[(particiones * repeticiones) + 1]] <- sample.int(1000, 1)

# DEFINICIÓN DEL ENTRENAMIENTO
#===============================================================================
control_train <- trainControl(method = "repeatedcv", number = particiones,
                              repeats = repeticiones, seeds = seeds,
                              returnResamp = "final", verboseIter = FALSE,
                              allowParallel = TRUE)

# AJUSTE DEL MODELO
# ==============================================================================
set.seed(342)
modelo_nb <- train(Survived ~ ., data = datos_train_prep,
                   method = "nb",
                   tuneGrid = hiperparametros,
                   metric = "Accuracy",
                   trControl = control_train)
modelo_nb

## Naive Bayes 
## 
## 714 samples
##  10 predictor
##   2 classes: 'No', 'Si' 
## 
## No pre-processing
## Resampling: Cross-Validated (10 fold, repeated 5 times) 
## Summary of sample sizes: 642, 643, 643, 643, 643, 642, ... 
## Resampling results:
## 
##   Accuracy   Kappa    
##   0.7877387  0.5568725
## 
## Tuning parameter 'fL' was held constant at a value of 0
## Tuning
##  parameter 'usekernel' was held constant at a value of FALSE
## Tuning
##  parameter 'adjust' was held constant at a value of 0

Empleando un modelo Naive Bayes con usekernel = FALSE, fL = 0, adjust = 0, se consigue un accuracy promedio de validación del 79%.

Regresión logística

Información detallada sobre regresión logística en Regresión logística simple y múltiple.

Ajuste, optimización y validación del modelo

El método glm de caret emplea la función glm() del paquete básico de R. Este algoritmo no tiene ningún hiperparámetro pero, para que efectúe una regresión logística, hay que indicar family = "binomial".

# PARALELIZACIÓN DE PROCESO
#===============================================================================
library(doMC)
registerDoMC(cores = 4)

# HIPERPARÁMETROS, NÚMERO DE REPETICIONES Y SEMILLAS PARA CADA REPETICIÓN
#===============================================================================
particiones  <- 10
repeticiones <- 5

# Hiperparámetros
hiperparametros <- data.frame(parameter = "none")

set.seed(123)
seeds <- vector(mode = "list", length = (particiones * repeticiones) + 1)
for (i in 1:(particiones * repeticiones)) {
  seeds[[i]] <- sample.int(1000, nrow(hiperparametros))
}
seeds[[(particiones * repeticiones) + 1]] <- sample.int(1000, 1)

# DEFINICIÓN DEL ENTRENAMIENTO
#===============================================================================
control_train <- trainControl(method = "repeatedcv", number = particiones,
                              repeats = repeticiones, seeds = seeds,
                              returnResamp = "final", verboseIter = FALSE,
                              allowParallel = TRUE)

# AJUSTE DEL MODELO
# ==============================================================================
set.seed(342)
modelo_logistic <- train(Survived ~ ., data = datos_train_prep,
                         method = "glm",
                         tuneGrid = hiperparametros,
                         metric = "Accuracy",
                         trControl = control_train,
                         family = "binomial")
modelo_logistic

## Generalized Linear Model 
## 
## 714 samples
##  10 predictor
##   2 classes: 'No', 'Si' 
## 
## No pre-processing
## Resampling: Cross-Validated (10 fold, repeated 5 times) 
## Summary of sample sizes: 642, 643, 643, 643, 643, 642, ... 
## Resampling results:
## 
##   Accuracy   Kappa    
##   0.7955829  0.5594967

summary(modelo_logistic$finalModel)

## 
## Call:
## NULL
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.6675  -0.6374  -0.4264   0.5978   2.4311  
## 
## Coefficients:
##                   Estimate Std. Error z value Pr(>|z|)    
## (Intercept)        2.53371    0.35004   7.238 4.55e-13 ***
## SibSp             -0.34097    0.15045  -2.266   0.0234 *  
## Parch             -0.19906    0.11099  -1.793   0.0729 .  
## Fare               0.09335    0.15828   0.590   0.5554    
## Pclass_X2         -0.71592    0.33390  -2.144   0.0320 *  
## Pclass_X3         -1.87473    0.32953  -5.689 1.28e-08 ***
## Sex_male          -2.78207    0.22850 -12.176  < 2e-16 ***
## Embarked_Q        -0.13009    0.43357  -0.300   0.7642    
## Embarked_S        -0.43086    0.26422  -1.631   0.1030    
## Age_grupo_anciano -0.85097    0.62372  -1.364   0.1725    
## Age_grupo_niño     1.74454    0.43994   3.965 7.33e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 950.86  on 713  degrees of freedom
## Residual deviance: 623.80  on 703  degrees of freedom
## AIC: 645.8
## 
## Number of Fisher Scoring iterations: 5

Empleando un modelo de regresión logística se consigue un accuracy promedio de validación del 79%.

LDA

Información detallada sobre LDA en Análisis discriminante lineal (LDA).

Ajuste, optimización y validación del modelo

El método lda de caret emplea la función lda() del paquete MASS. Este algoritmo no tiene ningún hiperparámetro.

# PARALELIZACIÓN DE PROCESO
#===============================================================================
library(doMC)
registerDoMC(cores = 4)

# HIPERPARÁMETROS, NÚMERO DE REPETICIONES Y SEMILLAS PARA CADA REPETICIÓN
#===============================================================================
particiones  <- 10
repeticiones <- 5

# Hiperparámetros
hiperparametros <- data.frame(parameter = "none")

set.seed(123)
seeds <- vector(mode = "list", length = (particiones * repeticiones) + 1)
for (i in 1:(particiones * repeticiones)) {
  seeds[[i]] <- sample.int(1000, nrow(hiperparametros))
}
seeds[[(particiones * repeticiones) + 1]] <- sample.int(1000, 1)

# DEFINICIÓN DEL ENTRENAMIENTO
#===============================================================================
control_train <- trainControl(method = "repeatedcv", number = particiones,
                              repeats = repeticiones, seeds = seeds,
                              returnResamp = "final", verboseIter = FALSE,
                              allowParallel = TRUE)

# AJUSTE DEL MODELO
# ==============================================================================
set.seed(342)
modelo_lda <- train(Survived ~ ., data = datos_train_prep,
                    method = "lda",
                    tuneGrid = hiperparametros,
                    metric = "Accuracy",
                    trControl = control_train)
modelo_lda

## Linear Discriminant Analysis 
## 
## 714 samples
##  10 predictor
##   2 classes: 'No', 'Si' 
## 
## No pre-processing
## Resampling: Cross-Validated (10 fold, repeated 5 times) 
## Summary of sample sizes: 642, 643, 643, 643, 643, 642, ... 
## Resampling results:
## 
##   Accuracy   Kappa    
##   0.8020305  0.5722749

modelo_lda$finalModel

## Call:
## lda(x, grouping = y)
## 
## Prior probabilities of groups:
##        No        Si 
## 0.6162465 0.3837535 
## 
## Group means:
##          SibSp       Parch       Fare Pclass_X2 Pclass_X3  Sex_male Embarked_Q
## No  0.01588641 -0.06492175 -0.2192080 0.1727273 0.6840909 0.8545455 0.09318182
## Si -0.02551103  0.10425390  0.3520128 0.2518248 0.3430657 0.3065693 0.08029197
##    Embarked_S Age_grupo_anciano Age_grupo_niño
## No  0.7590909        0.05454545     0.05227273
## Si  0.6386861        0.01459854     0.10948905
## 
## Coefficients of linear discriminants:
##                           LD1
## SibSp             -0.13809732
## Parch             -0.14008437
## Fare               0.03405313
## Pclass_X2         -0.53108942
## Pclass_X3         -1.27746486
## Sex_male          -2.16552490
## Embarked_Q        -0.02698187
## Embarked_S        -0.26103341
## Age_grupo_anciano -0.48861010
## Age_grupo_niño     1.08303593

Empleando como modelo un LDA, se consigue un accuracy promedio de validación del 80.0%.

Árbol de clasificación simple

Información detallada sobre árboles de clasificación en Árboles de predicción: bagging, random forest, boosting y C5.0.

Ajuste, optimización y validación del modelo

El método C5.0Tree de caret emplea la función C5.0.default() del paquete C50 para crear un árbol de clasificación simple. Este algoritmo no tiene ningún hiperparámetro.

# PARALELIZACIÓN DE PROCESO
#===============================================================================
library(doMC)
registerDoMC(cores = 4)

# HIPERPARÁMETROS, NÚMERO DE REPETICIONES Y SEMILLAS PARA CADA REPETICIÓN
#===============================================================================
particiones  <- 10
repeticiones <- 5

# Hiperparámetros
hiperparametros <- data.frame(parameter = "none")

set.seed(123)
seeds <- vector(mode = "list", length = (particiones * repeticiones) + 1)
for (i in 1:(particiones * repeticiones)) {
  seeds[[i]] <- sample.int(1000, nrow(hiperparametros))
}
seeds[[(particiones * repeticiones) + 1]] <- sample.int(1000, 1)

# DEFINICIÓN DEL ENTRENAMIENTO
#===============================================================================
control_train <- trainControl(method = "repeatedcv", number = particiones,
                              repeats = repeticiones, seeds = seeds,
                              returnResamp = "final", verboseIter = FALSE,
                              allowParallel = TRUE)

# AJUSTE DEL MODELO
# ==============================================================================
set.seed(342)
modelo_C50Tree <- train(Survived ~ ., data = datos_train_prep,
                    method = "C5.0Tree",
                    tuneGrid = hiperparametros,
                    metric = "Accuracy",
                    trControl = control_train)
modelo_C50Tree

## Single C5.0 Tree 
## 
## 714 samples
##  10 predictor
##   2 classes: 'No', 'Si' 
## 
## No pre-processing
## Resampling: Cross-Validated (10 fold, repeated 5 times) 
## Summary of sample sizes: 642, 643, 643, 643, 643, 642, ... 
## Resampling results:
## 
##   Accuracy   Kappa    
##   0.8118114  0.5909209

summary(modelo_C50Tree$finalModel)

## 
## Call:
## C50:::C5.0.default(x = x, y = y, weights = wts)
## 
## 
## C5.0 [Release 2.07 GPL Edition]      Sat Nov  7 23:05:58 2020
## -------------------------------
## 
## Class specified by attribute `outcome'
## 
## Read 714 cases (11 attributes) from undefined.data
## 
## Decision tree:
## 
## Sex_male <= 0:
## :...Pclass_X3 <= 0: Si (139/7)
## :   Pclass_X3 > 0:
## :   :...Fare > -0.1773338: No (19/2)
## :       Fare <= -0.1773338:
## :       :...Embarked_S <= 0:
## :           :...Parch <= -0.4618939: Si (31/9)
## :           :   Parch > -0.4618939:
## :           :   :...Fare <= -0.3244605: No (6/1)
## :           :       Fare > -0.3244605: Si (6)
## :           Embarked_S > 0:
## :           :...Age_grupo_niño > 0:
## :               :...SibSp <= 1.503186: Si (6/1)
## :               :   SibSp > 1.503186: No (2)
## :               Age_grupo_niño <= 0:
## :               :...SibSp > -0.467934: No (17/7)
## :                   SibSp <= -0.467934:
## :                   :...Fare <= -0.2952598: Si (26/11)
## :                       Fare > -0.2952598: No (2)
## Sex_male > 0:
## :...Age_grupo_niño > 0:
##     :...SibSp <= 1.503186: Si (15/1)
##     :   SibSp > 1.503186: No (11/1)
##     Age_grupo_niño <= 0:
##     :...Fare <= -0.1099475: No (326/33)
##         Fare > -0.1099475:
##         :...Pclass_X2 > 0: No (10)
##             Pclass_X2 <= 0:
##             :...Age_grupo_anciano > 0: No (5)
##                 Age_grupo_anciano <= 0:
##                 :...SibSp > 1.503186: No (8)
##                     SibSp <= 1.503186:
##                     :...SibSp > -0.467934:
##                         :...Pclass_X3 <= 0: Si (21/9)
##                         :   Pclass_X3 > 0: No (2)
##                         SibSp <= -0.467934:
##                         :...Pclass_X3 > 0: Si (5/1)
##                             Pclass_X3 <= 0:
##                             :...Fare <= -0.09310098: Si (13/3)
##                                 Fare > -0.09310098: No (44/10)
## 
## 
## Evaluation on training data (714 cases):
## 
##      Decision Tree   
##    ----------------  
##    Size      Errors  
## 
##      21   96(13.4%)   <<
## 
## 
##     (a)   (b)    <-classified as
##    ----  ----
##     398    42    (a): class No
##      54   220    (b): class Si
## 
## 
##  Attribute usage:
## 
##  100.00% Sex_male
##   76.89% Fare
##   71.85% Age_grupo_niño
##   47.48% Pclass_X3
##   24.09% SibSp
##   15.13% Pclass_X2
##   13.73% Age_grupo_anciano
##   13.45% Embarked_S
##    6.02% Parch
## 
## 
## Time: 0.0 secs

Empleando como modelo un árbol simple C5.0, se consigue un accuracy promedio de validación del 81%.

RandomForest

Información detallada sobre random forest en Árboles de predicción: bagging, random forest, boosting y C5.0

Ajuste, optimización y validación del modelo

El método ranger de caret emplea la función ranger() del paquete ranger. Este algoritmo tiene 3 hiperparámetros:

• mtry: número predictores seleccionados aleatoriamente en cada árbol.

• min.node.size: tamaño mínimo que tiene que tener un nodo para poder ser dividido.

• splitrule: criterio de división.

Aunque caret también incluye el método rf con la función rf() del paquete randomForest, este último solo permite optimizar el hiperparámetro mtry.

# PARALELIZACIÓN DE PROCESO
#===============================================================================
library(doMC)
registerDoMC(cores = 4)

# HIPERPARÁMETROS, NÚMERO DE REPETICIONES Y SEMILLAS PARA CADA REPETICIÓN
#===============================================================================
particiones  <- 10
repeticiones <- 5

# Hiperparámetros
hiperparametros <- expand.grid(mtry = c(3, 4, 5, 7),
                               min.node.size = c(2, 3, 4, 5, 10, 15, 20, 30),
                               splitrule = "gini")

set.seed(123)
seeds <- vector(mode = "list", length = (particiones * repeticiones) + 1)
for (i in 1:(particiones * repeticiones)) {
  seeds[[i]] <- sample.int(1000, nrow(hiperparametros))
}
seeds[[(particiones * repeticiones) + 1]] <- sample.int(1000, 1)

# DEFINICIÓN DEL ENTRENAMIENTO
#===============================================================================
control_train <- trainControl(method = "repeatedcv", number = particiones,
                              repeats = repeticiones, seeds = seeds,
                              returnResamp = "final", verboseIter = FALSE,
                              allowParallel = TRUE)

# AJUSTE DEL MODELO
# ==============================================================================
set.seed(342)
modelo_rf <- train(Survived ~ ., data = datos_train_prep,
                   method = "ranger",
                   tuneGrid = hiperparametros,
                   metric = "Accuracy",
                   trControl = control_train,
                   # Número de árboles ajustados
                   num.trees = 500)
modelo_rf

## Random Forest 
## 
## 714 samples
##  10 predictor
##   2 classes: 'No', 'Si' 
## 
## No pre-processing
## Resampling: Cross-Validated (10 fold, repeated 5 times) 
## Summary of sample sizes: 642, 643, 643, 643, 643, 642, ... 
## Resampling results across tuning parameters:
## 
##   mtry  min.node.size  Accuracy   Kappa    
##   3      2             0.8314554  0.6312850
##   3      3             0.8272613  0.6220921
##   3      4             0.8266980  0.6209059
##   3      5             0.8255829  0.6174576
##   3     10             0.8258294  0.6177823
##   3     15             0.8210720  0.6060377
##   3     20             0.8188380  0.6008616
##   3     30             0.8208020  0.6026653
##   4      2             0.8381651  0.6500802
##   4      3             0.8362246  0.6451766
##   4      4             0.8393036  0.6518719
##   4      5             0.8401369  0.6540397
##   4     10             0.8334272  0.6387734
##   4     15             0.8336972  0.6382629
##   4     20             0.8322926  0.6348692
##   4     30             0.8261228  0.6195223
##   5      2             0.8370618  0.6496388
##   5      3             0.8390102  0.6533939
##   5      4             0.8395814  0.6545875
##   5      5             0.8409781  0.6573115
##   5     10             0.8418036  0.6588622
##   5     15             0.8392919  0.6517671
##   5     20             0.8359194  0.6437111
##   5     30             0.8300509  0.6297057
##   7      2             0.8294679  0.6345876
##   7      3             0.8300391  0.6360452
##   7      4             0.8311659  0.6375872
##   7      5             0.8314593  0.6379587
##   7     10             0.8398318  0.6550971
##   7     15             0.8437598  0.6622692
##   7     20             0.8409624  0.6558266
##   7     30             0.8317097  0.6353084
## 
## Tuning parameter 'splitrule' was held constant at a value of gini
## Accuracy was used to select the optimal model using the largest value.
## The final values used for the model were mtry = 7, splitrule = gini
##  and min.node.size = 15.

modelo_rf$finalModel

## Ranger result
## 
## Call:
##  ranger::ranger(dependent.variable.name = ".outcome", data = x,      mtry = min(param$mtry, ncol(x)), min.node.size = param$min.node.size,      splitrule = as.character(param$splitrule), write.forest = TRUE,      probability = classProbs, ...) 
## 
## Type:                             Classification 
## Number of trees:                  500 
## Sample size:                      714 
## Number of independent variables:  10 
## Mtry:                             7 
## Target node size:                 15 
## Variable importance mode:         none 
## Splitrule:                        gini 
## OOB prediction error:             15.97 %

# REPRESENTACIÓN GRÁFICA
# ==============================================================================
ggplot(modelo_rf, highlight = TRUE) +
  scale_x_continuous(breaks = 1:30) +
  labs(title = "Evolución del accuracy del modelo Random Forest") +
  guides(color = guide_legend(title = "mtry"),
         shape = guide_legend(title = "mtry")) +
  theme_bw()

Empleando un modelo random forest con mtry = 7, min.node.size = 15 y splitrule = "gini", se consigue un accuracy promedio de validación del 84%.

Gradient Boosting

Información detallada sobre boosting en Árboles de predicción: bagging, random forest, boosting y C5.0

Ajuste, optimización y validación del modelo

El método gbm de caret emplea la función gbm() del paquete gbm. Este algoritmo tiene 4 hiperparámetros:

• n.trees: número de iteraciones del algoritmo de boosting, es decir, número de modelos que forman el ensemble. Cuanto mayor es este valor, más se reduce el error de entrenamiento, pudiendo llegar generarse overfitting.

• interaction.depth: complejidad de los árboles empleados como weak learner, en concreto, el número total de divisiones que tiene el árbol. Emplear árboles con ente 1 y 6 nodos suele dar buenos resultados.

• shrinkage: este parámetro, también conocido como learning rate, controla la influencia que tiene cada modelo sobre el conjunto del ensemble. Es preferible mejorar un modelo mediante muchos pasos pequeños que mediante unos pocos grandes. Por esta razón, se recomienda emplear un valor de shrinkage tan pequeño como sea posible, teniendo en cuenta que, cuanto menor sea, mayor el número de iteraciones necesarias. Por defecto es de 0.001.

• n.minobsinnode: número mínimo de observaciones que debe tener un nodo para poder ser dividido. Al igual que interaction.depth, permite controlar la complejidad de los weak learners basados en árboles.

Además de los hiperparámetros que permite controlar caret, la función gbm() tiene otros dos más que hay que tener en cuenta:

• distribution: determina la función de coste (loss function). Algunas de las más utilizadas son: gaussian (squared loss) para regresión, bernoulli para respuestas binarias, multinomial para variables respuesta con más de dos clases y adaboost para respuestas binarias y que emplea la función exponencial del algoritmo original AdaBoost.

• bag.fraction (subsampling fraction): fracción de observaciones del set de entrenamiento seleccionadas de forma aleatoria para ajustar cada weak learner. Si su valor es de 1, se emplea el algoritmo de Gradient Boosting, si es menor que 1, se emplea Stochastic Gradient Boosting. Por defecto su valor es de 0.5.

# PARALELIZACIÓN DE PROCESO
#===============================================================================
library(doMC)
registerDoMC(cores = 4)

# HIPERPARÁMETROS, NÚMERO DE REPETICIONES Y SEMILLAS PARA CADA REPETICIÓN
#===============================================================================
particiones  <- 10
repeticiones <- 5

# Hiperparámetros
hiperparametros <- expand.grid(interaction.depth = c(1, 2),
                               n.trees = c(500, 1000, 2000),
                               shrinkage = c(0.001, 0.01, 0.1),
                               n.minobsinnode = c(2, 5, 15))

set.seed(123)
seeds <- vector(mode = "list", length = (particiones * repeticiones) + 1)
for (i in 1:(particiones * repeticiones)) {
  seeds[[i]] <- sample.int(1000, nrow(hiperparametros))
}
seeds[[(particiones * repeticiones) + 1]] <- sample.int(1000, 1)

# DEFINICIÓN DEL ENTRENAMIENTO
#===============================================================================
control_train <- trainControl(method = "repeatedcv", number = particiones,
                              repeats = repeticiones, seeds = seeds,
                              returnResamp = "final", verboseIter = FALSE,
                              allowParallel = TRUE)

# AJUSTE DEL MODELO
# ==============================================================================
set.seed(342)
modelo_boost <- train(Survived ~ ., data = datos_train_prep,
                   method = "gbm",
                   tuneGrid = hiperparametros,
                   metric = "Accuracy",
                   trControl = control_train,
                   # Número de árboles ajustados
                   distribution = "adaboost",
                   verbose = FALSE)
modelo_boost

## Stochastic Gradient Boosting 
## 
## 714 samples
##  10 predictor
##   2 classes: 'No', 'Si' 
## 
## No pre-processing
## Resampling: Cross-Validated (10 fold, repeated 5 times) 
## Summary of sample sizes: 642, 643, 643, 643, 643, 642, ... 
## Resampling results across tuning parameters:
## 
##   shrinkage  interaction.depth  n.minobsinnode  n.trees  Accuracy   Kappa    
##   0.001      1                   2               500     0.7899609  0.5485232
##   0.001      1                   2              1000     0.7927660  0.5554307
##   0.001      1                   2              2000     0.7927660  0.5554307
##   0.001      1                   5               500     0.7913576  0.5519772
##   0.001      1                   5              1000     0.7927660  0.5554307
##   0.001      1                   5              2000     0.7927660  0.5554307
##   0.001      1                  15               500     0.7910759  0.5511284
##   0.001      1                  15              1000     0.7927660  0.5554307
##   0.001      1                  15              2000     0.7927660  0.5554307
##   0.001      2                   2               500     0.7914124  0.5122589
##   0.001      2                   2              1000     0.7788028  0.5013913
##   0.001      2                   2              2000     0.7939085  0.5505586
##   0.001      2                   5               500     0.7914124  0.5122589
##   0.001      2                   5              1000     0.7804734  0.5046415
##   0.001      2                   5              2000     0.7905516  0.5426988
##   0.001      2                  15               500     0.7914124  0.5122589
##   0.001      2                  15              1000     0.7821674  0.4999266
##   0.001      2                  15              2000     0.7902387  0.5415659
##   0.010      1                   2               500     0.7860329  0.5401062
##   0.010      1                   2              1000     0.7927387  0.5550821
##   0.010      1                   2              2000     0.7930282  0.5561399
##   0.010      1                   5               500     0.7846440  0.5368845
##   0.010      1                   5              1000     0.7941432  0.5579070
##   0.010      1                   5              2000     0.7935915  0.5575119
##   0.010      1                  15               500     0.7804186  0.5285964
##   0.010      1                  15              1000     0.7857355  0.5424127
##   0.010      1                  15              2000     0.7905008  0.5515624
##   0.010      2                   2               500     0.8070618  0.5806507
##   0.010      2                   2              1000     0.8073200  0.5835410
##   0.010      2                   2              2000     0.8120931  0.5927468
##   0.010      2                   5               500     0.8017449  0.5701443
##   0.010      2                   5              1000     0.8070344  0.5823929
##   0.010      2                   5              2000     0.8109859  0.5904235
##   0.010      2                  15               500     0.7975274  0.5584389
##   0.010      2                  15              1000     0.8047887  0.5752401
##   0.010      2                  15              2000     0.8081768  0.5808617
##   0.100      1                   2               500     0.7949883  0.5597297
##   0.100      1                   2              1000     0.7963967  0.5615612
##   0.100      1                   2              2000     0.8064828  0.5852633
##   0.100      1                   5               500     0.7938928  0.5575167
##   0.100      1                   5              1000     0.8008803  0.5710390
##   0.100      1                   5              2000     0.8062050  0.5826722
##   0.100      1                  15               500     0.8022770  0.5749268
##   0.100      1                  15              1000     0.8025822  0.5753210
##   0.100      1                  15              2000     0.8140728  0.6001542
##   0.100      2                   2               500     0.8331729  0.6390100
##   0.100      2                   2              1000     0.8236307  0.6189589
##   0.100      2                   2              2000     0.8216471  0.6170511
##   0.100      2                   5               500     0.8314710  0.6346001
##   0.100      2                   5              1000     0.8250235  0.6230662
##   0.100      2                   5              2000     0.8233333  0.6208026
##   0.100      2                  15               500     0.8230556  0.6146507
##   0.100      2                  15              1000     0.8252934  0.6205955
##   0.100      2                  15              2000     0.8258372  0.6231925
## 
## Accuracy was used to select the optimal model using the largest value.
## The final values used for the model were n.trees = 500, interaction.depth =
##  2, shrinkage = 0.1 and n.minobsinnode = 2.

# REPRESENTACIÓN GRÁFICA
# ==============================================================================
ggplot(modelo_boost, highlight = TRUE) +
  labs(title = "Evolución del accuracy del modelo Gradient Boosting") +
    guides(color = guide_legend(title = "shrinkage"),
           shape = guide_legend(title = "shrinkage")) +
  theme_bw() +
  theme(legend.position = "bottom")

Empleando un modelo boosting con n.trees = 500, minobsinnode = 2, interaction.depth = 2 y shrinkage = 0.1 se consigue un accuracy promedio de validación del 83%.

SVM

Información detallada sobre SVM en Máquinas de Vector Soporte (Support Vector Machines, SVMs)

Ajuste, optimización y validación del modelo

El método svmRadial de caret emplea la función ksvm() del paquete kernlab. Este algoritmo tiene 2 hiperparámetros:

• sigma: coeficiente del kernel radial.

• C: penalización por violaciones del margen del hiperplano.

# PARALELIZACIÓN DE PROCESO
#===============================================================================
library(doMC)
registerDoMC(cores = 4)

# HIPERPARÁMETROS, NÚMERO DE REPETICIONES Y SEMILLAS PARA CADA REPETICIÓN
#===============================================================================
particiones  <- 10
repeticiones <- 5

# Hiperparámetros
hiperparametros <- expand.grid(sigma = c(0.001, 0.01, 0.1, 0.5, 1),
                               C = c(1 , 20, 50, 100, 200, 500, 700))

set.seed(123)
seeds <- vector(mode = "list", length = (particiones * repeticiones) + 1)
for (i in 1:(particiones * repeticiones)) {
  seeds[[i]] <- sample.int(1000, nrow(hiperparametros))
}
seeds[[(particiones * repeticiones) + 1]] <- sample.int(1000, 1)

# DEFINICIÓN DEL ENTRENAMIENTO
#===============================================================================
control_train <- trainControl(method = "repeatedcv", number = particiones,
                              repeats = repeticiones, seeds = seeds,
                              returnResamp = "final", verboseIter = FALSE,
                              allowParallel = TRUE)

# AJUSTE DEL MODELO
# ==============================================================================
set.seed(342)
modelo_svmrad <- train(Survived ~ ., data = datos_train_prep,
                   method = "svmRadial",
                   tuneGrid = hiperparametros,
                   metric = "Accuracy",
                   trControl = control_train)
modelo_svmrad

## Support Vector Machines with Radial Basis Function Kernel 
## 
## 714 samples
##  10 predictor
##   2 classes: 'No', 'Si' 
## 
## No pre-processing
## Resampling: Cross-Validated (10 fold, repeated 5 times) 
## Summary of sample sizes: 642, 643, 643, 643, 643, 642, ... 
## Resampling results across tuning parameters:
## 
##   sigma  C    Accuracy   Kappa    
##   0.001    1  0.7749061  0.4872109
##   0.001   20  0.7927660  0.5554307
##   0.001   50  0.7938928  0.5570198
##   0.001  100  0.8031377  0.5763943
##   0.001  200  0.8219053  0.6170163
##   0.001  500  0.8277973  0.6289610
##   0.001  700  0.8277973  0.6290988
##   0.010    1  0.7972692  0.5631272
##   0.010   20  0.8272418  0.6272218
##   0.010   50  0.8210994  0.6094250
##   0.010  100  0.8183099  0.5996508
##   0.010  200  0.8157786  0.5926998
##   0.010  500  0.8104538  0.5832212
##   0.010  700  0.8068036  0.5765461
##   0.100    1  0.8160563  0.5949255
##   0.100   20  0.8115571  0.5908364
##   0.100   50  0.8051017  0.5765383
##   0.100  100  0.8045149  0.5750284
##   0.100  200  0.8042254  0.5745622
##   0.100  500  0.7972496  0.5602523
##   0.100  700  0.7950313  0.5568971
##   0.500    1  0.8112950  0.5862754
##   0.500   20  0.8070501  0.5849130
##   0.500   50  0.8045462  0.5807402
##   0.500  100  0.8014632  0.5742273
##   0.500  200  0.8034233  0.5776361
##   0.500  500  0.8023122  0.5741991
##   0.500  700  0.8011972  0.5711657
##   1.000    1  0.8073435  0.5778381
##   1.000   20  0.8039906  0.5772357
##   1.000   50  0.8025978  0.5746756
##   1.000  100  0.8009116  0.5701290
##   1.000  200  0.8008998  0.5698152
##   1.000  500  0.7950274  0.5558312
##   1.000  700  0.7939045  0.5531533
## 
## Accuracy was used to select the optimal model using the largest value.
## The final values used for the model were sigma = 0.001 and C = 500.

modelo_svmrad$finalModel

## Support Vector Machine object of class "ksvm" 
## 
## SV type: C-svc  (classification) 
##  parameter : cost C = 500 
## 
## Gaussian Radial Basis kernel function. 
##  Hyperparameter : sigma =  0.001 
## 
## Number of Support Vectors : 296 
## 
## Objective Function Value : -131506.7 
## Training error : 0.168067

# REPRESENTACIÓN GRÁFICA
# ==============================================================================
ggplot(modelo_svmrad, highlight = TRUE) +
  labs(title = "Evolución del accuracy del modelo SVM Radial") +
  theme_bw()

Empleando un modelo SVM Radial con sigma = 0.001 y C = 500, se consigue un accuracy promedio de validación del 83%.

Redes neuronales (NNET)

Ajuste, optimización y validación del modelo

El método nnet de caret emplea la función nnet() del paquete nnet para crear redes neuronales con una capa oculta. Este algoritmo tiene 2 hiperparámetros:

• size: número de neuronas en la capa oculta.

• decay: controla la regularización durante el entrenamiento de la red.

Además de estos hiperparámetros, la función nnet() tiene muchos otros argumentos que controlan el proceso de aprendizaje de la red. Entre ellos cabe destacar MaxNWts, que determina el número máximo de pesos. Por defecto, su valor es 1000, pero puede ocurrir que, si la red tiene muchas neuronas, se necesiten más pesos.

# PARALELIZACIÓN DE PROCESO
#===============================================================================
library(doMC)
registerDoMC(cores = 4)

# HIPERPARÁMETROS, NÚMERO DE REPETICIONES Y SEMILLAS PARA CADA REPETICIÓN
#===============================================================================
particiones  <- 10
repeticiones <- 5

# Hiperparámetros
hiperparametros <- expand.grid(size = c(10, 20, 50, 80, 100, 120),
                               decay = c(0.0001, 0.1, 0.5))

set.seed(123)
seeds <- vector(mode = "list", length = (particiones * repeticiones) + 1)
for (i in 1:(particiones * repeticiones)) {
  seeds[[i]] <- sample.int(1000, nrow(hiperparametros))
}
seeds[[(particiones * repeticiones) + 1]] <- sample.int(1000, 1)

# DEFINICIÓN DEL ENTRENAMIENTO
#===============================================================================
control_train <- trainControl(method = "repeatedcv", number = particiones,
                              repeats = repeticiones, seeds = seeds,
                              returnResamp = "final", verboseIter = FALSE,
                              allowParallel = TRUE)

# AJUSTE DEL MODELO
# ==============================================================================
set.seed(342)
modelo_nnet <- train(Survived ~ ., data = datos_train_prep,
                   method = "nnet",
                   tuneGrid = hiperparametros,
                   metric = "Accuracy",
                   trControl = control_train,
                   # Rango de inicialización de los pesos
                   rang = c(-0.7, 0.7),
                   # Número máximo de pesos
                   # se aumenta para poder incluir más meuronas
                   MaxNWts = 2000,
                   # Para que no se muestre cada iteración por pantalla
                   trace = FALSE)
modelo_nnet