Resumen

Este experimento pretende probar varias formas de modelar una predicción categórica, bien sea ordinal o no, a través de distintas técnicas como ordenar, balancear, transformar los valores categóricos a numéricos para hacer una regresión.

Para ello, he escogido el dataset de Titanic, en el que tenemos quién sobrevivió y quién no, la edad, género, edad, cuánto pagó, cuántos familiares viajaron con cada pasajero, clase de pasaje. Tengamos en cuenta que este dataset no es muy grande por lo que un par de casos clasificados bien o mal pueden variar notablemente algunos resultados.

Para nuestro experimento, queremos predecir a qué clase pertenecía cada pasajero (1era, 2da ó 3ra). Como se darán cuenta, estos valores para Pclass parecen numéricos pero realmente deberían ser categóricos, ordinales. Categóricos porque tenemos apenas 3 valores posibles en todos los pasajeros; fíjate que podrían haber sido perfectamente clases A, B y C. Ordinales porque hay un orden lógico, creciente o decreciente: no es lo mismo confundir 1era clase con 3era clase que 2da con 3ra; hay una referencia importante de posición o aproximación que un valor categórico normal no tiene.

¿Cómo podemos ayudar al modelo a entender las relaciones entre variables categóricas ordinales?

LIBRERÍAS

library(dplyr)
library(ggplot2)
library(lares)

DATOS

Cargamos los datos de Titanic. Este comando importará los datos y creará un data.frame llamado dft.

data(dft)

ESTRUCTURA: Chequeamos rápidamente en qué consisten los datos, cuántas columnas y filas tenemos, frecuencias, etc.

# Cómo está compuesto dft?
df_str(dft, quiet = TRUE)

# Veamos las primeras 10 filas
head(dft, 10)

# Veamos la frecuencia de todos los valores de cada columna
freqs(dft)

## 1 variables with more than 0.9 variance exluded: 'PassengerId'

TRANSFORMACIONES: Ahora transformamos Cabin por su alta varianza y nos deshacemos de PassengerId y Ticket ya que no tiene mayor sentido usarlos para este ejercicio predictivo.

dft <- dft %>% 
  select(-Ticket, -PassengerId) %>%
  mutate(Cabin = ifelse(Cabin == "", "NON", "WITH"))

ENTENDIMIENTO: Veamos si existen correlaciones que tengan sentido entre todas las variables:

corr_cross(dft)

De aquí, sin haber visto la película de Titanic, podemos entender qué pasó: 1era clase cuenta con cabinas (Cabin) definidas, pagaron (Fare) más y suelen ser mayores de edad (Age) que las otras clases. Los hombres (Sex) sobrevivieron (Survived) en menores proporciones. Los de 3era clase no tienen cabina asignada y pagaron menos. Los que tienen muchos hijos y padres abordo (parch) también suelen tener hermanos y cónyuges (SibSp), es decir, viajan en familias numerosas.

VALORES FALTANTES: Estudiemos los valores nulos o vacíos de nuestros datos:

missingness(dft, full = TRUE, plot = TRUE)

Podríamos imputar estos valores vacíos de Age (tiene 20% de sus valores faltantes), pero por simplicidad y para considerar en los siguientes experimentos, no lo haremos ahora. Si se quisiera hacer, se puede de la siguiente manera: dft <- impute(dft). Te dejamos ese experimento para que compares resultados ;)

VARIABLE A MODELAR: Analicemos la variable que vamos a modelar. Fíjate que la clase más popular y frecuente es la 3era, seguido por 1era. Un poco más de la mitad de los pasajeros se encuentran en 3era clase.

freqs(dft, Pclass)

EL EXPERIMENTO

Comparemos ahora distintas opciones de modelos bajo una misma métrica global: exactitud (accuracy). La exactitud es la cantidad de casos bien “adivinados” vs la cantidad total de casos.

NOTA: Todos los mensajes automáticos (muy útiles e informativos) que aparecen cuando corremos el h2o_automl serán excluidos de este reporte para tener una mejor y comprensible lectura. No recomiendo dejar quiet = TRUE si deseas replicar los ejemplos.

Modelo 1: Modelo balanceado (91.79%)

En este primer experimento, vamos a modelar balanceando todo el dataset para que hayan igual cantidad de casos para cada clase. Es decir, como 2da clase tiene menos pasajeros, vamos a hacer “under-sampling” de las otras clases para que tengan la misma cantidad de pasajeros en cada clase y no exista una sobre representación.

model1 <- h2o_automl(dft, "Pclass", max_models = 4, balance = TRUE, quiet = TRUE)

## Model type: Classifier

##       AUC     ACC
## 1 0.98663 0.91791

# mplot_full(model1$scores_test$tag, model1$scores_test$score,
#            select(model1$scores_test, -tag, -score))

model1$metrics$plots$conf_matrix

model1$scores_test %>% 
  mutate(correct = tag == score) %>% 
  freqs(score, correct, rel = TRUE) %>% 
  filter(correct == TRUE)

El resultado es bastante aceptable pero aún no tenemos modelos para comparar. Sigamos…

Modelo 2: Modelo no balanceado (93.28%)

Ahora corramos exactamente el mismo modelo pero sin balancear. Tomaremos todos los datos que tenemos para cada una de las 3 clases.

model2 <- h2o_automl(dft, "Pclass", max_models = 4, quiet = TRUE)

## Model type: Classifier

##       AUC     ACC
## 1 0.99058 0.93284

# mplot_full(model2$scores_test$tag, model2$scores_test$score,
#            select(model2$scores_test, -tag, -score))

model2$metrics$plots$conf_matrix

model2$scores_test %>% mutate(correct = tag == score) %>% 
  freqs(score, correct, rel = TRUE) %>% filter(correct == TRUE)

Vemos que los resultados son mejores esta vez. Comparando la exactitud, mejoramos 1.5 puntos porcentuales. Viendo clase por clase, empeoramos un poco en 3era clase (confundiéndolos un poco más con 2da pero menos en 1era), nos mantenemos muy bien con 1era clase y mejoramos 2da clase.

Modelo 3: Modelo como regresión (93.28% - 94.03%)

Ahora, modelemos las clases pero como valores numéricos ordenados. Tal y como podrían haber sido clases A, B y C, transformamos las clases a 1, 2 y 3. Luego, podríamos separar con umbrales 1.5 y 2.5 cada una de las 3 etiquetas. Por último, probamos si busando otros valores distintos a 1.5 conseguimos mayor exactitud.

model3 <-  dft %>% mutate(Pclass = as.numeric(Pclass)) %>%
  h2o_automl("Pclass", max_models = 4, thresh = 0, quiet = TRUE)

## Model type: Regression

##        rmse       mae       mape        mse    rsq  rsqa
## 1 0.2374868 0.1189057 0.04452558 0.05639998 0.9183 0.918

# mplot_full(model3$scores_test$tag, model3$scores_test$score, thresh = 0)

conf1 <- model3$scores_test %>% 
  mutate(label = case_when(
    score < 1.5 ~ "1", 
    score > 2.5 ~ "3", 
    TRUE ~ "2")) %>% 
  mutate(correct = tag == label)
mplot_conf(conf1$tag, conf1$label)

conf1 %>% freqs(tag, correct, rel = TRUE) %>% arrange(tag, desc(correct))

Interesante: vemos que estos resultados son idénticos que con el modelo anterior.

Ahora, podremos mejorar los resultados si logramos conseguir la combinación de umbrales que maximizan la exactitud global de las predicciones.

a <- seq(from = 1, to = 2, length.out = 10)
b <- seq(from = 2, to = 3, length.out = 50)
c <- expand.grid(a,b)
for (i in 1:nrow(c)) {
  if (i == 1) out <- data.frame(a = 0, b = 0, acc = 0)
  acc <- model3$scores_test %>% 
    mutate(label = case_when(
      score < c$Var1[i] ~ "1", 
      score > c$Var2[i] ~ "3", 
      TRUE ~ "2")) %>% 
    mutate(correct = tag == label) %>%
    freqs(correct) %>% filter(correct == TRUE) %>% .$p
  v <- data.frame(a = c$Var1[i], b = c$Var2[i], acc = acc)
  if (v$acc > max(out$acc)) out <- rbind(out, v)
}
print(tail(out, 1))

##           a        b   acc
## 20 1.444444 2.244898 94.03

conf2 <- model3$scores_test %>% 
  mutate(label = case_when(
    score < out$a[nrow(out)] ~ "1", 
    score > out$b[nrow(out)] ~ "3", 
    TRUE ~ "2")) %>% 
  mutate(correct = tag == label)
mplot_conf(conf2$tag, conf2$label)

conf2 %>% freqs(tag, correct, rel = TRUE) %>% arrange(tag, desc(correct))

Con estos nuevos cortes, vemos que 2da clase empeora (3 pasajeros) pero 3era mejora (5 pasajeros). Se sacrifica un poco la calidad de predicción pero se tiene mayor exactitud global. 1era clase queda intacta.

Veámoslo gráficamente:

Modelo 4: Otro algoritmo, categóricos ordenados (84,70%)

dft2 <- impute(dft, quiet = TRUE) # MASS::polr no admite valores nulos
dft2$Pclass <- factor(dft2$Pclass, ordered = TRUE)
split <- msplit(dft2, print = FALSE)
model4 <- MASS::polr(Pclass ~ ., data = split$train, Hess = TRUE)
summary(model4)

## Call:
## MASS::polr(formula = Pclass ~ ., data = split$train, Hess = TRUE)
## 
## Coefficients:
##                 Value Std. Error t value
## SurvivedTRUE -1.05394   0.299526  -3.519
## Sexmale      -0.45979   0.309646  -1.485
## Age          -0.04891   0.009202  -5.315
## SibSp         1.00288   0.227004   4.418
## Parch         0.51841   0.179020   2.896
## Fare         -0.12581   0.013795  -9.120
## CabinWITH    -2.81847   0.409023  -6.891
## EmbarkedC     6.87947   0.373512  18.418
## EmbarkedQ    10.01670   0.811359  12.346
## EmbarkedS     6.83830   0.326874  20.920
## 
## Intercepts:
##     Value   Std. Error t value
## 1|2 -0.6988  0.5453    -1.2817
## 2|3  2.6331  0.4284     6.1462
## 
## Residual Deviance: 544.967 
## AIC: 568.967

split$test$pred <- predict(model4, split$test)
100*sum(as.character(split$test$Pclass) == as.character(split$test$pred))/nrow(split$test)

## [1] 84.70149

No vale la pena continuar con este algoritmo por su mal desempeño para este caso en particular.

RESULTADOS

Modelo	Nombre del modelo	Exactitud
#1	Balanceado	91.79
#2	Sin balancear	93.27
#3.1	Regresión	93.30
#3.2	Regresión cortes óptimos	94.03
#4	Ordered Logistic con MASS	84.70

RECOMENDACIONES Y CONCLUSIONES

En todo momento nos tenemos que preguntar: en qué queremos ser mejores? Separando 2da de 3era clase? O mejor separar 3ra y 2da clase de la 1era? Queremos clasificar todos lo mejor posible pero sacrificando una clase en particular?
Cada caso es único así que no hay una respuesta correcta; se recomienda probar distintas combinaciones antes de llegar al modelo final.
Si tenemos dos categorías muy similares, consideremos unirlas si es irrelevante hacerlo. Si no, se podría hacer un segundo modelo para poder diferenciarlos dado cierto umbral o volumen. Así, cuando ciertas observaciones pasen ese umbral, podemos pasarlos por un segundo modelo que es “mejor” clasificando o separandos esas otras categorías.
Estos valores de exactitud dependen totalmente de los datos, el volumen y las proporciones entre cada una de las categorías. Si tenemos 5% de un caso, 5% de otro y 90% de otro, sin duda tendremos que balancear de alguna forma para que esos 5% no parezcan ruido ante el modelo y los trate por igual.
Del mismo modo, y lo digo con pinzas, si queremos que el modelo sea “mejor” prediciendo una categoría, podemos mostrarle más observaciones de esa cateogría.

Modelos Categóricos no Binarios

Bernardo Lares

2019-07-10